XANADOO: hzz_4e_analysis.cpp Source File

00001 
00002 #include <XANADOO/domain/interface/XANAEsdEvent.h>
00003 #include <XANADOO/XANAMcInfo/interface/XANAGeneratorEvent.h>
00004 #include <XANADOO/XANAMcInfo/interface/XANAGeantEvent.h>
00005 #include <XANADOO/XANAMcInfo/interface/XANAGeneratorParticle.h>
00006 #include <XANADOO/XANAElectronCandidate/interface/XANAElectronCandidate.h>
00007 #include <XANADOO/XANAElectronCandidate/interface/XANAElectronTrack.h>
00008 #include <XANADOO/XANAElectronCandidate/interface/XANAElectronGSFTrack.h>
00009 #include <XANADOO/XANAVertex/interface/XANAVertex.h>
00010 #include <XANADOO/XANATracks/interface/XANATrack.h>
00011 #include <XANADOO/XANAnalysisTools/interface/XANAmbiguityResolve.h>
00012 
00013 #include <XANADOO/XANAnalysisZZ/interface/XANASelectHZZ.h>
00014 #include <XANADOO/XANAnalysisZZ/interface/XANASelectMcHZZ.h>
00015 #include <XANADOO/XANAnalysisTools/interface/XANAEvtIsolation.h>
00016 #include <XANADOO/XANAnalysisTools/interface/XANAEvtRegionIsolation.h>
00017 #include <XANADOO/XANAnalysisTools/interface/XANAIsolation.h>
00018 #include <XANADOO/XANAnalysisTools/interface/XANAnalysisTrigger.h>
00019 #include <XANADOO/XANATriggerInfo/interface/XANATriggerInfo.h>
00020 #include <XANADOO/XANAnalysisTools/interface/XANAElectronRecGen.h>
00021 #include <XANADOO/XANAnalysisTools/interface/XANAEventSelect.h>
00022 #include <XANADOO/XANAnalysisTools/interface/XANAElectronClassification.h>
00023 #include <XANADOO/XANAnalysisTools/interface/XANAElectronSelect.h>
00024 #include <XANADOO/XANAnalysisTools/interface/XANAElectronCorrection.h>
00025 #include <XANADOO/XANAnalysisTools/interface/XANAElectronMomentum.h>
00026 
00027 #include "CLHEP/Geometry/Point3D.h"
00028 
00029 #include "TTree.h"
00030 #include "TFile.h"
00031 #include "TBranch.h"
00032 #include "TObjArray.h"
00033 #include "TH1.h"
00034 #include "TH2.h"
00035 #include "TProfile.h"
00036 #include "TStopwatch.h"
00037 #include "TTree.h"
00038 #include "TChain.h"
00039 #include <TBits.h>
00040 
00041 #include <string>
00042 
00043 
00044 using namespace std;
00045 
00046 void initSettings(TChain *T, XANAEsdEvent **eventaddress, XANAGeneratorEvent **generatorEventaddress, XANAGeantEvent **geantEventaddress)
00047 {
00048   // example of how to suppress reading of branches
00049   // has to be put before setting the branch addresses!
00050   
00051   // To improve the speed set status for all branches to 0 and then include those needed
00052 
00053   //  T->SetBranchStatus("Event.*",0);
00054   //  T->SetBranchStatus("GeneratorEvent.*",0); 
00055   T->SetBranchStatus("GeantEvent.*",0);
00056 
00057   //T->SetBranchStatus("Event.caloRecHits_",0);
00058   T->SetBranchStatus("Event.allTrackHits_",0); 
00059   // T->SetBranchStatus("Event.emClusters_",0);  
00060   T->SetBranchStatus("Event.preshClusters_",0);  
00061   T->SetBranchStatus("Event.hadClusters_",0);
00062   T->SetBranchStatus("Event.vertices_",0);        
00063   T->SetBranchStatus("Event.muonCandidates_",0); 
00064   T->SetBranchStatus("Event.muonTracks_",0);
00065   T->SetBranchStatus("Event.jets_",0);
00066   T->SetBranchStatus("Event.met_",0);
00067   
00068   /*
00069     T->SetBranchStatus("Event.electronCandidates_",1);
00070     T->SetBranchStatus("Event.numberOfElectronCandidates_",1);
00071     T->SetBranchStatus("Event.tracks_",1);
00072     T->SetBranchStatus("Event.triggerInfo_",1);  
00073   */ 
00074   
00075   // Set the branch addresses: first reconstructed event and then MC information (generator and GEANT event)
00076   T->SetBranchAddress("Event.",eventaddress);  
00077   T->SetBranchAddress("GeneratorEvent.",generatorEventaddress); 
00078   T->SetBranchAddress("GeantEvent.",geantEventaddress);  
00079 }
00080 
00081 void clean(XANAEsdEvent *event, XANAGeneratorEvent *generatorEvent, XANAGeantEvent *geantEvent){
00082   event->clear(); 
00083   generatorEvent->clear();
00084   geantEvent->clear();
00085   
00086 }
00087 
00088 int main(int argc, char** argv)
00089 {
00090   
00091   if(!argv[1]) {
00092     cout<<">>> No argument, please give an argument if you want me to do something <<<"<<endl;
00093     cout<<">>> list of possible arguments for the moment: hzz120 hzz150 and tt4e   <<<"<<endl;
00094     return 0;
00095   }
00096   cout<< ">>> Hello, welcome in hzz_4e_analysis program <<<"<<endl;
00097   
00098   TStopwatch stopWatch;
00099   stopWatch.Start();
00100 
00101   // Steering flags
00102   bool momentumCorrection = true;
00103   bool federicoEscaleCorrection = true;
00104   bool egammaEscaleCorrection = true;
00105   
00106   
00107   // output file definition 
00108   char name[80];  
00109   if ( federicoEscaleCorrection && momentumCorrection ) 
00110     sprintf(name,"HZZ_corr_momentum_%s.root",argv[1]);
00111   else if ( federicoEscaleCorrection && !momentumCorrection  ) 
00112     sprintf(name,"HZZ_corr_no_momentum_%s.root",argv[1]);
00113   else if ( !federicoEscaleCorrection && !momentumCorrection  ) 
00114     sprintf(name,"HZZ_no_corr_no_momentum_%s.root",argv[1]);
00115   else if ( !federicoEscaleCorrection && !egammaEscaleCorrection && !momentumCorrection  ) 
00116     sprintf(name,"HZZ_no_corr_at_all_%s.root",argv[1]);
00117   
00118   
00119   TFile f(name,"RECREATE");
00120   f.cd();
00121   
00122   // Histogram booking 
00123   
00124   // general info
00125   TH1F hElectronNumber("hElectronNumber", "Number of ElectronCandidates", 15, 0., 15.); 
00126   TH1F hElectronNumberAfterAmb("hElectronNumberAfterAmb", "Number of ElectronCandidates after ambiguity resolving", 15, 0., 15.);   
00127   TH1F hElectronNumberBeforeAfter("hElectronNumberBeforeAfter", "Number of ElectronCandidates before - after ambiguity resolving", 15, 0., 15.);   
00128   
00129   // isolation
00130   TH1F hIsoNumberTracks("hIsoNumberTracks", "Number of tracks inside cone", 40, 0., 40.);
00131   TH1F hIsoPtTracks("hIsoPtTracks", "Pt summ of all tracks inside cone", 100, 0., 100.);
00132   TH1F hIsoNumberECinsideCone("hIsoNumberECinsideCone", "Number of e candidates inside cone", 10, 0., 10.);
00133   TH1F hIsoPtMax("hIsoPtMax", "Max pT of tracks inside cone", 100, 0., 100.);
00134   
00135   TH1F hIsoRegionNumberTracks("hIsoRegionNumberTracks", "Number of tracks inside region", 40, 0., 40.);
00136   TH1F hIsoRegionPtTracks("hIsoRegionPtTracks", "Pt summ of all tracks inside region", 100, 0., 100.);
00137   TH1F hIsoRegionPtMax("hIsoRegionPtMax", "Max pT of tracks inside region", 100, 0., 100.);
00138   
00139   // Mc information
00140   TH1F hElectronRecGenPt("hElectronRecGenPt","ElectronCandidate: Rec/Gen pT", 500, 0., 5.);
00141   TH1F hElectronRecGenDr("hElectronRecGenDr","ElectronCandidate: Rec/Gen Delta R", 500, 0., 5.);
00142   
00143   // classification
00144   TH1F hElectronClasses("hElectronClasses","Electron Candidates Classification",23, -2., 21.);
00145   TH1F hElectronClass0("hElectronClass0","Electron Candidates Class: 0, Barrel-Golden",100, 0., 1.5);
00146   TH1F hElectronClass1("hElectronClass1","Electron Candidates Class: 1, Barrel-Big Brem",100, 0., 1.5);
00147   TH1F hElectronClass2("hElectronClass2","Electron Candidates Class: 2, Barrel-Narrow",100, 0., 1.5);
00148   TH1F hElectronClass3("hElectronClass3","Electron Candidates Class: 3, Barrel-Showering",100, 0., 1.5);
00149   TH1F hElectronClass4("hElectronClass4","Electron Candidates Class: 4, Barrel-Cracks",100, 0., 1.5);
00150   TH1F hElectronClass10("hElectronClass10","Electron Candidates Class: 0, Endcap-Golden",100, 0., 1.5);
00151   TH1F hElectronClass11("hElectronClass11","Electron Candidates Class: 1, Endcap-Big Brem",100, 0., 1.5);
00152   TH1F hElectronClass12("hElectronClass12","Electron Candidates Class: 2, Endcap-Narrow",100, 0., 1.5);
00153   TH1F hElectronClass13("hElectronClass13","Electron Candidates Class: 3, Endcap-Showering",100, 0., 1.5);
00154   
00155   TProfile hElectronClassProfile0("hElectronClassProfile0","Electron Candidates ClassProfile: 0, Barrel-Golden",100, 0., 2.7);
00156   TProfile hElectronClassProfile1("hElectronClassProfile1","Electron Candidates ClassProfile: 1, Barrel-Big Brem",100, 0., 2.7);
00157   TProfile hElectronClassProfile2("hElectronClassProfile2","Electron Candidates ClassProfile: 2, Barrel-Narrow",100, 0., 2.7);
00158   TProfile hElectronClassProfile3("hElectronClassProfile3","Electron Candidates ClassProfile: 3, Barrel-Showering",100, 0., 2.7);
00159   TProfile hElectronClassProfile4("hElectronClassProfile4","Electron Candidates ClassProfile: 4, Barrel-Cracks",100, 0., 2.7);
00160   TProfile hElectronClassProfile10("hElectronClassProfile10","Electron Candidates ClassProfile: 0, Endcap-Golden",100, 0., 2.7);
00161   TProfile hElectronClassProfile11("hElectronClassProfile11","Electron Candidates ClassProfile: 1, Endcap-Big Brem",100, 0., 2.7);
00162   TProfile hElectronClassProfile12("hElectronClassProfile12","Electron Candidates ClassProfile: 2, Endcap-Narrow",100, 0., 2.7);
00163   TProfile hElectronClassProfile13("hElectronClassProfile13","Electron Candidates ClassProfile: 3, Endcap-Showering",100, 0., 2.7);
00164   
00165   TH1F hElectronEta("hElectronEta","Electron Candidates Eta",100, 0., 2.7);
00166   TH1F hElectronEtaClass0("hElectronEtaClass0","Electron Candidates Eta Class: 0, Barrel-Golden",100, 0., 2.7);
00167   TH1F hElectronEtaClass1("hElectronEtaClass1","Electron Candidates Eta Class: 1, Barrel-Big Brem",100, 0., 2.7);
00168   TH1F hElectronEtaClass2("hElectronEtaClass2","Electron Candidates Eta Class: 2, Barrel-Narrow",100, 0., 2.7);
00169   TH1F hElectronEtaClass3("hElectronEtaClass3","Electron Candidates Eta Class: 3, Barrel-Showering",100, 0., 2.7);
00170   TH1F hElectronEtaClass4("hElectronEtaClass4","Electron Candidates Eta Class: 4, Barrel-Cracks",100, 0., 2.7);
00171   TH1F hElectronEtaClass10("hElectronEtaClass10","Electron Candidates Eta Class: 0, Endcap-Golden",100, 0., 2.7);
00172   TH1F hElectronEtaClass11("hElectronEtaClass11","Electron Candidates Eta Class: 1, Endcap-Big Brem",100, 0., 2.7);
00173   TH1F hElectronEtaClass12("hElectronEtaClass12","Electron Candidates Eta Class: 2, Endcap-Narrow",100, 0., 2.7);
00174   TH1F hElectronEtaClass13("hElectronEtaClass13","Electron Candidates Eta Class: 3, Endcap-Showering",100, 0., 2.7);
00175   
00176   // after corrections
00177   TH1F hElectronClassCorr0("hElectronClassCorr0","Electron Candidates Class After Corr: 0, Barrel-Golden",100, 0., 1.5);
00178   TH1F hElectronClassCorr1("hElectronClassCorr1","Electron Candidates Class After Corr: 1, Barrel-Big Brem",100, 0., 1.5);
00179   TH1F hElectronClassCorr2("hElectronClassCorr2","Electron Candidates Class After Corr: 2, Barrel-Narrow",100, 0., 1.5);
00180   TH1F hElectronClassCorr3("hElectronClassCorr3","Electron Candidates Class After Corr: 3, Barrel-Showering",100, 0., 1.5);
00181   TH1F hElectronClassCorr4("hElectronClassCorr4","Electron Candidates Class After Corr: 4, Barrel-Cracks",100, 0., 1.5);
00182   TH1F hElectronClassCorr10("hElectronClassCorr10","Electron Candidates Class After Corr: 0, Endcap-Golden",100, 0., 1.5);
00183   TH1F hElectronClassCorr11("hElectronClassCorr11","Electron Candidates Class After Corr: 1, Endcap-Big Brem",100, 0., 1.5);
00184   TH1F hElectronClassCorr12("hElectronClassCorr12","Electron Candidates Class After Corr: 2, Endcap-Narrow",100, 0., 1.5);
00185   TH1F hElectronClassCorr13("hElectronClassCorr13","Electron Candidates Class After Corr: 3, Endcap-Showering",100, 0., 1.5);
00186   
00187   TProfile hElectronClassProfileCorr0("hElectronClassProfileCorr0","Electron Candidates ClassProfileCorr: 0, Barrel-Golden",100, 0., 2.7);
00188   TProfile hElectronClassProfileCorr1("hElectronClassProfileCorr1","Electron Candidates ClassProfileCorr: 1, Barrel-Big Brem",100, 0., 2.7);
00189   TProfile hElectronClassProfileCorr2("hElectronClassProfileCorr2","Electron Candidates ClassProfileCorr: 2, Barrel-Narrow",100, 0., 2.7);
00190   TProfile hElectronClassProfileCorr3("hElectronClassProfileCorr3","Electron Candidates ClassProfileCorr: 3, Barrel-Showering",100, 0., 2.7);
00191   TProfile hElectronClassProfileCorr4("hElectronClassProfileCorr4","Electron Candidates ClassProfileCorr: 4, Barrel-Cracks",100, 0., 2.7);
00192   TProfile hElectronClassProfileCorr10("hElectronClassProfileCorr10","Electron Candidates ClassProfileCorr: 0, Endcap-Golden",100, 0., 2.7);
00193   TProfile hElectronClassProfileCorr11("hElectronClassProfileCorr11","Electron Candidates ClassProfileCorr: 1, Endcap-Big Brem",100, 0., 2.7);
00194   TProfile hElectronClassProfileCorr12("hElectronClassProfileCorr12","Electron Candidates ClassProfileCorr: 2, Endcap-Narrow",100, 0., 2.7);
00195   TProfile hElectronClassProfileCorr13("hElectronClassProfileCorr13","Electron Candidates ClassProfileCorr: 3, Endcap-Showering",100, 0., 2.7);
00196   
00197 
00198   // Momentum calculations
00199   TH1F hElectronMomentumEOverP("hElectronMomentumEOverP","E/p",200,0.,2.);
00200   TH1F hElectronMomentumEnergy("hElectronMomentumEnergy","E/Etrue",200,0.,2.);
00201   TH1F hElectronMomentumTrack("hElectronMomentumTrack","p/Etrue",200,0.,2.);
00202   TH1F hElectronMomentumFinal("hElectronMomentumFinal","Efinal/Etrue",200,0.,2.);  
00203   
00204   TH2F hElectronMomentumEnergy2("hElectronMomentumEnergy2","E/Etrue vs E/p",200,0.,2.,200,0.,2.);
00205   TH2F hElectronMomentumTrack2("hElectronMomentumTrack2","p/Etrue vs E/p",200,0.,2.,200,0.,2.5);
00206   TH2F hElectronMomentumFinal2("hElectronMomentumFinal2","Efinal/Etrue vs E/p",200,0.,2.,200,0.,2.);   
00207   
00208   // Higgs analysis   
00209   TH1F hRealZMass("hRealZMass","Higgs: Real Z mass", 150, 0., 150.);
00210   TH1F hVirtualZMass("hVirtualZMass","Higgs: Virtual Z mass", 150, 0., 150.);
00211   TH1F hHiggsMass("hHiggsMass","Higgs: Higgs mass", 180, 0., 180.);   
00212   TH1F hElectronPt("hElectronPt","XANASelectHZZ: Electron pT", 100, 0., 100.);
00213   TH1F hElectronRecGenPtHiggs("hElectronRecGenPtHiggs","ElectronCandidate for Higgs: Rec/Gen pT", 500, 0., 5.);
00214   TH1F hElectronRecGenDrHiggs("hElectronRecGenDrHiggs","ElectronCandidate for Higgs: Rec/Gen Delta R", 500, 0., 5.); 
00215   
00216   TH1F hRealZMassCuts("hRealZMassCuts","Higgs: Real Z mass, after cuts", 150, 0., 150.);
00217   TH1F hVirtualZMassCuts("hVirtualZMassCuts","Higgs: Virtual Z mass, after cuts", 150, 0., 150.);
00218   TH1F hHiggsMassCuts("hHiggsMassCuts","Higgs: Higgs mass, after cuts", 180, 0., 180.);   
00219   TH1F hElectronPtCuts("hElectronPtCuts","XANASelectHZZ: Electron pT, after cuts", 100, 0., 100.);
00220   TH1F hElectronRecGenPtHiggsCuts("hElectronRecGenPtHiggsCuts","ElectronCandidate for Higgs: Rec/Gen pT, after cuts", 500, 0., 5.);
00221   TH1F hElectronRecGenDrHiggsCuts("hElectronRecGenDrHiggsCuts","ElectronCandidate for Higgs: Rec/Gen Delta R, after cuts", 500, 0., 5.); 
00222   
00223   TH1F hElectronRecGenPtHiggsOrig("hElectronRecGenPtHiggsOrig","ElectronCandidate for Higgs: Rec/Gen (Original) pT", 500, 0., 5.);
00224   TH1F hElectronRecGenDrHiggsOrig("hElectronRecGenDrHiggsOrig","ElectronCandidate for Higgs: Rec/Gen (Original) Delta R", 500, 0., 5.); 
00225   TH1F hElectronRecGenPtHiggsOrigAftPh("hElectronRecGenPtHiggsOrigAftPh","ElectronCandidate for Higgs: Rec/Gen (Original) pT, After Photos", 500, 0., 5.);
00226   TH1F hElectronRecGenDrHiggsOrigAftPh("hElectronRecGenDrHiggsOrigAftPh","ElectronCandidate for Higgs: Rec/Gen (Original) Delta R, After Photos", 500, 0., 5.); 
00227 
00228   TH1F hElectronPhotosAftBefPt("hElectronPhotosAftBefPt","Photos: After/Before pT", 200, 0., 2.);
00229   TH1F hElectronPhotosAftBefE("hElectronPhotosAftBefE","Photos: After/Before E", 200, 0., 2.);
00230 
00231   TH1F hRealZMassRecoOrig("hRealZMassRecoOrig","Higgs: Real Z mass (Reconstructed/Original)", 200, 0., 2.);
00232   TH1F hVirtualZMassRecoOrig("hVirtualZMassRecoOrig","Higgs: Virtual Z mass (Reconstructed/Original)", 200, 0., 2.);
00233   TH1F hHiggsMassRecoOrig("hHiggsMassRecoOrig","Higgs: Higgs mass (Reconstructed/Original)", 200, 0., 2.);   
00234   
00235   // Electrons from Higgs
00236   
00237   TH1F hElectronClassesHiggs("hElectronClassesHiggs","Electron Candidates from Higgs - Classification",23, -2., 21.);
00238   
00239   TH1F hElectronClass0Higgs("hElectronClass0Higgs","Electron Candidates from Higgs Class: 0, Barrel-Golden",100, 0., 1.5);
00240   TH1F hElectronClass1Higgs("hElectronClass1Higgs","Electron Candidates from Higgs Class: 1, Barrel-Big Brem",100, 0., 1.5);
00241   TH1F hElectronClass2Higgs("hElectronClass2Higgs","Electron Candidates from Higgs Class: 2, Barrel-Narrow",100, 0., 1.5);
00242   TH1F hElectronClass3Higgs("hElectronClass3Higgs","Electron Candidates from Higgs Class: 3, Barrel-Showering",100, 0., 1.5);
00243   TH1F hElectronClass4Higgs("hElectronClass4Higgs","Electron Candidates from Higgs Class: 4, Barrel-Cracks",100, 0., 1.5);
00244   TH1F hElectronClass10Higgs("hElectronClass10Higgs","Electron Candidates from Higgs Class: 0, Endcap-Golden",100, 0., 1.5);
00245   TH1F hElectronClass11Higgs("hElectronClass11Higgs","Electron Candidates from Higgs Class: 1, Endcap-Big Brem",100, 0., 1.5);
00246   TH1F hElectronClass12Higgs("hElectronClass12Higgs","Electron Candidates from Higgs Class: 2, Endcap-Narrow",100, 0., 1.5);
00247   TH1F hElectronClass13Higgs("hElectronClass13Higgs","Electron Candidates from Higgs Class: 3, Endcap-Showering",100, 0., 1.5);
00248   
00249   TProfile hElectronClassProfile0Higgs("hElectronClassProfile0Higgs","Electron Candidates from Higgs ClassProfile: 0, Barrel-Golden",100, 0., 2.7);
00250   TProfile hElectronClassProfile1Higgs("hElectronClassProfile1Higgs","Electron Candidates from Higgs ClassProfile: 1, Barrel-Big Brem",100, 0., 2.7);
00251   TProfile hElectronClassProfile2Higgs("hElectronClassProfile2Higgs","Electron Candidates from Higgs ClassProfile: 2, Barrel-Narrow",100, 0., 2.7);
00252   TProfile hElectronClassProfile3Higgs("hElectronClassProfile3Higgs","Electron Candidates from Higgs ClassProfile: 3, Barrel-Showering",100, 0., 2.7);
00253   TProfile hElectronClassProfile4Higgs("hElectronClassProfile4Higgs","Electron Candidates from Higgs ClassProfile: 4, Barrel-Cracks",100, 0., 2.7);
00254   TProfile hElectronClassProfile10Higgs("hElectronClassProfile10Higgs","Electron Candidates from Higgs ClassProfile: 0, Endcap-Golden",100, 0., 2.7);
00255   TProfile hElectronClassProfile11Higgs("hElectronClassProfile11Higgs","Electron Candidates from Higgs ClassProfile: 1, Endcap-Big Brem",100, 0., 2.7);
00256   TProfile hElectronClassProfile12Higgs("hElectronClassProfile12Higgs","Electron Candidates from Higgs ClassProfile: 2, Endcap-Narrow",100, 0., 2.7);
00257   TProfile hElectronClassProfile13Higgs("hElectronClassProfile13Higgs","Electron Candidates from Higgs ClassProfile: 3, Endcap-Showering",100, 0., 2.7);
00258   
00259   TH1F hElectronEtaHiggs("hElectronEtaHiggs","Electron Candidates from Higgs Eta",100, 0., 2.7);
00260   TH1F hElectronEtaClass0Higgs("hElectronEtaClass0Higgs","Electron Candidates from Higgs Eta Class: 0, Barrel-Golden",100, 0., 2.7);
00261   TH1F hElectronEtaClass1Higgs("hElectronEtaClass1Higgs","Electron Candidates from Higgs Eta Class: 1, Barrel-Big Brem",100, 0., 2.7);
00262   TH1F hElectronEtaClass2Higgs("hElectronEtaClass2Higgs","Electron Candidates from Higgs Eta Class: 2, Barrel-Narrow",100, 0., 2.7);
00263   TH1F hElectronEtaClass3Higgs("hElectronEtaClass3Higgs","Electron Candidates from Higgs Eta Class: 3, Barrel-Showering",100, 0., 2.7);
00264   TH1F hElectronEtaClass4Higgs("hElectronEtaClass4Higgs","Electron Candidates from Higgs Eta Class: 4, Barrel-Cracks",100, 0., 2.7);
00265   TH1F hElectronEtaClass10Higgs("hElectronEtaClass10Higgs","Electron Candidates from Higgs Eta Class: 0, Endcap-Golden",100, 0., 2.7);
00266   TH1F hElectronEtaClass11Higgs("hElectronEtaClass11Higgs","Electron Candidates from Higgs Eta Class: 1, Endcap-Big Brem",100, 0., 2.7);
00267   TH1F hElectronEtaClass12Higgs("hElectronEtaClass12Higgs","Electron Candidates from Higgs Eta Class: 2, Endcap-Narrow",100, 0., 2.7);
00268   TH1F hElectronEtaClass13Higgs("hElectronEtaClass13Higgs","Electron Candidates from Higgs Eta Class: 3, Endcap-Showering",100, 0., 2.7);
00269   
00270   TH1F hElectronMomentumEOverPHiggs("hElectronMomentumEOverPHiggs","Electron Candidates from Higgs - E/p",200,0.,2.);
00271   TH1F hElectronMomentumEnergyHiggs("hElectronMomentumEnergyHiggs","Electron Candidates from Higgs - E/Etrue",200,0.,2.);
00272   TH1F hElectronMomentumTrackHiggs("hElectronMomentumTrackHiggs","Electron Candidates from Higgs - p/Etrue",200,0.,2.);
00273   TH1F hElectronMomentumFinalHiggs("hElectronMomentumFinalHiggs","Electron Candidates from Higgs - Efinal/Etrue",200,0.,2.);  
00274   
00275   TH2F hElectronMomentumEnergy2Higgs("hElectronMomentumEnergy2Higgs","Electron Candidates from Higgs - E/Etrue vs E/p",200,0.,2.,200,0.,2.);
00276   TH2F hElectronMomentumTrack2Higgs("hElectronMomentumTrack2Higgs","Electron Candidates from Higgs - p/Etrue vs E/p",200,0.,2.,200,0.,2.5);
00277   TH2F hElectronMomentumFinal2Higgs("hElectronMomentumFinal2Higgs","Electron Candidates from Higgs - Efinal/Etrue vs E/p",200,0.,2.,200,0.,2.);
00278     
00279   // Initialisation of some variables (mainly counters and temporary)  
00280   Int_t nb = 0; 
00281   Int_t nevent = 0, ievent = 0;
00282   int ne = 0, nTrigAccL1=0, nTrigAccHLT=0, nTrigAccHLTelec=0,nPresAcc = 0, nIsoAcc = 0, nAmbBef = 0, nAmbAft = 0, n2e2pAcc = 0;
00283   int nCutsIso = 0, nCutsPt = 0, nCutsRealZMass = 0, nCutsVirtualZMass = 0, nCutsHiggsMass = 0;
00284   int mc4e = 0, mc4t = 0, mc2e2t = 0, mc4eCuts = 0, mc4tCuts = 0, mc2e2tCuts = 0;
00285   int nconsidered=0,niso=0;
00286   TObjArray *electrons = 0, *superClusters = 0; 
00287   TObjArray *tracks=0;
00288   HepLorentzVector realZ, virtualZ, higgs, realZOrig, virtualZOrig, higgsOrig;
00289   HepLorentzVector tmpElRec, tmpElGen, tmpElGenOrigB, tmpElGenOrigA;
00290   XANAElectronCandidate *tmpElectronCandidate;
00291   
00292   // define XANA event
00293   XANAEsdEvent *event = new XANAEsdEvent();  
00294   XANAGeneratorEvent *generatorEvent = new XANAGeneratorEvent();
00295   XANAGeantEvent *geantEvent = new XANAGeantEvent();
00296   
00297   // version TChain
00298   Short_t numberOfFiles = 0;
00299   TChain *T=new TChain("ESD");
00300   char tmp[512];
00301   
00302   // files to load, dependant on the key word given at runtime
00303   // to be completed by other datasets
00304   if (strcmp(argv[1], "hzz150") == 0) {
00305     cout<<"hzz4e 150"<<endl;
00306     numberOfFiles = 18;
00307     for (Int_t m_fi = 0; m_fi < numberOfFiles ; m_fi++) {
00308       sprintf(tmp,"/sps/cms/baffioni/XanadooTrees/hg03_hzz_4e_150_XANA_153_pix/hg03_hzz_4e_150.root.%03d",m_fi);
00309       TFile F(tmp);
00310       if (!(F.GetNkeys())) { continue; }
00311       T->Add(tmp);
00312     }   
00313   } 
00314   else if (strcmp(argv[1], "hzz120") == 0) {
00315     cout<<"hzz4e 120"<<endl;
00316     numberOfFiles = 20;
00317     for (Int_t m_fi = 0; m_fi < numberOfFiles ; m_fi++) {
00318       sprintf(tmp,"/sps/cms/baffioni/XanadooTrees/hg03_hzz_4e_120_XANA_153_pix/hg03_hzz_4e_120.root.%03d",m_fi);
00319       TFile F(tmp);
00320       if (!(F.GetNkeys())) { continue; }
00321       T->Add(tmp);
00322     }   
00323   } 
00324   else if(strcmp(argv[1], "tt4e") == 0){       
00325     cout<<"ttbar->4e"<<endl;
00326     numberOfFiles = 160;
00327     for (Int_t m_fi = 0; m_fi < numberOfFiles ; m_fi++) {
00328       if (m_fi == 59 || m_fi == 106 ||m_fi == 158 || m_fi == 159 || m_fi == 149 || m_fi == 19|| m_fi == 31) {continue;}
00329       sprintf(tmp,"/sps/cms/baffioni/XanadooTrees/eg03_tt_4e_XANA_153_pix/eg03_tt_4e.root.%03d",m_fi); 
00330       TFile F(tmp);
00331       if (!(F.GetNkeys())) {continue; }
00332       T->Add(tmp);      
00333     }
00334   } 
00335   else{
00336     cout << "no such data yet"<<endl;
00337     return 0;
00338   }
00339 
00340   // init settings
00341   initSettings(T, &event, &generatorEvent, &geantEvent);
00342 
00343   // end of initialization phase
00344   
00345   //************************************************************************************************************
00346   
00347   //   Start main loop on all events
00348   nevent = (Int_t)T->GetEntries();
00349   cout << "================== Total number of events: " << nevent << " ================" << endl;  
00350   
00351   // for (Int_t i=0; i<T->GetEntries(); i++) {
00352   for (Int_t i=0; i<100; i++) {
00353     nb += T->GetEntry(i); 
00354     ievent++;
00355     
00356     if (i%100 == 0) cout << "================ Event: " << i << " ================" << endl;    
00357     //if (i%1 == 0) cout << "================ Event: " << i << " ================" << endl;    
00358 
00359     ne  = event->getNumberOfElectronCandidates(); electrons = event->getElectronCandidates();
00360     tracks= event->getTracks();
00361     hElectronNumber.Fill((Axis_t) ne);    
00362     
00363     /*
00364     // cross check
00365     for (int i3 = 0; i3 < ne; i3++) {
00366     tmpElectronCandidate = (XANAElectronCandidate *) electrons->At(i3);
00367     cout << "Super Cluster Pointer = " << tmpElectronCandidate->getSeedCluster() << endl;
00368     }
00369     */
00370     
00371     //**********************************************************************************************************
00372 
00373     // FIXME: problem with Super Clusters in endcap
00374     // this part of code should do nothing with data from XANADOO_1_5_3
00375     // to be deleted as soon as we don't use old data anymore
00376     
00377     //cout<<"fix superclusters"<<endl;
00378     superClusters = event->getSuperClusters();
00379 
00380     for (int i3 = 0; i3 < ne; i3++) {
00381       tmpElectronCandidate = (XANAElectronCandidate *) electrons->At(i3);
00382       
00383       // cout << "Super Cluster Pointer = " << tmpElectronCandidate->getSuperCluster() << endl;
00384       for (int isc = 0 ; isc < superClusters->GetLast() + 1; isc ++) {
00385         XANASuperCluster * sc = (XANASuperCluster *) superClusters->At(isc);
00386         if (tmpElectronCandidate->getSuperCluster()==0 &&
00387             tmpElectronCandidate->getSuperClusterPosition() == sc->getPosition()) {
00388           cout<<" !!! entering the super cluster bug correction "<<endl;
00389           cout<<" !!! should not be the case if data >= XANA1_5_3 are used"<<endl;
00390           // cout << " SC: Cand en, SC en, SS pos, brems: " << tmpElectronCandidate->getSuperClusterEnergy() << ", " <<
00391           // sc->getEnergy() << ", " << sc->getPosition().eta() << " " 
00392           // << sc->getNumberOfBrems() << endl;
00393           float oldEnergy = tmpElectronCandidate->getSuperClusterEnergy();
00394           tmpElectronCandidate->setSuperClusterEnergy(sc->getEnergy());
00395           //cout << "new supercluster endcap: old energy " << oldEnergy << " new energy " << sc->getEnergy() << endl;
00396           tmpElectronCandidate->setSuperCluster(sc);                    
00397           //CC also update MomentumAtVertex 
00398           Hep3Vector ptrack;
00399           if (tmpElectronCandidate->getElectronTrack()->getAlgoName() == "GSF") {
00400             XANAElectronGSFTrack *gsfetr = (XANAElectronGSFTrack *)(tmpElectronCandidate->getElectronTrack());
00401             ptrack = gsfetr->getGSFMomentumAtVertexMode();                                                
00402           } else {
00403             ptrack = tmpElectronCandidate->getElectronTrack()->getMomentumAtVertex();
00404           }
00405           HepLorentzVector momentum(ptrack.x(), ptrack.y(), ptrack.z(), 0.);
00406           momentum *= sc->getEnergy() / ptrack.mag();
00407           momentum.setE(sc->getEnergy());
00408           //cout << "new electron endcap: old momentum " << tmpElectronCandidate->getMomentumAtVertex() << " new momentum " << momentum << endl;
00409           tmpElectronCandidate->setMomentumAtVertex(momentum);
00410           //CC also update Eoverp
00411           //cout << "new electron endcap: old E/p " << tmpElectronCandidate->getESuperClusterOverP() << " new E/p " << sc->getEnergy() / ptrack.mag() << endl;
00412           tmpElectronCandidate->setMomentumAtVertex(momentum);
00413           tmpElectronCandidate->setESuperClusterOverP(sc->getEnergy()/ptrack.mag());
00414           //CC also update had/E
00415           //cout << "new electron endcap: old h/e " << tmpElectronCandidate->getHadronicOverEm() << " new h/e " << sc->getHadronicOverEm()*oldEnergy/sc->getEnergy() << endl;
00416           tmpElectronCandidate->setHadronicOverEm(sc->getHadronicOverEm()*oldEnergy/sc->getEnergy());
00417         }                       
00418       }
00419     }
00420     
00421     
00422     //****************************************************************************************
00423     // part of code to unescale electron candidates
00424     //  CC in case we do want to remove all corrections
00425 
00426     if (!egammaEscaleCorrection) {    
00427       //cout<< "unEscale electron candidates"<<endl;    
00428      for (int i3 = 0; i3 < ne; i3++) {
00429         tmpElectronCandidate = (XANAElectronCandidate *) electrons->At(i3);
00430         // from federico
00431         float newEnergy = tmpElectronCandidate->getSuperCluster()->getEnergy()/
00432           tmpElectronCandidate->getSuperCluster()->getEnergyScaleFactor();
00433         HepLorentzVector oldMomentum = tmpElectronCandidate->getMomentumAtVertex();
00434         HepLorentzVector newMomentum(oldMomentum.x()*newEnergy/oldMomentum.t(),
00435                                      oldMomentum.y()*newEnergy/oldMomentum.t(),
00436                                      oldMomentum.z()*newEnergy/oldMomentum.t(),
00437                                      newEnergy);
00438         tmpElectronCandidate->setMomentumAtVertex(newMomentum);         
00439       }
00440     }
00441   
00442   
00443     
00444     //**********************************************************************************
00445 
00446     // select only 4e decay for signal
00447     // -------------------------------------------------------
00448 
00449     //cout<<"XANASelectMc"<<endl;
00450 
00451     bool  issignal=false;
00452     XANASelectMcHZZ selectorMc(generatorEvent);
00453 
00454     //only for signal be careful to add a new signal if there's one
00455    if (strcmp(argv[1], "hzz150") == 0 || strcmp(argv[1], "hzz120") == 0) {
00456       issignal=true;
00457       // Analysing H-->ZZ decay tree
00458       if (selectorMc.isFourElectrons()) mc4e++;
00459       if (selectorMc.isFourTaus()) mc4t++;
00460       if (selectorMc.isTwoElectronsTwoTaus()) mc2e2t++;
00461       
00462       if (!selectorMc.isFourElectrons()) {clean(event, generatorEvent, geantEvent);continue;}    
00463     }
00464     nconsidered++;
00465 
00466     //**********************************************************************************
00467 
00469     // Trigger //
00471 
00472     // cout<<"XANAnalysisTrigger"<<endl;
00473     XANATriggerInfo triggerInfo = event->getTriggerInfo();
00474     XANAnalysisTrigger trigger(&triggerInfo); 
00475     
00476     if (trigger.getL1Decision() == true) nTrigAccL1++;
00477     else {clean(event, generatorEvent, geantEvent);continue; }  
00478 
00479     if (trigger.getHLTDecision() == true) nTrigAccHLT++;
00480     else {clean(event, generatorEvent, geantEvent);continue; }  
00481     
00482     if (trigger.getHLTElecDecision() == true) nTrigAccHLTelec++;
00483     else {clean(event, generatorEvent, geantEvent);continue; }     
00484     
00485     //**********************************************************************************
00486     
00488     // Preselection //
00490     
00491     // 1. Simple Cuts
00492     // --------------
00493     
00494     // Ask for 2 electrons and 2 positrons all with pT > 3 GeV with no isolation requirement    
00495     //cout<<"XANAEventSelec"<<endl;
00496     XANAEventSelect eventSelect(event);
00497     vector<float> ptCut;
00498     ptCut.push_back(3.); 
00499     ptCut.push_back(3.); 
00500     if ((eventSelect.NumberOfElCandSelect( 2, -1, 0, ptCut)) && 
00501         (eventSelect.NumberOfElCandSelect( 2, 1, 0, ptCut))) {
00502       nPresAcc++;
00503     }  else {clean(event, generatorEvent, geantEvent);continue;}
00504     
00505     
00506     // 2. Loose isolation
00507     // ------------------
00508     
00509     // Calculation isolation on all Electron candidate in the event   
00510     //cout<<"XANAEvtIsolation"<<endl;
00511     XANAEvtIsolation isolation(event);
00512     std::vector<XANAIsolation> isolationVector = isolation.getIsolationVector();
00513     std::vector<XANAIsolation>::const_iterator isoIter;
00514     
00515     // Fill some histograms
00516     for (isoIter = isolationVector.begin(); isoIter != isolationVector.end(); isoIter++){
00517       hIsoNumberTracks.Fill((Axis_t) isoIter->getNumberTracks());
00518       hIsoPtTracks.Fill((Axis_t) isoIter->getPtTracks());
00519       hIsoNumberECinsideCone.Fill((Axis_t) isoIter->getNumberECinsideCone());
00520       hIsoPtMax.Fill((Axis_t) isoIter->getPtMax());
00521     }  
00522     
00523 
00524     // Very loose isolation (CASE 1): Ask for four isolated electrons with pTsum < 10 GeV inside a cone 
00525     if (isolation.getIsolationPtSum(10.,4)) nIsoAcc++;  
00526     else {clean(event, generatorEvent, geantEvent);continue;}
00527     
00528     // Very loose isolation (CASE 2): Ask for four isolated electrons with pTmax < 5 GeV inside a cone 
00529     //if (isolation.getIsolationPtMax(5.,4)) nIsoAcc++;    
00530     
00531 
00532 
00533     //**********************************************************************************
00534 
00536     // Electron reconstruction //
00538 
00539     // Associate reconstructed and generated electrons and make some analysis
00540     //cout<<"XANAElectronRecGen"<<endl;
00541     XANAElectronRecGen anaMcEvent(event,generatorEvent);
00542     map<XANAElectronCandidate *, XANAGeneratorParticle *> electronsRecGen = anaMcEvent.getElectronsRecGen();
00543     map<XANAElectronCandidate *, XANAGeneratorParticle *>::const_iterator iterEleRecGen;
00544     
00545     for (iterEleRecGen = electronsRecGen.begin(); iterEleRecGen != electronsRecGen.end(); iterEleRecGen++){
00546       tmpElRec = (iterEleRecGen->first)->getMomentumAtVertex();
00547       tmpElGen = (iterEleRecGen->second)->getMomentum();
00548       hElectronRecGenPt.Fill(tmpElRec.perp()/tmpElGen.perp());
00549       hElectronRecGenDr.Fill(tmpElRec.deltaR(tmpElGen)); 
00550     }
00551     
00552 
00553     //**********************************************************************************
00554 
00556     // Electron preselection //
00558     
00559     //cout<<"XANAElectronSelect"<<endl;
00560     XANAElectronSelect electronSelector(electrons);
00561     electronSelector.ElSelect_dphi(0.1);
00562     TObjArray *selectedElectrons = electronSelector.getSelectedElectrons();
00563     
00564     
00565 
00567     // Ambiguity resolving //
00569     
00570     // single resolving    
00571     //cout<<"XANAmbiguityResolve"<<endl;
00572     XANAmbiguityResolve ambResolve(selectedElectrons,selectedElectrons->GetLast() + 1);
00573     ambResolve.ResolveByEoverP();
00574     
00575     // double resolving
00576 
00577     //XANAmbiguityResolve ambResolve1(selectedElectrons,selectedElectrons->GetLast() + 1);
00578     //ambResolve1.ResolveByEoverP();
00579         
00580     //  XANAmbiguityResolve ambResolve(ambResolve1.getResolvedElcands(),ambResolve1.getResolvedElcands()->GetLast() + 1);
00581     // ambResolve.ResolveByEoverP();
00582     
00583     // Prepare for calculation of average number of Electron Candidates before and after ambiguity resolving
00584     nAmbBef+=ne;
00585     nAmbAft+=ambResolve.getResolvedElcands()->GetLast()+1;
00586     // Fill some histograms
00587     hElectronNumberAfterAmb.Fill((Axis_t) ambResolve.getResolvedElcands()->GetLast() + 1);
00588     hElectronNumberBeforeAfter.Fill((Axis_t) ne - ambResolve.getResolvedElcands()->GetLast() - 1);
00589     
00590 
00591     //**********************************************************************************
00592    
00594     // Electron classification //
00596     
00597     //cout<<" XANAElectronClassification"<<endl;
00598     XANAElectronClassification electronClasses(ambResolve.getResolvedElcands()); 
00599     map<XANAElectronCandidate *, int>  elClassMap = electronClasses.getElectronClass();
00600     map<XANAElectronCandidate *, int>::const_iterator elClassMapIter;
00601     float uncorrEnergy = 0;
00602     
00603     // some checks
00604     // cout << "------- Electron classification ------- " << ambResolve.getResolvedElcands()->GetLast() + 1 <<  endl;
00605     for (elClassMapIter = elClassMap.begin(); elClassMapIter != elClassMap.end(); elClassMapIter++) {
00606       // cout << "El. candidate pointer and class : " << elClassMapIter->first << " " << elClassMapIter->second << endl;   
00607       hElectronClasses.Fill((Axis_t) elClassMapIter->second);
00608       tmpElectronCandidate = elClassMapIter->first;
00609       float myEta = tmpElectronCandidate->getSuperClusterPosition().eta();
00610       uncorrEnergy = tmpElectronCandidate->getSuperClusterEnergy();
00611       if (elClassMapIter->second >= 0 && elClassMapIter->second < 5 && tmpElectronCandidate->getSuperCluster()) 
00612         uncorrEnergy = uncorrEnergy/tmpElectronCandidate->getSuperCluster()->getEnergyScaleFactor();
00613       tmpElGen = ((electronsRecGen.find(tmpElectronCandidate))->second)->getMomentum();
00614       Axis_t ratio = uncorrEnergy/tmpElGen.e(); 
00615       hElectronEta.Fill(myEta);
00616       if (elClassMapIter->second == 0) {
00617         hElectronClass0.Fill(ratio); if (ratio > 0.8 && ratio < 2.) hElectronClassProfile0.Fill(myEta, ratio);
00618         hElectronEtaClass0.Fill(myEta);
00619       }
00620       else if (elClassMapIter->second == 1) {
00621         hElectronClass1.Fill(ratio); if (ratio > 0.8 && ratio  < 2.) hElectronClassProfile1.Fill(myEta, ratio);
00622         hElectronEtaClass1.Fill(myEta);
00623       }
00624       else if (elClassMapIter->second == 2) {
00625         hElectronClass2.Fill(ratio); if (ratio > 0.8 && ratio < 2.) hElectronClassProfile2.Fill(myEta, ratio);
00626         hElectronEtaClass2.Fill(myEta);
00627       }
00628       else if (elClassMapIter->second == 3)  {
00629         hElectronClass3.Fill(ratio); if (ratio > 0.8 && ratio < 2.) hElectronClassProfile3.Fill(myEta, ratio);
00630         hElectronEtaClass3.Fill(myEta);
00631       }
00632       else if (elClassMapIter->second == 4)  {
00633         hElectronClass4.Fill(ratio); if (ratio > 0.8 && ratio < 2.) hElectronClassProfile4.Fill(myEta, ratio);
00634         hElectronEtaClass4.Fill(myEta);
00635       }
00636       else if (elClassMapIter->second == 10) {
00637         hElectronClass10.Fill(ratio); if (ratio > 0.8 && ratio < 2.) hElectronClassProfile10.Fill(myEta, ratio);
00638         hElectronEtaClass10.Fill(myEta);
00639       }
00640       else if (elClassMapIter->second == 11) {
00641         hElectronClass11.Fill(ratio); if (ratio > 0.8 && ratio < 2.) hElectronClassProfile11.Fill(myEta, ratio);
00642         hElectronEtaClass11.Fill(myEta);
00643       }
00644       else if (elClassMapIter->second == 12) {
00645         hElectronClass12.Fill(ratio); if (ratio > 0.8 && ratio < 2.) hElectronClassProfile12.Fill(myEta, ratio);
00646         hElectronEtaClass12.Fill(myEta);
00647       }
00648       else if (elClassMapIter->second == 13) {
00649         hElectronClass13.Fill(ratio); if (ratio > 0.8 && ratio < 2.)hElectronClassProfile13.Fill(myEta, ratio);
00650         hElectronEtaClass13.Fill(myEta);
00651       }
00652     }
00653     
00654     
00655     //**********************************************************************************
00656     
00658     // Electron corrections //
00660     
00661     float corrEnergy;  
00662     
00663     if (federicoEscaleCorrection) {
00664       
00665       XANAElectronCorrection electronCorrection(elClassMap);
00666       corrEnergy = 0.;  
00667       
00668       // cross checks
00669       
00670       for (elClassMapIter = elClassMap.begin(); elClassMapIter != elClassMap.end(); elClassMapIter++) {
00671         tmpElectronCandidate = elClassMapIter->first;
00672         float myEta = tmpElectronCandidate->getSuperClusterPosition().eta();
00673         tmpElGen = ((electronsRecGen.find(tmpElectronCandidate))->second)->getMomentum();
00674         corrEnergy = tmpElectronCandidate->getSuperClusterEnergy();
00675         Axis_t ratio = corrEnergy/tmpElGen.e(); 
00676         if (elClassMapIter->second == 0)      { hElectronClassCorr0.Fill(ratio); if (ratio > 0.8 && ratio < 2.) hElectronClassProfileCorr0.Fill(myEta, ratio); }
00677         else if (elClassMapIter->second == 1)  { hElectronClassCorr1.Fill(ratio); if (ratio > 0.8 && ratio < 2.) hElectronClassProfileCorr1.Fill(myEta, ratio); }
00678         else if (elClassMapIter->second == 2)  { hElectronClassCorr2.Fill(ratio); if (ratio > 0.8 && ratio < 2.) hElectronClassProfileCorr2.Fill(myEta, ratio); }
00679         else if (elClassMapIter->second == 3)  { hElectronClassCorr3.Fill(ratio); if (ratio > 0.8 && ratio < 2.) hElectronClassProfileCorr3.Fill(myEta, ratio); }
00680         else if (elClassMapIter->second == 4)  { hElectronClassCorr4.Fill(ratio); if (ratio > 0.8 && ratio < 2.) hElectronClassProfileCorr4.Fill(myEta, ratio); }
00681         else if (elClassMapIter->second == 10) { hElectronClassCorr10.Fill(ratio); if (ratio > 0.8 && ratio < 2.) hElectronClassProfileCorr10.Fill(myEta, ratio); }
00682         else if (elClassMapIter->second == 11) { hElectronClassCorr11.Fill(ratio); if (ratio > 0.8 && ratio < 2.) hElectronClassProfileCorr11.Fill(myEta, ratio); }
00683         else if (elClassMapIter->second == 12) { hElectronClassCorr12.Fill(ratio); if (ratio > 0.8 && ratio < 2.) hElectronClassProfileCorr12.Fill(myEta, ratio); }
00684         else if (elClassMapIter->second == 13) { hElectronClassCorr13.Fill(ratio); if (ratio > 0.8 && ratio < 2.) hElectronClassProfileCorr13.Fill(myEta, ratio); }
00685       }   
00686       
00687     }
00688     
00689     //**********************************************************************************
00690       
00692     // Electron momentum calculation //
00694     if (momentumCorrection) {
00695 
00696       XANAElectronMomentum electronMomentum(elClassMap);
00697       
00698       // cross checks
00699       for (elClassMapIter = elClassMap.begin(); elClassMapIter != elClassMap.end(); elClassMapIter++) {
00700         tmpElectronCandidate = elClassMapIter->first;    
00701         XANAElectronGSFTrack *gsfTrack = (XANAElectronGSFTrack *) tmpElectronCandidate->getElectronTrack();
00702         HepVector3D tmpTrackRec = gsfTrack->getGSFMomentumAtVertexMode();       
00703         tmpElGen = ((electronsRecGen.find(tmpElectronCandidate))->second)->getMomentum();
00704         Axis_t eOverP = tmpElectronCandidate->getSuperClusterEnergy()/tmpTrackRec.mag();
00705         Axis_t eOverETrue = tmpElectronCandidate->getSuperClusterEnergy()/tmpElGen.t();
00706         Axis_t pOverETrue = tmpTrackRec.mag()/tmpElGen.t();
00707         HepLorentzVector final = tmpElectronCandidate->getMomentumAtVertex();
00708         Axis_t finalOverETrue = final.t()/tmpElGen.t();
00709         
00710         hElectronMomentumEOverP.Fill(eOverP);
00711         hElectronMomentumEnergy.Fill(eOverETrue);    hElectronMomentumEnergy2.Fill(eOverP, eOverETrue);
00712         hElectronMomentumTrack.Fill(pOverETrue);     hElectronMomentumTrack2.Fill(eOverP, pOverETrue);
00713         hElectronMomentumFinal.Fill(finalOverETrue); hElectronMomentumFinal2.Fill(eOverP, finalOverETrue);       
00714       }     
00715       
00716     }
00717                         
00718     //********************************************************************
00719 
00721     // Higgs Analysis //
00723 
00724     // Isolation 
00725     // select events with >=2e+ and >=2e- isolated electrons with Pt>5GeV
00726     
00727     
00728     XANAEventSelect eventSelect2(ambResolve.getResolvedElcands(),tracks);
00729     vector<float> ptCut2;
00730     ptCut2.push_back(5.); 
00731     ptCut2.push_back(5.); 
00732 
00733     if ((eventSelect2.NumberOfElCandSelect( 2, -1, 1, ptCut2)) && 
00734         (eventSelect2.NumberOfElCandSelect( 2, 1, 1, ptCut2))) {
00735       niso++;
00736     }  else {clean(event, generatorEvent, geantEvent);continue;}
00737     TObjArray  * isoel=eventSelect2.getGoodElCands(0,1,5);
00738 
00739  
00740     // Reconstruct Z, Z*, Higgs (using unambigous isolated electrons only)
00741     // ----------------------------------------------------------
00742     // XANASelectHZZ selector(electrons,ne);
00743     XANASelectHZZ selector(isoel,isoel->GetLast()+1);  
00744     
00745     // internal consisteny, asks for 2e- and 2e+
00746     if (selector.isValid()) n2e2pAcc++;
00747     else {clean(event, generatorEvent, geantEvent);continue; }
00748     
00749     // Fill temporary array for easier navigation
00750     XANAElectronCandidate *electronsHiggs[4] = {selector.getElectronFromRealZ(), 
00751                                                 selector.getPositronFromRealZ(), 
00752                                                 selector.getElectronFromVirtualZ(), 
00753                                                 selector.getPositronFromVirtualZ()};    
00754     
00755     
00756     
00757     realZ = selector.getRealZ(); hRealZMass.Fill(realZ.m());
00758     virtualZ = selector.getVirtualZ(); hVirtualZMass.Fill(virtualZ.m());
00759     higgs = selector.getHiggs(); hHiggsMass.Fill(higgs.m());
00760     
00761 
00762 
00763     // Fill temporary array for easier navigation
00764     XANAGeneratorParticle *electronsHiggsOrig[4];
00765     XANAGeneratorParticle *electronsHiggsOrigAftPh[4];
00766     if(issignal){
00767       // original electrons before PHOTOS
00768 
00769       electronsHiggsOrig [0]= selectorMc.getElectronFromRealZ();
00770       electronsHiggsOrig [1]= selectorMc.getPositronFromRealZ();
00771       electronsHiggsOrig [2]= selectorMc.getElectronFromVirtualZ();
00772       electronsHiggsOrig [3]= selectorMc.getPositronFromVirtualZ() ;
00773       
00774       // original electrons after PHOTOS
00775       electronsHiggsOrigAftPh[0] =  selectorMc.getElectronFromRealZAftPh();
00776       electronsHiggsOrigAftPh[1] =  selectorMc.getPositronFromRealZAftPh();
00777       electronsHiggsOrigAftPh[2] =  selectorMc.getElectronFromVirtualZAftPh();
00778       electronsHiggsOrigAftPh[3] =  selectorMc.getPositronFromVirtualZAftPh() ;
00779       
00780       // Some analysis
00781       realZOrig = selectorMc.getRealZ()->getMomentum(); hRealZMassRecoOrig.Fill(realZ.m()/realZOrig.m());
00782       virtualZOrig = selectorMc.getVirtualZ()->getMomentum(); hVirtualZMassRecoOrig.Fill(virtualZ.m()/virtualZOrig.m());
00783       higgsOrig = selectorMc.getHiggs()->getMomentum(); hHiggsMassRecoOrig.Fill(higgs.m()/higgsOrig.m());
00784     }
00785     
00786 
00787     // Electrons from Higgs - Analysis
00788     // Reconstructed versus generated (for Higgs electrons) and reconstructed versus original electrons
00789     for (int ieh = 0; ieh < 4; ieh++){
00790       hElectronPt.Fill(electronsHiggs[ieh]->getPt());
00791       tmpElRec = electronsHiggs[ieh]->getMomentumAtVertex();
00792       tmpElGen = ((electronsRecGen.find(electronsHiggs[ieh]))->second)->getMomentum();
00793       hElectronRecGenPtHiggs.Fill(tmpElRec.perp()/tmpElGen.perp()); 
00794       hElectronRecGenDrHiggs.Fill(tmpElRec.deltaR(tmpElGen));
00795 
00796       if (issignal && selectorMc.isFourElectrons()) {
00797         tmpElGenOrigB = electronsHiggsOrig[ieh]->getMomentum();
00798         hElectronRecGenPtHiggsOrig.Fill(tmpElRec.perp()/tmpElGenOrigB.perp());
00799         hElectronRecGenDrHiggsOrig.Fill(tmpElRec.deltaR(tmpElGenOrigB));
00800         tmpElGenOrigA = electronsHiggsOrigAftPh[ieh]->getMomentum();              
00801 
00802         hElectronRecGenPtHiggsOrigAftPh.Fill(tmpElRec.perp()/tmpElGenOrigA.perp());
00803         hElectronRecGenDrHiggsOrigAftPh.Fill(tmpElRec.deltaR(tmpElGenOrigA));
00804         hElectronPhotosAftBefPt.Fill(tmpElGenOrigA.perp()/tmpElGenOrigB.perp());
00805         hElectronPhotosAftBefE.Fill(tmpElGenOrigA.t()/tmpElGenOrigB.t());
00806         // if (tmpElGenOrigA.perp()/tmpElGenOrigB.perp()> 1)
00807         //         cout << tmpElGenOrigA.perp()/tmpElGenOrigB.perp() << " | " << tmpElGenOrigA.t()/tmpElGenOrigB.t() << endl;
00808       }
00809       
00810 
00811       // find to which classes electrons belongs
00812       int elHiggsClass = (elClassMap.find(electronsHiggs[ieh]))->second;
00813       hElectronClassesHiggs.Fill((Axis_t) elHiggsClass);
00814       
00815       // use the same analysis as above
00816       tmpElectronCandidate = electronsHiggs[ieh];        
00817       XANAElectronGSFTrack *gsfTrack = (XANAElectronGSFTrack *) tmpElectronCandidate->getElectronTrack();
00818       HepVector3D tmpTrackRec = gsfTrack->getGSFMomentumAtVertexMode();
00819       Axis_t eOverP = tmpElectronCandidate->getSuperClusterEnergy()/tmpTrackRec.mag();
00820       Axis_t eOverETrue = tmpElectronCandidate->getSuperClusterEnergy()/tmpElGen.t();
00821       Axis_t pOverETrue = tmpTrackRec.mag()/tmpElGen.t();
00822       HepLorentzVector final = tmpElectronCandidate->getMomentumAtVertex();
00823       //std::cout << "tmpElGen " << tmpElGen << std::endl;
00824       //std::cout << "tmpElRec " << tmpElRec << std::endl;
00825       //std::cout << "final    " << final << std::endl;
00826       Axis_t finalOverETrue = final.t()/tmpElGen.t();
00827       
00828       hElectronMomentumEOverPHiggs.Fill(eOverP);
00829       hElectronMomentumEnergyHiggs.Fill(eOverETrue);    hElectronMomentumEnergy2Higgs.Fill(eOverP, eOverETrue);
00830       hElectronMomentumTrackHiggs.Fill(pOverETrue);     hElectronMomentumTrack2Higgs.Fill(eOverP, pOverETrue);
00831       hElectronMomentumFinalHiggs.Fill(finalOverETrue); hElectronMomentumFinal2Higgs.Fill(eOverP, finalOverETrue);       
00832       
00833       // Classes analysis
00834       int electronClass = (elClassMap.find(electronsHiggs[ieh]))->second;
00835       float myEta = electronsHiggs[ieh]->getSuperClusterPosition().eta();
00836       tmpElGen = ((electronsRecGen.find(electronsHiggs[ieh]))->second)->getMomentum();
00837       corrEnergy = electronsHiggs[ieh]->getSuperClusterEnergy();
00838       Axis_t ratio = corrEnergy/tmpElGen.e(); 
00839       hElectronEtaHiggs.Fill(myEta);
00840       if (electronClass == 0) { 
00841         hElectronClass0Higgs.Fill(ratio); if (ratio > 0.8 && ratio < 2.) hElectronClassProfile0Higgs.Fill(myEta, ratio);
00842        hElectronEtaClass0Higgs.Fill(myEta);
00843       }
00844       else if (electronClass == 1) { 
00845         hElectronClass1Higgs.Fill(ratio); if (ratio > 0.8 && ratio < 2.) hElectronClassProfile1Higgs.Fill(myEta, ratio); 
00846         hElectronEtaClass1Higgs.Fill(myEta);
00847       }
00848       else if (electronClass == 2)  {
00849         hElectronClass2Higgs.Fill(ratio); if (ratio > 0.8 && ratio < 2.) hElectronClassProfile2Higgs.Fill(myEta, ratio);
00850         hElectronEtaClass2Higgs.Fill(myEta);
00851       }
00852       else if (electronClass == 3)  {
00853         hElectronClass3Higgs.Fill(ratio); if (ratio > 0.8 && ratio < 2.) hElectronClassProfile3Higgs.Fill(myEta, ratio);
00854         hElectronEtaClass3Higgs.Fill(myEta);
00855       }
00856       else if (electronClass == 4)  {
00857         hElectronClass4Higgs.Fill(ratio); if (ratio > 0.8 && ratio < 2.) hElectronClassProfile4Higgs.Fill(myEta, ratio);
00858         hElectronEtaClass4Higgs.Fill(myEta);
00859       }
00860       else if (electronClass == 10) {
00861         hElectronClass10Higgs.Fill(ratio); if (ratio > 0.8 && ratio < 2.) hElectronClassProfile10Higgs.Fill(myEta, ratio);
00862         hElectronEtaClass10Higgs.Fill(myEta);
00863       }
00864       else if (electronClass == 11) {
00865         hElectronClass11Higgs.Fill(ratio); if (ratio > 0.8 && ratio < 2.) hElectronClassProfile11Higgs.Fill(myEta, ratio);
00866         hElectronEtaClass11Higgs.Fill(myEta);
00867       }
00868       else if (electronClass == 12) {
00869         hElectronClass12Higgs.Fill(ratio); if (ratio > 0.8 && ratio < 2.) hElectronClassProfile12Higgs.Fill(myEta, ratio);
00870         hElectronEtaClass12Higgs.Fill(myEta);
00871       }
00872       else if (electronClass == 13) {
00873         hElectronClass13Higgs.Fill(ratio); if (ratio > 0.8 && ratio < 2.) hElectronClassProfile13Higgs.Fill(myEta, ratio);
00874         hElectronEtaClass13Higgs.Fill(myEta);
00875       }
00876     }   
00877 
00878     //***********************************************************************
00879 
00880 
00881     //this part is commented because the isolation is now done before for the moment
00882 
00883     // // Calculate isolation in the region around four selected electrons
00884     //     // -----------------------------------------------------------------
00885     //     map<XANAElectronCandidate *,bool> mapSelectedElectrons;
00886     //     for (int ieh = 0; ieh < 4; ieh++) mapSelectedElectrons[electronsHiggs[ieh]] = true;
00887     //     XANAEvtRegionIsolation regionIsolation(event,mapSelectedElectrons);
00888 
00889     //     // Fill some histograms
00890     //     hIsoRegionNumberTracks.Fill((Axis_t) regionIsolation.getNumberTracks());
00891     //     hIsoRegionPtTracks.Fill((Axis_t) regionIsolation.getPtTracks());
00892     //     hIsoRegionPtMax.Fill((Axis_t) regionIsolation.getPtMax()); 
00893     
00894 
00895 
00896     //**********************************************************************
00897     
00898     // Cuts
00899     // ----
00900     
00901     // 1. Isolation cut
00902     // ................
00903     //  already done
00904     //  if (regionIsolation.getPtTracks() < 20.) nCutsIso++; // will be replaced by isIsolated() method    
00905     //     else {clean(event, generatorEvent, geantEvent);continue;}
00906 
00907                 
00908 
00909     // 2. pT electrons cut
00910     // ...................
00911     if (selector.isPtCut(5.,10.,15.,20.)) nCutsPt++;
00912     else {clean(event, generatorEvent, geantEvent);continue;}
00913     
00914     // 3. Real Z mass cut
00915     // ..................
00916     if (selector.isRealZMassCut(XANASelectHZZ::ZMASSPDG-15.,XANASelectHZZ::ZMASSPDG+15.)) nCutsRealZMass++;
00917     //if (selector.isRealZMassCut(40.,120.)) nCutsRealZMass++;
00918     else {clean(event, generatorEvent, geantEvent);continue; } 
00919 
00920     // 4. Virtual Z mass cut
00921     // .....................
00922     if (selector.isVirtualZMassCut(10.,60.)) nCutsVirtualZMass++;
00923     //if (selector.isVirtualZMassCut(0.,100.)) nCutsVirtualZMass++;
00924     else {clean(event, generatorEvent, geantEvent);continue; }   
00925     
00926     // 5. Higgs mass cut
00927     // .....................
00928     if (selector.isHiggsMassCut(140.,160.)) nCutsHiggsMass++;
00929     //if (selector.isHiggsMassCut(80.,180.)) nCutsHiggsMass++;
00930     else {clean(event, generatorEvent, geantEvent);continue; }
00931     
00932     
00933     // some analysis after all cuts
00934     
00935     hRealZMassCuts.Fill(realZ.m());
00936     hVirtualZMassCuts.Fill(virtualZ.m());
00937     hHiggsMassCuts.Fill(higgs.m());    
00938     
00939     for (int ieh = 0; ieh < 4; ieh++){
00940       hElectronPtCuts.Fill(electronsHiggs[ieh]->getPt());
00941       tmpElRec = electronsHiggs[ieh]->getMomentumAtVertex();
00942       tmpElGen = ((electronsRecGen.find(electronsHiggs[ieh]))->second)->getMomentum();
00943       hElectronRecGenPtHiggsCuts.Fill(tmpElRec.perp()/tmpElGen.perp()); 
00944       hElectronRecGenDrHiggsCuts.Fill(tmpElRec.deltaR(tmpElGen));        
00945     }
00946     
00947     // Surviving probability for different classes 
00948     if (issignal&&selectorMc.isFourElectrons()) mc4eCuts++;
00949     if (issignal&&selectorMc.isFourTaus()) mc4tCuts++;
00950     if (issignal&&selectorMc.isTwoElectronsTwoTaus()) mc2e2tCuts++;   
00951     
00952     
00953     event->clear(); 
00954     generatorEvent->clear();
00955     geantEvent->clear();
00956     
00957   }
00958         
00959   //        cout << "XANASelectMcHZZ: " << " 4E = " << 100.*mc4e/nevent << "% | 4T = " <<
00960   //                                         100.*mc4t/nevent << "% | 2E2T = " << 100.*mc2e2t/nevent << "% | OTHER(2E2M?) = " <<
00961   //                                              100.*(nevent-mc4e-mc4t-mc2e2t)/nevent << "%" << endl << endl;
00962   cout << " -------------------------------------------------------------------------------" << endl;
00963   cout << "                                 ACCEPTANCES                                    " << endl;
00964   cout << "                           TOTAL NUMBER OF EVENTS       = " << nevent << endl;
00965   cout << "                           TOTAL NUMBER OF EVENTS read  = " << ievent << endl;
00966   cout << "                           NUMBER of 4e events (signal) = " <<nconsidered<<endl;
00967   cout << " -------------------------------------------------------------------------------" << endl;
00968   cout << "                             #EVENTS       RELATIVE ACC. (%)      APSOLUTE ACC. (%)        " << endl;
00969   cout << " -------------------------------------------------------------------------------" << endl;
00970   cout << " 1. TRIGGER (L1)       = \t" << nTrigAccL1    << "\t\t"   << 100.*nTrigAccL1/nconsidered 
00971        << "\t\t\t" << 100.*nTrigAccL1/nconsidered << endl; 
00972   cout << " 1. TRIGGER (HLT)      = \t" << nTrigAccHLT       << "\t\t"   << 100.*nTrigAccHLT/nTrigAccL1
00973        << "\t\t\t" << 100.*nTrigAccHLT/nconsidered << endl;
00974   cout << " 1. TRIGGER (HLT elec) = \t" << nTrigAccHLTelec       << "\t\t"   << 100.*nTrigAccHLTelec/nTrigAccHLT 
00975        << "\t\t\t" << 100.*nTrigAccHLTelec/nconsidered << endl;
00976   cout << " 2. SIMPLE CUTS        = \t" << nPresAcc          << "\t\t"   << 100.*nPresAcc/nTrigAccHLTelec 
00977        << "\t\t\t" << 100.*nPresAcc/nconsidered << endl;
00978   cout << " 3. LOOSE ISOLATION    = \t" << nIsoAcc           << "\t\t"   << 100.*nIsoAcc/nPresAcc 
00979        << "\t\t\t" << 100.*nIsoAcc/nconsidered << endl;
00980   cout << " 3. TIGH ISOLATION     = \t" << niso              << "\t\t"   << 100.*niso/nIsoAcc
00981        << "\t\t\t" << 100.*niso/nconsidered << endl;
00982 
00983 
00984   cout << " 4. HZZ RECO           = \t" << n2e2pAcc          << "\t\t"   << 100.*n2e2pAcc/niso 
00985        << "\t\t\t" << 100.*n2e2pAcc/nconsidered << endl;
00986  //  cout << " 5. ISOLATION          = \t" << nCutsIso          << "\t\t"   << 100.*nCutsIso/n2e2pAcc 
00987 //        << "\t\t\t" << 100.*nCutsIso/nconsidered << endl;
00988   cout << " 6. ELECTRON PT        = \t" << nCutsPt           << "\t\t"   << 100.*nCutsPt/n2e2pAcc
00989        << "\t\t\t" << 100.*nCutsPt/nconsidered << endl;
00990   cout << " 7. REAL Z MASS        = \t" << nCutsRealZMass    << "\t\t"   << 100.*nCutsRealZMass/nCutsPt 
00991        << "\t\t\t" << 100.*nCutsRealZMass/nconsidered << endl;
00992   cout << " 8. VIRTUAL Z MASS     = \t" << nCutsVirtualZMass << "\t\t"   << 100.*nCutsVirtualZMass/nCutsRealZMass 
00993        << "\t\t\t" << 100.*nCutsVirtualZMass/nconsidered << endl;
00994   cout << " 9. HIGGS MASS         = \t" << nCutsHiggsMass    << "\t\t"   << 100.*nCutsHiggsMass/nCutsVirtualZMass 
00995        << "\t\t\t" << 100.*nCutsHiggsMass/nconsidered << endl;
00996   cout << " -------------------------------------------------------------------------------" << endl;  
00997   cout << "                              TOTAL ACCEPTANCE = " << 100.*nCutsHiggsMass/nconsidered << "%" << endl;
00998   cout << " -------------------------------------------------------------------------------" << endl; 
00999   
01000   cout << nevent << " events and " << nb << " bytes read " << endl;  
01001   cout << endl;
01002   
01003   // Exporting histograms to file
01004   
01005   f.Write();
01006   f.Close();
01007  
01008   event->Delete();
01009   generatorEvent->Delete();
01010   geantEvent->Delete();
01011   stopWatch.Stop();
01012   cout << "Time elapsed - CPU: " << stopWatch.CpuTime() <<
01013     " real: " << stopWatch.RealTime() << endl;
01014   
01015 }