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>
40 #include "AliEMCALFast.h"
42 #include "AliEMCALHit.h"
43 #include "AliEMCALGeometry.h"
44 #include "AliEMCALGetter.h"
45 #include "AliEMCALLoader.h"
46 #include "AliGenEventHeader.h"
47 #include "AliGenPythiaEventHeader.h"
48 #include "AliGenerator.h"
49 #include "AliHeader.h"
54 ClassImp(AliEMCALJetFinderInputSimPrep)
56 AliEMCALJetFinderInputSimPrep::AliEMCALJetFinderInputSimPrep()
58 // Default constructor
59 if (fDebug > 0) Info("AliEMCALJetFinderInputSimPrep","Beginning Constructor");
62 fSmearType = kSmearEffic;
64 fTrackType = kCharged;
68 AliEMCALJetFinderInputSimPrep::~AliEMCALJetFinderInputSimPrep()
70 if (fDebug > 0) Info("~AliEMCALJetFinderInputSimPrep","Beginning Destructor");
74 void AliEMCALJetFinderInputSimPrep::Reset(AliEMCALJetFinderResetType_t resettype)
76 // Method to reset data
77 if (fDebug > 1) Info("Reset","Beginning Reset");
81 fInputObject.Reset(resettype);
84 fInputObject.Reset(resettype);
87 fInputObject.Reset(resettype);
89 case kResetParameters:
90 fSmearType = kSmearEffic;
92 fTrackType = kCharged;
95 fSmearType = kSmearEffic;
97 fTrackType = kCharged;
98 fInputObject.Reset(resettype);
101 Warning("FillFromFile", "kResetPartons not implemented") ;
103 case kResetParticles:
104 Warning("FillFromFile", "kResetParticles not implemented") ;
107 Warning("FillFromFile", "kResetJets not implemented") ;
113 Int_t AliEMCALJetFinderInputSimPrep::FillFromFile(TString *filename, AliEMCALJetFinderFileType_t filetype,Int_t EventNumber)
115 if (fDebug > 1) Info("FillFromFile","Beginning FillFromFile");
116 fFileType = filetype;
117 AliEMCALGetter *gime = AliEMCALGetter::Instance(filename->Data());
118 if (fDebug > 1) Info("FillFromFile","Instantiated Getter with %s as the file",filename->Data());
119 gime->Event(EventNumber,"XH") ;
120 if (fDebug > 1) Info("FillFromFile","Got event %i with option \"XH\"",EventNumber);
122 if ( fEMCALType == kHits ||
123 fEMCALType == kTimeCut ) FillHits();
124 if ( fTrackType != kNoTracks ) FillTracks();
125 if ( fFileType != kData){
128 /* gime->EmcalLoader()->UnloadRecParticles();
129 gime->EmcalLoader()->UnloadTracks();
134 void AliEMCALJetFinderInputSimPrep::FillHits()
136 // Fill from the hits to input object from simulation
137 if (fDebug > 1) Info("FillHits","Beginning FillHits");
139 // Access hit information
141 // AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
142 // Info("AliEMCALJetFinderInputSimPrep","Title of the geometry is %s",pEMCAL->GetTitle());
143 AliEMCALGeometry* geom = AliEMCALGetter::Instance()->EMCALGeometry();
144 AliEMCALGetter *gime = AliEMCALGetter::Instance();
146 // TTree *treeH = AliEMCALGetter::Instance()->TreeH();
147 // Int_t ntracks = (Int_t) treeH->GetEntries();
153 /* for (Int_t track=0; track<ntracks;track++) {
155 nbytes += treeH->GetEvent(track);*/
159 for(Int_t i =0; i < gime->Hits()->GetEntries() ;i++)
162 if (fEMCALType == kTimeCut &&
163 (mHit->GetTime()>fTimeCut) ) continue;
164 Float_t x = mHit->X(); // x-pos of hit
165 Float_t y = mHit->Y(); // y-pos
166 Float_t z = mHit->Z(); // z-pos
167 Float_t eloss = mHit->GetEnergy(); // deposited energy
169 Float_t r = TMath::Sqrt(x*x+y*y);
170 Float_t theta = TMath::ATan2(r,z);
171 Float_t eta = -TMath::Log(TMath::Tan(theta/2.));
172 Float_t phi = TMath::ATan2(y,x);
173 etH = eloss*TMath::Sin(theta);
174 if (fDebug > 10) Info("FillHits","Values of hit energy %i",Int_t(1e7*etH));
175 if (geom->TowerIndexFromEtaPhi(eta,180.0/TMath::Pi()*phi) == -1)
179 Error("FillHits","Hit not inside EMCAL!");
184 fInputObject.AddEnergyToDigit(geom->TowerIndexFromEtaPhi(eta,180.0/TMath::Pi()*phi)-1,Int_t(1e7*etH));
190 void AliEMCALJetFinderInputSimPrep::FillTracks()
192 // Fill from particles simulating a TPC to input object from simulation
194 if (fDebug > 1) Info("FillTracks","Beginning FillTracks");
196 TParticlePDG* pdgP = 0;
198 Int_t npart = (gAlice->GetHeader())->GetNprimary();
199 if (fDebug > 1) Info("FillTracks","Header says there are %i primaries",npart);
200 //AliEMCALGetter *gime = AliEMCALGetter::Instance();
201 Float_t bfield,rEMCAL;
202 /* for (Int_t i =0; i<gime->NPrimaries(); i++)
204 cout<<gime->Primary(i)->GetFirstMother()<<endl;
205 if (gime->Primary(i)->GetFirstMother() == -1) {
207 gime->Primary(i)->Dump();
211 if (fDebug > 1) Info("FillTracks","Defining the geometry");
213 AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
214 AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
215 fEtaMax = geom->GetArm1EtaMax();
216 fEtaMin = geom->GetArm1EtaMin();
217 fPhiMax = TMath::Pi()/180.0*geom->GetArm1PhiMax();
218 fPhiMin = TMath::Pi()/180.0*geom->GetArm1PhiMin();
220 if (fDebug > 1) Info("FillTracks","Starting particle loop");
222 if (fDebug > 1) Info("FillTracks","Particle loop of %i",npart);
223 for (Int_t part = 0; part < npart; part++) {
224 mPart = gAlice->GetMCApp()->Particle(part);
225 //if (part%10) gObjectTable->Print();
226 pdgP = mPart->GetPDG();
228 if (fDebug > 5) Info("FillTracks","Checking if track is a primary");
230 if (fFileType == kPythia) {
231 if (mPart->GetStatusCode() != 1) continue;
232 } else if (fFileType == kHijing) {
233 if (mPart->GetFirstDaughter() >= 0 && mPart->GetFirstDaughter() < npart) continue;
236 if (fDebug > 15) Info("FillTracks","Checking if track (eta - %f, phi - %f) is in acceptance",mPart->Eta(),mPart->Phi());
237 if (fDebug > 10) Info("FillTracks","Checking if EMCAL acceptance ( %f < eta < %f, %f < phi < %f) is in acceptance",fEtaMin,fEtaMax,fPhiMin,fPhiMax);
239 if ((!fPythiaComparison)&&(mPart->Eta() > fEtaMax || mPart->Eta() < fEtaMin)) continue;
240 if ((!fPythiaComparison)&&(mPart->Phi() > fPhiMax || mPart->Phi() < fPhiMin)) continue;
243 {kProton, kProtonBar, kElectron, kPositron,
244 kNuE, kNuEBar, kGamma, kNeutron, kNeutronBar,
245 kMuonPlus, kMuonMinus, kK0Long , kPiPlus, kPiMinus,
246 kKPlus, kKMinus, kLambda0, kLambda0Bar, kK0Short,
247 kSigmaMinus, kSigmaPlus, kSigma0, kPi0, kK0, kK0Bar,
251 if (fDebug > 15) Info("FillTracks","Checking if track is rejected");
253 if ((fSmearType == kEfficiency ||
254 fSmearType == kSmearEffic) &&
256 if (AliEMCALFast::RandomReject(fEfficiency)) {
261 bfield = gAlice->Field()->SolenoidField();
262 rEMCAL = AliEMCALGeometry::GetInstance()->GetIPDistance();
263 Float_t rB = 3335.6 * mPart->Pt() / bfield; // [cm] (case of |charge|=1)
264 if (2.*rB < rEMCAL) continue; // track curls away
266 //if (part%10) gObjectTable->Print();
270 case kAll: // All Stable particles to be included
271 if (fDebug > 5) Info("FillTracks","Storing track");
272 if (fSmearType == kSmear ||
273 fSmearType == kSmearEffic ){
275 TParticle *tmp = Smear(MPart);
276 fInputObject.AddTrack(Smear(MPart));
279 fInputObject.AddTrack(*mPart);
282 case kEM: // All Electromagnetic particles
283 if (mPart->GetPdgCode() == kElectron ||
284 mPart->GetPdgCode() == kMuonPlus ||
285 mPart->GetPdgCode() == kMuonMinus ||
286 mPart->GetPdgCode() == kPositron ){
287 if (fDebug > 5) Info("FillTracks","Storing electron or positron");
288 if (fSmearType == kSmear ||
289 fSmearType == kSmearEffic ){
291 TParticle *tmp = Smear(MPart);
292 fInputObject.AddTrack(tmp);
295 fInputObject.AddTrack(*mPart);
298 if ( mPart->GetPdgCode() == kGamma ){
299 fInputObject.AddTrack(*mPart);
300 if (fDebug > 5) Info("FillTracks","Storing gamma");
304 case kCharged: // All Particles with non-zero charge
305 if (pdgP->Charge() != 0) {
306 if (fDebug > 5) Info("FillTracks","Storing charged track");
307 if (fSmearType == kSmear ||
308 fSmearType == kSmearEffic ){
310 TParticle *tmp = Smear(MPart);
311 fInputObject.AddTrack(tmp);
314 fInputObject.AddTrack(*mPart);
318 case kNeutral: // All particles with no charge
319 if (pdgP->Charge() == 0){
320 fInputObject.AddTrack(*mPart);
321 if (fDebug > 5) Info("FillTracks","Storing neutral");
324 case kHadron: //All hadrons
325 if (mPart->GetPdgCode() != kElectron &&
326 mPart->GetPdgCode() != kPositron &&
327 mPart->GetPdgCode() != kMuonPlus &&
328 mPart->GetPdgCode() != kMuonMinus &&
329 mPart->GetPdgCode() != kGamma )
331 if (fDebug > 5) Info("FillTracks","Storing hadron");
332 if (pdgP->Charge() == 0){
333 fInputObject.AddTrack(*mPart);
335 if (fSmearType == kSmear ||
336 fSmearType == kSmearEffic ){
338 TParticle *tmp = Smear(MPart);
339 fInputObject.AddTrack(tmp);
342 fInputObject.AddTrack(*mPart);
347 case kChargedHadron: // only charged hadrons
348 if (mPart->GetPdgCode() != kElectron &&
349 mPart->GetPdgCode() != kPositron &&
350 mPart->GetPdgCode() != kGamma &&
351 mPart->GetPdgCode() != kMuonPlus &&
352 mPart->GetPdgCode() != kMuonMinus &&
353 pdgP->Charge() != 0 ){
354 if (fDebug > 5) Info("FillTracks","Storing charged hadron");
355 if (fSmearType == kSmear ||
356 fSmearType == kSmearEffic ){
358 TParticle *tmp = Smear(MPart);
359 fInputObject.AddTrack(tmp);
362 fInputObject.AddTrack(*mPart);
369 if (pdgP->Charge() != 0 || mPart->GetPdgCode() == kPi0 ||
370 mPart->GetPdgCode() == kGamma )
372 if (fDebug > 5) Info("FillTracks","Storing charged track");
373 if (fSmearType == kSmear ||
374 fSmearType == kSmearEffic ){
376 TParticle *tmp = Smear(MPart);
377 fInputObject.AddTrack(tmp);
380 fInputObject.AddTrack(*mPart);
388 // Info("FillTracks","After Particle Storage");
389 //if (part%10) gObjectTable->Print();
391 } //end of particle loop
392 //gObjectTable->Print();
395 void AliEMCALJetFinderInputSimPrep::FillPartons()
398 // input object from simulation
399 if (fDebug > 1) Info("FillPartons","Beginning FillPartons");
401 AliGenEventHeader* evHeader = ((AliHeader*)(gAlice->GetHeader()))->GenEventHeader();
402 Float_t triggerJetValues[4];
403 AliEMCALParton tempParton;
405 if (fFileType == kPythia)
407 Int_t ntriggerjets = ((AliGenPythiaEventHeader*)evHeader)->NTriggerJets();
408 if (fDebug > 1) Info("FillPartons","Number of trigger jets --> %i",ntriggerjets);
409 for (Int_t jetcounter = 0 ; jetcounter < ntriggerjets; jetcounter++){
410 ((AliGenPythiaEventHeader*)evHeader)->TriggerJet(jetcounter,triggerJetValues);
411 TLorentzVector tempLParton;
412 tempLParton.SetPxPyPzE(triggerJetValues[0],triggerJetValues[1],triggerJetValues[2],triggerJetValues[3]);
413 tempParton.SetEnergy(tempLParton.Energy());
414 tempParton.SetEta(tempLParton.Eta());
415 tempParton.SetPhi(tempLParton.Phi());
416 FillPartonTracks(&tempParton);
417 fInputObject.AddParton(&tempParton);
423 void AliEMCALJetFinderInputSimPrep::FillParticles()
425 // Fill particles to input object from simulation
426 if (fDebug > 1) Info("FillParticles","Beginning FillParticles");
428 Int_t npart = (gAlice->GetHeader())->GetNprimary();
429 TParticlePDG* pdgP = 0;
431 for (Int_t part = 0; part < npart; part++) {
432 TParticle *mPart = gAlice->GetMCApp()->Particle(part);
433 pdgP = mPart->GetPDG();
435 if (fDebug > 10) Info("FillParticles","Checking if particle is a primary");
437 if (fFileType == kPythia) {
438 if (mPart->GetStatusCode() != 1) continue;
439 } else if (fFileType == kHijing) {
440 if (mPart->GetFirstDaughter() >= 0 && mPart->GetFirstDaughter() < npart) continue;
443 if (fDebug > 10) Info("FillParticles","Checking if particle is in acceptance");
445 if ((!fPythiaComparison)&&(mPart->Eta() > fEtaMax || mPart->Eta() < fEtaMin)) continue;
446 if ((!fPythiaComparison)&&(mPart->Phi() > fPhiMax || mPart->Phi() < fPhiMin)) continue;
450 {kProton, kProtonBar, kElectron, kPositron,
451 kNuE, kNuEBar, kGamma, kNeutron, kNeutronBar,
452 kMuonPlus, kMuonMinus, kK0Long , kPiPlus, kPiMinus,
453 kKPlus, kKMinus, kLambda0, kLambda0Bar, kK0Short,
454 kSigmaMinus, kSigmaPlus, kSigma0, kPi0, kK0, kK0Bar,
460 case kAll: // All Stable particles to be included
461 if (fDebug > 5) Info("FillParticles","Storing particle");
462 fInputObject.AddParticle(mPart);
464 case kEM: // All Electromagnetic particles
465 if (mPart->GetPdgCode() == kElectron ||
466 mPart->GetPdgCode() == kPositron ||
467 mPart->GetPdgCode() == kGamma){
468 if (fDebug > 5) Info("FillParticles","Storing electromagnetic particle");
469 fInputObject.AddParticle(mPart);
472 case kCharged: // All Particles with non-zero charge
473 if (pdgP->Charge() != 0) {
474 if (fDebug > 5) Info("FillParticles","Storing charged particle");
475 fInputObject.AddParticle(mPart);
478 case kNeutral: // All particles with no charge
479 if (pdgP->Charge() == 0){
480 if (fDebug > 5) Info("FillParticles","Storing neutral particle");
481 fInputObject.AddParticle(mPart);
484 case kHadron: //All hadrons
485 if (mPart->GetPdgCode() != kElectron &&
486 mPart->GetPdgCode() != kPositron &&
487 mPart->GetPdgCode() != kGamma )
490 if (fDebug > 5) Info("FillParticles","Storing hadron");
491 fInputObject.AddParticle(mPart);
494 case kChargedHadron: // only charged hadrons
495 if (mPart->GetPdgCode() != kElectron &&
496 mPart->GetPdgCode() != kPositron &&
497 mPart->GetPdgCode() != kGamma &&
498 pdgP->Charge() != 0 ){
499 if (fDebug > 5) Info("FillParticles","Storing charged hadron");
500 fInputObject.AddParticle(mPart);
506 if (pdgP->Charge() != 0 || mPart->GetPdgCode() == kPi0 ||
507 mPart->GetPdgCode() == kGamma )
509 if (fDebug > 5) Info("FillTracks","Storing charged track");
510 if (fSmearType == kSmear ||
511 fSmearType == kSmearEffic ){
513 TParticle *tmp = Smear(MPart);
514 fInputObject.AddTrack(tmp);
517 fInputObject.AddTrack(*mPart);
524 }// end of loop over particles
527 void AliEMCALJetFinderInputSimPrep::FillDigits()
529 // Fill digits to input object
533 void AliEMCALJetFinderInputSimPrep::Smear(TParticle *particle)
535 // Smear particle momentum
536 if (fDebug > 5) Info("Smear","Beginning Smear");
538 Float_t tmp = AliEMCALFast::SmearMomentum(1,particle->P());
539 Float_t tmpt = (tmp/particle->P()) * particle->Pt();
540 if (fDebug > 10) Info("Smear","Smearing changed P from %f to %f",particle->Pt(),tmpt);
541 TLorentzVector tmpvector;
542 tmpvector.SetPtEtaPhiM(tmpt,particle->Eta(), particle->Phi(), particle->GetMass());
543 // create a new particle with smearing momentum - and no daughter or parent information and the same
546 TParticle tmparticle(particle->GetPdgCode(), 1, 0, 0,0,0,
547 tmpvector.Px() , tmpvector.Py(), tmpvector.Pz(), tmpvector.Energy(),
548 particle->Vx(), particle->Vy(), particle->Vz(), particle->T());
549 fInputObject.AddTrack(tmparticle);
553 void AliEMCALJetFinderInputSimPrep::FillPartonTracks(AliEMCALParton *parton)
555 // Populate parton tracks for input distributions
556 if (fDebug>1) Info("FillPartonTracks","Beginning FillPartonTracks");
557 Int_t npart = (gAlice->GetHeader())->GetNprimary();
559 TParticlePDG *getpdg;
561 for (Int_t part = 0; part < npart; part++)
563 tempPart = gAlice->GetMCApp()->Particle(part);
564 if (tempPart->GetStatusCode() != 1) continue;
565 if ((!fPythiaComparison)&&(tempPart->Eta() > fEtaMax || tempPart->Eta() < fEtaMin ||
566 tempPart->Phi() > fPhiMax || tempPart->Phi() < fPhiMin) ){
567 if (fDebug>10) Info("FillPartonTracks","Excluding particle not pointing at the EMCAL");
570 getpdg = tempPart->GetPDG();
571 if (getpdg->Charge() == 0.0 ) {
572 if (fDebug>10) Info("FillPartonTracks","Excluding a neutral particle");
575 if ( (parton->Eta() - tempPart->Eta())*(parton->Eta() - tempPart->Eta()) +
576 (parton->Phi() - tempPart->Phi())*(parton->Phi() - tempPart->Phi()) < 1.0 ) ntracks++;
578 Float_t *energy = new Float_t[ntracks];
579 Float_t *eta = new Float_t[ntracks];
580 Float_t *phi = new Float_t[ntracks];
581 Int_t *pdg = new Int_t[ntracks];
583 for (Int_t part = 0; part < npart; part++)
585 tempPart = gAlice->GetMCApp()->Particle(part);
586 if (tempPart->GetStatusCode() != 1) continue;
587 if ((!fPythiaComparison)&&(tempPart->Eta() > fEtaMax || tempPart->Eta() < fEtaMin ||
588 tempPart->Phi() > fPhiMax || tempPart->Phi() < fPhiMin) ){
589 if (fDebug>10) Info("FillPartonTracks","Excluding particle not pointing at the EMCAL");
592 if (tempPart->GetStatusCode() != 1) continue;
593 getpdg = tempPart->GetPDG();
594 if (getpdg->Charge() == 0.0 ) {
595 if (fDebug>10) Info("FillPartonTracks","Excluding a neutral particle");
598 if ( (parton->Eta() - tempPart->Eta())*(parton->Eta() - tempPart->Eta()) +
599 (parton->Phi() - tempPart->Phi())*(parton->Phi() - tempPart->Phi()) < 1.0 )
601 energy[ntracks] = tempPart->Pt();
602 eta[ntracks] = tempPart->Eta();
603 phi[ntracks] = tempPart->Phi();
604 pdg[ntracks] = tempPart->GetPdgCode();
608 parton->SetTrackList(ntracks,energy,eta,phi,pdg);