Changeset 258 in svn for trunk/Resolutions.cpp

Timestamp:

Feb 9, 2009, 11:31:56 AM (16 years ago)

Author:

severine ovyn

Message:

resolution from delphes output

File:

• trunk/Resolutions.cpp (modified) (13 diffs)

trunk/Resolutions.cpp

@@ -60 +60 @@
        if(abs(gen1->PID)==15)
          {
+cout<<"au moins on a un tau "<<endl;
             int d1=gen1->D1;
             int d2=gen1->D2;
-
+cout<<"il a des filles? "<<endl;
             if((d1 < array.GetEntries()) && (d1 > 0) && (d2 < array.GetEntries()) && (d2 > 0))
                {
                  tauhad=true;
                  for(int d=d1; d < d2+1; d++)
-                    {
+                    {cout<<abs(array[d]->PID)<<"  "<<endl;
                       if(abs(array[d]->PID)== pE || abs(array[d]->PID)== pMU)tauhad=false;
                     }
@@ -77 +78 @@
  }
 
-
-
-void PairingJet(TLorentzVector &JETSm, const TLorentzVector &JET,  vector<fastjet::PseudoJet> sorted_jetsS)
+double EnergySmallCone(const vector<TLorentzVector> &towers, const float eta, const float phi,float energy_scone,float JET_seed) {
+  double Energie=0;
+  for(unsigned int i=0; i < towers.size(); i++) {
+      if(towers[i].Pt() < JET_seed) continue;
+      if((DeltaR(phi,eta,towers[i].Phi(),towers[i].Eta()) < energy_scone)) {
+          Energie += towers[i].E();
+      }
+  }
+  return Energie;
+}
+
+
+void PairingJet(TLorentzVector &JETSm, const TLorentzVector &JET, const TClonesArray *branchJet)
 {
   JETSm.SetPxPyPzE(0,0,0,0);
   float deltaRtest=5000;
-  for (unsigned int i = 0; i < sorted_jetsS.size(); i++) {
+  TIter itJet((TCollection*)branchJet);
+  TRootJet *jet;
+  itJet.Reset();
+  while( (jet = (TRootJet*) itJet.Next()) )
+   {
       TLorentzVector Att;
-      Att.SetPxPyPzE(sorted_jetsS[i].px(),sorted_jetsS[i].py(),sorted_jetsS[i].pz(),sorted_jetsS[i].E());
+      Att.SetPxPyPzE(jet->Px,jet->Py,jet->Pz,jet->E);
       if(DeltaR(JET.Phi(),JET.Eta(),Att.Phi(),Att.Eta()) < deltaRtest)
         {
@@ -97 +112 @@
 }
 
+void PairingElec(TLorentzVector &ELECSm, const TLorentzVector &ELEC, const TClonesArray *branchElec)
+{
+  ELECSm.SetPxPyPzE(0,0,0,0);
+  float deltaRtest=5000;
+  TIter itElec((TCollection*)branchElec);
+  TRootElectron *elec;
+  itElec.Reset();
+  while( (elec = (TRootElectron*) itElec.Next()) )
+   {
+      TLorentzVector Att;
+      Att.SetPxPyPzE(elec->Px,elec->Py,elec->Pz,elec->E);
+      if(DeltaR(ELEC.Phi(),ELEC.Eta(),Att.Phi(),Att.Eta()) < deltaRtest)
+        {
+          deltaRtest = DeltaR(ELEC.Phi(),ELEC.Eta(),Att.Phi(),Att.Eta());
+          if(deltaRtest < 0.025)
+             {
+               ELECSm = Att;
+             }
+        }
+    }
+}
+
+void PairingMuon(TLorentzVector &MUONSm, const TLorentzVector &MUON, const TClonesArray *branchMuon)
+{
+  MUONSm.SetPxPyPzE(0,0,0,0);
+  float deltaRtest=5000;
+  TIter itMuon((TCollection*)branchMuon);
+  TRootMuon *muon;
+  itMuon.Reset();
+  while( (muon = (TRootMuon*) itMuon.Next()) )
+   {
+      TLorentzVector Att;
+      Att.SetPxPyPzE(muon->Px,muon->Py,muon->Pz,muon->E);
+      if(DeltaR(MUON.Phi(),MUON.Eta(),Att.Phi(),Att.Eta()) < deltaRtest)
+        {
+          deltaRtest = DeltaR(MUON.Phi(),MUON.Eta(),Att.Phi(),Att.Eta());
+          if(deltaRtest < 0.025)
+             {
+               MUONSm = Att;
+             }
+        }
+    }
+}
+
+unsigned int NumTracks(const TClonesArray *branchTracks, const float pt_track, const float eta, const float phi,float track_scone) {
+  unsigned int numtrack=0;
+  TIter itTrack((TCollection*)branchTracks);
+  TRootTracks *track;
+  itTrack.Reset();
+  while( (track = (TRootTracks*) itTrack.Next()) )
+   {
+      if((track->PT < pt_track )||
+         (DeltaR(phi,eta,track->Phi,track->Eta) > track_scone)
+        )continue;
+      numtrack++;
+  }
+  return numtrack;
+}
+
+
+
 int main(int argc, char *argv[])
 {
   int appargc = 2;
-  char *appName = "Smearing";
+  char *appName = "Resolution";
   char *appargv[] = {appName, "-b"};
   TApplication app(appName, &appargc, appargv);
   
-  if(argc != 4 && argc != 3) {
-      cout << " Usage: " << argv[0] << " input_file" << " output_file" << " data_card " << endl;
-      cout << " input_list - list of files in Ntpl, StdHep of LHEF format," << endl;
+  if(argc != 3) {
+      cout << " Usage: " << argv[0] << " input_file" << " output_file" << endl;
+      cout << " input_list - list of files in root format," << endl;
       cout << " output_file - output file." << endl;
-      cout << " data_card - Datacard containing resolution variables for the detector simulation (optional) "<<endl;
       exit(1);
   }
@@ -115 +190 @@
   
   //read the input TROOT file
-  string inputFileList(argv[1]), outputfilename(argv[2]);
+  string inputfilename(argv[1]), outputfilename(argv[2]);
+
   if(outputfilename.find(".root") > outputfilename.length() ) {
     cout << "output_file should be a .root file!\n";
@@ -121 +197 @@
   }
 
+
+
   TFile *outputFile = TFile::Open(outputfilename.c_str(), "RECREATE"); // Creates the file, but should be closed just after
   outputFile->Close();
   
-  string line;
-  ifstream infile(inputFileList.c_str());
-  infile >> line; // the first line determines the type of input files
-  
-  DataConverter *converter=0;
-  
-  if(strstr(line.c_str(),".hep"))
-    {                           
-      cout<<"*************************************************************************"<<endl;
-      cout<<"************         StdHEP file format detected           **************"<<endl;
-      cout<<"************     Starting convertion to TRoot format       **************"<<endl;
-      cout<<"*************************************************************************"<<endl;
-      converter = new STDHEPConverter(inputFileList,outputfilename);//case ntpl file in input list
-    }
-  else if(strstr(line.c_str(),".lhe"))
-    {
-      cout<<"*************************************************************************"<<endl;
-      cout<<"************          LHEF file format detected            **************"<<endl;
-      cout<<"************     Starting convertion to TRoot format       **************"<<endl;
-      cout<<"*************************************************************************"<<endl;
-      converter = new LHEFConverter(inputFileList,outputfilename);//case ntpl file in input list
-    }
-  else if(strstr(line.c_str(),".root"))
-    {
-      cout<<"*************************************************************************"<<endl;
-      cout<<"************         h2root file format detected           **************"<<endl;
-      cout<<"************     Starting convertion to TRoot format       **************"<<endl;
-      cout<<"*************************************************************************"<<endl;
-      converter = new HEPEVTConverter(inputFileList,outputfilename);//case ntpl file in input list
-    }
-  else { cout << "***  " << line.c_str() << "\n***  file format not identified\n***  Exiting\n"; return -1;};
-  
-  TChain chain("GEN");
-  chain.Add(outputfilename.c_str());
+  TChain chainGEN("GEN");
+  chainGEN.Add(inputfilename.c_str());
+  ExRootTreeReader *treeReaderGEN = new ExRootTreeReader(&chainGEN);
+  TChain chain("Analysis");
+  chain.Add(inputfilename.c_str());
   ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
-  const TClonesArray *branchGen = treeReader->UseBranch("Particle");
+  const TClonesArray *branchJet = treeReader->UseBranch("Jet");
+  const TClonesArray *branchElec = treeReader->UseBranch("Electron");
+  const TClonesArray *branchMuon = treeReader->UseBranch("Muon");
+  const TClonesArray *branchTracks = treeReader->UseBranch("Tracks");
+  const TClonesArray *branchTowers = treeReader->UseBranch("CaloTower");
+  const TClonesArray *branchGen = treeReaderGEN->UseBranch("Particle");
   TIter itGen((TCollection*)branchGen);
   
@@ -178 +232 @@
   ETMIS *elementEtmis;
   
-int numTau=0;
-int numTauRec=0;
-  
-  //read the datacard input file
-  string DetDatacard("data/DataCardDet.dat");
-  if(argc==4)  DetDatacard =argv[3];
-  RESOLution *DET = new RESOLution();   
-  DET->ReadDataCard(DetDatacard);
+  int numTau=0;
+  int numTauRec=0;
+  
+  RESOLution *DET = new RESOLution();
+  DET->ReadDataCard("data/Datacard_CMS.dat");
+
 
   //Jet information
-  JetsUtil *JETRUN = new JetsUtil(DetDatacard);
-
-  //Propagation of tracks in the B field
-  //TrackPropagation *TRACP = new TrackPropagation(DetDatacard);
-   
+  JetsUtil *JETRUN = new JetsUtil("data/Datacard_CMS.dat");
+
   TLorentzVector genMomentum(0,0,0,0);//TLorentzVector containing generator level information
-  TLorentzVector recoMomentumCalo(0,0,0,0);
-  TLorentzVector recoMomentum(0,0,0,0);//TLorentzVector containing Reco level information
+  TLorentzVector recoMomentum(0,0,0,0);//TLorentzVector containing generator level information
   LorentzVector jetMomentum;
-  vector<TLorentzVector> TrackCentral;
   
   vector<fastjet::PseudoJet> input_particlesGEN;//for FastJet algorithm
   vector<fastjet::PseudoJet> sorted_jetsGEN;
   
-  vector<fastjet::PseudoJet> input_particlesReco;//for FastJet algorithm
-  vector<fastjet::PseudoJet> sorted_jetsReco;
-
-  vector<PhysicsTower> towers;
- 
- float iPhi=0,iEta=0;
+  vector<TLorentzVector> towers;
  
   // Loop over all events
@@ -217 +259 @@
       TLorentzVector PTmisGEN(0,0,0,0);
       treeReader->ReadEntry(entry);
+      treeReaderGEN->ReadEntry(entry);
       treeWriter->Clear();
       if((entry % 100) == 0 && entry > 0 )  cout << "** Processing element # " << entry << endl;
@@ -222 +265 @@
       TSimpleArray<TRootGenParticle> bGen;
       itGen.Reset();
-      TrackCentral.clear();
       TSimpleArray<TRootGenParticle> NFCentralQ;
       
-      input_particlesReco.clear();
       input_particlesGEN.clear();
       towers.clear();
@@ -245 +286 @@
           //Calculate ETMIS from generated particles
           if((pid == pNU1) || (pid == pNU2) || (pid == pNU3))PTmisGEN = PTmisGEN + genMomentum;
-          
-          if( (particle->Status == 1)    &&
-              ((pid != pNU1) && (pid != pNU2) && (pid != pNU3)) &&
-              (fabs(particle->Eta) < DET->CEN_max_calo_fwd)
-              )
+         
+          if( (particle->Status == 1)    &&
+            ((pid != pNU1) && (pid != pNU2) && (pid != pNU3)) &&
+             (fabs(particle->Eta) < DET->CEN_max_calo_fwd)
+            )
             {
-              recoMomentum = genMomentum;          
-              //use of the magnetic field propagation
-              //TRACP->Propagation(particle,recoMomentum);
-              //              cout<<"eta "<<eta<<endl;
-              eta=fabs(recoMomentum.Eta());
-              //            cout<<"eta apres"<<eta<<endl;
-              
-              switch(pid) {
-                
-              case pE: // all electrons with eta < DET->MAX_CALO_FWD
-                DET->SmearElectron(recoMomentum);
-                if(recoMomentum.E() !=0 && eta < DET->CEN_max_tracker){
-                  elementElec=(RESOLELEC*) branchelec->NewEntry();
-                  elementElec->E = genMomentum.E();
-                  elementElec->SmearedE = recoMomentum.E();}
-                break; // case pE
-              case pGAMMA: // all photons with eta < DET->MAX_CALO_FWD
-                DET->SmearElectron(recoMomentum);
-                break; // case pGAMMA
-              case pMU: // all muons with eta < DET->MAX_MU
-                DET->SmearMu(recoMomentum);
-                if(recoMomentum.E() !=0){
-                elementMuon = (RESOLMUON*) branchmuon->NewEntry();
-                elementMuon->OverPT = (1/genMomentum.Pt());
-                elementMuon->OverSmearedPT = (1/recoMomentum.Pt());}
-                break; // case pMU
-              case pLAMBDA: // all lambdas with eta < DET->MAX_CALO_FWD
-              case pK0S:    // all K0s with eta < DET->MAX_CALO_FWD
-                DET->SmearHadron(recoMomentum, 0.7);
-                break; // case hadron
-              default:   // all other final particles with eta < DET->MAX_CALO_FWD
-                DET->SmearHadron(recoMomentum, 1.0);
-                break;
-              } // switch (pid)
-              
-              //information to reconstruct jets from reconstructed towers
-               int charge=Charge(pid);
-                if(recoMomentum.E() !=0 && pid != pMU) {
-                        // in case the Bfield is not simulated, checks that charged particles have enough pt to reach the calos
-                  if ( !DET->FLAG_bfield && charge!=0 && genMomentum.Pt() <= DET->TRACK_ptmin ) { /* particules do not reach calos */ }
-                  else { // particles reach calos
-                  DET->BinEtaPhi(recoMomentum.Phi(), recoMomentum.Eta(), iPhi, iEta);
-                  if(iEta != -100 && iPhi != -100)
-                    {
-                      recoMomentumCalo.SetPtEtaPhiE(recoMomentum.Pt(),iEta,iPhi,recoMomentum.E());
-
-                      PhysicsTower Tower(LorentzVector(recoMomentumCalo.Px(),recoMomentumCalo.Py(),recoMomentumCalo.Pz(), recoMomentumCalo.E()));
-                      towers.push_back(Tower);
-                    }
-                }
-              }
-             
-              // all final charged particles
-              if(
-                 (recoMomentum.E()!=0) &&
-                 (fabs(recoMomentum.Eta()) < DET->CEN_max_tracker) &&
-                 (DET->FLAG_bfield || ( !DET->FLAG_bfield && genMomentum.Pt() > DET->TRACK_ptmin )) &&
-                        // if bfield not simulated, pt should be high enough to be taken into account
-                 ((rand()%100) < DET->TRACK_eff)  &&
-                 (charge!=0)
-                 )
-                {TrackCentral.push_back(recoMomentum);}
-              
-            } // switch
+              eta=fabs(genMomentum.Eta());
+
+              switch(pid) {
+
+                 case pE: // all electrons with eta < DET->MAX_CALO_FWD
+                    PairingElec(recoMomentum,genMomentum,branchElec);
+                    if(recoMomentum.Pt()>1){
+                     elementElec=(RESOLELEC*) branchelec->NewEntry();
+                     elementElec->E = genMomentum.E();
+                     elementElec->SmearedE = recoMomentum.E();}
+                   break; // case pE
+                 case pMU: // all muons with eta < DET->MAX_MU
+                    PairingMuon(recoMomentum,genMomentum,branchMuon);
+                    if(recoMomentum.E() !=0){
+                     elementMuon = (RESOLMUON*) branchmuon->NewEntry();
+                     elementMuon->OverPT = (1/genMomentum.Pt());
+                     elementMuon->OverSmearedPT = (1/recoMomentum.Pt());}
+                   break; // case pMU
+                 default:
+                   break;
+              } // switch (pid)
+           }
+
         } // while
       
       //compute missing transverse energy from calo towers    
+      TIter itCalo((TCollection*)branchTowers);
+      TRootCalo *calo;
+      itCalo.Reset();
       TLorentzVector Att(0.,0.,0.,0.);
       float ScalarEt=0;
-      for(unsigned int i=0; i < towers.size(); i++)
-        {
-          Att.SetPxPyPzE(towers[i].fourVector.px, towers[i].fourVector.py, towers[i].fourVector.pz, towers[i].fourVector.E);
+      while( (calo = (TRootCalo*) itCalo.Next()) )
+           {
+          //Att.SetPxPyPzE(towers[i].fourVector.px, towers[i].fourVector.py, towers[i].fourVector.pz, towers[i].fourVector.E);
+          Att.SetPtEtaPhiE(calo->PT,calo->Eta,calo->Phi,calo->E);
+          towers.push_back(Att);
           if(fabs(Att.Eta()) < DET->CEN_max_calo_fwd)
             {
                ScalarEt = ScalarEt + Att.Et();
                PTmisReco = PTmisReco + Att;
-               // create a fastjet::PseudoJet with these components and put it onto
-               // back of the input_particles vector
-               input_particlesReco.push_back(fastjet::PseudoJet(towers[i].fourVector.px,towers[i].fourVector.py,towers[i].fourVector.pz,towers[i].fourVector.E));
             }
         }
@@ -344 +346 @@
       
       sorted_jetsGEN=JETRUN->RunJets(input_particlesGEN);
-      sorted_jetsReco=JETRUN->RunJets(input_particlesReco);
 
       TSimpleArray<TRootGenParticle> TausHadr = TauHadr(branchGen);
+cout<<"nombre de tau-jets "<<TausHadr.GetEntries()<<endl;
 
       TLorentzVector JETreco(0,0,0,0);
@@ -352 +354 @@
         TLorentzVector JETgen(0,0,0,0);
         JETgen.SetPxPyPzE(sorted_jetsGEN[i].px(),sorted_jetsGEN[i].py(),sorted_jetsGEN[i].pz(),sorted_jetsGEN[i].E());
-        PairingJet(JETreco,JETgen,sorted_jetsReco);
+        PairingJet(JETreco,JETgen,branchJet);
         if(JETreco.Pt()>1)
           {
@@ -361 +363 @@
       }
       numTau = numTau+TausHadr.GetEntries();
-      for (unsigned int i = 0; i < sorted_jetsReco.size(); i++) {
-        TLorentzVector JETT(0,0,0,0);
-        JETT.SetPxPyPzE(sorted_jetsReco[i].px(),sorted_jetsReco[i].py(),sorted_jetsReco[i].pz(),sorted_jetsReco[i].E());
-        if(fabs(JETT.Eta()) < (DET->CEN_max_tracker - DET->TAU_track_scone)) 
+
+      TIter itJet((TCollection*)branchJet);
+      TRootJet *jet;
+      itJet.Reset();
+//cout<<"a"<<endl;
+      while( (jet = (TRootJet*) itJet.Next()) )
+        {
+//cout<<"b"<<endl;
+          TLorentzVector JETT(0,0,0,0);
+          JETT.SetPxPyPzE(jet->Px,jet->Py,jet->Pz,jet->E);
+          if(fabs(JETT.Eta()) < (DET->CEN_max_tracker - DET->TAU_track_scone)) 
            {
+//cout<<"c"<<endl;
             for(Int_t i=0; i<TausHadr.GetEntries();i++)
               {
+//cout<<"d"<<endl;
                 if(DeltaR(TausHadr[i]->Phi,TausHadr[i]->Eta,JETT.Phi(),JETT.Eta())<0.1)
                   {
+//cout<<"e"<<endl;
                     elementTaujet= (TAUHAD*) branchtaujet->NewEntry();
-                    elementTaujet->EnergieCen = (DET->EnergySmallCone(towers,JETT.Eta(),JETT.Phi())/JETT.E());
-                    elementTaujet->NumTrack = DET->NumTracks(TrackCentral,DET->TAU_track_pt,JETT.Eta(),JETT.Phi());
-                    if( (DET->EnergySmallCone(towers,JETT.Eta(),JETT.Phi())/JETT.E()) > 0.95 
-                     &&  (DET->NumTracks(TrackCentral,DET->TAU_track_pt,JETT.Eta(),JETT.Phi()))==1)numTauRec++;
+                    elementTaujet->EnergieCen = (EnergySmallCone(towers,JETT.Eta(),JETT.Phi(),DET->TAU_energy_scone,DET->JET_seed)/JETT.E());
+                    elementTaujet->NumTrack = NumTracks(branchTracks,DET->TAU_track_pt,JETT.Eta(),JETT.Phi(),DET->TAU_track_scone);
+                    if( (EnergySmallCone(towers,JETT.Eta(),JETT.Phi(),DET->TAU_energy_scone,DET->JET_seed)/JETT.E()) > 0.95 
+                     &&  (NumTracks(branchTracks,DET->TAU_track_pt,JETT.Eta(),JETT.Phi(),DET->TAU_track_scone))==1)numTauRec++;
+//cout<<"f"<<endl;
                     
                  }
+//cout<<"g"<<endl;
               }
            }
+//cout<<"h"<<endl;
          
         
       } // for itJet : loop on all jets
+//cout<<"i"<<endl;
       
       treeWriter->Fill();
@@ -397 +413 @@
   delete treeReader;
   delete DET;
-  if(converter) delete converter;
-}
-
+}
+