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 "AliGenEventHeader.h"
46 #include "AliGenPythiaEventHeader.h"
47 #include "AliGenerator.h"
48 #include "AliHeader.h"
52 ClassImp(AliEMCALJetFinderInputSimPrep)
54 AliEMCALJetFinderInputSimPrep::AliEMCALJetFinderInputSimPrep()
56 // Default constructor
57 if (fDebug > 0) Info("AliEMCALJetFinderInputSimPrep","Beginning Constructor");
60 fSmearType = kSmearEffic;
62 fTrackType = kCharged;
66 AliEMCALJetFinderInputSimPrep::~AliEMCALJetFinderInputSimPrep()
68 if (fDebug > 0) Info("~AliEMCALJetFinderInputSimPrep","Beginning Destructor");
72 void AliEMCALJetFinderInputSimPrep::Reset(AliEMCALJetFinderResetType_t resettype)
74 // Method to reset data
75 if (fDebug > 1) Info("Reset","Beginning Reset");
79 fInputObject.Reset(resettype);
82 fInputObject.Reset(resettype);
85 fInputObject.Reset(resettype);
87 case kResetParameters:
88 fSmearType = kSmearEffic;
90 fTrackType = kCharged;
93 fSmearType = kSmearEffic;
95 fTrackType = kCharged;
96 fInputObject.Reset(resettype);
99 Warning("FillFromFile", "kResetPartons not implemented") ;
101 case kResetParticles:
102 Warning("FillFromFile", "kResetParticles not implemented") ;
105 Warning("FillFromFile", "kResetJets not implemented") ;
111 Int_t AliEMCALJetFinderInputSimPrep::FillFromFile(TString *filename, AliEMCALJetFinderFileType_t filetype,Int_t EventNumber)
113 // gObjectTable->Print();
114 // Test that file exists - Getter doesn't like bogus filenames
115 if (fDebug > 1) Info("FillFromFile","Beginning FillFromFile");
116 fFileType = filetype;
117 TFile file(filename->Data());
119 Error("FillFromFile","Could not open file in FillFromFile");
124 gAlice = static_cast<AliRun *>(file.Get("gAlice"));
125 AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
126 if (gAlice->GetEvent(EventNumber) < 0){
127 Error("FillFromFile","Could not open event in FillFromFile");
131 AliEMCALGetter *gime = AliEMCALGetter::Instance(filename->Data());
132 gime->Event(EventNumber) ;
134 if ( fEMCALType == kHits ||
135 fEMCALType == kTimeCut ) FillHits();
136 if ( fTrackType != kNoTracks ) FillTracks();
137 if ( fFileType != kData){
142 // gObjectTable->Print();
147 void AliEMCALJetFinderInputSimPrep::FillHits()
149 // Fill from the hits to input object from simulation
150 if (fDebug > 1) Info("FillHits","Beginning FillHits");
152 // Access hit information
154 AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
155 Info("AliEMCALJetFinderInputSimPrep","Title of the geometry is %s",pEMCAL->GetTitle());
156 AliEMCALGeometry* geom = AliEMCALGetter::Instance()->EMCALGeometry();
157 // Float_t samplingF = geom->GetSampling();
158 Float_t samplingF = 11.6;
159 Info("AliEMCALJetFinderInputSimPrep","Sampling fraction is %f",samplingF);
160 TTree *treeH = AliEMCALGetter::Instance()->TreeH();
161 Int_t ntracks = (Int_t) treeH->GetEntries();
168 for (Int_t track=0; track<ntracks;track++) {
170 nbytes += treeH->GetEvent(track);
174 for(mHit=(AliEMCALHit*) pEMCAL->FirstHit(-1);
176 mHit=(AliEMCALHit*) pEMCAL->NextHit())
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 = samplingF*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
210 if (fDebug > 1) Info("FillTracks","Beginning FillTracks");
212 TParticlePDG* pdgP = 0;
214 Int_t npart = (gAlice->GetHeader())->GetNprimary();
215 Float_t bfield,rEMCAL;
217 if (fDebug > 1) Info("FillTracks","Defining the geometry");
219 AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
220 AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
221 fEtaMax = geom->GetArm1EtaMax();
222 fEtaMin = geom->GetArm1EtaMin();
223 fPhiMax = TMath::Pi()/180.0*geom->GetArm1PhiMax();
224 fPhiMin = TMath::Pi()/180.0*geom->GetArm1PhiMin();
226 if (fDebug > 1) Info("FillTracks","Starting particle loop");
228 for (Int_t part = 0; part < npart; part++) {
229 mPart = gAlice->GetMCApp()->Particle(part);
230 //if (part%10) gObjectTable->Print();
231 pdgP = mPart->GetPDG();
233 if (fDebug > 5) Info("FillTracks","Checking if track is a primary");
235 if (fFileType == kPythia) {
236 if (mPart->GetStatusCode() != 1) continue;
237 } else if (fFileType == kHijing) {
238 if (mPart->GetFirstDaughter() >= 0 && mPart->GetFirstDaughter() < npart) continue;
241 if (fDebug > 15) Info("FillTracks","Checking if track (eta - %f, phi - %f) is in acceptance",mPart->Eta(),mPart->Phi());
242 if (fDebug > 10) Info("FillTracks","Checking if EMCAL acceptance ( %f < eta < %f, %f < phi < %f) is in acceptance",fEtaMin,fEtaMax,fPhiMin,fPhiMax);
244 if (mPart->Eta() > fEtaMax || mPart->Eta() < fEtaMin) continue;
245 if (mPart->Phi() > fPhiMax || mPart->Phi() < fPhiMin ) continue;
248 {kProton, kProtonBar, kElectron, kPositron,
249 kNuE, kNuEBar, kGamma, kNeutron, kNeutronBar,
250 kMuonPlus, kMuonMinus, kK0Long , kPiPlus, kPiMinus,
251 kKPlus, kKMinus, kLambda0, kLambda0Bar, kK0Short,
252 kSigmaMinus, kSigmaPlus, kSigma0, kPi0, kK0, kK0Bar,
256 if (fDebug > 15) Info("FillTracks","Checking if track is rejected");
258 if ((fSmearType == kEfficiency ||
259 fSmearType == kSmearEffic) &&
261 if (AliEMCALFast::RandomReject(fEfficiency)) {
266 bfield = gAlice->Field()->SolenoidField();
267 rEMCAL = AliEMCALGeometry::GetInstance()->GetIPDistance();
268 Float_t rB = 3335.6 * mPart->Pt() / bfield; // [cm] (case of |charge|=1)
269 if (2.*rB < rEMCAL) continue; // track curls away
271 //if (part%10) gObjectTable->Print();
275 case kAll: // All Stable particles to be included
276 if (fDebug > 5) Info("FillTracks","Storing track");
277 if (fSmearType == kSmear ||
278 fSmearType == kSmearEffic ){
280 TParticle *tmp = Smear(MPart);
281 fInputObject.AddTrack(Smear(MPart));
284 fInputObject.AddTrack(*mPart);
287 case kEM: // All Electromagnetic particles
288 if (mPart->GetPdgCode() == kElectron ||
289 mPart->GetPdgCode() == kMuonPlus ||
290 mPart->GetPdgCode() == kMuonMinus ||
291 mPart->GetPdgCode() == kPositron ){
292 if (fDebug > 5) Info("FillTracks","Storing electron or positron");
293 if (fSmearType == kSmear ||
294 fSmearType == kSmearEffic ){
296 TParticle *tmp = Smear(MPart);
297 fInputObject.AddTrack(tmp);
300 fInputObject.AddTrack(*mPart);
303 if ( mPart->GetPdgCode() == kGamma ){
304 fInputObject.AddTrack(*mPart);
305 if (fDebug > 5) Info("FillTracks","Storing gamma");
309 case kCharged: // All Particles with non-zero charge
310 if (pdgP->Charge() != 0) {
311 if (fDebug > 5) Info("FillTracks","Storing charged track");
312 if (fSmearType == kSmear ||
313 fSmearType == kSmearEffic ){
315 TParticle *tmp = Smear(MPart);
316 fInputObject.AddTrack(tmp);
319 fInputObject.AddTrack(*mPart);
323 case kNeutral: // All particles with no charge
324 if (pdgP->Charge() == 0){
325 fInputObject.AddTrack(*mPart);
326 if (fDebug > 5) Info("FillTracks","Storing neutral");
329 case kHadron: //All hadrons
330 if (mPart->GetPdgCode() != kElectron &&
331 mPart->GetPdgCode() != kPositron &&
332 mPart->GetPdgCode() != kMuonPlus &&
333 mPart->GetPdgCode() != kMuonMinus &&
334 mPart->GetPdgCode() != kGamma )
336 if (fDebug > 5) Info("FillTracks","Storing hadron");
337 if (pdgP->Charge() == 0){
338 fInputObject.AddTrack(*mPart);
340 if (fSmearType == kSmear ||
341 fSmearType == kSmearEffic ){
343 TParticle *tmp = Smear(MPart);
344 fInputObject.AddTrack(tmp);
347 fInputObject.AddTrack(*mPart);
352 case kChargedHadron: // only charged hadrons
353 if (mPart->GetPdgCode() != kElectron &&
354 mPart->GetPdgCode() != kPositron &&
355 mPart->GetPdgCode() != kGamma &&
356 mPart->GetPdgCode() != kMuonPlus &&
357 mPart->GetPdgCode() != kMuonMinus &&
358 pdgP->Charge() != 0 ){
359 if (fDebug > 5) Info("FillTracks","Storing charged hadron");
360 if (fSmearType == kSmear ||
361 fSmearType == kSmearEffic ){
363 TParticle *tmp = Smear(MPart);
364 fInputObject.AddTrack(tmp);
367 fInputObject.AddTrack(*mPart);
377 // Info("FillTracks","After Particle Storage");
378 //if (part%10) gObjectTable->Print();
380 } //end of particle loop
381 //gObjectTable->Print();
384 void AliEMCALJetFinderInputSimPrep::FillPartons()
387 // input object from simulation
388 if (fDebug > 1) Info("FillPartons","Beginning FillPartons");
390 AliGenEventHeader* evHeader = ((AliHeader*)(gAlice->GetHeader()))->GenEventHeader();
391 Float_t triggerJetValues[4];
392 AliEMCALParton tempParton;
394 if (fFileType == kPythia)
396 Int_t ntriggerjets = ((AliGenPythiaEventHeader*)evHeader)->NTriggerJets();
397 if (fDebug > 1) Info("FillPartons","Number of trigger jets --> %i",ntriggerjets);
398 for (Int_t jetcounter = 0 ; jetcounter < ntriggerjets; jetcounter++){
399 ((AliGenPythiaEventHeader*)evHeader)->TriggerJet(jetcounter,triggerJetValues);
400 TLorentzVector tempLParton;
401 tempLParton.SetPxPyPzE(triggerJetValues[0],triggerJetValues[1],triggerJetValues[2],triggerJetValues[3]);
402 tempParton.SetEnergy(tempLParton.Energy());
403 tempParton.SetEta(tempLParton.Eta());
404 tempParton.SetPhi(tempLParton.Phi());
405 FillPartonTracks(&tempParton);
406 fInputObject.AddParton(&tempParton);
412 void AliEMCALJetFinderInputSimPrep::FillParticles()
414 // Fill particles to input object from simulation
415 if (fDebug > 1) Info("FillParticles","Beginning FillParticles");
417 Int_t npart = (gAlice->GetHeader())->GetNprimary();
418 TParticlePDG* pdgP = 0;
420 for (Int_t part = 0; part < npart; part++) {
421 TParticle *mPart = gAlice->GetMCApp()->Particle(part);
422 pdgP = mPart->GetPDG();
424 if (fDebug > 10) Info("FillParticles","Checking if particle is a primary");
426 if (fFileType == kPythia) {
427 if (mPart->GetStatusCode() != 1) continue;
428 } else if (fFileType == kHijing) {
429 if (mPart->GetFirstDaughter() >= 0 && mPart->GetFirstDaughter() < npart) continue;
432 if (fDebug > 10) Info("FillParticles","Checking if particle is in acceptance");
434 if (mPart->Eta() > fEtaMax || mPart->Eta() < fEtaMin) continue;
435 if (mPart->Phi() > fPhiMax || mPart->Phi() < fPhiMin ) continue;
439 {kProton, kProtonBar, kElectron, kPositron,
440 kNuE, kNuEBar, kGamma, kNeutron, kNeutronBar,
441 kMuonPlus, kMuonMinus, kK0Long , kPiPlus, kPiMinus,
442 kKPlus, kKMinus, kLambda0, kLambda0Bar, kK0Short,
443 kSigmaMinus, kSigmaPlus, kSigma0, kPi0, kK0, kK0Bar,
449 case kAll: // All Stable particles to be included
450 if (fDebug > 5) Info("FillParticles","Storing particle");
451 fInputObject.AddParticle(mPart);
453 case kEM: // All Electromagnetic particles
454 if (mPart->GetPdgCode() == kElectron ||
455 mPart->GetPdgCode() == kPositron ||
456 mPart->GetPdgCode() == kGamma){
457 if (fDebug > 5) Info("FillParticles","Storing electromagnetic particle");
458 fInputObject.AddParticle(mPart);
461 case kCharged: // All Particles with non-zero charge
462 if (pdgP->Charge() != 0) {
463 if (fDebug > 5) Info("FillParticles","Storing charged particle");
464 fInputObject.AddParticle(mPart);
467 case kNeutral: // All particles with no charge
468 if (pdgP->Charge() == 0){
469 if (fDebug > 5) Info("FillParticles","Storing neutral particle");
470 fInputObject.AddParticle(mPart);
473 case kHadron: //All hadrons
474 if (mPart->GetPdgCode() != kElectron &&
475 mPart->GetPdgCode() != kPositron &&
476 mPart->GetPdgCode() != kGamma )
479 if (fDebug > 5) Info("FillParticles","Storing hadron");
480 fInputObject.AddParticle(mPart);
483 case kChargedHadron: // only charged hadrons
484 if (mPart->GetPdgCode() != kElectron &&
485 mPart->GetPdgCode() != kPositron &&
486 mPart->GetPdgCode() != kGamma &&
487 pdgP->Charge() != 0 ){
488 if (fDebug > 5) Info("FillParticles","Storing charged hadron");
489 fInputObject.AddParticle(mPart);
497 }// end of loop over particles
500 void AliEMCALJetFinderInputSimPrep::FillDigits()
502 // Fill digits to input object
506 void AliEMCALJetFinderInputSimPrep::Smear(TParticle *particle)
508 // Smear particle momentum
509 if (fDebug > 5) Info("Smear","Beginning Smear");
511 Float_t tmp = AliEMCALFast::SmearMomentum(1,particle->P());
512 Float_t tmpt = (tmp/particle->P()) * particle->Pt();
513 if (fDebug > 10) Info("Smear","Smearing changed P from %f to %f",particle->Pt(),tmpt);
514 TLorentzVector tmpvector;
515 tmpvector.SetPtEtaPhiM(tmpt,particle->Eta(), particle->Phi(), particle->GetMass());
516 // create a new particle with smearing momentum - and no daughter or parent information and the same
519 TParticle tmparticle(particle->GetPdgCode(), 1, 0, 0,0,0,
520 tmpvector.Px() , tmpvector.Py(), tmpvector.Pz(), tmpvector.Energy(),
521 particle->Vx(), particle->Vy(), particle->Vz(), particle->T());
522 fInputObject.AddTrack(tmparticle);
526 void AliEMCALJetFinderInputSimPrep::FillPartonTracks(AliEMCALParton *parton)
528 // Populate parton tracks for input distributions
529 if (fDebug>1) Info("FillPartonTracks","Beginning FillPartonTracks");
530 Int_t npart = (gAlice->GetHeader())->GetNprimary();
532 TParticlePDG *getpdg;
534 for (Int_t part = 0; part < npart; part++)
536 tempPart = gAlice->GetMCApp()->Particle(part);
537 if (tempPart->GetStatusCode() != 1) continue;
538 if (tempPart->Eta() > fEtaMax || tempPart->Eta() < fEtaMin ||
539 tempPart->Phi() > fPhiMax || tempPart->Phi() < fPhiMin ){
540 if (fDebug>10) Info("FillPartonTracks","Excluding particle not pointing at the EMCAL");
543 getpdg = tempPart->GetPDG();
544 if (getpdg->Charge() == 0.0 ) {
545 if (fDebug>10) Info("FillPartonTracks","Excluding a neutral particle");
548 if ( (parton->Eta() - tempPart->Eta())*(parton->Eta() - tempPart->Eta()) +
549 (parton->Phi() - tempPart->Phi())*(parton->Phi() - tempPart->Phi()) < 1.0 ) ntracks++;
551 Float_t *energy = new Float_t[ntracks];
552 Float_t *eta = new Float_t[ntracks];
553 Float_t *phi = new Float_t[ntracks];
554 Int_t *pdg = new Int_t[ntracks];
556 for (Int_t part = 0; part < npart; part++)
558 tempPart = gAlice->GetMCApp()->Particle(part);
559 if (tempPart->GetStatusCode() != 1) continue;
560 if (tempPart->Eta() > fEtaMax || tempPart->Eta() < fEtaMin ||
561 tempPart->Phi() > fPhiMax || tempPart->Phi() < fPhiMin ){
562 if (fDebug>10) Info("FillPartonTracks","Excluding particle not pointing at the EMCAL");
565 if (tempPart->GetStatusCode() != 1) continue;
566 getpdg = tempPart->GetPDG();
567 if (getpdg->Charge() == 0.0 ) {
568 if (fDebug>10) Info("FillPartonTracks","Excluding a neutral particle");
571 if ( (parton->Eta() - tempPart->Eta())*(parton->Eta() - tempPart->Eta()) +
572 (parton->Phi() - tempPart->Phi())*(parton->Phi() - tempPart->Phi()) < 1.0 )
574 energy[ntracks] = tempPart->Pt();
575 eta[ntracks] = tempPart->Eta();
576 phi[ntracks] = tempPart->Phi();
577 pdg[ntracks] = tempPart->GetPdgCode();
581 parton->SetTrackList(ntracks,energy,eta,phi,pdg);