Fork me on GitHub

Changeset 30 in svn


Ignore:
Timestamp:
Nov 13, 2008, 5:10:52 PM (16 years ago)
Author:
severine ovyn
Message:

remove bug

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Delphes.cpp

    r29 r30  
    123123    }
    124124  else { cout << "***  " << line.c_str() << "\n***  file format not identified\n***  Exiting\n"; return -1;};
    125 
     125 
    126126  TChain chain("GEN");
    127127  chain.Add(outputfilename.c_str());
     
    143143  ExRootTreeBranch *branchRP220 = treeWriter->NewBranch("RP220hits", TRootRomanPotHits::Class());
    144144  ExRootTreeBranch *branchFP420 = treeWriter->NewBranch("FP420hits", TRootRomanPotHits::Class());
    145 
    146  
     145 
     146  
    147147  TRootGenParticle *particle;
    148148  TRootETmis *elementEtmis;
     
    163163  DET->ReadDataCard(DetDatacard);
    164164 
    165  
    166165  TLorentzVector genMomentum(0,0,0,0);
    167166  LorentzVector jetMomentum;
    168167  vector<TLorentzVector> TrackCentral; 
    169168  vector<PhysicsTower> towers;
    170 
     169 
    171170  vector<fastjet::PseudoJet> input_particles;//for FastJet algorithm
    172171  vector<fastjet::PseudoJet> inclusive_jets;
    173172  vector<fastjet::PseudoJet> sorted_jets;
     173 
     174  vector<TLorentzVector> electron;
     175  vector<int> elecPID;
     176  vector<TLorentzVector> muon;
     177  vector<int> muonPID;
     178  TSimpleArray<TRootGenParticle> NFCentralQ;
    174179 
    175180  //Initialisation of Hector
     
    196201  beamline2->add(rp220_2);
    197202  beamline2->add(rp420_2);
    198 
     203 
    199204  // we will have four jet definitions, and the first three will be
    200205  // plugins
    201206  fastjet::JetDefinition jet_def;
    202207  fastjet::JetDefinition::Plugin * plugins;
    203 
     208 
    204209  switch(DET->JETALGO) {
    205210  default:
    206211  case 1: {
    207 
    208    // set up a CDF midpoint jet definition
    209    #ifdef ENABLE_PLUGIN_CDFCONES
    210    plugins = new fastjet::CDFJetCluPlugin(DET->C_SEEDTHRESHOLD,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->C_OVERLAPTHRESHOLD);
    211    jet_def = fastjet::JetDefinition(plugins);
    212    #else
    213    plugins = NULL;
    214    #endif
    215    }
    216    break;
    217 
    218    case 2: {
    219 
    220    // set up a CDF midpoint jet definition
    221    #ifdef ENABLE_PLUGIN_CDFCONES
    222    plugins = new fastjet::CDFMidPointPlugin (DET->M_SEEDTHRESHOLD,DET->CONERADIUS,DET->M_CONEAREAFRACTION,DET->M_MAXPAIRSIZE,DET->M_MAXPAIRSIZE,DET->C_OVERLAPTHRESHOLD);
    223    jet_def = fastjet::JetDefinition(plugins);
    224    #else
    225    plugins = NULL;
    226    #endif
    227    }
    228    break;
    229    case 3: {
    230    // set up a siscone jet definition
    231    #ifdef ENABLE_PLUGIN_SISCONE
    232    int npass = 0;               // do infinite number of passes
    233    double protojet_ptmin = 0.0; // use all protojets
    234    plugins = new fastjet::SISConePlugin (DET->CONERADIUS,DET->C_OVERLAPTHRESHOLD,npass, protojet_ptmin);
    235    jet_def = fastjet::JetDefinition(plugins);
    236    #else
    237    plugins = NULL;
    238    #endif
    239    }
    240    break;
    241 
     212   
     213    // set up a CDF midpoint jet definition
     214#ifdef ENABLE_PLUGIN_CDFCONES
     215    plugins = new fastjet::CDFJetCluPlugin(DET->C_SEEDTHRESHOLD,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->C_OVERLAPTHRESHOLD);
     216    jet_def = fastjet::JetDefinition(plugins);
     217#else
     218    plugins = NULL;
     219#endif
     220  }
     221    break;
     222   
     223  case 2: {
     224   
     225    // set up a CDF midpoint jet definition
     226#ifdef ENABLE_PLUGIN_CDFCONES
     227    plugins = new fastjet::CDFMidPointPlugin (DET->M_SEEDTHRESHOLD,DET->CONERADIUS,DET->M_CONEAREAFRACTION,DET->M_MAXPAIRSIZE,DET->M_MAXPAIRSIZE,DET->C_OVERLAPTHRESHOLD);
     228    jet_def = fastjet::JetDefinition(plugins);
     229#else
     230    plugins = NULL;
     231#endif
     232  }
     233    break;
     234  case 3: {
     235    // set up a siscone jet definition
     236#ifdef ENABLE_PLUGIN_SISCONE
     237    int npass = 0;               // do infinite number of passes
     238    double protojet_ptmin = 0.0; // use all protojets
     239    plugins = new fastjet::SISConePlugin (DET->CONERADIUS,DET->C_OVERLAPTHRESHOLD,npass, protojet_ptmin);
     240    jet_def = fastjet::JetDefinition(plugins);
     241#else
     242    plugins = NULL;
     243#endif
     244  }
     245    break;
     246   
    242247  case 4: {
    243   jet_def = fastjet::JetDefinition(fastjet::kt_algorithm, DET->CONERADIUS);
    244   //jet_defs[4] = fastjet::JetDefinition(fastjet::cambridge_algorithm,jet_radius);
    245   //jet_defs[5] = fastjet::JetDefinition(fastjet::antikt_algorithm,jet_radius);
    246   }
    247   break;
    248   }
    249  
    250    // Loop over all events
     248    jet_def = fastjet::JetDefinition(fastjet::kt_algorithm, DET->CONERADIUS);
     249    //jet_defs[4] = fastjet::JetDefinition(fastjet::cambridge_algorithm,jet_radius);
     250    //jet_defs[5] = fastjet::JetDefinition(fastjet::antikt_algorithm,jet_radius);
     251  }
     252    break;
     253  }
     254  
     255  // Loop over all events
    251256  Long64_t entry, allEntries = treeReader->GetEntries();
    252257  cout << "** Chain contains " << allEntries << " events" << endl;
     
    259264      if((entry % 100) == 0 && entry > 0 )  cout << "** Processing element # " << entry << endl;
    260265     
    261       TSimpleArray<TRootGenParticle> electron;
    262       TSimpleArray<TRootGenParticle> muon;
     266      electron.clear();
     267      muon.clear();
     268      elecPID.clear();
     269      muonPID.clear();
     270      NFCentralQ.Clear();
     271     
    263272      itGen.Reset();
    264273      TrackCentral.clear();
    265274      towers.clear();
    266275      input_particles.clear();
    267       TSimpleArray<TRootGenParticle> NFCentralQ;
    268276      inclusive_jets.clear();
    269277      sorted_jets.clear();
    270 
     278     
    271279      // Loop over all particles in event
    272280      while( (particle = (TRootGenParticle*) itGen.Next()) )
    273281        {
    274282          genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
    275 
     283         
    276284          int pid = abs(particle->PID);
    277285          float eta = fabs(particle->Eta);
    278 
     286         
    279287          //subarray of the quarks (i.e. not final particles, i.e status not equal to 1)
    280288          // in the tracker (i.e. for those we know the tracks)
     
    288296            NFCentralQ.Add(particle);
    289297          }
    290 
    291 
     298         
     299         
    292300          // keeps only final particles, visible by the central detector, including the fiducial volume
    293301          // the ordering of conditions have been optimised for speed : put first the STATUS condition
    294302          if( (particle->Status == 1)    &&
    295303              (
    296                 (pid == pMU  && eta < DET->MAX_MU) ||
    297                 (pid != pMU  && (pid != pNU1) && (pid != pNU2) && (pid != pNU3) && eta < DET->MAX_CALO_FWD)
    298               )
    299             ) {
    300                 switch(pid) {
    301                    
    302                     case pE: // all electrons with eta < DET->MAX_CALO_FWD
    303                       DET->SmearElectron(genMomentum);
    304                       electron.Add(particle);
    305                       break; // case pE
    306                     case pGAMMA: // all photons with eta < DET->MAX_CALO_FWD
    307                       DET->SmearElectron(genMomentum);
    308                       if(genMomentum.E()!=0 && eta < DET->MAX_TRACKER) {
    309                           elementPhoton = (TRootPhoton*) branchPhoton->NewEntry();
    310                           elementPhoton->Set(genMomentum);
    311                       }
    312                       break; // case pGAMMA
    313                     case pMU: // all muons with eta < DET->MAX_MU
    314                       DET->SmearMu(genMomentum);
    315                       muon.Add(particle);
    316                       break; // case pMU
    317                     case pLAMBDA: // all lambdas with eta < DET->MAX_CALO_FWD
    318                     case pK0S:    // all K0s with eta < DET->MAX_CALO_FWD
    319                       DET->SmearHadron(genMomentum, 0.7);
    320                       break; // case hadron
    321                     default:   // all other final particles with eta < DET->MAX_CALO_FWD
    322                       DET->SmearHadron(genMomentum, 1.0);
    323                       break;
    324                 } // switch (pid)
    325                
    326                 // all final particles but muons and neutrinos 
    327                 // for calorimetric towers and mission PT
    328                 if(genMomentum.E()!=0) {
    329                   PTmis = PTmis + genMomentum;//ptmis
    330                   if(pid !=pMU) {
    331                     towers.push_back(PhysicsTower(LorentzVector(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())));
    332                     // create a fastjet::PseudoJet with these components and put it onto
    333                     // back of the input_particles vector
    334                     input_particles.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
    335                     elementCalo = (TRootCalo*) branchCalo->NewEntry();
    336                     elementCalo->Set(genMomentum);
    337                   }
    338                 }
    339                
    340                
    341                 // all final charged particles
    342                 if(
    343                    ((rand()%100) < DET->TRACKING_EFF)  &&
    344                    (genMomentum.E()!=0) &&
    345                    (fabs(particle->Eta) < DET->MAX_TRACKER) &&
    346                    (genMomentum.Pt() > DET->PT_TRACKS_MIN ) &&     // pt too small to be taken into account
    347                    (pid != pGAMMA)    &&
    348                    (pid != pPI0)      &&
    349                    (pid != pK0L)      &&
    350                    (pid != pN)        &&
    351                    (pid != pSIGMA0)   &&
    352                    (pid != pDELTA0)   &&
    353                    (pid != pK0S)   // not charged particles : invisible by tracker
    354                    )
    355                   {
    356                     elementTracks = (TRootTracks*) branchTracks->NewEntry();
    357                     elementTracks->Set(genMomentum);
    358                     TrackCentral.push_back(genMomentum);
    359                   }
     304               (pid == pMU  && eta < DET->MAX_MU) ||
     305               (pid != pMU  && (pid != pNU1) && (pid != pNU2) && (pid != pNU3) && eta < DET->MAX_CALO_FWD)
     306                )
     307              ) {
     308            switch(pid) {
     309             
     310            case pE: // all electrons with eta < DET->MAX_CALO_FWD
     311              DET->SmearElectron(genMomentum);
     312              electron.push_back(genMomentum);
     313              elecPID.push_back(particle->PID);
     314              break; // case pE
     315            case pGAMMA: // all photons with eta < DET->MAX_CALO_FWD
     316              DET->SmearElectron(genMomentum);
     317              if(genMomentum.E()!=0 && eta < DET->MAX_TRACKER) {
     318                elementPhoton = (TRootPhoton*) branchPhoton->NewEntry();
     319                elementPhoton->Set(genMomentum);
     320              }
     321              break; // case pGAMMA
     322            case pMU: // all muons with eta < DET->MAX_MU
     323              DET->SmearMu(genMomentum);
     324              muonPID.push_back(particle->PID);
     325              muon.push_back(genMomentum);
     326              break; // case pMU
     327            case pLAMBDA: // all lambdas with eta < DET->MAX_CALO_FWD
     328            case pK0S:    // all K0s with eta < DET->MAX_CALO_FWD
     329              DET->SmearHadron(genMomentum, 0.7);
     330              break; // case hadron
     331            default:   // all other final particles with eta < DET->MAX_CALO_FWD
     332              DET->SmearHadron(genMomentum, 1.0);
     333              break;
     334            } // switch (pid)
     335           
     336            // all final particles but muons and neutrinos     
     337            // for calorimetric towers and mission PT
     338            if(genMomentum.E()!=0) {
     339              PTmis = PTmis + genMomentum;//ptmis
     340              if(pid !=pMU) {
     341                towers.push_back(PhysicsTower(LorentzVector(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())));
     342                // create a fastjet::PseudoJet with these components and put it onto
     343                // back of the input_particles vector
     344                input_particles.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
     345                elementCalo = (TRootCalo*) branchCalo->NewEntry();
     346                elementCalo->Set(genMomentum);
     347              }
     348            }
     349           
     350           
     351            // all final charged particles
     352            if(
     353               ((rand()%100) < DET->TRACKING_EFF)  &&
     354               (genMomentum.E()!=0) &&
     355               (fabs(particle->Eta) < DET->MAX_TRACKER) &&
     356               (genMomentum.Pt() > DET->PT_TRACKS_MIN ) &&     // pt too small to be taken into account
     357               (pid != pGAMMA)    &&
     358               (pid != pPI0)      &&
     359               (pid != pK0L)      &&
     360               (pid != pN)        &&
     361               (pid != pSIGMA0)   &&
     362               (pid != pDELTA0)   &&
     363               (pid != pK0S)   // not charged particles : invisible by tracker
     364               )
     365              {
     366                elementTracks = (TRootTracks*) branchTracks->NewEntry();
     367                elementTracks->Set(genMomentum);
     368                TrackCentral.push_back(genMomentum);
     369              }
    360370          } // switch
    361 
    362           //branchElectron->Reset();
    363          
    364           if(electron.GetEntries()>0) {
    365           for(int i=0; i < electron.GetEntries();i++) {
    366                genMomentum.SetPxPyPzE(electron[i]->Px,electron[i]->Py,electron[i]->Pz,electron[i]->E);
    367                if(genMomentum.E()!=0 && fabs(genMomentum.Eta()) < DET->MAX_TRACKER) {
    368                elementElec = (TRootElectron*) branchElectron->NewEntry();
    369                elementElec->Set(genMomentum);
    370                elementElec->Charge = sign(electron[i]->PID);
    371                if(fabs(genMomentum.Eta()))elementElec->IsolFlag = DET->Isolation(genMomentum.Phi(),genMomentum.Eta(),TrackCentral,2.0);
    372                }
    373           }}
    374           if(muon.GetEntries()>0) {
    375           for(int i=0; i < muon.GetEntries();i++) {
    376                genMomentum.SetPxPyPzE(muon[i]->Px,muon[i]->Py,muon[i]->Pz,muon[i]->E);
    377                if(genMomentum.E()!=0 && fabs(genMomentum.Eta()) < DET->MAX_MU) {
    378                  elementMu = (TRootMuon*) branchMuon->NewEntry();
    379                  elementMu->Charge = sign(particle->PID);
    380                  elementMu->Set(genMomentum);
    381                  elementMu->IsolFlag = DET->Isolation(genMomentum.Phi(),genMomentum.Eta(),TrackCentral,2.0);
    382                }
    383           }}
    384 
     371         
    385372          // Forward particles in CASTOR ?
    386373          /*      if (particle->Status == 1 && (fabs(particle->Eta) > DET->MIN_CALO_VFWD)
     
    461448        } // while
    462449     
     450      for(unsigned int i=0; i < electron.size(); i++) {
     451        if(electron[i].E()!=0 && fabs(electron[i].Eta()) < DET->MAX_TRACKER)
     452          {
     453            elementElec = (TRootElectron*) branchElectron->NewEntry();
     454            elementElec->Set(electron[i]);
     455            elementElec->Charge = sign(elecPID[i]);
     456            elementElec->IsolFlag = DET->Isolation(electron[i].Phi(),electron[i].Eta(),TrackCentral,2.0);
     457          }
     458      }
     459      for(unsigned int i=0; i < muon.size(); i++) {
     460        if(muon[i].E()!=0 && fabs(muon[i].Eta()) < DET->MAX_MU)
     461          {
     462            elementMu = (TRootMuon*) branchMuon->NewEntry();
     463            elementMu->Charge = sign(muonPID[i]);
     464            elementMu->Set(muon[i]);
     465            elementMu->IsolFlag = DET->Isolation(muon[i].Phi(),muon[i].Eta(),TrackCentral,2.0);
     466          }
     467      }
     468     
     469     
    463470      // computes the Missing Transverse Momentum
    464471      elementEtmis = (TRootETmis*) branchETmis->NewEntry();
     
    469476     
    470477      //*****************************
    471 
     478     
    472479      // run the jet clustering with the above jet definition
    473480      if(input_particles.size()!=0)
Note: See TracChangeset for help on using the changeset viewer.