Changeset 30 in svn

Timestamp:

Nov 13, 2008, 5:10:52 PM (16 years ago)

Author:

severine ovyn

Message:

remove bug

File:

• trunk/Delphes.cpp (modified) (8 diffs)

trunk/Delphes.cpp

@@ -123 +123 @@
     }
   else { cout << "***  " << line.c_str() << "\n***  file format not identified\n***  Exiting\n"; return -1;};
-
+  
   TChain chain("GEN");
   chain.Add(outputfilename.c_str());
@@ -143 +143 @@
   ExRootTreeBranch *branchRP220 = treeWriter->NewBranch("RP220hits", TRootRomanPotHits::Class());
   ExRootTreeBranch *branchFP420 = treeWriter->NewBranch("FP420hits", TRootRomanPotHits::Class());
-
- 
+  
+  
   TRootGenParticle *particle; 
   TRootETmis *elementEtmis;
@@ -163 +163 @@
   DET->ReadDataCard(DetDatacard);
   
-  
   TLorentzVector genMomentum(0,0,0,0);
   LorentzVector jetMomentum;
   vector<TLorentzVector> TrackCentral;  
   vector<PhysicsTower> towers;
-
+  
   vector<fastjet::PseudoJet> input_particles;//for FastJet algorithm
   vector<fastjet::PseudoJet> inclusive_jets;
   vector<fastjet::PseudoJet> sorted_jets;
+  
+  vector<TLorentzVector> electron;
+  vector<int> elecPID;
+  vector<TLorentzVector> muon;
+  vector<int> muonPID;
+  TSimpleArray<TRootGenParticle> NFCentralQ;
   
   //Initialisation of Hector
@@ -196 +201 @@
   beamline2->add(rp220_2);
   beamline2->add(rp420_2);
-
+  
   // we will have four jet definitions, and the first three will be
   // plugins
   fastjet::JetDefinition jet_def;
   fastjet::JetDefinition::Plugin * plugins;
-
+  
   switch(DET->JETALGO) {
   default:
   case 1: {
-
-   // set up a CDF midpoint jet definition
-   #ifdef ENABLE_PLUGIN_CDFCONES
-   plugins = new fastjet::CDFJetCluPlugin(DET->C_SEEDTHRESHOLD,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->C_OVERLAPTHRESHOLD);
-   jet_def = fastjet::JetDefinition(plugins);
-   #else
-   plugins = NULL;
-   #endif
-   }
-   break;
-
-   case 2: {
-
-   // set up a CDF midpoint jet definition
-   #ifdef ENABLE_PLUGIN_CDFCONES
-   plugins = new fastjet::CDFMidPointPlugin (DET->M_SEEDTHRESHOLD,DET->CONERADIUS,DET->M_CONEAREAFRACTION,DET->M_MAXPAIRSIZE,DET->M_MAXPAIRSIZE,DET->C_OVERLAPTHRESHOLD);
-   jet_def = fastjet::JetDefinition(plugins);
-   #else
-   plugins = NULL;
-   #endif
-   }
-   break;
-   case 3: {
-   // set up a siscone jet definition
-   #ifdef ENABLE_PLUGIN_SISCONE
-   int npass = 0;               // do infinite number of passes
-   double protojet_ptmin = 0.0; // use all protojets
-   plugins = new fastjet::SISConePlugin (DET->CONERADIUS,DET->C_OVERLAPTHRESHOLD,npass, protojet_ptmin);
-   jet_def = fastjet::JetDefinition(plugins);
-   #else
-   plugins = NULL;
-   #endif
-   }
-   break;
-
+    
+    // set up a CDF midpoint jet definition
+#ifdef ENABLE_PLUGIN_CDFCONES
+    plugins = new fastjet::CDFJetCluPlugin(DET->C_SEEDTHRESHOLD,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->C_OVERLAPTHRESHOLD);
+    jet_def = fastjet::JetDefinition(plugins);
+#else
+    plugins = NULL;
+#endif
+  }
+    break;
+    
+  case 2: {
+    
+    // set up a CDF midpoint jet definition
+#ifdef ENABLE_PLUGIN_CDFCONES
+    plugins = new fastjet::CDFMidPointPlugin (DET->M_SEEDTHRESHOLD,DET->CONERADIUS,DET->M_CONEAREAFRACTION,DET->M_MAXPAIRSIZE,DET->M_MAXPAIRSIZE,DET->C_OVERLAPTHRESHOLD);
+    jet_def = fastjet::JetDefinition(plugins);
+#else
+    plugins = NULL;
+#endif
+  }
+    break;
+  case 3: {
+    // set up a siscone jet definition
+#ifdef ENABLE_PLUGIN_SISCONE
+    int npass = 0;               // do infinite number of passes
+    double protojet_ptmin = 0.0; // use all protojets
+    plugins = new fastjet::SISConePlugin (DET->CONERADIUS,DET->C_OVERLAPTHRESHOLD,npass, protojet_ptmin);
+    jet_def = fastjet::JetDefinition(plugins);
+#else
+    plugins = NULL;
+#endif
+  }
+    break;
+    
   case 4: {
-  jet_def = fastjet::JetDefinition(fastjet::kt_algorithm, DET->CONERADIUS);
-  //jet_defs[4] = fastjet::JetDefinition(fastjet::cambridge_algorithm,jet_radius);
-  //jet_defs[5] = fastjet::JetDefinition(fastjet::antikt_algorithm,jet_radius);
-  }
-  break;
-  }
- 
-   // Loop over all events
+    jet_def = fastjet::JetDefinition(fastjet::kt_algorithm, DET->CONERADIUS);
+    //jet_defs[4] = fastjet::JetDefinition(fastjet::cambridge_algorithm,jet_radius);
+    //jet_defs[5] = fastjet::JetDefinition(fastjet::antikt_algorithm,jet_radius);
+  }
+    break;
+  }
+  
+  // Loop over all events
   Long64_t entry, allEntries = treeReader->GetEntries();
   cout << "** Chain contains " << allEntries << " events" << endl;
@@ -259 +264 @@
       if((entry % 100) == 0 && entry > 0 )  cout << "** Processing element # " << entry << endl;
       
-      TSimpleArray<TRootGenParticle> electron;
-      TSimpleArray<TRootGenParticle> muon;
+      electron.clear();
+      muon.clear();
+      elecPID.clear();
+      muonPID.clear();
+      NFCentralQ.Clear();
+      
       itGen.Reset();
       TrackCentral.clear();
       towers.clear();
       input_particles.clear();
-      TSimpleArray<TRootGenParticle> NFCentralQ;
       inclusive_jets.clear();
       sorted_jets.clear();
-
+      
       // Loop over all particles in event
       while( (particle = (TRootGenParticle*) itGen.Next()) )
         {
           genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
-
+          
           int pid = abs(particle->PID);
           float eta = fabs(particle->Eta);
-
+          
           //subarray of the quarks (i.e. not final particles, i.e status not equal to 1) 
           // in the tracker (i.e. for those we know the tracks)
@@ -288 +296 @@
             NFCentralQ.Add(particle);
           }
-
-
+          
+          
           // keeps only final particles, visible by the central detector, including the fiducial volume
           // the ordering of conditions have been optimised for speed : put first the STATUS condition
           if( (particle->Status == 1)    &&
               ( 
-                (pid == pMU  && eta < DET->MAX_MU) ||
-                (pid != pMU  && (pid != pNU1) && (pid != pNU2) && (pid != pNU3) && eta < DET->MAX_CALO_FWD) 
-              )
-            ) { 
-                switch(pid) {
-                    
-                    case pE: // all electrons with eta < DET->MAX_CALO_FWD
-                      DET->SmearElectron(genMomentum);
-                      electron.Add(particle);
-                      break; // case pE
-                    case pGAMMA: // all photons with eta < DET->MAX_CALO_FWD
-                      DET->SmearElectron(genMomentum);
-                      if(genMomentum.E()!=0 && eta < DET->MAX_TRACKER) {
-                          elementPhoton = (TRootPhoton*) branchPhoton->NewEntry();
-                          elementPhoton->Set(genMomentum);
-                      }
-                      break; // case pGAMMA
-                    case pMU: // all muons with eta < DET->MAX_MU
-                      DET->SmearMu(genMomentum);
-                      muon.Add(particle);
-                      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(genMomentum, 0.7);
-                      break; // case hadron
-                    default:   // all other final particles with eta < DET->MAX_CALO_FWD
-                      DET->SmearHadron(genMomentum, 1.0);
-                      break;
-                } // switch (pid)
-                
-                // all final particles but muons and neutrinos  
-                // for calorimetric towers and mission PT
-                if(genMomentum.E()!=0) {
-                  PTmis = PTmis + genMomentum;//ptmis
-                  if(pid !=pMU) {
-                    towers.push_back(PhysicsTower(LorentzVector(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())));
-                    // create a fastjet::PseudoJet with these components and put it onto
-                    // back of the input_particles vector
-                    input_particles.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())); 
-                    elementCalo = (TRootCalo*) branchCalo->NewEntry();
-                    elementCalo->Set(genMomentum);
-                  }
-                }
-                
-                
-                // all final charged particles
-                if(
-                   ((rand()%100) < DET->TRACKING_EFF)  &&
-                   (genMomentum.E()!=0) &&
-                   (fabs(particle->Eta) < DET->MAX_TRACKER) && 
-                   (genMomentum.Pt() > DET->PT_TRACKS_MIN ) &&     // pt too small to be taken into account
-                   (pid != pGAMMA)    && 
-                   (pid != pPI0)      &&
-                   (pid != pK0L)      &&
-                   (pid != pN)        &&
-                   (pid != pSIGMA0)   &&
-                   (pid != pDELTA0)   &&
-                   (pid != pK0S)   // not charged particles : invisible by tracker
-                   )
-                  {
-                    elementTracks = (TRootTracks*) branchTracks->NewEntry();
-                    elementTracks->Set(genMomentum);
-                    TrackCentral.push_back(genMomentum);
-                  }
+               (pid == pMU  && eta < DET->MAX_MU) ||
+               (pid != pMU  && (pid != pNU1) && (pid != pNU2) && (pid != pNU3) && eta < DET->MAX_CALO_FWD) 
+                )
+              ) { 
+            switch(pid) {
+              
+            case pE: // all electrons with eta < DET->MAX_CALO_FWD
+              DET->SmearElectron(genMomentum);
+              electron.push_back(genMomentum);
+              elecPID.push_back(particle->PID);
+              break; // case pE
+            case pGAMMA: // all photons with eta < DET->MAX_CALO_FWD
+              DET->SmearElectron(genMomentum);
+              if(genMomentum.E()!=0 && eta < DET->MAX_TRACKER) {
+                elementPhoton = (TRootPhoton*) branchPhoton->NewEntry();
+                elementPhoton->Set(genMomentum);
+              }
+              break; // case pGAMMA
+            case pMU: // all muons with eta < DET->MAX_MU
+              DET->SmearMu(genMomentum);
+              muonPID.push_back(particle->PID);
+              muon.push_back(genMomentum);
+              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(genMomentum, 0.7);
+              break; // case hadron
+            default:   // all other final particles with eta < DET->MAX_CALO_FWD
+              DET->SmearHadron(genMomentum, 1.0);
+              break;
+            } // switch (pid)
+            
+            // all final particles but muons and neutrinos      
+            // for calorimetric towers and mission PT
+            if(genMomentum.E()!=0) {
+              PTmis = PTmis + genMomentum;//ptmis
+              if(pid !=pMU) {
+                towers.push_back(PhysicsTower(LorentzVector(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())));
+                // create a fastjet::PseudoJet with these components and put it onto
+                // back of the input_particles vector
+                input_particles.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())); 
+                elementCalo = (TRootCalo*) branchCalo->NewEntry();
+                elementCalo->Set(genMomentum);
+              }
+            }
+            
+            
+            // all final charged particles
+            if(
+               ((rand()%100) < DET->TRACKING_EFF)  &&
+               (genMomentum.E()!=0) &&
+               (fabs(particle->Eta) < DET->MAX_TRACKER) && 
+               (genMomentum.Pt() > DET->PT_TRACKS_MIN ) &&     // pt too small to be taken into account
+               (pid != pGAMMA)    && 
+               (pid != pPI0)      &&
+               (pid != pK0L)      &&
+               (pid != pN)        &&
+               (pid != pSIGMA0)   &&
+               (pid != pDELTA0)   &&
+               (pid != pK0S)   // not charged particles : invisible by tracker
+               )
+              {
+                elementTracks = (TRootTracks*) branchTracks->NewEntry();
+                elementTracks->Set(genMomentum);
+                TrackCentral.push_back(genMomentum);
+              }
           } // switch
-
-          //branchElectron->Reset();
-          
-          if(electron.GetEntries()>0) {
-          for(int i=0; i < electron.GetEntries();i++) { 
-               genMomentum.SetPxPyPzE(electron[i]->Px,electron[i]->Py,electron[i]->Pz,electron[i]->E);
-               if(genMomentum.E()!=0 && fabs(genMomentum.Eta()) < DET->MAX_TRACKER) {
-               elementElec = (TRootElectron*) branchElectron->NewEntry();
-               elementElec->Set(genMomentum);
-               elementElec->Charge = sign(electron[i]->PID);
-               if(fabs(genMomentum.Eta()))elementElec->IsolFlag = DET->Isolation(genMomentum.Phi(),genMomentum.Eta(),TrackCentral,2.0);
-               }
-          }}
-          if(muon.GetEntries()>0) {
-          for(int i=0; i < muon.GetEntries();i++) {
-               genMomentum.SetPxPyPzE(muon[i]->Px,muon[i]->Py,muon[i]->Pz,muon[i]->E);
-               if(genMomentum.E()!=0 && fabs(genMomentum.Eta()) < DET->MAX_MU) {
-                 elementMu = (TRootMuon*) branchMuon->NewEntry();
-                 elementMu->Charge = sign(particle->PID);
-                 elementMu->Set(genMomentum);
-                 elementMu->IsolFlag = DET->Isolation(genMomentum.Phi(),genMomentum.Eta(),TrackCentral,2.0);
-               }
-          }}
-
+          
           // Forward particles in CASTOR ?
           /*      if (particle->Status == 1 && (fabs(particle->Eta) > DET->MIN_CALO_VFWD) 
@@ -461 +448 @@
         } // while
       
+      for(unsigned int i=0; i < electron.size(); i++) {
+        if(electron[i].E()!=0 && fabs(electron[i].Eta()) < DET->MAX_TRACKER)
+          {
+            elementElec = (TRootElectron*) branchElectron->NewEntry();
+            elementElec->Set(electron[i]);
+            elementElec->Charge = sign(elecPID[i]);
+            elementElec->IsolFlag = DET->Isolation(electron[i].Phi(),electron[i].Eta(),TrackCentral,2.0);
+          }
+      }
+      for(unsigned int i=0; i < muon.size(); i++) {
+        if(muon[i].E()!=0 && fabs(muon[i].Eta()) < DET->MAX_MU) 
+          {
+            elementMu = (TRootMuon*) branchMuon->NewEntry();
+            elementMu->Charge = sign(muonPID[i]);
+            elementMu->Set(muon[i]);
+            elementMu->IsolFlag = DET->Isolation(muon[i].Phi(),muon[i].Eta(),TrackCentral,2.0);
+          }
+      }
+      
+      
       // computes the Missing Transverse Momentum
       elementEtmis = (TRootETmis*) branchETmis->NewEntry();
@@ -469 +476 @@
       
       //*****************************
-
+      
       // run the jet clustering with the above jet definition
       if(input_particles.size()!=0)