1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
19 //_________________________________________________________________________
20 // Class for JetFinder Input preparation from simulated data
22 //*-- Author: Mark Horner (LBL/UCT)
27 #include "AliEMCALJetFinderInputSimPrep.h"
29 #include <TParticle.h>
30 #include <TParticlePDG.h>
34 #include <TObjectTable.h>
39 #include "AliRunLoader.h"
42 #include "AliEMCALFast.h"
44 #include "AliEMCALHit.h"
45 #include "AliEMCALGeometry.h"
46 #include "AliEMCALLoader.h"
47 #include "AliGenEventHeader.h"
48 #include "AliGenPythiaEventHeader.h"
49 #include "AliGenerator.h"
50 #include "AliHeader.h"
55 ClassImp(AliEMCALJetFinderInputSimPrep)
57 AliEMCALJetFinderInputSimPrep::AliEMCALJetFinderInputSimPrep()
59 // Default constructor
60 if (fDebug > 0) Info("AliEMCALJetFinderInputSimPrep","Beginning Constructor");
63 fSmearType = kSmearEffic;
65 fTrackType = kCharged;
69 AliEMCALJetFinderInputSimPrep::~AliEMCALJetFinderInputSimPrep()
71 if (fDebug > 0) Info("~AliEMCALJetFinderInputSimPrep","Beginning Destructor");
75 void AliEMCALJetFinderInputSimPrep::Reset(AliEMCALJetFinderResetType_t resettype)
77 // Method to reset data
78 if (fDebug > 1) Info("Reset","Beginning Reset");
82 fInputObject.Reset(resettype);
85 fInputObject.Reset(resettype);
88 fInputObject.Reset(resettype);
90 case kResetParameters:
91 fSmearType = kSmearEffic;
93 fTrackType = kCharged;
96 fSmearType = kSmearEffic;
98 fTrackType = kCharged;
99 fInputObject.Reset(resettype);
102 Warning("FillFromFile", "kResetPartons not implemented") ;
104 case kResetParticles:
105 Warning("FillFromFile", "kResetParticles not implemented") ;
108 Warning("FillFromFile", "kResetJets not implemented") ;
114 Int_t AliEMCALJetFinderInputSimPrep::FillFromFile(TString *filename, AliEMCALJetFinderFileType_t filetype,Int_t EventNumber,TString data)
116 if (fDebug > 1) Info("FillFromFile","Beginning FillFromFile");
117 fFileType = filetype;
118 AliRunLoader *rl=AliRunLoader::Open(filename->Data());
119 if (fDebug > 1) Info("FillFromFile","Opened file %s",filename->Data());
120 if (data.Contains("S"))
122 rl->LoadKinematics();
124 rl->GetEvent(EventNumber);
127 rl->LoadKinematics();
129 rl->GetEvent(EventNumber);
131 if (fDebug > 1) Info("FillFromFile","Got event %i with option \"XH\"",EventNumber);
133 if ( fEMCALType == kHits ||
134 fEMCALType == kTimeCut )
136 if (data.Contains("S"))
144 if ( fTrackType != kNoTracks ) FillTracks();
145 if ( fFileType != kData){
151 void AliEMCALJetFinderInputSimPrep::FillHits()
153 // Fill from the hits to input object from simulation
154 if (fDebug > 1) Info("FillHits","Beginning FillHits");
156 // Access hit information
157 // AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
158 // Info("AliEMCALJetFinderInputSimPrep","Title of the geometry is %s",pEMCAL->GetTitle());
159 AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
160 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
162 // TTree *treeH = AliEMCALGetter::Instance()->TreeH();
163 // Int_t ntracks = (Int_t) treeH->GetEntries();
169 /* for (Int_t track=0; track<ntracks;track++) {
171 nbytes += treeH->GetEvent(track);*/
175 for(Int_t i =0; i < emcalLoader->Hits()->GetEntries() ;i++)
177 const AliEMCALHit *mHit = emcalLoader->Hit(i);
178 if (fEMCALType == kTimeCut &&
179 (mHit->GetTime()>fTimeCut) ) continue;
180 Float_t x = mHit->X(); // x-pos of hit
181 Float_t y = mHit->Y(); // y-pos
182 Float_t z = mHit->Z(); // z-pos
183 Float_t eloss = mHit->GetEnergy(); // deposited energy
185 Float_t r = TMath::Sqrt(x*x+y*y);
186 Float_t theta = TMath::ATan2(r,z);
187 Float_t eta = -TMath::Log(TMath::Tan(theta/2.));
188 Float_t phi = TMath::ATan2(y,x);
189 etH = eloss*TMath::Sin(theta);
190 if (fDebug > 10) Info("FillHits","Values of hit energy %i",Int_t(1e7*etH));
191 if (geom->TowerIndexFromEtaPhi(eta,180.0/TMath::Pi()*phi) == -1)
195 Error("FillHits","Hit not inside EMCAL!");
200 fInputObject.AddEnergyToDigit(geom->TowerIndexFromEtaPhi(eta,180.0/TMath::Pi()*phi)-1,Int_t(1e7*etH));
206 void AliEMCALJetFinderInputSimPrep::FillTracks()
208 // Fill from particles simulating a TPC to input object from simulation
209 if (fDebug > 1) Info("FillTracks","Beginning FillTracks");
211 AliRunLoader *rl = AliRunLoader::GetRunLoader();
212 TParticlePDG* pdgP = 0;
214 Int_t npart = rl->Stack()->GetNprimary();
215 if (fDebug > 1) Info("FillTracks","Header says there are %i primaries",npart);
216 Float_t bfield,rEMCAL;
218 if (fDebug > 1) Info("FillTracks","Defining the geometry");
220 AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
221 AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
222 if (geom == 0 && pEMCAL)
223 geom = AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
225 geom = AliEMCALGeometry::GetInstance("EMCAL_55_25","");//","SHISH_77_TRD1_2X2_FINAL_110DEG");
226 fEtaMax = geom->GetArm1EtaMax();
227 fEtaMin = geom->GetArm1EtaMin();
228 fPhiMax = TMath::Pi()/180.0*geom->GetArm1PhiMax();
229 fPhiMin = TMath::Pi()/180.0*geom->GetArm1PhiMin();
231 if (fDebug > 1) Info("FillTracks","Starting particle loop");
233 if (fDebug > 1) Info("FillTracks","Particle loop of %i",npart);
234 for (Int_t part = 0; part < npart; part++) {
235 mPart = rl->Stack()->Particle(part);
236 pdgP = mPart->GetPDG();
238 if (fDebug > 5) Info("FillTracks","Checking if track is a primary");
240 if (fFileType == kPythia) {
241 if (mPart->GetStatusCode() != 1) continue;
242 } else if (fFileType == kHijing) {
243 if (mPart->GetFirstDaughter() >= 0 && mPart->GetFirstDaughter() < npart) continue;
246 if (fDebug > 15) Info("FillTracks","Checking if track (eta - %f, phi - %f) is in acceptance",mPart->Eta(),mPart->Phi());
247 if (fDebug > 10) Info("FillTracks","Checking if EMCAL acceptance ( %f < eta < %f, %f < phi < %f) is in acceptance",fEtaMin,fEtaMax,fPhiMin,fPhiMax);
249 if ((!fPythiaComparison)&&(mPart->Eta() > fEtaMax || mPart->Eta() < fEtaMin)) continue;
250 if ((!fPythiaComparison)&&(mPart->Phi() > fPhiMax || mPart->Phi() < fPhiMin)) continue;
253 {kProton, kProtonBar, kElectron, kPositron,
254 kNuE, kNuEBar, kGamma, kNeutron, kNeutronBar,
255 kMuonPlus, kMuonMinus, kK0Long , kPiPlus, kPiMinus,
256 kKPlus, kKMinus, kLambda0, kLambda0Bar, kK0Short,
257 kSigmaMinus, kSigmaPlus, kSigma0, kPi0, kK0, kK0Bar,
261 if (fDebug > 15) Info("FillTracks","Checking if track is rejected");
263 if ((fSmearType == kEfficiency ||
264 fSmearType == kSmearEffic) &&
266 if (AliEMCALFast::RandomReject(fEfficiency)) {
271 if (gAlice && gAlice->Field())
272 bfield = gAlice->Field()->SolenoidField();
275 rEMCAL = AliEMCALGeometry::GetInstance()->GetIPDistance();
276 Float_t rB = 3335.6 * mPart->Pt() / bfield; // [cm] (case of |charge|=1)
277 if (2.*rB < rEMCAL) continue; // track curls away
279 //if (part%10) gObjectTable->Print();
283 case kAllP: // All Stable particles to be included
284 if (fDebug > 5) Info("FillTracks","Storing track");
285 if (fSmearType == kSmear ||
286 fSmearType == kSmearEffic ){
288 TParticle *tmp = Smear(MPart);
289 fInputObject.AddTrack(Smear(MPart));
292 fInputObject.AddTrack(*mPart);
295 case kEM: // All Electromagnetic particles
296 if (mPart->GetPdgCode() == kElectron ||
297 mPart->GetPdgCode() == kMuonPlus ||
298 mPart->GetPdgCode() == kMuonMinus ||
299 mPart->GetPdgCode() == kPositron ){
300 if (fDebug > 5) Info("FillTracks","Storing electron or positron");
301 if (fSmearType == kSmear ||
302 fSmearType == kSmearEffic ){
304 TParticle *tmp = Smear(MPart);
305 fInputObject.AddTrack(tmp);
308 fInputObject.AddTrack(*mPart);
311 if ( mPart->GetPdgCode() == kGamma ){
312 fInputObject.AddTrack(*mPart);
313 if (fDebug > 5) Info("FillTracks","Storing gamma");
317 case kCharged: // All Particles with non-zero charge
318 if (pdgP->Charge() != 0) {
319 if (fDebug > 5) Info("FillTracks","Storing charged track");
320 if (fSmearType == kSmear ||
321 fSmearType == kSmearEffic ){
323 TParticle *tmp = Smear(MPart);
324 fInputObject.AddTrack(tmp);
327 fInputObject.AddTrack(*mPart);
331 case kNeutral: // All particles with no charge
332 if (pdgP->Charge() == 0){
333 fInputObject.AddTrack(*mPart);
334 if (fDebug > 5) Info("FillTracks","Storing neutral");
337 case kHadron: //All hadrons
338 if (mPart->GetPdgCode() != kElectron &&
339 mPart->GetPdgCode() != kPositron &&
340 mPart->GetPdgCode() != kMuonPlus &&
341 mPart->GetPdgCode() != kMuonMinus &&
342 mPart->GetPdgCode() != kGamma )
344 if (fDebug > 5) Info("FillTracks","Storing hadron");
345 if (pdgP->Charge() == 0){
346 fInputObject.AddTrack(*mPart);
348 if (fSmearType == kSmear ||
349 fSmearType == kSmearEffic ){
351 TParticle *tmp = Smear(MPart);
352 fInputObject.AddTrack(tmp);
355 fInputObject.AddTrack(*mPart);
360 case kChargedHadron: // only charged hadrons
361 if (mPart->GetPdgCode() != kElectron &&
362 mPart->GetPdgCode() != kPositron &&
363 mPart->GetPdgCode() != kGamma &&
364 mPart->GetPdgCode() != kMuonPlus &&
365 mPart->GetPdgCode() != kMuonMinus &&
366 pdgP->Charge() != 0 ){
367 if (fDebug > 5) Info("FillTracks","Storing charged hadron");
368 if (fSmearType == kSmear ||
369 fSmearType == kSmearEffic ){
371 TParticle *tmp = Smear(MPart);
372 fInputObject.AddTrack(tmp);
375 fInputObject.AddTrack(*mPart);
382 if (pdgP->Charge() != 0 || mPart->GetPdgCode() == kPi0 ||
383 mPart->GetPdgCode() == kGamma )
385 if (fDebug > 5) Info("FillTracks","Storing charged track");
386 if (fSmearType == kSmear ||
387 fSmearType == kSmearEffic ){
389 TParticle *tmp = Smear(MPart);
390 fInputObject.AddTrack(tmp);
393 fInputObject.AddTrack(*mPart);
397 case kNoNeutronNeutrinoKlong:
398 if ( mPart->GetPdgCode() != kNeutron &&
399 mPart->GetPdgCode() != kNeutronBar &&
400 mPart->GetPdgCode() != kK0Long &&
401 mPart->GetPdgCode() != kNuE &&
402 mPart->GetPdgCode() != kNuEBar &&
403 mPart->GetPdgCode() != kNuMu &&
404 mPart->GetPdgCode() != kNuMuBar &&
405 mPart->GetPdgCode() != kNuTau &&
406 mPart->GetPdgCode() != kNuTauBar )
408 if (fDebug > 5) Info("FillTracks","Storing charged track");
409 if (fSmearType == kSmear ||
410 fSmearType == kSmearEffic ){
412 TParticle *tmp = Smear(MPart);
413 fInputObject.AddTrack(tmp);
416 fInputObject.AddTrack(*mPart);
424 // Info("FillTracks","After Particle Storage");
425 //if (part%10) gObjectTable->Print();
427 } //end of particle loop
428 //gObjectTable->Print();
431 void AliEMCALJetFinderInputSimPrep::FillPartons()
434 // input object from simulation
435 if (fDebug > 1) Info("FillPartons","Beginning FillPartons");
437 AliRunLoader *rl = AliRunLoader::GetRunLoader();
438 AliGenEventHeader* evHeader = rl->GetHeader()->GenEventHeader();
439 Float_t triggerJetValues[4];
440 AliEMCALParton tempParton;
442 if (fFileType == kPythia)
444 Int_t ntriggerjets = ((AliGenPythiaEventHeader*)evHeader)->NTriggerJets();
445 if (fDebug > 1) Info("FillPartons","Number of trigger jets --> %i",ntriggerjets);
446 for (Int_t jetcounter = 0 ; jetcounter < ntriggerjets; jetcounter++){
447 ((AliGenPythiaEventHeader*)evHeader)->TriggerJet(jetcounter,triggerJetValues);
448 TLorentzVector tempLParton;
449 tempLParton.SetPxPyPzE(triggerJetValues[0],triggerJetValues[1],triggerJetValues[2],triggerJetValues[3]);
450 tempParton.SetEnergy(tempLParton.Energy());
451 tempParton.SetEta(tempLParton.Eta());
452 tempParton.SetPhi(tempLParton.Phi());
453 FillPartonTracks(&tempParton);
454 fInputObject.AddParton(&tempParton);
460 void AliEMCALJetFinderInputSimPrep::FillParticles()
462 // Fill particles to input object from simulation
463 if (fDebug > 1) Info("FillParticles","Beginning FillParticles");
465 Int_t npart = (gAlice->GetHeader())->GetNprimary();
466 TParticlePDG* pdgP = 0;
468 for (Int_t part = 0; part < npart; part++) {
469 TParticle *mPart = gAlice->GetMCApp()->Particle(part);
470 pdgP = mPart->GetPDG();
472 if (fDebug > 10) Info("FillParticles","Checking if particle is a primary");
474 if (fFileType == kPythia) {
475 if (mPart->GetStatusCode() != 1) continue;
476 } else if (fFileType == kHijing) {
477 if (mPart->GetFirstDaughter() >= 0 && mPart->GetFirstDaughter() < npart) continue;
480 if (fDebug > 10) Info("FillParticles","Checking if particle is in acceptance");
482 if ((!fPythiaComparison)&&(mPart->Eta() > fEtaMax || mPart->Eta() < fEtaMin)) continue;
483 if ((!fPythiaComparison)&&(mPart->Phi() > fPhiMax || mPart->Phi() < fPhiMin)) continue;
487 {kProton, kProtonBar, kElectron, kPositron,
488 kNuE, kNuEBar, kGamma, kNeutron, kNeutronBar,
489 kMuonPlus, kMuonMinus, kK0Long , kPiPlus, kPiMinus,
490 kKPlus, kKMinus, kLambda0, kLambda0Bar, kK0Short,
491 kSigmaMinus, kSigmaPlus, kSigma0, kPi0, kK0, kK0Bar,
497 case kAllP: // All Stable particles to be included
498 if (fDebug > 5) Info("FillParticles","Storing particle");
499 fInputObject.AddParticle(mPart);
501 case kEM: // All Electromagnetic particles
502 if (mPart->GetPdgCode() == kElectron ||
503 mPart->GetPdgCode() == kPositron ||
504 mPart->GetPdgCode() == kGamma){
505 if (fDebug > 5) Info("FillParticles","Storing electromagnetic particle");
506 fInputObject.AddParticle(mPart);
509 case kCharged: // All Particles with non-zero charge
510 if (pdgP->Charge() != 0) {
511 if (fDebug > 5) Info("FillParticles","Storing charged particle");
512 fInputObject.AddParticle(mPart);
515 case kNeutral: // All particles with no charge
516 if (pdgP->Charge() == 0){
517 if (fDebug > 5) Info("FillParticles","Storing neutral particle");
518 fInputObject.AddParticle(mPart);
521 case kHadron: //All hadrons
522 if (mPart->GetPdgCode() != kElectron &&
523 mPart->GetPdgCode() != kPositron &&
524 mPart->GetPdgCode() != kGamma )
527 if (fDebug > 5) Info("FillParticles","Storing hadron");
528 fInputObject.AddParticle(mPart);
531 case kChargedHadron: // only charged hadrons
532 if (mPart->GetPdgCode() != kElectron &&
533 mPart->GetPdgCode() != kPositron &&
534 mPart->GetPdgCode() != kGamma &&
535 pdgP->Charge() != 0 ){
536 if (fDebug > 5) Info("FillParticles","Storing charged hadron");
537 fInputObject.AddParticle(mPart);
543 if (pdgP->Charge() != 0 || mPart->GetPdgCode() == kPi0 ||
544 mPart->GetPdgCode() == kGamma )
546 if (fDebug > 5) Info("FillTracks","Storing kEMChargedPi0 track");
547 if (fSmearType == kSmear ||
548 fSmearType == kSmearEffic ){
550 TParticle *tmp = Smear(MPart);
551 fInputObject.AddTrack(tmp);
554 fInputObject.AddTrack(*mPart);
558 case kNoNeutronNeutrinoKlong:
559 if ( mPart->GetPdgCode() != kNeutron &&
560 mPart->GetPdgCode() != kNeutronBar &&
561 mPart->GetPdgCode() != kK0Long &&
562 mPart->GetPdgCode() != kNuE &&
563 mPart->GetPdgCode() != kNuEBar &&
564 mPart->GetPdgCode() != kNuMu &&
565 mPart->GetPdgCode() != kNuMuBar &&
566 mPart->GetPdgCode() != kNuTau &&
567 mPart->GetPdgCode() != kNuTauBar )
569 if (fDebug > 5) Info("FillTracks","Storing kNoNeutronNeutrinoKlong track");
570 if (fSmearType == kSmear ||
571 fSmearType == kSmearEffic ){
573 TParticle *tmp = Smear(MPart);
574 fInputObject.AddTrack(tmp);
577 fInputObject.AddTrack(*mPart);
584 }// end of loop over particles
587 void AliEMCALJetFinderInputSimPrep::FillDigits()
589 // Fill digits to input object
592 void AliEMCALJetFinderInputSimPrep::FillSDigits()
594 Info("FillSDigits","Beginning FillSDigits");
595 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
597 // Fill digits to input object
598 for (Int_t towerid=0; towerid < emcalLoader->SDigits()->GetEntries(); towerid++)
600 fInputObject.AddEnergyToDigit(emcalLoader->SDigit(towerid)->GetId(), emcalLoader->SDigit(towerid)->GetAmp());
605 void AliEMCALJetFinderInputSimPrep::Smear(TParticle *particle)
607 // Smear particle momentum
608 if (fDebug > 5) Info("Smear","Beginning Smear");
610 Float_t tmp = AliEMCALFast::SmearMomentum(1,particle->P());
611 Float_t tmpt = (tmp/particle->P()) * particle->Pt();
612 if (fDebug > 10) Info("Smear","Smearing changed P from %f to %f",particle->Pt(),tmpt);
613 TLorentzVector tmpvector;
614 tmpvector.SetPtEtaPhiM(tmpt,particle->Eta(), particle->Phi(), particle->GetMass());
615 // create a new particle with smearing momentum - and no daughter or parent information and the same
618 TParticle tmparticle(particle->GetPdgCode(), 1, 0, 0,0,0,
619 tmpvector.Px() , tmpvector.Py(), tmpvector.Pz(), tmpvector.Energy(),
620 particle->Vx(), particle->Vy(), particle->Vz(), particle->T());
621 fInputObject.AddTrack(tmparticle);
625 void AliEMCALJetFinderInputSimPrep::FillPartonTracks(AliEMCALParton *parton)
627 // Populate parton tracks for input distributions
628 if (fDebug>1) Info("FillPartonTracks","Beginning FillPartonTracks");
629 AliRunLoader *rl = AliRunLoader::GetRunLoader();
631 Int_t npart = rl->Stack()->GetNprimary();
633 TParticlePDG *getpdg;
635 for (Int_t part = 0; part < npart; part++)
637 tempPart = rl->Stack()->Particle(part);
638 if (tempPart->GetStatusCode() != 1) continue;
639 if ((!fPythiaComparison)&&(tempPart->Eta() > fEtaMax || tempPart->Eta() < fEtaMin ||
640 tempPart->Phi() > fPhiMax || tempPart->Phi() < fPhiMin) ){
641 if (fDebug>10) Info("FillPartonTracks","Excluding particle not pointing at the EMCAL");
644 getpdg = tempPart->GetPDG();
645 if (getpdg->Charge() == 0.0 ) {
646 if (fDebug>10) Info("FillPartonTracks","Excluding a neutral particle");
649 if ( (parton->Eta() - tempPart->Eta())*(parton->Eta() - tempPart->Eta()) +
650 (parton->Phi() - tempPart->Phi())*(parton->Phi() - tempPart->Phi()) < 1.0 ) ntracks++;
652 Float_t *energy = new Float_t[ntracks];
653 Float_t *eta = new Float_t[ntracks];
654 Float_t *phi = new Float_t[ntracks];
655 Int_t *pdg = new Int_t[ntracks];
657 for (Int_t part = 0; part < npart; part++)
659 tempPart = rl->Stack()->Particle(part);
660 if (tempPart->GetStatusCode() != 1) continue;
661 if ((!fPythiaComparison)&&(tempPart->Eta() > fEtaMax || tempPart->Eta() < fEtaMin ||
662 tempPart->Phi() > fPhiMax || tempPart->Phi() < fPhiMin) ){
663 if (fDebug>10) Info("FillPartonTracks","Excluding particle not pointing at the EMCAL");
666 if (tempPart->GetStatusCode() != 1) continue;
667 getpdg = tempPart->GetPDG();
668 if (getpdg->Charge() == 0.0 ) {
669 if (fDebug>10) Info("FillPartonTracks","Excluding a neutral particle");
672 if ( (parton->Eta() - tempPart->Eta())*(parton->Eta() - tempPart->Eta()) +
673 (parton->Phi() - tempPart->Phi())*(parton->Phi() - tempPart->Phi()) < 1.0 )
675 energy[ntracks] = tempPart->Pt();
676 eta[ntracks] = tempPart->Eta();
677 phi[ntracks] = tempPart->Phi();
678 pdg[ntracks] = tempPart->GetPdgCode();
682 parton->SetTrackList(ntracks,energy,eta,phi,pdg);