Changeset 30 in svn for trunk/Delphes.cpp
- Timestamp:
- Nov 13, 2008, 5:10:52 PM (16 years ago)
- File:
-
- 1 edited
-
trunk/Delphes.cpp (modified) (8 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/Delphes.cpp
r29 r30 123 123 } 124 124 else { cout << "*** " << line.c_str() << "\n*** file format not identified\n*** Exiting\n"; return -1;}; 125 125 126 126 TChain chain("GEN"); 127 127 chain.Add(outputfilename.c_str()); … … 143 143 ExRootTreeBranch *branchRP220 = treeWriter->NewBranch("RP220hits", TRootRomanPotHits::Class()); 144 144 ExRootTreeBranch *branchFP420 = treeWriter->NewBranch("FP420hits", TRootRomanPotHits::Class()); 145 146 145 146 147 147 TRootGenParticle *particle; 148 148 TRootETmis *elementEtmis; … … 163 163 DET->ReadDataCard(DetDatacard); 164 164 165 166 165 TLorentzVector genMomentum(0,0,0,0); 167 166 LorentzVector jetMomentum; 168 167 vector<TLorentzVector> TrackCentral; 169 168 vector<PhysicsTower> towers; 170 169 171 170 vector<fastjet::PseudoJet> input_particles;//for FastJet algorithm 172 171 vector<fastjet::PseudoJet> inclusive_jets; 173 172 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; 174 179 175 180 //Initialisation of Hector … … 196 201 beamline2->add(rp220_2); 197 202 beamline2->add(rp420_2); 198 203 199 204 // we will have four jet definitions, and the first three will be 200 205 // plugins 201 206 fastjet::JetDefinition jet_def; 202 207 fastjet::JetDefinition::Plugin * plugins; 203 208 204 209 switch(DET->JETALGO) { 205 210 default: 206 211 case 1: { 207 208 // set up a CDF midpoint jet definition209 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 213 plugins = NULL;214 215 216 break;217 218 219 220 // set up a CDF midpoint jet definition221 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 225 plugins = NULL;226 227 228 break;229 230 // set up a siscone jet definition231 232 int npass = 0; // do infinite number of passes233 double protojet_ptmin = 0.0; // use all protojets234 plugins = new fastjet::SISConePlugin (DET->CONERADIUS,DET->C_OVERLAPTHRESHOLD,npass, protojet_ptmin);235 jet_def = fastjet::JetDefinition(plugins);236 237 plugins = NULL;238 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 242 247 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 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 251 256 Long64_t entry, allEntries = treeReader->GetEntries(); 252 257 cout << "** Chain contains " << allEntries << " events" << endl; … … 259 264 if((entry % 100) == 0 && entry > 0 ) cout << "** Processing element # " << entry << endl; 260 265 261 TSimpleArray<TRootGenParticle> electron; 262 TSimpleArray<TRootGenParticle> muon; 266 electron.clear(); 267 muon.clear(); 268 elecPID.clear(); 269 muonPID.clear(); 270 NFCentralQ.Clear(); 271 263 272 itGen.Reset(); 264 273 TrackCentral.clear(); 265 274 towers.clear(); 266 275 input_particles.clear(); 267 TSimpleArray<TRootGenParticle> NFCentralQ;268 276 inclusive_jets.clear(); 269 277 sorted_jets.clear(); 270 278 271 279 // Loop over all particles in event 272 280 while( (particle = (TRootGenParticle*) itGen.Next()) ) 273 281 { 274 282 genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E); 275 283 276 284 int pid = abs(particle->PID); 277 285 float eta = fabs(particle->Eta); 278 286 279 287 //subarray of the quarks (i.e. not final particles, i.e status not equal to 1) 280 288 // in the tracker (i.e. for those we know the tracks) … … 288 296 NFCentralQ.Add(particle); 289 297 } 290 291 298 299 292 300 // keeps only final particles, visible by the central detector, including the fiducial volume 293 301 // the ordering of conditions have been optimised for speed : put first the STATUS condition 294 302 if( (particle->Status == 1) && 295 303 ( 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 } 360 370 } // 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 385 372 // Forward particles in CASTOR ? 386 373 /* if (particle->Status == 1 && (fabs(particle->Eta) > DET->MIN_CALO_VFWD) … … 461 448 } // while 462 449 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 463 470 // computes the Missing Transverse Momentum 464 471 elementEtmis = (TRootETmis*) branchETmis->NewEntry(); … … 469 476 470 477 //***************************** 471 478 472 479 // run the jet clustering with the above jet definition 473 480 if(input_particles.size()!=0)
Note:
See TracChangeset
for help on using the changeset viewer.
