Reference Manual

Go to the documentation of this file.
00001 // -*- C++ -*
00002 //
00003 // Package:    IsolatedTracksNxN
00004 // Class:      IsolatedTracksNxN
00005 // 
00013 //
00014 // Original Author:  Seema Sharma
00015 //         Created:  Mon Aug 10 15:30:40 CST 2009
00016 //
00017 //
00018 
00019 #include "DataFormats/TrackReco/interface/TrackBase.h"
00020 #include "Calibration/IsolatedParticles/plugins/IsolatedTracksNxN.h"
00021 
00022 #include "Calibration/IsolatedParticles/interface/CaloPropagateTrack.h"
00023 #include "Calibration/IsolatedParticles/interface/ChargeIsolation.h"
00024 #include "Calibration/IsolatedParticles/interface/eHCALMatrixExtra.h"
00025 #include "MagneticField/Engine/interface/MagneticField.h"
00026 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00027 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
00028 #include "Geometry/CaloTopology/interface/EcalTrigTowerConstituentsMap.h"
00029 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
00030 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
00031 #include "CondFormats/L1TObjects/interface/L1GtTriggerMask.h"
00032 #include "CondFormats/DataRecord/interface/L1GtTriggerMaskAlgoTrigRcd.h"
00033 #include "CondFormats/DataRecord/interface/L1GtTriggerMaskTechTrigRcd.h"
00034 
00035 IsolatedTracksNxN::IsolatedTracksNxN(const edm::ParameterSet& iConfig) {
00036 
00037   //now do what ever initialization is needed
00038   doMC                   = iConfig.getUntrackedParameter<bool>  ("DoMC", false); 
00039   writeAllTracks         = iConfig.getUntrackedParameter<bool>  ("WriteAllTracks", false ); 
00040   myverbose_             = iConfig.getUntrackedParameter<int>   ("Verbosity", 5          );
00041   pvTracksPtMin_         = iConfig.getUntrackedParameter<double>("PVTracksPtMin", 1.0    );
00042   debugTrks_             = iConfig.getUntrackedParameter<int>   ("DebugTracks"           );
00043   printTrkHitPattern_    = iConfig.getUntrackedParameter<bool>  ("PrintTrkHitPattern"    );
00044   minTrackP_             = iConfig.getUntrackedParameter<double>("minTrackP", 1.0        );
00045   maxTrackEta_           = iConfig.getUntrackedParameter<double>("maxTrackEta", 5.0      );
00046   debugL1Info_           = iConfig.getUntrackedParameter<bool>  ("DebugL1Info", false    );
00047   L1TriggerAlgoInfo_     = iConfig.getUntrackedParameter<bool>  ("L1TriggerAlgoInfo");
00048   L1extraTauJetSource_   = iConfig.getParameter<edm::InputTag>  ("L1extraTauJetSource"   );
00049   L1extraCenJetSource_   = iConfig.getParameter<edm::InputTag>  ("L1extraCenJetSource"   );
00050   L1extraFwdJetSource_   = iConfig.getParameter<edm::InputTag>  ("L1extraFwdJetSource"   );
00051   L1extraMuonSource_     = iConfig.getParameter<edm::InputTag>  ("L1extraMuonSource"     );
00052   L1extraIsoEmSource_    = iConfig.getParameter<edm::InputTag>  ("L1extraIsoEmSource"    );
00053   L1extraNonIsoEmSource_ = iConfig.getParameter<edm::InputTag>  ("L1extraNonIsoEmSource" );
00054   L1GTReadoutRcdSource_  = iConfig.getParameter<edm::InputTag>  ("L1GTReadoutRcdSource"  );
00055   L1GTObjectMapRcdSource_= iConfig.getParameter<edm::InputTag>  ("L1GTObjectMapRcdSource");
00056   JetSrc_                = iConfig.getParameter<edm::InputTag>  ("JetSource");
00057   JetExtender_           = iConfig.getParameter<edm::InputTag>  ("JetExtender");
00058   HBHERecHitSource_      = iConfig.getParameter<edm::InputTag>  ("HBHERecHitSource");
00059   tMinE_                 = iConfig.getUntrackedParameter<double>("TimeMinCutECAL", -500.);
00060   tMaxE_                 = iConfig.getUntrackedParameter<double>("TimeMaxCutECAL",  500.);
00061   tMinH_                 = iConfig.getUntrackedParameter<double>("TimeMinCutHCAL", -500.);
00062   tMaxH_                 = iConfig.getUntrackedParameter<double>("TimeMaxCutHCAL",  500.);
00063   nbad                   = 0;
00064 
00065   if(myverbose_>=0) {
00066     std::cout <<"Parameters read from config file \n" 
00067               <<" doMC         "               << doMC
00068               <<"\t myverbose_ "               << myverbose_        
00069               <<"\t minTrackP_ "               << minTrackP_     
00070               << "\t maxTrackEta_ "            << maxTrackEta_
00071               << "\t tMinE_ "                  << tMinE_
00072               << "\t tMaxE_ "                  << tMaxE_
00073               << "\t tMinH_ "                  << tMinH_
00074               << "\t tMaxH_ "                  << tMaxH_
00075               << "\n debugL1Info_ "            << debugL1Info_  
00076               << "\t L1TriggerAlgoInfo_ "      << L1TriggerAlgoInfo_
00077               << "\t L1extraTauJetSource_ "    << L1extraTauJetSource_
00078               << "\t L1extraCenJetSource_ "    << L1extraCenJetSource_
00079               << "\t L1extraFwdJetSource_ "    << L1extraFwdJetSource_
00080               << "\t L1extraMuonSource_ "      << L1extraMuonSource_
00081               << "\t L1extraIsoEmSource_ "     << L1extraIsoEmSource_
00082               << "\t L1extraNonIsoEmSource_ "  << L1extraNonIsoEmSource_
00083               << "\t L1GTReadoutRcdSource_ "   << L1GTReadoutRcdSource_
00084               << "\t L1GTObjectMapRcdSource_ " << L1GTObjectMapRcdSource_
00085               << "\t JetSrc_  "                << JetSrc_
00086               << "\t JetExtender_ "            << JetExtender_
00087               << "\t HBHERecHitSource_ "       << HBHERecHitSource_ 
00088               << "\n" << std::endl;
00089   }
00090 
00091   initL1 = false;
00092 
00093 }
00094 
00095 IsolatedTracksNxN::~IsolatedTracksNxN() {
00096 }
00097 
00098 void IsolatedTracksNxN::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
00099 
00100   bool haveIsoTrack=false;  
00101 
00102   edm::ESHandle<MagneticField> bFieldH;
00103   iSetup.get<IdealMagneticFieldRecord>().get(bFieldH);
00104   bField = bFieldH.product();
00105 
00106   clearTreeVectors();
00107 
00108   t_RunNo = iEvent.id().run();
00109   t_EvtNo = iEvent.id().event();
00110   t_Lumi  = iEvent.luminosityBlock();
00111   t_Bunch = iEvent.bunchCrossing();
00112 
00113   nEventProc++;
00114 
00115   edm::Handle<reco::TrackCollection> trkCollection;
00116   iEvent.getByLabel("generalTracks", trkCollection);
00117   reco::TrackCollection::const_iterator trkItr;
00118   if(debugTrks_>1){
00119     std::cout << "Track Collection: " << std::endl;
00120     std::cout << "Number of Tracks " << trkCollection->size() << std::endl;
00121   }
00122   std::string theTrackQuality = "highPurity";
00123   reco::TrackBase::TrackQuality trackQuality_=reco::TrackBase::qualityByName(theTrackQuality);
00124 
00125   //===================== save L1 Trigger information =======================
00126   if( L1TriggerAlgoInfo_ ) {
00127 
00128     bool useL1EventSetup = true;
00129     bool useL1GtTriggerMenuLite = true;
00130 
00131     m_l1GtUtils.getL1GtRunCache(iEvent, iSetup, useL1EventSetup, useL1GtTriggerMenuLite);
00132 
00133     int iErrorCode = -1;
00134     int l1ConfCode = -1;
00135     const bool l1Conf = m_l1GtUtils.availableL1Configuration(iErrorCode, l1ConfCode);
00136     if( !l1Conf) { 
00137       std::cout << "\nL1 configuration code:" << l1ConfCode
00138                 << "\nNo valid L1 trigger configuration available."
00139                 << "\nSee text above for error code interpretation"
00140                 << "\nNo return here, in order to test each method, protected against configuration error."
00141                 << std::endl;
00142     }
00143 
00144     const L1GtTriggerMenu* m_l1GtMenu   = m_l1GtUtils.ptrL1TriggerMenuEventSetup(iErrorCode);
00145     const AlgorithmMap&    algorithmMap = m_l1GtMenu->gtAlgorithmMap();
00146     const std::string&     menuName     = m_l1GtMenu->gtTriggerMenuName();
00147 
00148     if (!initL1) {
00149       initL1=true;
00150       std::cout << "menuName " << menuName << std::endl;
00151       for (CItAlgo itAlgo = algorithmMap.begin(); itAlgo != algorithmMap.end(); itAlgo++) {
00152         std::string algName = itAlgo->first;
00153         int algBitNumber = ( itAlgo->second ).algoBitNumber();
00154         l1AlgoMap.insert( std::pair<std::pair<unsigned int,std::string>,int>( std::pair<unsigned int,std::string>(algBitNumber, algName) , 0) ) ;
00155       }
00156       std::map< std::pair<unsigned int,std::string>, int>::iterator itr;
00157       for(itr=l1AlgoMap.begin(); itr!=l1AlgoMap.end(); itr++) {
00158         std::cout << " ********** " << (itr->first).first <<" "<<(itr->first).second <<" "<<itr->second << std::endl;
00159       }
00160     }
00161 
00162     std::vector<int> algbits;
00163     for (CItAlgo itAlgo        = algorithmMap.begin(); itAlgo != algorithmMap.end(); itAlgo++) {
00164       std::string algName      = itAlgo->first;
00165       int algBitNumber         = ( itAlgo->second ).algoBitNumber();
00166       //bool algResultBeforeMask = m_l1GtUtils.decisionBeforeMask(iEvent, itAlgo->first, iErrorCode);
00167       //bool algResultAfterMask  = m_l1GtUtils.decisionAfterMask (iEvent, itAlgo->first, iErrorCode);
00168       //int  triggerMask         = m_l1GtUtils.triggerMask       (iEvent, itAlgo->first, iErrorCode);
00169       bool decision            = m_l1GtUtils.decision          (iEvent, itAlgo->first, iErrorCode);
00170       int  preScale            = m_l1GtUtils.prescaleFactor    (iEvent, itAlgo->first, iErrorCode);
00171     
00172       // save the algo names which fired
00173       if( decision ){
00174         l1AlgoMap[std::pair<unsigned int,std::string>(algBitNumber, algName)] += 1;
00175         h_L1AlgoNames->Fill( algBitNumber );
00176         t_L1AlgoNames->push_back(itAlgo->first); 
00177         t_L1PreScale->push_back(preScale);
00178         t_L1Decision[algBitNumber] = 1;
00179         algbits.push_back(algBitNumber);
00180       }
00181     }
00182 
00183     if(debugL1Info_) {
00184       for(unsigned int ii=0; ii<t_L1AlgoNames->size(); ii++){
00185         std::cout<<ii<<" "<<(*t_L1AlgoNames)[ii]<<" "<<(*t_L1PreScale)[ii]<<" "<<algbits[ii]<<std::endl;
00186       }
00187       std::cout << "L1Decision: ";
00188       for (int i=0; i<128; ++i) {
00189         std::cout << " " << i << ":" << t_L1Decision[i];
00190         if ((i+1)%32 == 0) std::cout << std::endl;
00191       }
00192     }
00193 
00194     // L1Taus 
00195     edm::Handle<l1extra::L1JetParticleCollection> l1TauHandle;
00196     iEvent.getByLabel(L1extraTauJetSource_,l1TauHandle);
00197     l1extra::L1JetParticleCollection::const_iterator itr;
00198     int iL1Obj=0;
00199     for(itr = l1TauHandle->begin(),iL1Obj=0; itr != l1TauHandle->end(); ++itr,iL1Obj++) {    
00200       if(iL1Obj<1) {
00201         t_L1TauJetPt      ->push_back( itr->pt() );
00202         t_L1TauJetEta     ->push_back( itr->eta() );
00203         t_L1TauJetPhi     ->push_back( itr->phi() );
00204       }
00205       if(debugL1Info_) {
00206         std::cout << "tauJ p/pt  " << itr->momentum() << " " << itr->pt() 
00207                   << "  eta/phi " << itr->eta() << " " << itr->phi()
00208                   << std::endl;
00209       }
00210     }
00211 
00212     // L1 Central Jets
00213     edm::Handle<l1extra::L1JetParticleCollection> l1CenJetHandle;
00214     iEvent.getByLabel(L1extraCenJetSource_,l1CenJetHandle);
00215     for( itr = l1CenJetHandle->begin(),iL1Obj=0;  itr != l1CenJetHandle->end(); ++itr,iL1Obj++ ) {
00216       if(iL1Obj<1) {
00217         t_L1CenJetPt    ->push_back( itr->pt() );
00218         t_L1CenJetEta   ->push_back( itr->eta() );
00219         t_L1CenJetPhi   ->push_back( itr->phi() );
00220       }
00221       if(debugL1Info_) {
00222         std::cout << "cenJ p/pt     " << itr->momentum() << " " << itr->pt() 
00223                   << "  eta/phi " << itr->eta() << " " << itr->phi()
00224                   << std::endl;
00225       }
00226     }
00227   
00228     // L1 Forward Jets
00229     edm::Handle<l1extra::L1JetParticleCollection> l1FwdJetHandle;
00230     iEvent.getByLabel(L1extraFwdJetSource_,l1FwdJetHandle);
00231     for( itr = l1FwdJetHandle->begin(),iL1Obj=0;  itr != l1FwdJetHandle->end(); ++itr,iL1Obj++ ) {
00232       if(iL1Obj<1) {
00233         t_L1FwdJetPt    ->push_back( itr->pt() );
00234         t_L1FwdJetEta   ->push_back( itr->eta() );
00235         t_L1FwdJetPhi   ->push_back( itr->phi() );
00236       }
00237       if(debugL1Info_) {
00238         std::cout << "fwdJ p/pt     " << itr->momentum() << " " << itr->pt() 
00239                   << "  eta/phi " << itr->eta() << " " << itr->phi()
00240                   << std::endl;
00241       }
00242     }
00243 
00244     // L1 Isolated EM onjects
00245     l1extra::L1EmParticleCollection::const_iterator itrEm;
00246     edm::Handle<l1extra::L1EmParticleCollection> l1IsoEmHandle ;
00247     iEvent.getByLabel(L1extraIsoEmSource_, l1IsoEmHandle);
00248     for( itrEm = l1IsoEmHandle->begin(),iL1Obj=0;  itrEm != l1IsoEmHandle->end(); ++itrEm,iL1Obj++ ) {
00249       if(iL1Obj<1) {
00250         t_L1IsoEMPt     ->push_back(  itrEm->pt() );
00251         t_L1IsoEMEta    ->push_back(  itrEm->eta() );
00252         t_L1IsoEMPhi    ->push_back(  itrEm->phi() );
00253       }
00254       if(debugL1Info_) {
00255         std::cout << "isoEm p/pt    " << itrEm->momentum() << " " << itrEm->pt()
00256                   << "  eta/phi " << itrEm->eta() << " " << itrEm->phi()
00257                   << std::endl;
00258       }
00259     }
00260   
00261     // L1 Non-Isolated EM onjects
00262     edm::Handle<l1extra::L1EmParticleCollection> l1NonIsoEmHandle ;
00263     iEvent.getByLabel(L1extraNonIsoEmSource_, l1NonIsoEmHandle);
00264     for( itrEm = l1NonIsoEmHandle->begin(),iL1Obj=0;  itrEm != l1NonIsoEmHandle->end(); ++itrEm,iL1Obj++ ) {
00265       if(iL1Obj<1) {
00266         t_L1NonIsoEMPt  ->push_back( itrEm->pt() );
00267         t_L1NonIsoEMEta ->push_back( itrEm->eta() );
00268         t_L1NonIsoEMPhi ->push_back( itrEm->phi() );
00269       }
00270       if(debugL1Info_) {
00271         std::cout << "nonIsoEm p/pt " << itrEm->momentum() << " " << itrEm->pt()
00272                   << "  eta/phi " << itrEm->eta() << " " << itrEm->phi()
00273                   << std::endl;
00274       }
00275     }
00276   
00277     // L1 Muons
00278     l1extra::L1MuonParticleCollection::const_iterator itrMu;
00279     edm::Handle<l1extra::L1MuonParticleCollection> l1MuHandle ;
00280     iEvent.getByLabel(L1extraMuonSource_, l1MuHandle);
00281     for( itrMu = l1MuHandle->begin(),iL1Obj=0;  itrMu != l1MuHandle->end(); ++itrMu,iL1Obj++ ) {
00282       if(iL1Obj<1) {
00283         t_L1MuonPt      ->push_back( itrMu->pt() );
00284         t_L1MuonEta     ->push_back( itrMu->eta() );
00285         t_L1MuonPhi     ->push_back( itrMu->phi() );
00286       }
00287       if(debugL1Info_) {
00288         std::cout << "l1muon p/pt   " << itrMu->momentum() << " " << itrMu->pt()
00289                   << "  eta/phi " << itrMu->eta() << " " << itrMu->phi()
00290                   << std::endl;
00291       }
00292     }
00293   }
00294 
00295   //============== store the information about all the Non-Fake vertices ===============
00296 
00297   edm::Handle<reco::VertexCollection> recVtxs;
00298   iEvent.getByLabel("offlinePrimaryVertices",recVtxs);
00299   
00300   std::vector<reco::Track> svTracks;
00301   math::XYZPoint leadPV(0,0,0);
00302   double sumPtMax = -1.0;
00303   for(unsigned int ind=0;ind<recVtxs->size();ind++) {
00304 
00305     if ( !((*recVtxs)[ind].isFake()) ) {
00306       double vtxTrkSumPt=0.0, vtxTrkSumPtWt=0.0;
00307       int    vtxTrkNWt =0;
00308       double vtxTrkSumPtHP=0.0, vtxTrkSumPtHPWt=0.0;
00309       int    vtxTrkNHP =0, vtxTrkNHPWt =0;
00310 
00311       reco::Vertex::trackRef_iterator vtxTrack = (*recVtxs)[ind].tracks_begin();
00312 
00313       for (vtxTrack = (*recVtxs)[ind].tracks_begin(); vtxTrack!=(*recVtxs)[ind].tracks_end(); vtxTrack++) {
00314 
00315         if((*vtxTrack)->pt()<pvTracksPtMin_) continue;
00316 
00317         bool trkQuality  = (*vtxTrack)->quality(trackQuality_);
00318         
00319         vtxTrkSumPt += (*vtxTrack)->pt();
00320         
00321         vtxTrkSumPt += (*vtxTrack)->pt();
00322         if( trkQuality ) {
00323           vtxTrkSumPtHP += (*vtxTrack)->pt();
00324           vtxTrkNHP++;
00325         }
00326 
00327         double weight = (*recVtxs)[ind].trackWeight(*vtxTrack);
00328         h_PVTracksWt ->Fill(weight);
00329         if( weight>0.5) {
00330           vtxTrkSumPtWt += (*vtxTrack)->pt();
00331           vtxTrkNWt++;
00332           if( trkQuality ) {
00333             vtxTrkSumPtHPWt += (*vtxTrack)->pt();
00334             vtxTrkNHPWt++;
00335           }
00336         }
00337       }
00338 
00339       if(vtxTrkSumPt>sumPtMax) {
00340         sumPtMax = vtxTrkSumPt;
00341         leadPV = math::XYZPoint( (*recVtxs)[ind].x(),(*recVtxs)[ind].y(), (*recVtxs)[ind].z() );
00342       } 
00343 
00344       t_PVx            ->push_back( (*recVtxs)[ind].x() );
00345       t_PVy            ->push_back( (*recVtxs)[ind].y() );
00346       t_PVz            ->push_back( (*recVtxs)[ind].z() );
00347       t_PVisValid      ->push_back( (*recVtxs)[ind].isValid() );
00348       t_PVNTracks      ->push_back( (*recVtxs)[ind].tracksSize() );
00349       t_PVndof         ->push_back( (*recVtxs)[ind].ndof() );
00350       t_PVNTracksWt    ->push_back( vtxTrkNWt );
00351       t_PVTracksSumPt  ->push_back( vtxTrkSumPt );
00352       t_PVTracksSumPtWt->push_back( vtxTrkSumPtWt );
00353 
00354       t_PVNTracksHP      ->push_back( vtxTrkNHP );
00355       t_PVNTracksHPWt    ->push_back( vtxTrkNHPWt );
00356       t_PVTracksSumPtHP  ->push_back( vtxTrkSumPtHP );
00357       t_PVTracksSumPtHPWt->push_back( vtxTrkSumPtHPWt );
00358 
00359       if(myverbose_==4) {
00360         std::cout<<"PV "<<ind<<" isValid "<<(*recVtxs)[ind].isValid()<<" isFake "<<(*recVtxs)[ind].isFake()
00361                  <<" hasRefittedTracks() "<<ind<<" "<<(*recVtxs)[ind].hasRefittedTracks()
00362                  <<" refittedTrksSize "<<(*recVtxs)[ind].refittedTracks().size()
00363                  <<"  tracksSize() "<<(*recVtxs)[ind].tracksSize()<<" sumPt "<<vtxTrkSumPt
00364                  <<std::endl;
00365       }
00366     } // if vtx is not Fake
00367   } // loop over PVs
00368   //===================================================================================
00369 
00370   // Get the beamspot
00371   edm::Handle<reco::BeamSpot> beamSpotH;
00372   iEvent.getByLabel("offlineBeamSpot", beamSpotH);
00373   math::XYZPoint bspot;
00374   bspot = ( beamSpotH.isValid() ) ? beamSpotH->position() : math::XYZPoint(0, 0, 0);
00375 
00376   //=====================================================================
00377 
00378   edm::Handle<reco::CaloJetCollection> jets;
00379   iEvent.getByLabel(JetSrc_,jets);
00380   //  edm::Handle<reco::JetExtendedAssociation::Container> jetExtender;
00381   //  iEvent.getByLabel(JetExtender_,jetExtender);
00382 
00383   for(unsigned int ijet=0;ijet<(*jets).size();ijet++) {
00384     t_jetPt       ->push_back( (*jets)[ijet].pt()     );
00385     t_jetEta      ->push_back( (*jets)[ijet].eta()    );
00386     t_jetPhi      ->push_back( (*jets)[ijet].phi()    );
00387     t_nTrksJetVtx ->push_back( -1.0);
00388     t_nTrksJetCalo->push_back( -1.0);
00389   }
00390 
00391   //===================================================================================
00392   
00393   // get handles to calogeometry and calotopology
00394   edm::ESHandle<CaloGeometry> pG;
00395   iSetup.get<CaloGeometryRecord>().get(pG);
00396   const CaloGeometry* geo = pG.product();
00397   
00398   edm::ESHandle<CaloTopology> theCaloTopology;
00399   iSetup.get<CaloTopologyRecord>().get(theCaloTopology); 
00400   const CaloTopology *caloTopology = theCaloTopology.product();
00401   
00402   edm::ESHandle<HcalTopology> htopo;
00403   iSetup.get<IdealGeometryRecord>().get(htopo);
00404   const HcalTopology* theHBHETopology = htopo.product();
00405   
00406   edm::Handle<EcalRecHitCollection> barrelRecHitsHandle;
00407   edm::Handle<EcalRecHitCollection> endcapRecHitsHandle;
00408   iEvent.getByLabel("ecalRecHit","EcalRecHitsEB",barrelRecHitsHandle);
00409   iEvent.getByLabel("ecalRecHit","EcalRecHitsEE",endcapRecHitsHandle);
00410 
00411   // Retrieve the good/bad ECAL channels from the DB
00412   edm::ESHandle<EcalChannelStatus> ecalChStatus;
00413   iSetup.get<EcalChannelStatusRcd>().get(ecalChStatus);
00414   const EcalChannelStatus* theEcalChStatus = ecalChStatus.product();
00415 
00416   // Retrieve trigger tower map
00417   edm::ESHandle<EcalTrigTowerConstituentsMap> hTtmap;
00418   iSetup.get<IdealGeometryRecord>().get(hTtmap);
00419   const EcalTrigTowerConstituentsMap& ttMap = *hTtmap;
00420 
00421   edm::Handle<HBHERecHitCollection> hbhe;
00422   iEvent.getByLabel(HBHERecHitSource_, hbhe);
00423   if (!hbhe.isValid()) {
00424     nbad++;
00425     if (nbad < 10) std::cout << "No HBHE rechit collection\n";
00426     return;
00427   }
00428   const HBHERecHitCollection Hithbhe = *(hbhe.product());
00429 
00430   //get Handles to SimTracks and SimHits
00431   edm::Handle<edm::SimTrackContainer> SimTk;
00432   if (doMC) iEvent.getByLabel("g4SimHits",SimTk);
00433   edm::SimTrackContainer::const_iterator simTrkItr;
00434 
00435   edm::Handle<edm::SimVertexContainer> SimVtx;
00436   if (doMC) iEvent.getByLabel("g4SimHits",SimVtx);
00437 
00438   //get Handles to PCaloHitContainers of eb/ee/hbhe
00439   edm::Handle<edm::PCaloHitContainer> pcaloeb;
00440   if (doMC) iEvent.getByLabel("g4SimHits", "EcalHitsEB", pcaloeb);
00441 
00442   edm::Handle<edm::PCaloHitContainer> pcaloee;
00443   if (doMC) iEvent.getByLabel("g4SimHits", "EcalHitsEE", pcaloee);
00444 
00445   edm::Handle<edm::PCaloHitContainer> pcalohh;
00446   if (doMC) iEvent.getByLabel("g4SimHits", "HcalHits", pcalohh);
00447   
00448   //associates tracker rechits/simhits to a track
00449   TrackerHitAssociator* associate=0;
00450   if (doMC) associate = new TrackerHitAssociator(iEvent);
00451 
00452   //===================================================================================
00453   
00454   h_nTracks->Fill(trkCollection->size());
00455   
00456   int nTracks        = 0;
00457 
00458   t_nTracks = trkCollection->size();
00459 
00460   // get the list of DetIds closest to the impact point of track on surface calorimeters
00461   std::vector<spr::propagatedTrackID> trkCaloDets;
00462   spr::propagateCALO(trkCollection, geo, bField, theTrackQuality, trkCaloDets, false);
00463   std::vector<spr::propagatedTrackID>::const_iterator trkDetItr;
00464 
00465   if(myverbose_>2) {
00466     for(trkDetItr = trkCaloDets.begin(); trkDetItr != trkCaloDets.end(); trkDetItr++){
00467       std::cout<<trkDetItr->trkItr->p()<<" "<<trkDetItr->trkItr->eta()<<" "<<trkDetItr->okECAL<<" ";
00468       if(trkDetItr->detIdECAL.subdetId() == EcalBarrel) std::cout << (EBDetId)trkDetItr->detIdECAL <<" ";
00469       else                                              std::cout << (EEDetId)trkDetItr->detIdECAL <<" ";
00470       std::cout<<trkDetItr->okHCAL<<" ";
00471       if(trkDetItr->okHCAL) std::cout<<(HcalDetId)trkDetItr->detIdHCAL;
00472       std::cout << std::endl;
00473     }
00474   }
00475   
00476   int nvtxTracks=0;
00477   for(trkDetItr = trkCaloDets.begin(),nTracks=0; trkDetItr != trkCaloDets.end(); trkDetItr++,nTracks++){
00478     
00479     const reco::Track* pTrack = &(*(trkDetItr->trkItr));
00480     
00481     // find vertex index the track is associated with
00482     int pVtxTkId = -1;
00483     for(unsigned int ind=0; ind<recVtxs->size(); ind++) {
00484       if (!((*recVtxs)[ind].isFake())) {
00485         reco::Vertex::trackRef_iterator vtxTrack = (*recVtxs)[ind].tracks_begin();
00486         for (vtxTrack = (*recVtxs)[ind].tracks_begin(); vtxTrack!=(*recVtxs)[ind].tracks_end(); vtxTrack++) {
00487           
00488           const edm::RefToBase<reco::Track> pvtxTrack = (*vtxTrack);
00489           if ( pTrack == pvtxTrack.get() ) { 
00490             pVtxTkId = ind;
00491             break;
00492             if(myverbose_>2) {
00493               if( pTrack->pt()>1.0) {
00494                 std::cout<<"Debug the track association with vertex "<<std::endl;
00495                 std::cout<< pTrack << " "<< pvtxTrack.get() << std::endl;
00496                 std::cout<<" trkVtxIndex "<<nvtxTracks<<" vtx "<<ind<<" pt "<<pTrack->pt()
00497                          <<" eta "<<pTrack->eta()<<" "<<pTrack->pt()-pvtxTrack->pt()
00498                          <<" "<< pTrack->eta()-pvtxTrack->eta()
00499                          << std::endl;
00500                 nvtxTracks++;
00501               }
00502             }
00503           }
00504         }
00505       }
00506     }
00507     
00508     const reco::HitPattern& hitp    = pTrack->hitPattern();
00509     const reco::HitPattern& hitpIn  = pTrack->trackerExpectedHitsInner();
00510     const reco::HitPattern& hitpOut = pTrack->trackerExpectedHitsOuter();
00511 
00512     int nLayersCrossed = hitp.trackerLayersWithMeasurement() ;        
00513     int nOuterHits     = hitp.stripTOBLayersWithMeasurement()+hitp.stripTECLayersWithMeasurement() ;
00514     
00515     bool   ifGood      = pTrack->quality(trackQuality_);
00516     double pt1         = pTrack->pt();
00517     double p1          = pTrack->p();
00518     double eta1        = pTrack->momentum().eta();
00519     double phi1        = pTrack->momentum().phi();
00520     double etaEcal1    = trkDetItr->etaECAL;
00521     double phiEcal1    = trkDetItr->phiECAL;
00522     double etaHcal1    = trkDetItr->etaHCAL;
00523     double phiHcal1    = trkDetItr->phiHCAL;
00524     double dxy1        = pTrack->dxy();
00525     double dz1         = pTrack->dz();
00526     double chisq1      = pTrack->normalizedChi2();
00527     double dxybs1      = beamSpotH.isValid() ? pTrack->dxy(bspot) : pTrack->dxy();      
00528     double dzbs1       = beamSpotH.isValid() ? pTrack->dz(bspot)  : pTrack->dz();       
00529     double dxypv1      = pTrack->dxy(); 
00530     double dzpv1       = pTrack->dz();  
00531     if(pVtxTkId>=0) {
00532       math::XYZPoint thisTkPV = math::XYZPoint( (*recVtxs)[pVtxTkId].x(),(*recVtxs)[pVtxTkId].y(), (*recVtxs)[pVtxTkId].z() );
00533       dxypv1 = pTrack->dxy(thisTkPV);
00534       dzpv1  = pTrack->dz (thisTkPV);
00535     }
00536 
00537     h_recEtaPt_0->Fill(eta1, pt1);
00538     h_recEtaP_0 ->Fill(eta1, p1);
00539     h_recPt_0   ->Fill(pt1);
00540     h_recP_0    ->Fill(p1);
00541     h_recEta_0  ->Fill(eta1);
00542     h_recPhi_0  ->Fill(phi1);
00543     
00544     if(ifGood && nLayersCrossed>7 ) {
00545       h_recEtaPt_1->Fill(eta1, pt1);
00546       h_recEtaP_1 ->Fill(eta1, p1);
00547       h_recPt_1   ->Fill(pt1);
00548       h_recP_1    ->Fill(p1);
00549       h_recEta_1  ->Fill(eta1);
00550       h_recPhi_1  ->Fill(phi1);
00551     }
00552     
00553     if( ! ifGood ) continue;
00554 
00555     if( writeAllTracks && p1>2.0 && nLayersCrossed>7) {
00556       t_trackPAll             ->push_back( p1 );
00557       t_trackEtaAll           ->push_back( eta1 );
00558       t_trackPhiAll           ->push_back( phi1 );
00559       t_trackPtAll            ->push_back( pt1 );
00560       t_trackDxyAll           ->push_back( dxy1 );      
00561       t_trackDzAll            ->push_back( dz1 );       
00562       t_trackDxyPVAll         ->push_back( dxypv1 );    
00563       t_trackDzPVAll          ->push_back( dzpv1 );     
00564       t_trackChiSqAll         ->push_back( chisq1 );    
00565     }
00566     if (doMC) {
00567       edm::SimTrackContainer::const_iterator matchedSimTrkAll = spr::matchedSimTrack(iEvent, SimTk, SimVtx, pTrack, *associate, false); 
00568       if( writeAllTracks && matchedSimTrkAll != SimTk->end())     t_trackPdgIdAll->push_back( matchedSimTrkAll->type() );
00569     }
00570     
00571     if( pt1>minTrackP_ && std::abs(eta1)<maxTrackEta_ && trkDetItr->okECAL) { 
00572       
00573       double maxNearP31x31=999.0, maxNearP25x25=999.0, maxNearP21x21=999.0, maxNearP15x15=999.0;
00574       maxNearP31x31 = spr::chargeIsolationEcal(nTracks, trkCaloDets, geo, caloTopology, 15,15);
00575       maxNearP25x25 = spr::chargeIsolationEcal(nTracks, trkCaloDets, geo, caloTopology, 12,12);
00576       maxNearP21x21 = spr::chargeIsolationEcal(nTracks, trkCaloDets, geo, caloTopology, 10,10);
00577       maxNearP15x15 = spr::chargeIsolationEcal(nTracks, trkCaloDets, geo, caloTopology,  7, 7);
00578 
00579       int iTrkEtaBin=-1, iTrkMomBin=-1;
00580       for(unsigned int ieta=0; ieta<NEtaBins; ieta++) {
00581         if(std::abs(eta1)>genPartEtaBins[ieta] && std::abs(eta1)<genPartEtaBins[ieta+1] ) iTrkEtaBin = ieta;
00582       }
00583       for(unsigned int ipt=0;  ipt<NPBins;   ipt++) {
00584         if( p1>genPartPBins[ipt] &&  p1<genPartPBins[ipt+1] )  iTrkMomBin = ipt;
00585       }
00586       if( iTrkMomBin>=0 && iTrkEtaBin>=0 ) {
00587         h_maxNearP31x31[iTrkMomBin][iTrkEtaBin]->Fill( maxNearP31x31 );
00588         h_maxNearP25x25[iTrkMomBin][iTrkEtaBin]->Fill( maxNearP25x25 );
00589         h_maxNearP21x21[iTrkMomBin][iTrkEtaBin]->Fill( maxNearP21x21 );
00590         h_maxNearP15x15[iTrkMomBin][iTrkEtaBin]->Fill( maxNearP15x15 );
00591       }
00592       if( maxNearP31x31<0.0 && nLayersCrossed>7 && nOuterHits>4) {
00593         h_recEtaPt_2->Fill(eta1, pt1);
00594         h_recEtaP_2 ->Fill(eta1, p1);
00595         h_recPt_2   ->Fill(pt1);
00596         h_recP_2    ->Fill(p1);
00597         h_recEta_2  ->Fill(eta1);
00598         h_recPhi_2  ->Fill(phi1);
00599       }
00600       
00601       // if charge isolated in ECAL, store the further quantities
00602       if( maxNearP31x31<1.0) {
00603 
00604         haveIsoTrack = true;
00605         
00606         // get the matching simTrack
00607         double simTrackP = -1;
00608         if (doMC) {
00609           edm::SimTrackContainer::const_iterator matchedSimTrk = spr::matchedSimTrack(iEvent, SimTk, SimVtx, pTrack, *associate, false);
00610           if( matchedSimTrk != SimTk->end() )simTrackP = matchedSimTrk->momentum().P();
00611         }
00612 
00613         edm::ESHandle<EcalSeverityLevelAlgo> sevlv;
00614         iSetup.get<EcalSeverityLevelAlgoRcd>().get(sevlv);
00615 
00616         // get ECal Tranverse Profile
00617         std::pair<double, bool>  e7x7P, e9x9P, e11x11P, e15x15P;
00618         std::pair<double, bool>  e7x7_10SigP, e9x9_10SigP, e11x11_10SigP, e15x15_10SigP;
00619         std::pair<double, bool>  e7x7_15SigP, e9x9_15SigP, e11x11_15SigP, e15x15_15SigP;
00620         std::pair<double, bool>  e7x7_20SigP, e9x9_20SigP, e11x11_20SigP, e15x15_20SigP;
00621         std::pair<double, bool>  e7x7_25SigP, e9x9_25SigP, e11x11_25SigP, e15x15_25SigP;
00622         std::pair<double, bool>  e7x7_30SigP, e9x9_30SigP, e11x11_30SigP, e15x15_30SigP;
00623 
00624         spr::caloSimInfo simInfo3x3,   simInfo5x5,   simInfo7x7,   simInfo9x9;
00625         spr::caloSimInfo simInfo11x11, simInfo13x13, simInfo15x15, simInfo21x21, simInfo25x25, simInfo31x31;
00626         double trkEcalEne=0;
00627 
00628         const DetId isoCell = trkDetItr->detIdECAL;
00629         e7x7P         = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),3,3,   -100.0, -100.0, tMinE_,tMaxE_);
00630         e9x9P         = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),4,4,   -100.0, -100.0, tMinE_,tMaxE_);
00631         e11x11P       = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),5,5,   -100.0, -100.0, tMinE_,tMaxE_);
00632         e15x15P       = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),7,7,   -100.0, -100.0, tMinE_,tMaxE_);
00633         
00634         e7x7_10SigP   = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),3,3,   0.030,  0.150, tMinE_,tMaxE_);
00635         e9x9_10SigP   = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),4,4,   0.030,  0.150, tMinE_,tMaxE_);
00636         e11x11_10SigP = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),5,5,   0.030,  0.150, tMinE_,tMaxE_);
00637         e15x15_10SigP = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),7,7,   0.030,  0.150, tMinE_,tMaxE_);
00638 
00639         e7x7_15SigP   = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology, sevlv.product(),ttMap, 3,3, 0.20,0.45, tMinE_,tMaxE_);
00640         e9x9_15SigP   = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology, sevlv.product(),ttMap, 4,4, 0.20,0.45, tMinE_,tMaxE_);
00641         e11x11_15SigP = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(), ttMap, 5,5, 0.20,0.45, tMinE_,tMaxE_);
00642         e15x15_15SigP = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology, sevlv.product(),ttMap, 7,7, 0.20,0.45, tMinE_,tMaxE_, false);
00643         
00644         e7x7_20SigP   = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),3,3,   0.060,  0.300, tMinE_,tMaxE_);
00645         e9x9_20SigP   = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),4,4,   0.060,  0.300, tMinE_,tMaxE_);
00646         e11x11_20SigP = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),5,5,   0.060,  0.300, tMinE_,tMaxE_);
00647         e15x15_20SigP = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),7,7,   0.060,  0.300, tMinE_,tMaxE_);
00648         
00649         e7x7_25SigP   = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),3,3,   0.075,  0.375, tMinE_,tMaxE_);
00650         e9x9_25SigP   = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),4,4,   0.075,  0.375, tMinE_,tMaxE_);
00651         e11x11_25SigP = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),5,5,   0.075,  0.375, tMinE_,tMaxE_);
00652         e15x15_25SigP = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),7,7,   0.075,  0.375, tMinE_,tMaxE_);
00653         
00654         e7x7_30SigP   = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),3,3,   0.090,  0.450, tMinE_,tMaxE_);
00655         e9x9_30SigP   = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),4,4,   0.090,  0.450, tMinE_,tMaxE_);
00656         e11x11_30SigP = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),5,5,   0.090,  0.450, tMinE_,tMaxE_);
00657         e15x15_30SigP = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),7,7,   0.090,  0.450, tMinE_,tMaxE_);
00658         if(myverbose_ == 2) {
00659           std::cout << "clean  ecal rechit " << std::endl;
00660           std::cout<<"e7x7 "<<e7x7P.first<<" e9x9 "<<e9x9P.first<<" e11x11 " << e11x11P.first << " e15x15 "<<e15x15P.first<<std::endl;
00661           std::cout<<"e7x7_10Sig "<<e7x7_10SigP.first<<" e11x11_10Sig "<<e11x11_10SigP.first<<" e15x15_10Sig "<<e15x15_10SigP.first<<std::endl;
00662         }
00663 
00664         if (doMC) {
00665           // check the energy from SimHits
00666           spr::eECALSimInfo(iEvent,isoCell,geo,caloTopology,pcaloeb,pcaloee,SimTk,SimVtx,pTrack, *associate, 1,1, simInfo3x3);
00667           spr::eECALSimInfo(iEvent,isoCell,geo,caloTopology,pcaloeb,pcaloee,SimTk,SimVtx,pTrack, *associate, 2,2, simInfo5x5);
00668           spr::eECALSimInfo(iEvent,isoCell,geo,caloTopology,pcaloeb,pcaloee,SimTk,SimVtx,pTrack, *associate, 3,3, simInfo7x7);
00669           spr::eECALSimInfo(iEvent,isoCell,geo,caloTopology,pcaloeb,pcaloee,SimTk,SimVtx,pTrack, *associate, 4,4, simInfo9x9);
00670           spr::eECALSimInfo(iEvent,isoCell,geo,caloTopology,pcaloeb,pcaloee,SimTk,SimVtx,pTrack, *associate, 5,5, simInfo11x11);
00671           spr::eECALSimInfo(iEvent,isoCell,geo,caloTopology,pcaloeb,pcaloee,SimTk,SimVtx,pTrack, *associate, 6,6, simInfo13x13);
00672           spr::eECALSimInfo(iEvent,isoCell,geo,caloTopology,pcaloeb,pcaloee,SimTk,SimVtx,pTrack, *associate, 7,7, simInfo15x15, 150.0,false);
00673           spr::eECALSimInfo(iEvent,isoCell,geo,caloTopology,pcaloeb,pcaloee,SimTk,SimVtx,pTrack, *associate, 10,10, simInfo21x21);
00674           spr::eECALSimInfo(iEvent,isoCell,geo,caloTopology,pcaloeb,pcaloee,SimTk,SimVtx,pTrack, *associate, 12,12, simInfo25x25);
00675           spr::eECALSimInfo(iEvent,isoCell,geo,caloTopology,pcaloeb,pcaloee,SimTk,SimVtx,pTrack, *associate, 15,15, simInfo31x31);
00676   
00677           trkEcalEne   = spr::eCaloSimInfo(iEvent, geo, pcaloeb,pcaloee, SimTk, SimVtx, pTrack, *associate);
00678            if(myverbose_ == 1) {
00679             std::cout << "Track momentum " << pt1 << std::endl;
00680 
00681             std::cout << "ecal siminfo " << std::endl;
00682             std::cout << "simInfo3x3: " << "eTotal " << simInfo3x3.eTotal << " eMatched " << simInfo3x3.eMatched << " eRest " << simInfo3x3.eRest << " eGamma "<<simInfo3x3.eGamma << " eNeutralHad " << simInfo3x3.eNeutralHad << " eChargedHad " << simInfo3x3.eChargedHad << std::endl;
00683             std::cout << "simInfo5x5: " << "eTotal " << simInfo5x5.eTotal << " eMatched " << simInfo5x5.eMatched << " eRest " << simInfo5x5.eRest << " eGamma "<<simInfo5x5.eGamma << " eNeutralHad " << simInfo5x5.eNeutralHad << " eChargedHad " << simInfo5x5.eChargedHad << std::endl;
00684             std::cout << "simInfo7x7: " << "eTotal " << simInfo7x7.eTotal << " eMatched " << simInfo7x7.eMatched << " eRest " << simInfo7x7.eRest << " eGamma "<<simInfo7x7.eGamma << " eNeutralHad " << simInfo7x7.eNeutralHad << " eChargedHad " << simInfo7x7.eChargedHad << std::endl;
00685             std::cout << "simInfo9x9: " << "eTotal " << simInfo9x9.eTotal << " eMatched " << simInfo9x9.eMatched << " eRest " << simInfo9x9.eRest << " eGamma "<<simInfo9x9.eGamma << " eNeutralHad " << simInfo9x9.eNeutralHad << " eChargedHad " << simInfo9x9.eChargedHad << std::endl;
00686             std::cout << "simInfo11x11: " << "eTotal " << simInfo11x11.eTotal << " eMatched " << simInfo11x11.eMatched << " eRest " << simInfo11x11.eRest << " eGamma "<<simInfo11x11.eGamma << " eNeutralHad " << simInfo11x11.eNeutralHad << " eChargedHad " << simInfo11x11.eChargedHad << std::endl;
00687             std::cout << "simInfo15x15: " << "eTotal " << simInfo15x15.eTotal << " eMatched " << simInfo15x15.eMatched << " eRest " << simInfo15x15.eRest << " eGamma "<<simInfo15x15.eGamma << " eNeutralHad " << simInfo15x15.eNeutralHad << " eChargedHad " << simInfo15x15.eChargedHad << std::endl;
00688             std::cout << "simInfo31x31: " << "eTotal " << simInfo31x31.eTotal << " eMatched " << simInfo31x31.eMatched << " eRest " << simInfo31x31.eRest << " eGamma "<<simInfo31x31.eGamma << " eNeutralHad " << simInfo31x31.eNeutralHad << " eChargedHad " << simInfo31x31.eChargedHad << std::endl;
00689             std::cout << "trkEcalEne" << trkEcalEne << std::endl;
00690            }
00691         }
00692 
00693         // =======  Get HCAL information 
00694         double hcalScale=1.0;
00695         if( std::abs(pTrack->eta())<1.4 ) {
00696           hcalScale=120.0;
00697         } else {
00698           hcalScale=135.0;
00699         }
00700         
00701         double maxNearHcalP3x3=-1, maxNearHcalP5x5=-1, maxNearHcalP7x7=-1;
00702         maxNearHcalP3x3  = spr::chargeIsolationHcal(nTracks, trkCaloDets, theHBHETopology, 1,1);
00703         maxNearHcalP5x5  = spr::chargeIsolationHcal(nTracks, trkCaloDets, theHBHETopology, 2,2);
00704         maxNearHcalP7x7  = spr::chargeIsolationHcal(nTracks, trkCaloDets, theHBHETopology, 3,3);
00705 
00706         double h3x3=0,    h5x5=0,    h7x7=0;
00707         double h3x3Sig=0, h5x5Sig=0, h7x7Sig=0;
00708         double trkHcalEne = 0;
00709         spr::caloSimInfo hsimInfo3x3, hsimInfo5x5, hsimInfo7x7;
00710         
00711         if(trkDetItr->okHCAL) {
00712           const DetId ClosestCell(trkDetItr->detIdHCAL);
00713           // bool includeHO=false, bool algoNew=true, bool debug=false
00714           h3x3 = spr::eHCALmatrix(theHBHETopology, ClosestCell, hbhe,1,1, false, true, -100.0, -100.0, -100.0, -100.0, tMinH_,tMaxH_);  
00715           h5x5 = spr::eHCALmatrix(theHBHETopology, ClosestCell, hbhe,2,2, false, true, -100.0, -100.0, -100.0, -100.0, tMinH_,tMaxH_);  
00716           h7x7 = spr::eHCALmatrix(theHBHETopology, ClosestCell, hbhe,3,3, false, true, -100.0, -100.0, -100.0, -100.0, tMinH_,tMaxH_);  
00717           h3x3Sig = spr::eHCALmatrix(theHBHETopology, ClosestCell, hbhe,1,1, false, true, 0.7, 0.8, -100.0, -100.0, tMinH_,tMaxH_);  
00718           h5x5Sig = spr::eHCALmatrix(theHBHETopology, ClosestCell, hbhe,2,2, false, true, 0.7, 0.8, -100.0, -100.0, tMinH_,tMaxH_);  
00719           h7x7Sig = spr::eHCALmatrix(theHBHETopology, ClosestCell, hbhe,3,3, false, true, 0.7, 0.8, -100.0, -100.0, tMinH_,tMaxH_);  
00720 
00721           if(myverbose_==2) {
00722             std::cout << "HCAL 3x3 " << h3x3 << " " << h3x3Sig << " 5x5 " <<  h5x5 << " " << h5x5Sig << " 7x7 " << h7x7 << " " << h7x7Sig << std::endl;
00723           }
00724           
00725           if (doMC) {
00726             spr::eHCALSimInfo(iEvent, theHBHETopology, ClosestCell, geo,pcalohh, SimTk, SimVtx, pTrack, *associate, 1,1, hsimInfo3x3);
00727             spr::eHCALSimInfo(iEvent, theHBHETopology, ClosestCell, geo,pcalohh, SimTk, SimVtx, pTrack, *associate, 2,2, hsimInfo5x5);
00728             spr::eHCALSimInfo(iEvent, theHBHETopology, ClosestCell, geo,pcalohh, SimTk, SimVtx, pTrack, *associate, 3,3, hsimInfo7x7, 150.0, false,false);
00729             trkHcalEne  = spr::eCaloSimInfo(iEvent, geo,pcalohh, SimTk, SimVtx, pTrack, *associate);
00730             if(myverbose_ == 1) {
00731               std::cout << "Hcal siminfo " << std::endl;
00732               std::cout << "hsimInfo3x3: " << "eTotal " << hsimInfo3x3.eTotal << " eMatched " << hsimInfo3x3.eMatched << " eRest " << hsimInfo3x3.eRest << " eGamma "<<hsimInfo3x3.eGamma << " eNeutralHad " << hsimInfo3x3.eNeutralHad << " eChargedHad " << hsimInfo3x3.eChargedHad << std::endl;
00733               std::cout << "hsimInfo5x5: " << "eTotal " << hsimInfo5x5.eTotal << " eMatched " << hsimInfo5x5.eMatched << " eRest " << hsimInfo5x5.eRest << " eGamma "<<hsimInfo5x5.eGamma << " eNeutralHad " << hsimInfo5x5.eNeutralHad << " eChargedHad " << hsimInfo5x5.eChargedHad << std::endl;
00734               std::cout << "hsimInfo7x7: " << "eTotal " << hsimInfo7x7.eTotal << " eMatched " << hsimInfo7x7.eMatched << " eRest " << hsimInfo7x7.eRest << " eGamma "<<hsimInfo7x7.eGamma << " eNeutralHad " << hsimInfo7x7.eNeutralHad << " eChargedHad " << hsimInfo7x7.eChargedHad << std::endl;
00735               std::cout << "trkHcalEne " << trkHcalEne << std::endl;
00736             }
00737           }
00738 
00739           // debug the ecal and hcal matrix
00740           if(myverbose_==4) {
00741             std::cout<<"Run "<<iEvent.id().run()<<"  Event "<<iEvent.id().event()<<std::endl; 
00742             std::vector<std::pair<DetId,double> > v7x7 = spr::eHCALmatrixCell(theHBHETopology, ClosestCell, hbhe,3,3, false, false);
00743             double sumv=0.0;
00744             
00745             for(unsigned int iv=0; iv<v7x7.size(); iv++) { 
00746               sumv += v7x7[iv].second;
00747             }
00748             std::cout<<"h7x7 "<<h7x7<<" v7x7 "<<sumv << " in " << v7x7.size() <<std::endl;
00749             for(unsigned int iv=0; iv<v7x7.size(); iv++) { 
00750               HcalDetId id = v7x7[iv].first;
00751               std::cout << " Cell " << iv << " 0x" << std::hex << v7x7[iv].first() << std::dec << " " << id << " Energy " << v7x7[iv].second << std::endl;
00752             }
00753           }
00754           
00755         }
00756         if (doMC) {
00757           trkHcalEne  = spr::eCaloSimInfo(iEvent, geo,pcalohh, SimTk, SimVtx, pTrack, *associate);
00758         }
00759 
00760         // ====================================================================================================
00761         // get diff between track outermost hit position and the propagation point at outermost surface of tracker      
00762         std::pair<math::XYZPoint,double> point2_TK0 = spr::propagateTrackerEnd( pTrack, bField, false);
00763         math::XYZPoint diff(pTrack->outerPosition().X()-point2_TK0.first.X(), 
00764                             pTrack->outerPosition().Y()-point2_TK0.first.Y(), 
00765                             pTrack->outerPosition().Z()-point2_TK0.first.Z() );
00766         double trackOutPosOutHitDr = diff.R();
00767         double trackL              = point2_TK0.second;
00768         if (myverbose_==5) {
00769           std::cout<<" propagted "<<point2_TK0.first<<" "<< point2_TK0.first.eta()<<" "<<point2_TK0.first.phi()<<std::endl;
00770           std::cout<<" outerPosition() "<< pTrack->outerPosition() << " "<< pTrack->outerPosition().eta()<< " " << pTrack->outerPosition().phi()<< std::endl;
00771           std::cout<<"diff " << diff << " diffR " <<diff.R()<<" diffR/L "<<diff.R()/point2_TK0.second <<std::endl;
00772         }
00773 
00774         for(unsigned int ind=0;ind<recVtxs->size();ind++) {
00775           if (!((*recVtxs)[ind].isFake())) {
00776             reco::Vertex::trackRef_iterator vtxTrack = (*recVtxs)[ind].tracks_begin();
00777             if( DeltaR(eta1,phi1, (*vtxTrack)->eta(),(*vtxTrack)->phi()) < 0.01 ) t_trackPVIdx ->push_back( ind );
00778             else                                                                  t_trackPVIdx ->push_back( -1 );
00779           }
00780         }
00781 
00782         // Fill the tree Branches here 
00783         t_trackPVIdx            ->push_back( pVtxTkId );
00784         t_trackP                ->push_back( p1 );
00785         t_trackPt               ->push_back( pt1 );
00786         t_trackEta              ->push_back( eta1 );
00787         t_trackPhi              ->push_back( phi1 );
00788         t_trackEcalEta          ->push_back( etaEcal1 );
00789         t_trackEcalPhi          ->push_back( phiEcal1 );
00790         t_trackHcalEta          ->push_back( etaHcal1 );
00791         t_trackHcalPhi          ->push_back( phiHcal1 );
00792         t_trackDxy              ->push_back( dxy1 );    
00793         t_trackDz               ->push_back( dz1 );     
00794         t_trackDxyBS            ->push_back( dxybs1 );  
00795         t_trackDzBS             ->push_back( dzbs1 );   
00796         t_trackDxyPV            ->push_back( dxypv1 );  
00797         t_trackDzPV             ->push_back( dzpv1 );   
00798         t_trackChiSq            ->push_back( chisq1 );  
00799         t_trackNOuterHits       ->push_back( nOuterHits );
00800         t_NLayersCrossed        ->push_back( nLayersCrossed );
00801         
00802         t_trackHitsTOB          ->push_back( hitp.stripTOBLayersWithMeasurement()       ); 
00803         t_trackHitsTEC          ->push_back( hitp.stripTECLayersWithMeasurement()       );
00804         t_trackHitInMissTOB     ->push_back( hitpIn.stripTOBLayersWithoutMeasurement()  ); 
00805         t_trackHitInMissTEC     ->push_back( hitpIn.stripTECLayersWithoutMeasurement()  );  
00806         t_trackHitInMissTIB     ->push_back( hitpIn.stripTIBLayersWithoutMeasurement()  );  
00807         t_trackHitInMissTID     ->push_back( hitpIn.stripTIDLayersWithoutMeasurement()  );
00808         t_trackHitInMissTIBTID  ->push_back( hitpIn.stripTIBLayersWithoutMeasurement() + hitpIn.stripTIDLayersWithoutMeasurement() );  
00809 
00810         t_trackHitOutMissTOB    ->push_back( hitpOut.stripTOBLayersWithoutMeasurement() );
00811         t_trackHitOutMissTEC    ->push_back( hitpOut.stripTECLayersWithoutMeasurement() ); 
00812         t_trackHitOutMissTIB    ->push_back( hitpOut.stripTIBLayersWithoutMeasurement() );
00813         t_trackHitOutMissTID    ->push_back( hitpOut.stripTIDLayersWithoutMeasurement() );
00814         t_trackHitOutMissTOBTEC ->push_back( hitpOut.stripTOBLayersWithoutMeasurement() + hitpOut.stripTECLayersWithoutMeasurement() );
00815 
00816         t_trackHitInMeasTOB     ->push_back( hitpIn.stripTOBLayersWithMeasurement()  ); 
00817         t_trackHitInMeasTEC     ->push_back( hitpIn.stripTECLayersWithMeasurement()  );  
00818         t_trackHitInMeasTIB     ->push_back( hitpIn.stripTIBLayersWithMeasurement()  );  
00819         t_trackHitInMeasTID     ->push_back( hitpIn.stripTIDLayersWithMeasurement()  );
00820         t_trackHitOutMeasTOB    ->push_back( hitpOut.stripTOBLayersWithMeasurement() );
00821         t_trackHitOutMeasTEC    ->push_back( hitpOut.stripTECLayersWithMeasurement() ); 
00822         t_trackHitOutMeasTIB    ->push_back( hitpOut.stripTIBLayersWithMeasurement() );
00823         t_trackHitOutMeasTID    ->push_back( hitpOut.stripTIDLayersWithMeasurement() );
00824         t_trackOutPosOutHitDr   ->push_back( trackOutPosOutHitDr                     );
00825         t_trackL                ->push_back( trackL                                  );
00826 
00827         t_maxNearP31x31         ->push_back( maxNearP31x31 );
00828         t_maxNearP21x21         ->push_back( maxNearP21x21 );
00829         
00830         t_ecalSpike11x11        ->push_back( e11x11P.second );
00831         t_e7x7                  ->push_back( e7x7P.first );
00832         t_e9x9                  ->push_back( e9x9P.first );
00833         t_e11x11                ->push_back( e11x11P.first );
00834         t_e15x15                ->push_back( e15x15P.first );
00835 
00836         t_e7x7_10Sig            ->push_back( e7x7_10SigP.first   ); 
00837         t_e9x9_10Sig            ->push_back( e9x9_10SigP.first   ); 
00838         t_e11x11_10Sig          ->push_back( e11x11_10SigP.first ); 
00839         t_e15x15_10Sig          ->push_back( e15x15_10SigP.first );
00840         t_e7x7_15Sig            ->push_back( e7x7_15SigP.first   ); 
00841         t_e9x9_15Sig            ->push_back( e9x9_15SigP.first   ); 
00842         t_e11x11_15Sig          ->push_back( e11x11_15SigP.first ); 
00843         t_e15x15_15Sig          ->push_back( e15x15_15SigP.first );
00844         t_e7x7_20Sig            ->push_back( e7x7_20SigP.first   ); 
00845         t_e9x9_20Sig            ->push_back( e9x9_20SigP.first   ); 
00846         t_e11x11_20Sig          ->push_back( e11x11_20SigP.first ); 
00847         t_e15x15_20Sig          ->push_back( e15x15_20SigP.first );
00848         t_e7x7_25Sig            ->push_back( e7x7_25SigP.first   ); 
00849         t_e9x9_25Sig            ->push_back( e9x9_25SigP.first   ); 
00850         t_e11x11_25Sig          ->push_back( e11x11_25SigP.first ); 
00851         t_e15x15_25Sig          ->push_back( e15x15_25SigP.first );
00852         t_e7x7_30Sig            ->push_back( e7x7_30SigP.first   ); 
00853         t_e9x9_30Sig            ->push_back( e9x9_30SigP.first   ); 
00854         t_e11x11_30Sig          ->push_back( e11x11_30SigP.first ); 
00855         t_e15x15_30Sig          ->push_back( e15x15_30SigP.first );
00856 
00857         if (doMC) {
00858           t_esim7x7               ->push_back( simInfo7x7.eTotal );
00859           t_esim9x9               ->push_back( simInfo9x9.eTotal );
00860           t_esim11x11             ->push_back( simInfo11x11.eTotal );
00861           t_esim15x15             ->push_back( simInfo15x15.eTotal );
00862         
00863           t_esim7x7Matched        ->push_back( simInfo7x7.eMatched );
00864           t_esim9x9Matched        ->push_back( simInfo9x9.eMatched );
00865           t_esim11x11Matched      ->push_back( simInfo11x11.eMatched );
00866           t_esim15x15Matched      ->push_back( simInfo15x15.eMatched );
00867         
00868           t_esim7x7Rest           ->push_back( simInfo7x7.eRest );
00869           t_esim9x9Rest           ->push_back( simInfo9x9.eRest );
00870           t_esim11x11Rest         ->push_back( simInfo11x11.eRest );
00871           t_esim15x15Rest         ->push_back( simInfo15x15.eRest );
00872         
00873           t_esim7x7Photon         ->push_back( simInfo7x7.eGamma );
00874           t_esim9x9Photon         ->push_back( simInfo9x9.eGamma );
00875           t_esim11x11Photon       ->push_back( simInfo11x11.eGamma );
00876           t_esim15x15Photon       ->push_back( simInfo15x15.eGamma );
00877         
00878           t_esim7x7NeutHad        ->push_back( simInfo7x7.eNeutralHad );
00879           t_esim9x9NeutHad        ->push_back( simInfo9x9.eNeutralHad );
00880           t_esim11x11NeutHad      ->push_back( simInfo11x11.eNeutralHad );
00881           t_esim15x15NeutHad      ->push_back( simInfo15x15.eNeutralHad );
00882         
00883           t_esim7x7CharHad        ->push_back( simInfo7x7.eChargedHad );
00884           t_esim9x9CharHad        ->push_back( simInfo9x9.eChargedHad );
00885           t_esim11x11CharHad      ->push_back( simInfo11x11.eChargedHad );
00886           t_esim15x15CharHad      ->push_back( simInfo15x15.eChargedHad );
00887         
00888           t_trkEcalEne            ->push_back( trkEcalEne );
00889           t_simTrackP             ->push_back( simTrackP );
00890           t_esimPdgId             ->push_back( simInfo11x11.pdgMatched );
00891         }
00892 
00893         t_maxNearHcalP3x3       ->push_back( maxNearHcalP3x3 );
00894         t_maxNearHcalP5x5       ->push_back( maxNearHcalP5x5 );
00895         t_maxNearHcalP7x7       ->push_back( maxNearHcalP7x7 );
00896         
00897         t_h3x3                  ->push_back( h3x3 );
00898         t_h5x5                  ->push_back( h5x5 );
00899         t_h7x7                  ->push_back( h7x7 );
00900         t_h3x3Sig               ->push_back( h3x3Sig );
00901         t_h5x5Sig               ->push_back( h5x5Sig );
00902         t_h7x7Sig               ->push_back( h7x7Sig );
00903         
00904         t_infoHcal              ->push_back( trkDetItr->okHCAL );
00905         if (doMC) {
00906           t_trkHcalEne            ->push_back( hcalScale*trkHcalEne );
00907         
00908           t_hsim3x3               ->push_back( hcalScale*hsimInfo3x3.eTotal );
00909           t_hsim5x5               ->push_back( hcalScale*hsimInfo5x5.eTotal );
00910           t_hsim7x7               ->push_back( hcalScale*hsimInfo7x7.eTotal );
00911         
00912           t_hsim3x3Matched        ->push_back( hcalScale*hsimInfo3x3.eMatched );
00913           t_hsim5x5Matched        ->push_back( hcalScale*hsimInfo5x5.eMatched );
00914           t_hsim7x7Matched        ->push_back( hcalScale*hsimInfo7x7.eMatched );
00915         
00916           t_hsim3x3Rest           ->push_back( hcalScale*hsimInfo3x3.eRest );
00917           t_hsim5x5Rest           ->push_back( hcalScale*hsimInfo5x5.eRest );
00918           t_hsim7x7Rest           ->push_back( hcalScale*hsimInfo7x7.eRest );
00919         
00920           t_hsim3x3Photon         ->push_back( hcalScale*hsimInfo3x3.eGamma );
00921           t_hsim5x5Photon         ->push_back( hcalScale*hsimInfo5x5.eGamma );
00922           t_hsim7x7Photon         ->push_back( hcalScale*hsimInfo7x7.eGamma );
00923         
00924           t_hsim3x3NeutHad        ->push_back( hcalScale*hsimInfo3x3.eNeutralHad );
00925           t_hsim5x5NeutHad        ->push_back( hcalScale*hsimInfo5x5.eNeutralHad );
00926           t_hsim7x7NeutHad        ->push_back( hcalScale*hsimInfo7x7.eNeutralHad );
00927         
00928           t_hsim3x3CharHad        ->push_back( hcalScale*hsimInfo3x3.eChargedHad );
00929           t_hsim5x5CharHad        ->push_back( hcalScale*hsimInfo5x5.eChargedHad );
00930           t_hsim7x7CharHad        ->push_back( hcalScale*hsimInfo7x7.eChargedHad );
00931         }
00932         /*
00933         if(hcalScale*hsimInfo3x3["eTotal"] > 50.0) {
00934 
00935           std::cout << "Loosely Iso Track : eta " << eta1 << " Rec Mom " << p1 << " SimMom " << simTrackP << " h3x3 " << h3x3 << std::endl;
00936 
00937           std::cout <<"Closest cell Hcal (atHCAL) " << (HcalDetId)ClosestCell << std::endl;
00938 
00939           std::cout <<"trkHcalEne, etotal, matched, rest " <<hcalScale*trkHcalEne<<std::setw(15)<<hcalScale*hsimInfo3x3["eTotal"]
00940                     <<std::setw(15)<<hcalScale*hsimInfo3x3["eMatched"]<<std::setw(15)<<hcalScale*hsimInfo3x3["eRest"]
00941                     <<std::endl;
00942           unsigned int nn = t_trkHcalEne->size();
00943           std::cout <<"in Tree                           " << (*t_trkHcalEne)[nn-1] <<std::setw(15)<< (*t_hsim3x3)[nn-1]
00944                     <<std::setw(15)<< (*t_hsim3x3Matched)[nn-1] <<std::setw(15)<< (*t_hsim3x3Rest)[nn-1]
00945                     << std::endl;
00946 
00947           std::cout << "debug output \n" << std::endl;
00948           std::map<std::string, double> hsimInfo3x3_debug = spr::eHCALSimInfo(iEvent, theHBHETopology, ClosestCell, geo,pcalohh, SimTk, SimVtx, pTrack, *associate, 1,1, 150.0, true);
00949 
00950         }
00951         */
00952 
00953 
00954       } // if loosely isolated track
00955     } // check p1/eta
00956   } // loop over track collection
00957   
00958   if(haveIsoTrack) tree->Fill();
00959 }
00960 
00961 // ----- method called once each job just before starting event loop ----
00962 void IsolatedTracksNxN::beginJob() {
00963 
00964   nEventProc=0;
00965 
00966   //  double tempgen_TH[21] = { 1.0,  2.0,  3.0,  4.0,  5.0, 
00967  double tempgen_TH[16] = { 0.0,  1.0,  2.0,  3.0,  4.0,  
00968                            5.0,  6.0,  7.0,  9.0, 11.0, 
00969                           15.0, 20.0, 30.0, 50.0, 75.0, 100.0};
00970 
00971   for(int i=0; i<16; i++)  genPartPBins[i]  = tempgen_TH[i];
00972 
00973   double tempgen_Eta[4] = {0.0, 1.131, 1.653, 2.172};
00974 
00975   for(int i=0; i<4; i++) genPartEtaBins[i] = tempgen_Eta[i];
00976 
00977   BookHistograms();
00978 }
00979 
00980 // ----- method called once each job just after ending the event loop ----
00981 void IsolatedTracksNxN::endJob() {
00982 
00983   if(L1TriggerAlgoInfo_) {
00984     std::map< std::pair<unsigned int,std::string>, int>::iterator itr;
00985     for(itr=l1AlgoMap.begin(); itr!=l1AlgoMap.end(); itr++) {
00986       std::cout << " ****endjob**** " << (itr->first).first <<" "
00987                 <<(itr->first).second <<" "<<itr->second 
00988                 << std::endl;
00989       int ibin = (itr->first).first;
00990       TString name( (itr->first).second );
00991       h_L1AlgoNames->GetXaxis()->SetBinLabel(ibin+1,name);
00992     }
00993     std::cout << "Number of Events Processed " << nEventProc << std::endl;
00994   }
00995 
00996 }
00997 
00998 //========================================================================================================
00999 
01000 void IsolatedTracksNxN::clearTreeVectors() {
01001 
01002   t_PVx               ->clear();
01003   t_PVy               ->clear();
01004   t_PVz               ->clear();
01005   t_PVisValid         ->clear();
01006   t_PVndof            ->clear();
01007   t_PVNTracks         ->clear();
01008   t_PVNTracksWt       ->clear();
01009   t_PVTracksSumPt     ->clear();
01010   t_PVTracksSumPtWt   ->clear();
01011   t_PVNTracksHP       ->clear();
01012   t_PVNTracksHPWt     ->clear();
01013   t_PVTracksSumPtHP   ->clear();
01014   t_PVTracksSumPtHPWt ->clear();
01015 
01016   for(int i=0; i<128; i++) t_L1Decision[i]=0;
01017   t_L1AlgoNames       ->clear();
01018   t_L1PreScale        ->clear();
01019 
01020   t_L1CenJetPt        ->clear();
01021   t_L1CenJetEta       ->clear();    
01022   t_L1CenJetPhi       ->clear();
01023   t_L1FwdJetPt        ->clear();
01024   t_L1FwdJetEta       ->clear();
01025   t_L1FwdJetPhi       ->clear();
01026   t_L1TauJetPt        ->clear();
01027   t_L1TauJetEta       ->clear();     
01028   t_L1TauJetPhi       ->clear();
01029   t_L1MuonPt          ->clear();
01030   t_L1MuonEta         ->clear();     
01031   t_L1MuonPhi         ->clear();
01032   t_L1IsoEMPt         ->clear();
01033   t_L1IsoEMEta        ->clear();
01034   t_L1IsoEMPhi        ->clear();
01035   t_L1NonIsoEMPt      ->clear();
01036   t_L1NonIsoEMEta     ->clear();
01037   t_L1NonIsoEMPhi     ->clear();
01038   t_L1METPt           ->clear();
01039   t_L1METEta          ->clear();
01040   t_L1METPhi          ->clear();
01041 
01042   t_jetPt             ->clear();
01043   t_jetEta            ->clear();
01044   t_jetPhi            ->clear();
01045   t_nTrksJetCalo      ->clear();  
01046   t_nTrksJetVtx       ->clear();
01047 
01048   t_trackPAll         ->clear();
01049   t_trackEtaAll       ->clear();
01050   t_trackPhiAll       ->clear();
01051   t_trackPdgIdAll     ->clear();
01052   t_trackPtAll        ->clear();
01053   t_trackDxyAll       ->clear();
01054   t_trackDzAll        ->clear();
01055   t_trackDxyPVAll     ->clear();
01056   t_trackDzPVAll      ->clear();
01057   t_trackChiSqAll     ->clear();
01058 
01059   t_trackP            ->clear();
01060   t_trackPt           ->clear();
01061   t_trackEta          ->clear();
01062   t_trackPhi          ->clear();
01063   t_trackEcalEta      ->clear();
01064   t_trackEcalPhi      ->clear();
01065   t_trackHcalEta      ->clear();
01066   t_trackHcalPhi      ->clear();
01067   t_NLayersCrossed    ->clear();
01068   t_trackNOuterHits   ->clear();
01069   t_trackDxy          ->clear();
01070   t_trackDxyBS        ->clear();
01071   t_trackDz           ->clear();
01072   t_trackDzBS         ->clear();
01073   t_trackDxyPV        ->clear();
01074   t_trackDzPV         ->clear();
01075   t_trackChiSq        ->clear();
01076   t_trackPVIdx        ->clear();
01077   t_trackHitsTOB          ->clear(); 
01078   t_trackHitsTEC          ->clear();
01079   t_trackHitInMissTOB     ->clear(); 
01080   t_trackHitInMissTEC     ->clear();  
01081   t_trackHitInMissTIB     ->clear();  
01082   t_trackHitInMissTID     ->clear();
01083   t_trackHitInMissTIBTID  ->clear();  
01084   t_trackHitOutMissTOB    ->clear();
01085   t_trackHitOutMissTEC    ->clear(); 
01086   t_trackHitOutMissTIB    ->clear();
01087   t_trackHitOutMissTID    ->clear(); 
01088   t_trackHitOutMissTOBTEC ->clear();
01089   t_trackHitInMeasTOB     ->clear(); 
01090   t_trackHitInMeasTEC     ->clear();  
01091   t_trackHitInMeasTIB     ->clear();  
01092   t_trackHitInMeasTID     ->clear();
01093   t_trackHitOutMeasTOB    ->clear();
01094   t_trackHitOutMeasTEC    ->clear(); 
01095   t_trackHitOutMeasTIB    ->clear();
01096   t_trackHitOutMeasTID    ->clear();
01097   t_trackOutPosOutHitDr   ->clear();
01098   t_trackL                ->clear();
01099 
01100   t_maxNearP31x31     ->clear();
01101   t_maxNearP21x21     ->clear();
01102 
01103   t_ecalSpike11x11    ->clear();
01104   t_e7x7              ->clear();
01105   t_e9x9              ->clear();
01106   t_e11x11            ->clear();
01107   t_e15x15            ->clear();
01108 
01109   t_e7x7_10Sig        ->clear();
01110   t_e9x9_10Sig        ->clear();
01111   t_e11x11_10Sig      ->clear();
01112   t_e15x15_10Sig      ->clear();
01113   t_e7x7_15Sig        ->clear();
01114   t_e9x9_15Sig        ->clear();
01115   t_e11x11_15Sig      ->clear();
01116   t_e15x15_15Sig      ->clear();
01117   t_e7x7_20Sig        ->clear();
01118   t_e9x9_20Sig        ->clear();
01119   t_e11x11_20Sig      ->clear();
01120   t_e15x15_20Sig      ->clear();
01121   t_e7x7_25Sig        ->clear();
01122   t_e9x9_25Sig        ->clear();
01123   t_e11x11_25Sig      ->clear();
01124   t_e15x15_25Sig      ->clear();
01125   t_e7x7_30Sig        ->clear();
01126   t_e9x9_30Sig        ->clear();
01127   t_e11x11_30Sig      ->clear();
01128   t_e15x15_30Sig      ->clear();
01129 
01130   if (doMC) {
01131     t_simTrackP         ->clear();
01132     t_esimPdgId         ->clear();
01133     t_trkEcalEne        ->clear();
01134 
01135     t_esim7x7           ->clear();
01136     t_esim9x9           ->clear();
01137     t_esim11x11         ->clear();
01138     t_esim15x15         ->clear();
01139   
01140     t_esim7x7Matched    ->clear();
01141     t_esim9x9Matched    ->clear();
01142     t_esim11x11Matched  ->clear();
01143     t_esim15x15Matched  ->clear();
01144 
01145     t_esim7x7Rest       ->clear();
01146     t_esim9x9Rest       ->clear();
01147     t_esim11x11Rest     ->clear();
01148     t_esim15x15Rest     ->clear();
01149 
01150     t_esim7x7Photon     ->clear();
01151     t_esim9x9Photon     ->clear();
01152     t_esim11x11Photon   ->clear();
01153     t_esim15x15Photon   ->clear();
01154 
01155     t_esim7x7NeutHad    ->clear();
01156     t_esim9x9NeutHad    ->clear();
01157     t_esim11x11NeutHad  ->clear();
01158     t_esim15x15NeutHad  ->clear();
01159 
01160     t_esim7x7CharHad    ->clear();
01161     t_esim9x9CharHad    ->clear();
01162     t_esim11x11CharHad  ->clear();
01163     t_esim15x15CharHad  ->clear();
01164   }
01165 
01166   t_maxNearHcalP3x3   ->clear();
01167   t_maxNearHcalP5x5   ->clear();
01168   t_maxNearHcalP7x7   ->clear();
01169 
01170   t_h3x3              ->clear();
01171   t_h5x5              ->clear();
01172   t_h7x7              ->clear();
01173   t_h3x3Sig           ->clear();
01174   t_h5x5Sig           ->clear();
01175   t_h7x7Sig           ->clear();
01176 
01177   t_infoHcal          ->clear();
01178 
01179   if (doMC) {
01180     t_trkHcalEne        ->clear();
01181 
01182     t_hsim3x3           ->clear();
01183     t_hsim5x5           ->clear();
01184     t_hsim7x7           ->clear();
01185     t_hsim3x3Matched    ->clear();
01186     t_hsim5x5Matched    ->clear();
01187     t_hsim7x7Matched    ->clear();
01188     t_hsim3x3Rest       ->clear();
01189     t_hsim5x5Rest       ->clear();
01190     t_hsim7x7Rest       ->clear();
01191     t_hsim3x3Photon     ->clear();
01192     t_hsim5x5Photon     ->clear();
01193     t_hsim7x7Photon     ->clear();
01194     t_hsim3x3NeutHad    ->clear();
01195     t_hsim5x5NeutHad    ->clear();
01196     t_hsim7x7NeutHad    ->clear();
01197     t_hsim3x3CharHad    ->clear();
01198     t_hsim5x5CharHad    ->clear();
01199     t_hsim7x7CharHad    ->clear();
01200   }
01201 }
01202 
01203 void IsolatedTracksNxN::BookHistograms(){
01204 
01205   char hname[100], htit[100];
01206 
01207   TFileDirectory dir = fs->mkdir("nearMaxTrackP");
01208 
01209   for(unsigned int ieta=0; ieta<NEtaBins; ieta++) {
01210     double lowEta=-5.0, highEta= 5.0;
01211     lowEta  = genPartEtaBins[ieta];
01212     highEta = genPartEtaBins[ieta+1];
01213 
01214     for(unsigned int ipt=0; ipt<NPBins; ipt++) {
01215       double lowP=0.0, highP=300.0;
01216       lowP    = genPartPBins[ipt];
01217       highP   = genPartPBins[ipt+1];
01218       sprintf(hname, "h_maxNearP31x31_ptBin%i_etaBin%i",ipt, ieta);
01219       sprintf(htit,  "maxNearP in 31x31 (%3.2f<|#eta|<%3.2f), (%2.0f<trkP<%3.0f)", lowEta, highEta, lowP, highP );
01220       h_maxNearP31x31[ipt][ieta] = dir.make<TH1F>(hname, htit, 220, -2.0, 20.0);
01221       h_maxNearP31x31[ipt][ieta] ->Sumw2();
01222       sprintf(hname, "h_maxNearP25x25_ptBin%i_etaBin%i",ipt, ieta);
01223       sprintf(htit,  "maxNearP in 25x25 (%3.2f<|#eta|<%3.2f), (%2.0f<trkP<%3.0f)", lowEta, highEta, lowP, highP );
01224       h_maxNearP25x25[ipt][ieta] = dir.make<TH1F>(hname, htit, 220, -2.0, 20.0);
01225       h_maxNearP25x25[ipt][ieta] ->Sumw2();
01226       sprintf(hname, "h_maxNearP21x21_ptBin%i_etaBin%i",ipt, ieta);
01227       sprintf(htit,  "maxNearP in 21x21 (%3.2f<|#eta|<%3.2f), (%2.0f<trkP<%3.0f)", lowEta, highEta, lowP, highP );
01228       h_maxNearP21x21[ipt][ieta] = dir.make<TH1F>(hname, htit, 220, -2.0, 20.0);
01229       h_maxNearP21x21[ipt][ieta] ->Sumw2();
01230       sprintf(hname, "h_maxNearP15x15_ptBin%i_etaBin%i",ipt, ieta);
01231       sprintf(htit,  "maxNearP in 15x15 (%3.2f<|#eta|<%3.2f), (%2.0f<trkP<%3.0f)", lowEta, highEta, lowP, highP );
01232       h_maxNearP15x15[ipt][ieta] = dir.make<TH1F>(hname, htit, 220, -2.0, 20.0);
01233       h_maxNearP15x15[ipt][ieta] ->Sumw2();
01234      }
01235   }
01236 
01237   h_L1AlgoNames = fs->make<TH1I>("h_L1AlgoNames", "h_L1AlgoNames:Bin Labels", 128, -0.5, 127.5);  
01238 
01239   // Reconstructed Tracks
01240 
01241   h_PVTracksWt = fs->make<TH1F>("h_PVTracksWt", "h_PVTracksWt", 600, -0.1, 1.1);
01242 
01243   h_nTracks = fs->make<TH1F>("h_nTracks", "h_nTracks", 1000, -0.5, 999.5);
01244 
01245   sprintf(hname, "h_recEtaPt_0");
01246   sprintf(htit,  "h_recEtaPt (all tracks Eta vs pT)");
01247   h_recEtaPt_0 = fs->make<TH2F>(hname, htit, 30, -3.0,3.0, 15, genPartPBins);
01248 
01249   sprintf(hname, "h_recEtaP_0");
01250   sprintf(htit,  "h_recEtaP (all tracks Eta vs pT)");
01251   h_recEtaP_0 = fs->make<TH2F>(hname, htit, 30, -3.0,3.0, 15, genPartPBins);
01252 
01253   h_recPt_0  = fs->make<TH1F>("h_recPt_0",  "Pt (all tracks)",  15, genPartPBins);
01254   h_recP_0   = fs->make<TH1F>("h_recP_0",   "P  (all tracks)",  15, genPartPBins);
01255   h_recEta_0 = fs->make<TH1F>("h_recEta_0", "Eta (all tracks)", 60, -3.0,   3.0);
01256   h_recPhi_0 = fs->make<TH1F>("h_recPhi_0", "Phi (all tracks)", 100, -3.2,   3.2);
01257   //-------------------------
01258   sprintf(hname, "h_recEtaPt_1");
01259   sprintf(htit,  "h_recEtaPt (all good tracks Eta vs pT)");
01260   h_recEtaPt_1 = fs->make<TH2F>(hname, htit, 30, -3.0,3.0, 15, genPartPBins);
01261 
01262   sprintf(hname, "h_recEtaP_1");
01263   sprintf(htit,  "h_recEtaP (all good tracks Eta vs pT)");
01264   h_recEtaP_1 = fs->make<TH2F>(hname, htit, 30, -3.0, 3.0, 15, genPartPBins);
01265 
01266   h_recPt_1  = fs->make<TH1F>("h_recPt_1",  "Pt (all good tracks)",  15, genPartPBins);
01267   h_recP_1   = fs->make<TH1F>("h_recP_1",   "P  (all good tracks)",  15, genPartPBins);
01268   h_recEta_1 = fs->make<TH1F>("h_recEta_1", "Eta (all good tracks)", 60, -3.0,   3.0);
01269   h_recPhi_1 = fs->make<TH1F>("h_recPhi_1", "Phi (all good tracks)", 100, -3.2,   3.2);
01270   //-------------------------
01271   sprintf(hname, "h_recEtaPt_2");
01272   sprintf(htit,  "h_recEtaPt (charge isolation Eta vs pT)");
01273   h_recEtaPt_2 = fs->make<TH2F>(hname, htit, 30, -3.0, 3.0, 15, genPartPBins);
01274 
01275   sprintf(hname, "h_recEtaP_2");
01276   sprintf(htit,  "h_recEtaP (charge isolation Eta vs pT)");
01277   h_recEtaP_2 = fs->make<TH2F>(hname, htit, 30, -3.0, 3.0, 15, genPartPBins);
01278   
01279   h_recPt_2  = fs->make<TH1F>("h_recPt_2",  "Pt (charge isolation)",  15, genPartPBins);
01280   h_recP_2   = fs->make<TH1F>("h_recP_2",   "P  (charge isolation)",  15, genPartPBins);
01281   h_recEta_2 = fs->make<TH1F>("h_recEta_2", "Eta (charge isolation)", 60, -3.0,   3.0);
01282   h_recPhi_2 = fs->make<TH1F>("h_recPhi_2", "Phi (charge isolation)", 100, -3.2,   3.2);
01283 
01284 
01285   tree = fs->make<TTree>("tree", "tree");
01286   tree->SetAutoSave(10000);
01287 
01288 
01289   tree->Branch("t_EvtNo"              ,&t_EvtNo               ,"t_EvtNo/I");
01290   tree->Branch("t_RunNo"              ,&t_RunNo               ,"t_RunNo/I");
01291   tree->Branch("t_Lumi"               ,&t_Lumi                ,"t_Lumi/I");
01292   tree->Branch("t_Bunch"              ,&t_Bunch               ,"t_Bunch/I");
01293 
01294 
01295   t_PVx               = new std::vector<double>();
01296   t_PVy               = new std::vector<double>();
01297   t_PVz               = new std::vector<double>();
01298   t_PVisValid         = new std::vector<int>();
01299   t_PVndof            = new std::vector<int>(); 
01300   t_PVNTracks         = new std::vector<int>();
01301   t_PVNTracksWt       = new std::vector<int>();
01302   t_PVTracksSumPt     = new std::vector<double>();
01303   t_PVTracksSumPtWt   = new std::vector<double>();
01304   t_PVNTracksHP       = new std::vector<int>();
01305   t_PVNTracksHPWt     = new std::vector<int>();
01306   t_PVTracksSumPtHP   = new std::vector<double>();
01307   t_PVTracksSumPtHPWt = new std::vector<double>();
01308 
01309   tree->Branch("PVx"                  ,"vector<double>"      ,&t_PVx);
01310   tree->Branch("PVy"                  ,"vector<double>"      ,&t_PVy);
01311   tree->Branch("PVz"                  ,"vector<double>"      ,&t_PVz);
01312   tree->Branch("PVisValid"            ,"vector<int>"         ,&t_PVisValid);
01313   tree->Branch("PVndof"               ,"vector<int>"         ,&t_PVndof);
01314   tree->Branch("PVNTracks"            ,"vector<int>"         ,&t_PVNTracks);
01315   tree->Branch("PVNTracksWt"          ,"vector<int>"         ,&t_PVNTracksWt);
01316   tree->Branch("t_PVTracksSumPt"      ,"vector<double>"      ,&t_PVTracksSumPt);
01317   tree->Branch("t_PVTracksSumPtWt"    ,"vector<double>"      ,&t_PVTracksSumPtWt);
01318   tree->Branch("PVNTracksHP"          ,"vector<int>"         ,&t_PVNTracksHP);
01319   tree->Branch("PVNTracksHPWt"        ,"vector<int>"         ,&t_PVNTracksHPWt);
01320   tree->Branch("t_PVTracksSumPtHP"    ,"vector<double>"      ,&t_PVTracksSumPtHP);
01321   tree->Branch("t_PVTracksSumPtHPWt"  ,"vector<double>"      ,&t_PVTracksSumPtHPWt);
01322 
01323   //----- L1Trigger information
01324   for(int i=0; i<128; i++) t_L1Decision[i]=0;
01325   t_L1AlgoNames       = new std::vector<std::string>();
01326   t_L1PreScale        = new std::vector<int>();
01327   t_L1CenJetPt        = new std::vector<double>();
01328   t_L1CenJetEta       = new std::vector<double>();    
01329   t_L1CenJetPhi       = new std::vector<double>();
01330   t_L1FwdJetPt        = new std::vector<double>();
01331   t_L1FwdJetEta       = new std::vector<double>();
01332   t_L1FwdJetPhi       = new std::vector<double>();
01333   t_L1TauJetPt        = new std::vector<double>();
01334   t_L1TauJetEta       = new std::vector<double>();     
01335   t_L1TauJetPhi       = new std::vector<double>();
01336   t_L1MuonPt          = new std::vector<double>();
01337   t_L1MuonEta         = new std::vector<double>();     
01338   t_L1MuonPhi         = new std::vector<double>();
01339   t_L1IsoEMPt         = new std::vector<double>();
01340   t_L1IsoEMEta        = new std::vector<double>();
01341   t_L1IsoEMPhi        = new std::vector<double>();
01342   t_L1NonIsoEMPt      = new std::vector<double>();
01343   t_L1NonIsoEMEta     = new std::vector<double>();
01344   t_L1NonIsoEMPhi     = new std::vector<double>();
01345   t_L1METPt           = new std::vector<double>();
01346   t_L1METEta          = new std::vector<double>();
01347   t_L1METPhi          = new std::vector<double>();
01348   
01349   tree->Branch("t_L1Decision",        t_L1Decision,     "t_L1Decision[128]/I");
01350   tree->Branch("t_L1AlgoNames",       "vector<string>", &t_L1AlgoNames);
01351   tree->Branch("t_L1PreScale",        "vector<int>",    &t_L1PreScale);
01352   tree->Branch("t_L1CenJetPt",        "vector<double>", &t_L1CenJetPt);
01353   tree->Branch("t_L1CenJetEta",       "vector<double>", &t_L1CenJetEta);
01354   tree->Branch("t_L1CenJetPhi",       "vector<double>", &t_L1CenJetPhi);
01355   tree->Branch("t_L1FwdJetPt",        "vector<double>", &t_L1FwdJetPt);
01356   tree->Branch("t_L1FwdJetEta",       "vector<double>", &t_L1FwdJetEta);
01357   tree->Branch("t_L1FwdJetPhi",       "vector<double>", &t_L1FwdJetPhi);
01358   tree->Branch("t_L1TauJetPt",        "vector<double>", &t_L1TauJetPt);
01359   tree->Branch("t_L1TauJetEta",       "vector<double>", &t_L1TauJetEta);     
01360   tree->Branch("t_L1TauJetPhi",       "vector<double>", &t_L1TauJetPhi);
01361   tree->Branch("t_L1MuonPt",          "vector<double>", &t_L1MuonPt);
01362   tree->Branch("t_L1MuonEta",         "vector<double>", &t_L1MuonEta);
01363   tree->Branch("t_L1MuonPhi",         "vector<double>", &t_L1MuonPhi);
01364   tree->Branch("t_L1IsoEMPt",         "vector<double>", &t_L1IsoEMPt);
01365   tree->Branch("t_L1IsoEMEta",        "vector<double>", &t_L1IsoEMEta);
01366   tree->Branch("t_L1IsoEMPhi",        "vector<double>", &t_L1IsoEMPhi);
01367   tree->Branch("t_L1NonIsoEMPt",      "vector<double>", &t_L1NonIsoEMPt);
01368   tree->Branch("t_L1NonIsoEMEta",     "vector<double>", &t_L1NonIsoEMEta);
01369   tree->Branch("t_L1NonIsoEMPhi",     "vector<double>", &t_L1NonIsoEMPhi);
01370   tree->Branch("t_L1METPt",           "vector<double>", &t_L1METPt);
01371   tree->Branch("t_L1METEta",          "vector<double>", &t_L1METEta);
01372   tree->Branch("t_L1METPhi",          "vector<double>", &t_L1METPhi);
01373 
01374   t_jetPt            = new std::vector<double>();
01375   t_jetEta           = new std::vector<double>();
01376   t_jetPhi           = new std::vector<double>();
01377   t_nTrksJetCalo     = new std::vector<double>();  
01378   t_nTrksJetVtx      = new std::vector<double>();
01379   tree->Branch("t_jetPt",             "vector<double>",&t_jetPt);
01380   tree->Branch("t_jetEta",            "vector<double>",&t_jetEta);
01381   tree->Branch("t_jetPhi",            "vector<double>",&t_jetPhi);
01382   tree->Branch("t_nTrksJetCalo",      "vector<double>",&t_nTrksJetCalo);  
01383   tree->Branch("t_nTrksJetVtx",       "vector<double>",&t_nTrksJetVtx);
01384 
01385   t_trackPAll        = new std::vector<double>();
01386   t_trackEtaAll      = new std::vector<double>();
01387   t_trackPhiAll      = new std::vector<double>();
01388   t_trackPdgIdAll    = new std::vector<double>();
01389   t_trackPtAll       = new std::vector<double>();
01390   t_trackDxyAll      = new std::vector<double>();
01391   t_trackDzAll       = new std::vector<double>();
01392   t_trackDxyPVAll    = new std::vector<double>();
01393   t_trackDzPVAll     = new std::vector<double>();
01394   t_trackChiSqAll    = new std::vector<double>();
01395   tree->Branch("t_trackPAll",         "vector<double>", &t_trackPAll    );
01396   tree->Branch("t_trackPhiAll",       "vector<double>", &t_trackPhiAll  );
01397   tree->Branch("t_trackEtaAll",       "vector<double>", &t_trackEtaAll  );
01398   tree->Branch("t_trackPtAll",        "vector<double>", &t_trackPtAll    );
01399   tree->Branch("t_trackDxyAll",       "vector<double>", &t_trackDxyAll    );
01400   tree->Branch("t_trackDzAll",        "vector<double>", &t_trackDzAll    );
01401   tree->Branch("t_trackDxyPVAll",     "vector<double>", &t_trackDxyPVAll    );
01402   tree->Branch("t_trackDzPVAll",      "vector<double>", &t_trackDzPVAll    );
01403   tree->Branch("t_trackChiSqAll",     "vector<double>", &t_trackChiSqAll    );
01404   //tree->Branch("t_trackPdgIdAll",     "vector<double>", &t_trackPdgIdAll);
01405 
01406   t_trackP               = new std::vector<double>();
01407   t_trackPt              = new std::vector<double>();
01408   t_trackEta             = new std::vector<double>();
01409   t_trackPhi             = new std::vector<double>();
01410   t_trackEcalEta         = new std::vector<double>();
01411   t_trackEcalPhi         = new std::vector<double>();
01412   t_trackHcalEta         = new std::vector<double>();
01413   t_trackHcalPhi         = new std::vector<double>();
01414   t_trackNOuterHits      = new std::vector<int>();
01415   t_NLayersCrossed       = new std::vector<int>();
01416   t_trackDxy             = new std::vector<double>();
01417   t_trackDxyBS           = new std::vector<double>();
01418   t_trackDz              = new std::vector<double>();
01419   t_trackDzBS            = new std::vector<double>();
01420   t_trackDxyPV           = new std::vector<double>();
01421   t_trackDzPV            = new std::vector<double>();
01422   t_trackPVIdx           = new std::vector<int>();
01423   t_trackChiSq           = new std::vector<double>();
01424   t_trackHitsTOB         = new std::vector<int>(); 
01425   t_trackHitsTEC         = new std::vector<int>();
01426   t_trackHitInMissTOB    = new std::vector<int>(); 
01427   t_trackHitInMissTEC    = new std::vector<int>();  
01428   t_trackHitInMissTIB    = new std::vector<int>();  
01429   t_trackHitInMissTID    = new std::vector<int>();
01430   t_trackHitOutMissTOB   = new std::vector<int>();
01431   t_trackHitOutMissTEC   = new std::vector<int>(); 
01432   t_trackHitOutMissTIB   = new std::vector<int>();
01433   t_trackHitOutMissTID   = new std::vector<int>();
01434   t_trackHitInMissTIBTID = new std::vector<int>();  
01435   t_trackHitOutMissTOB   = new std::vector<int>();
01436   t_trackHitOutMissTEC   = new std::vector<int>(); 
01437   t_trackHitOutMissTIB   = new std::vector<int>();
01438   t_trackHitOutMissTID   = new std::vector<int>();
01439   t_trackHitOutMissTOBTEC= new std::vector<int>();
01440   t_trackHitInMeasTOB    = new std::vector<int>(); 
01441   t_trackHitInMeasTEC    = new std::vector<int>();  
01442   t_trackHitInMeasTIB    = new std::vector<int>();  
01443   t_trackHitInMeasTID    = new std::vector<int>();
01444   t_trackHitOutMeasTOB   = new std::vector<int>();
01445   t_trackHitOutMeasTEC   = new std::vector<int>(); 
01446   t_trackHitOutMeasTIB   = new std::vector<int>();
01447   t_trackHitOutMeasTID   = new std::vector<int>();
01448   t_trackOutPosOutHitDr  = new std::vector<double>();
01449   t_trackL               = new std::vector<double>();
01450 
01451   tree->Branch("t_trackP",            "vector<double>", &t_trackP            );
01452   tree->Branch("t_trackPt",           "vector<double>", &t_trackPt           );
01453   tree->Branch("t_trackEta",          "vector<double>", &t_trackEta          );
01454   tree->Branch("t_trackPhi",          "vector<double>", &t_trackPhi          );
01455   tree->Branch("t_trackEcalEta",      "vector<double>", &t_trackEcalEta      );
01456   tree->Branch("t_trackEcalPhi",      "vector<double>", &t_trackEcalPhi      );
01457   tree->Branch("t_trackHcalEta",      "vector<double>", &t_trackHcalEta      );
01458   tree->Branch("t_trackHcalPhi",      "vector<double>", &t_trackHcalPhi      );
01459 
01460   tree->Branch("t_trackNOuterHits",      "vector<int>",    &t_trackNOuterHits   );
01461   tree->Branch("t_NLayersCrossed",       "vector<int>",    &t_NLayersCrossed    );
01462   tree->Branch("t_trackHitsTOB",         "vector<int>",   &t_trackHitsTOB      ); 
01463   tree->Branch("t_trackHitsTEC",         "vector<int>",   &t_trackHitsTEC      );
01464   tree->Branch("t_trackHitInMissTOB",    "vector<int>",   &t_trackHitInMissTOB ); 
01465   tree->Branch("t_trackHitInMissTEC",    "vector<int>",   &t_trackHitInMissTEC );  
01466   tree->Branch("t_trackHitInMissTIB",    "vector<int>",   &t_trackHitInMissTIB );  
01467   tree->Branch("t_trackHitInMissTID",    "vector<int>",   &t_trackHitInMissTID );
01468   tree->Branch("t_trackHitInMissTIBTID", "vector<int>",   &t_trackHitInMissTIBTID );  
01469   tree->Branch("t_trackHitOutMissTOB",   "vector<int>",   &t_trackHitOutMissTOB);
01470   tree->Branch("t_trackHitOutMissTEC",   "vector<int>",   &t_trackHitOutMissTEC); 
01471   tree->Branch("t_trackHitOutMissTIB",   "vector<int>",   &t_trackHitOutMissTIB);
01472   tree->Branch("t_trackHitOutMissTID",   "vector<int>",   &t_trackHitOutMissTID);
01473   tree->Branch("t_trackHitOutMissTOBTEC","vector<int>",   &t_trackHitOutMissTOBTEC);
01474   tree->Branch("t_trackHitInMeasTOB",    "vector<int>",   &t_trackHitInMeasTOB ); 
01475   tree->Branch("t_trackHitInMeasTEC",    "vector<int>",   &t_trackHitInMeasTEC );  
01476   tree->Branch("t_trackHitInMeasTIB",    "vector<int>",   &t_trackHitInMeasTIB );  
01477   tree->Branch("t_trackHitInMeasTID",    "vector<int>",   &t_trackHitInMeasTID );
01478   tree->Branch("t_trackHitOutMeasTOB",   "vector<int>",   &t_trackHitOutMeasTOB);
01479   tree->Branch("t_trackHitOutMeasTEC",   "vector<int>",   &t_trackHitOutMeasTEC); 
01480   tree->Branch("t_trackHitOutMeasTIB",   "vector<int>",   &t_trackHitOutMeasTIB);
01481   tree->Branch("t_trackHitOutMeasTID",   "vector<int>",   &t_trackHitOutMeasTID);
01482   tree->Branch("t_trackOutPosOutHitDr",  "vector<double>", &t_trackOutPosOutHitDr);
01483   tree->Branch("t_trackL",               "vector<double>", &t_trackL);
01484 
01485   tree->Branch("t_trackDxy",          "vector<double>", &t_trackDxy     );
01486   tree->Branch("t_trackDxyBS",        "vector<double>", &t_trackDxyBS   );
01487   tree->Branch("t_trackDz",           "vector<double>", &t_trackDz      );
01488   tree->Branch("t_trackDzBS",         "vector<double>", &t_trackDzBS    );
01489   tree->Branch("t_trackDxyPV",        "vector<double>", &t_trackDxyPV   );
01490   tree->Branch("t_trackDzPV",         "vector<double>", &t_trackDzPV    );
01491   tree->Branch("t_trackChiSq",        "vector<double>", &t_trackChiSq   );
01492   tree->Branch("t_trackPVIdx",        "vector<int>",    &t_trackPVIdx   );
01493 
01494   t_maxNearP31x31     = new std::vector<double>();
01495   t_maxNearP21x21     = new std::vector<double>();
01496 
01497   tree->Branch("t_maxNearP31x31",     "vector<double>", &t_maxNearP31x31);
01498   tree->Branch("t_maxNearP21x21",     "vector<double>", &t_maxNearP21x21);
01499 
01500   t_ecalSpike11x11    = new std::vector<int>();
01501   t_e7x7              = new std::vector<double>();
01502   t_e9x9              = new std::vector<double>();
01503   t_e11x11            = new std::vector<double>();
01504   t_e15x15            = new std::vector<double>();
01505 
01506   tree->Branch("t_ecalSpike11x11",    "vector<int>",    &t_ecalSpike11x11);
01507   tree->Branch("t_e7x7",              "vector<double>", &t_e7x7);
01508   tree->Branch("t_e9x9",              "vector<double>", &t_e9x9);
01509   tree->Branch("t_e11x11",            "vector<double>", &t_e11x11);
01510   tree->Branch("t_e15x15",            "vector<double>", &t_e15x15);
01511 
01512   t_e7x7_10Sig        = new std::vector<double>();
01513   t_e9x9_10Sig        = new std::vector<double>();
01514   t_e11x11_10Sig      = new std::vector<double>();
01515   t_e15x15_10Sig      = new std::vector<double>();
01516   t_e7x7_15Sig        = new std::vector<double>();
01517   t_e9x9_15Sig        = new std::vector<double>();
01518   t_e11x11_15Sig      = new std::vector<double>();
01519   t_e15x15_15Sig      = new std::vector<double>();
01520   t_e7x7_20Sig        = new std::vector<double>();
01521   t_e9x9_20Sig        = new std::vector<double>();
01522   t_e11x11_20Sig      = new std::vector<double>();
01523   t_e15x15_20Sig      = new std::vector<double>();
01524   t_e7x7_25Sig        = new std::vector<double>();
01525   t_e9x9_25Sig        = new std::vector<double>();
01526   t_e11x11_25Sig      = new std::vector<double>();
01527   t_e15x15_25Sig      = new std::vector<double>();
01528   t_e7x7_30Sig        = new std::vector<double>();
01529   t_e9x9_30Sig        = new std::vector<double>();
01530   t_e11x11_30Sig      = new std::vector<double>();
01531   t_e15x15_30Sig      = new std::vector<double>();
01532 
01533   tree->Branch("t_e7x7_10Sig"        ,"vector<double>", &t_e7x7_10Sig);
01534   tree->Branch("t_e9x9_10Sig"        ,"vector<double>", &t_e9x9_10Sig);
01535   tree->Branch("t_e11x11_10Sig"      ,"vector<double>", &t_e11x11_10Sig);
01536   tree->Branch("t_e15x15_10Sig"      ,"vector<double>", &t_e15x15_10Sig);
01537   tree->Branch("t_e7x7_15Sig"        ,"vector<double>", &t_e7x7_15Sig);
01538   tree->Branch("t_e9x9_15Sig"        ,"vector<double>", &t_e9x9_15Sig);
01539   tree->Branch("t_e11x11_15Sig"      ,"vector<double>", &t_e11x11_15Sig);
01540   tree->Branch("t_e15x15_15Sig"      ,"vector<double>", &t_e15x15_15Sig);
01541   tree->Branch("t_e7x7_20Sig"        ,"vector<double>", &t_e7x7_20Sig);
01542   tree->Branch("t_e9x9_20Sig"        ,"vector<double>", &t_e9x9_20Sig);
01543   tree->Branch("t_e11x11_20Sig"      ,"vector<double>", &t_e11x11_20Sig);
01544   tree->Branch("t_e15x15_20Sig"      ,"vector<double>", &t_e15x15_20Sig);
01545   tree->Branch("t_e7x7_25Sig"        ,"vector<double>", &t_e7x7_25Sig);
01546   tree->Branch("t_e9x9_25Sig"        ,"vector<double>", &t_e9x9_25Sig);
01547   tree->Branch("t_e11x11_25Sig"      ,"vector<double>", &t_e11x11_25Sig);
01548   tree->Branch("t_e15x15_25Sig"      ,"vector<double>", &t_e15x15_25Sig);
01549   tree->Branch("t_e7x7_30Sig"        ,"vector<double>", &t_e7x7_30Sig);
01550   tree->Branch("t_e9x9_30Sig"        ,"vector<double>", &t_e9x9_30Sig);
01551   tree->Branch("t_e11x11_30Sig"      ,"vector<double>", &t_e11x11_30Sig);
01552   tree->Branch("t_e15x15_30Sig"      ,"vector<double>", &t_e15x15_30Sig);
01553 
01554   if (doMC) {
01555     t_esim7x7              = new std::vector<double>();
01556     t_esim9x9              = new std::vector<double>();
01557     t_esim11x11            = new std::vector<double>();
01558     t_esim15x15            = new std::vector<double>();
01559 
01560     t_esim7x7Matched       = new std::vector<double>();
01561     t_esim9x9Matched       = new std::vector<double>();
01562     t_esim11x11Matched     = new std::vector<double>();
01563     t_esim15x15Matched     = new std::vector<double>();
01564 
01565     t_esim7x7Rest          = new std::vector<double>();
01566     t_esim9x9Rest          = new std::vector<double>();
01567     t_esim11x11Rest        = new std::vector<double>();
01568     t_esim15x15Rest        = new std::vector<double>();
01569 
01570     t_esim7x7Photon        = new std::vector<double>();
01571     t_esim9x9Photon        = new std::vector<double>();
01572     t_esim11x11Photon      = new std::vector<double>();
01573     t_esim15x15Photon      = new std::vector<double>();
01574 
01575     t_esim7x7NeutHad       = new std::vector<double>();
01576     t_esim9x9NeutHad       = new std::vector<double>();
01577     t_esim11x11NeutHad     = new std::vector<double>();
01578     t_esim15x15NeutHad     = new std::vector<double>();
01579 
01580     t_esim7x7CharHad       = new std::vector<double>();
01581     t_esim9x9CharHad       = new std::vector<double>();
01582     t_esim11x11CharHad     = new std::vector<double>();
01583     t_esim15x15CharHad     = new std::vector<double>();
01584 
01585     t_trkEcalEne           = new std::vector<double>();
01586     t_simTrackP            = new std::vector<double>();
01587     t_esimPdgId            = new std::vector<double>();
01588 
01589     tree->Branch("t_esim7x7",             "vector<double>", &t_esim7x7);
01590     tree->Branch("t_esim9x9",             "vector<double>", &t_esim9x9);
01591     tree->Branch("t_esim11x11",           "vector<double>", &t_esim11x11);
01592     tree->Branch("t_esim15x15",           "vector<double>", &t_esim15x15);
01593 
01594     tree->Branch("t_esim7x7Matched",      "vector<double>", &t_esim7x7Matched);
01595     tree->Branch("t_esim9x9Matched",      "vector<double>", &t_esim9x9Matched);
01596     tree->Branch("t_esim11x11Matched",    "vector<double>", &t_esim11x11Matched);
01597     tree->Branch("t_esim15x15Matched",    "vector<double>", &t_esim15x15Matched);
01598 
01599     tree->Branch("t_esim7x7Rest",         "vector<double>", &t_esim7x7Rest);
01600     tree->Branch("t_esim9x9Rest",         "vector<double>", &t_esim9x9Rest);
01601     tree->Branch("t_esim11x11Rest",       "vector<double>", &t_esim11x11Rest);
01602     tree->Branch("t_esim15x15Rest",       "vector<double>", &t_esim15x15Rest);
01603 
01604     tree->Branch("t_esim7x7Photon",       "vector<double>", &t_esim7x7Photon);
01605     tree->Branch("t_esim9x9Photon",       "vector<double>", &t_esim9x9Photon);
01606     tree->Branch("t_esim11x11Photon",     "vector<double>", &t_esim11x11Photon);
01607     tree->Branch("t_esim15x15Photon",     "vector<double>", &t_esim15x15Photon);
01608 
01609 ;
01610     tree->Branch("t_esim7x7NeutHad",      "vector<double>", &t_esim7x7NeutHad);
01611     tree->Branch("t_esim9x9NeutHad",      "vector<double>", &t_esim9x9NeutHad);
01612     tree->Branch("t_esim11x11NeutHad",    "vector<double>", &t_esim11x11NeutHad);
01613     tree->Branch("t_esim15x15NeutHad",    "vector<double>", &t_esim15x15NeutHad);
01614 
01615     tree->Branch("t_esim7x7CharHad",      "vector<double>", &t_esim7x7CharHad);
01616     tree->Branch("t_esim9x9CharHad",      "vector<double>", &t_esim9x9CharHad);
01617     tree->Branch("t_esim11x11CharHad",    "vector<double>", &t_esim11x11CharHad);
01618     tree->Branch("t_esim15x15CharHad",    "vector<double>", &t_esim15x15CharHad);
01619 
01620     tree->Branch("t_trkEcalEne",          "vector<double>", &t_trkEcalEne);
01621     tree->Branch("t_simTrackP",           "vector<double>", &t_simTrackP);
01622     tree->Branch("t_esimPdgId",           "vector<double>", &t_esimPdgId);
01623   }
01624 
01625   t_maxNearHcalP3x3      = new std::vector<double>();
01626   t_maxNearHcalP5x5      = new std::vector<double>();
01627   t_maxNearHcalP7x7      = new std::vector<double>();
01628   t_h3x3                 = new std::vector<double>();
01629   t_h5x5                 = new std::vector<double>();
01630   t_h7x7                 = new std::vector<double>();
01631   t_h3x3Sig              = new std::vector<double>();
01632   t_h5x5Sig              = new std::vector<double>();
01633   t_h7x7Sig              = new std::vector<double>();
01634   t_infoHcal             = new std::vector<int>();
01635 
01636   if (doMC) {
01637     t_trkHcalEne           = new std::vector<double>();
01638     t_hsim3x3              = new std::vector<double>();
01639     t_hsim5x5              = new std::vector<double>();
01640     t_hsim7x7              = new std::vector<double>();
01641     t_hsim3x3Matched       = new std::vector<double>();
01642     t_hsim5x5Matched       = new std::vector<double>();
01643     t_hsim7x7Matched       = new std::vector<double>();
01644     t_hsim3x3Rest          = new std::vector<double>();
01645     t_hsim5x5Rest          = new std::vector<double>();
01646     t_hsim7x7Rest          = new std::vector<double>();
01647     t_hsim3x3Photon        = new std::vector<double>();
01648     t_hsim5x5Photon        = new std::vector<double>();
01649     t_hsim7x7Photon        = new std::vector<double>();
01650     t_hsim3x3NeutHad       = new std::vector<double>();
01651     t_hsim5x5NeutHad       = new std::vector<double>();
01652     t_hsim7x7NeutHad       = new std::vector<double>();
01653     t_hsim3x3CharHad       = new std::vector<double>();
01654     t_hsim5x5CharHad       = new std::vector<double>();
01655     t_hsim7x7CharHad       = new std::vector<double>();
01656   }
01657 
01658   tree->Branch("t_maxNearHcalP3x3",     "vector<double>", &t_maxNearHcalP3x3);
01659   tree->Branch("t_maxNearHcalP5x5",     "vector<double>", &t_maxNearHcalP5x5);
01660   tree->Branch("t_maxNearHcalP7x7",     "vector<double>", &t_maxNearHcalP7x7);
01661   tree->Branch("t_h3x3",                "vector<double>", &t_h3x3);
01662   tree->Branch("t_h5x5",                "vector<double>", &t_h5x5);
01663   tree->Branch("t_h7x7",                "vector<double>", &t_h7x7);
01664   tree->Branch("t_h3x3Sig",             "vector<double>", &t_h3x3Sig);
01665   tree->Branch("t_h5x5Sig",             "vector<double>", &t_h5x5Sig);
01666   tree->Branch("t_h7x7Sig",             "vector<double>", &t_h7x7Sig);
01667   tree->Branch("t_infoHcal",            "vector<int>",    &t_infoHcal);
01668 
01669   if (doMC) {
01670     tree->Branch("t_trkHcalEne",          "vector<double>", &t_trkHcalEne);
01671     tree->Branch("t_hsim3x3",             "vector<double>", &t_hsim3x3);
01672     tree->Branch("t_hsim5x5",             "vector<double>", &t_hsim5x5);
01673     tree->Branch("t_hsim7x7",             "vector<double>", &t_hsim7x7);
01674     tree->Branch("t_hsim3x3Matched",      "vector<double>", &t_hsim3x3Matched);
01675     tree->Branch("t_hsim5x5Matched",      "vector<double>", &t_hsim5x5Matched);
01676     tree->Branch("t_hsim7x7Matched",      "vector<double>", &t_hsim7x7Matched);
01677     tree->Branch("t_hsim3x3Rest",         "vector<double>", &t_hsim3x3Rest);
01678     tree->Branch("t_hsim5x5Rest",         "vector<double>", &t_hsim5x5Rest);
01679     tree->Branch("t_hsim7x7Rest",         "vector<double>", &t_hsim7x7Rest);
01680     tree->Branch("t_hsim3x3Photon",       "vector<double>", &t_hsim3x3Photon);
01681     tree->Branch("t_hsim5x5Photon",       "vector<double>", &t_hsim5x5Photon);
01682     tree->Branch("t_hsim7x7Photon",       "vector<double>", &t_hsim7x7Photon);
01683     tree->Branch("t_hsim3x3NeutHad",      "vector<double>", &t_hsim3x3NeutHad);
01684     tree->Branch("t_hsim5x5NeutHad",      "vector<double>", &t_hsim5x5NeutHad);
01685     tree->Branch("t_hsim7x7NeutHad",      "vector<double>", &t_hsim7x7NeutHad);
01686     tree->Branch("t_hsim3x3CharHad",      "vector<double>", &t_hsim3x3CharHad);
01687     tree->Branch("t_hsim5x5CharHad",      "vector<double>", &t_hsim5x5CharHad);
01688     tree->Branch("t_hsim7x7CharHad",      "vector<double>", &t_hsim7x7CharHad);
01689   }
01690   tree->Branch("t_nTracks",             &t_nTracks,       "t_nTracks/I");
01691 
01692 }
01693 
01694 
01695 double IsolatedTracksNxN::DeltaPhi(double v1, double v2) {
01696   // Computes the correctly normalized phi difference
01697   // v1, v2 = phi of object 1 and 2
01698   
01699   double pi    = 3.141592654;
01700   double twopi = 6.283185307;
01701   
01702   double diff = std::abs(v2 - v1);
01703   double corr = twopi - diff;
01704   if (diff < pi){ return diff;} else { return corr;} 
01705 }
01706 
01707 double IsolatedTracksNxN::DeltaR(double eta1, double phi1, double eta2, double phi2) {
01708   double deta = eta1 - eta2;
01709   double dphi = DeltaPhi(phi1, phi2);
01710   return std::sqrt(deta*deta + dphi*dphi);
01711 }
01712 
01713 void IsolatedTracksNxN::printTrack(const reco::Track* pTrack) {
01714   
01715   std::string theTrackQuality = "highPurity";
01716   reco::TrackBase::TrackQuality trackQuality_ = reco::TrackBase::qualityByName(theTrackQuality);
01717 
01718   std::cout << " Reference Point " << pTrack->referencePoint() <<"\n"
01719             << " TrackMmentum " << pTrack->momentum()
01720             << " (pt,eta,phi)(" << pTrack->pt()<<","<<pTrack->eta()<<","<<pTrack->phi()<<")"
01721             << " p " << pTrack->p() << "\n"
01722             << " Normalized chi2 " << pTrack->normalizedChi2() <<"  charge " << pTrack->charge()
01723             << " qoverp() " << pTrack->qoverp() <<"+-" << pTrack->qoverpError()
01724             << " d0 " << pTrack->d0() << "\n"
01725             << " NValidHits " << pTrack->numberOfValidHits() << "  NLostHits " << pTrack->numberOfLostHits()
01726             << " TrackQuality " << pTrack->qualityName(trackQuality_) << " " << pTrack->quality(trackQuality_) 
01727             << std::endl;
01728   
01729   if( printTrkHitPattern_ ) {
01730     const reco::HitPattern& p = pTrack->hitPattern();
01731     const reco::HitPattern& p1 = pTrack->trackerExpectedHitsInner();
01732     const reco::HitPattern& p2 = pTrack->trackerExpectedHitsOuter();
01733 
01734     std::cout<<"default " << std::endl;
01735     for (int i=0; i<p.numberOfHits(); i++) {
01736       p.printHitPattern(i, std::cout);
01737     }
01738     std::cout<<"trackerExpectedHitsInner() " << std::endl;
01739     for (int i=0; i<p1.numberOfHits(); i++) {
01740       p1.printHitPattern(i, std::cout);
01741     }
01742     std::cout<<"trackerExpectedHitsOuter() " << std::endl;
01743     for (int i=0; i<p2.numberOfHits(); i++) {
01744       p2.printHitPattern(i, std::cout);
01745     }
01746 
01747 
01748     std::cout << "\n \t trackerLayersWithMeasurement() "     << p.trackerLayersWithMeasurement() 
01749               << "\n \t pixelLayersWithMeasurement() "       << p.pixelLayersWithMeasurement() 
01750               << "\n \t stripLayersWithMeasurement() "       << p.stripLayersWithMeasurement()  
01751               << "\n \t pixelBarrelLayersWithMeasurement() " << p.pixelBarrelLayersWithMeasurement()
01752               << "\n \t pixelEndcapLayersWithMeasurement() " << p.pixelEndcapLayersWithMeasurement()
01753               << "\n \t stripTIBLayersWithMeasurement() "    << p.stripTIBLayersWithMeasurement()
01754               << "\n \t stripTIDLayersWithMeasurement() "    << p.stripTIDLayersWithMeasurement()
01755               << "\n \t stripTOBLayersWithMeasurement() "    << p.stripTOBLayersWithMeasurement()
01756               << "\n \t stripTECLayersWithMeasurement() "    << p.stripTECLayersWithMeasurement()
01757               << std::endl;
01758 
01759   }
01760 }
01761 
01762 //define this as a plug-in
01763 DEFINE_FWK_MODULE(IsolatedTracksNxN);
/data/refman/pasoursint/CMSSW_5_2_7_hltpatch2/src/Calibration/IsolatedParticles/plugins/IsolatedTracksNxN.cc
