CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_4_1_8_patch13/src/Calibration/IsolatedParticles/plugins/IsolatedTracksNxN.cc

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