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 **************************************************************************/
17 //======================================================================
19 // Fill Unit Array class
20 // Class used by AliJetESDReader to fill a UnitArray from the information extracted
21 // from the particle tracks
22 // Author: magali.estienne@ires.in2p3.fr
23 //======================================================================
26 // --- ROOT system ---
28 #include <TProcessID.h>
30 #include "AliJetUnitArray.h"
31 #include "AliJetHadronCorrectionv1.h"
32 #include "AliJetAODFillUnitArrayTracks.h"
34 // --- ROOT system ---
39 // --- AliRoot header files ---
42 ClassImp(AliJetAODFillUnitArrayTracks)
44 //_____________________________________________________________________________
45 AliJetAODFillUnitArrayTracks::AliJetAODFillUnitArrayTracks()
46 : AliJetFillUnitArray(),
49 fApplyMIPCorrection(kTRUE),
60 //_____________________________________________________________________________
61 AliJetAODFillUnitArrayTracks::AliJetAODFillUnitArrayTracks(AliAODEvent* aod)
62 : AliJetFillUnitArray(),
65 fApplyMIPCorrection(kTRUE),
76 //_____________________________________________________________________________
77 AliJetAODFillUnitArrayTracks::AliJetAODFillUnitArrayTracks(const AliJetAODFillUnitArrayTracks &det)
78 : AliJetFillUnitArray(det),
79 fNumUnits(det.fNumUnits),
80 fHadCorr(det.fHadCorr),
81 fApplyMIPCorrection(det.fApplyMIPCorrection),
92 //_____________________________________________________________________________
93 AliJetAODFillUnitArrayTracks& AliJetAODFillUnitArrayTracks::operator=(const AliJetAODFillUnitArrayTracks& other)
97 fNumUnits = other.fNumUnits;
98 fHadCorr = other.fHadCorr;
99 fApplyMIPCorrection = other.fApplyMIPCorrection;
101 fGrid0 = other.fGrid0;
102 fGrid1 = other.fGrid1;
103 fGrid2 = other.fGrid2;
104 fGrid3 = other.fGrid3;
105 fGrid4 = other.fGrid4;
110 //____________________________________________________________________________
111 void AliJetAODFillUnitArrayTracks::InitParameters()
113 // fHadCorr = 0; // For hadron correction
114 fNumUnits = fGeom->GetNCells(); // Number of towers in EMCAL
116 fTPCGrid->GetAccParam(fNphi,fNeta,fPhiMin,
117 fPhiMax,fEtaMin,fEtaMax);
118 fTPCGrid->GetBinParam(fPhiBinInTPCAcc,fEtaBinInTPCAcc,
119 fPhiBinInEMCalAcc,fEtaBinInEMCalAcc,fNbinPhi);
121 fIndex = fTPCGrid->GetIndexObject();
124 for(Int_t i=0; i<fNphi+1; i++)
125 for(Int_t j=0; j<fNeta+1; j++) {cout << "fIndex[" << i << "," << j << "] : " <<
126 (*fIndex)(i,j) << endl; }
128 if(fDebug>1) printf("\n Parameters initiated ! \n");
131 //_____________________________________________________________________________
132 AliJetAODFillUnitArrayTracks::~AliJetAODFillUnitArrayTracks()
137 //_____________________________________________________________________________
138 void AliJetAODFillUnitArrayTracks::Exec(Option_t* const /*option*/)
139 //void AliJetAODFillUnitArrayTracks::Exec(Option_t* option)
143 // Fill the unit array with the charged particle information in AOD
146 fDebug = fReaderHeader->GetDebug();
152 // get number of tracks in event (for the loop)
155 //(not used ?) Int_t nt2 = 0;
161 nt = fAOD->GetNumberOfTracks();
162 if(fDebug>1) cout << "Number of Tracks in AOD : " << nt << endl;
164 // temporary storage of signal and pt cut flag
165 Int_t* sflag = new Int_t[nt];
166 Int_t* cflag = new Int_t[nt];
168 // get cuts set by user
169 Float_t ptMin = fReaderHeader->GetPtCut();
170 Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
171 Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
172 fOpt = fReaderHeader->GetDetector();
173 fDZ = fReaderHeader->GetDZ();
174 UInt_t filterMask = ((AliJetAODReaderHeader*)fReaderHeader)->GetTestFilterMask();
176 Int_t nTracksEmcal = 0;
177 Int_t nTracksEmcalDZ = 0;
178 Int_t nTracksTpc = 0;
179 Int_t nTracksTpcOnly = 0;
180 Int_t nTracksEmcalCut = 0;
181 Int_t nTracksEmcalDZCut = 0;
182 Int_t nTracksTpcCut = 0;
183 Int_t nTracksTpcOnlyCut = 0;
185 fGrid = fTPCGrid->GetGridType();
189 for (Int_t it = 0; it < nmax; it++) {
191 track = fAOD->GetTrack(it);
192 UInt_t status = track->GetStatus();
195 track->GetPxPyPz(mom);
197 p3.SetXYZ(mom[0],mom[1],mom[2]);
201 phi = ( (p3.Phi()) < 0) ? (p3.Phi()) + 2. * TMath::Pi() : (p3.Phi());
202 if (status == 0) continue;
203 if((filterMask>0)&&!(track->TestFilterBit(filterMask)))continue;
205 if ( (eta > etaMax) || (eta < etaMin)) continue; // checking eta cut
207 if (goodTrack == 0) new(fRef) TRefArray(TProcessID::GetProcessWithUID(track));
210 if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
212 if (pt > ptMin) cflag[goodTrack]=1; // pt cut
217 // Only TPC filled from its grid in its total acceptance
219 Int_t idTPC = fTPCGrid->GetIndex(phi,eta);
223 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(idTPC-1);
224 TRefArray *reference = uArray->GetUnitTrackRef();
225 if (reference->GetEntries() == 0) {
226 new(reference) TRefArray(TProcessID::GetProcessWithUID(track));
229 reference->Add(track);
231 Float_t unitEnergy = 0.;
232 unitEnergy = uArray->GetUnitEnergy();
233 // nTracksTpcOnly is necessary to count the number of candidate cells
234 // but it doesn't give the real multiplicity -> it will be extracted
235 // in the jet finder through the reference to tracks
242 // Fill energy in TPC acceptance
243 uArray->SetUnitEnergy(unitEnergy + pt);
246 if(uArray->GetUnitEnergy()<ptMin){
247 uArray->SetUnitCutFlag(kPtSmaller);
250 uArray->SetUnitCutFlag(kPtHigher);
251 if(ok) nTracksTpcOnlyCut++;
255 if(sflag[goodTrack] == 1) {
256 uArray->SetUnitSignalFlag(kGood);
257 uArray->SetUnitSignalFlagC(kFALSE,kGood);
258 } else uArray->SetUnitSignalFlagC(kFALSE,kBad);
260 if(uArray->GetUnitEnergy()>0 && ref){
263 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
265 fRefArray->Add(uArray);
272 Int_t nElements = fTPCGrid->GetNEntries();
273 // Fill track information in EMCAL acceptance
274 if((eta >= fEtaMin && eta <= fEtaMax) &&
275 (phi >= fPhiMin && phi <= fPhiMax))// &&
277 // Include dead-zones
280 Double_t phimin0 = 0., phimin1 = 0., phimin2 = 0., phimin3 = 0., phimin4 = 0.;
281 Double_t phimax0 = 0., phimax1 = 0., phimax2 = 0., phimax3 = 0., phimax4 = 0.;
282 fGeom->GetPhiBoundariesOfSMGap(0,phimin0,phimax0);
283 fGeom->GetPhiBoundariesOfSMGap(1,phimin1,phimax1);
284 fGeom->GetPhiBoundariesOfSMGap(2,phimin2,phimax2);
285 fGeom->GetPhiBoundariesOfSMGap(3,phimin3,phimax3);
286 fGeom->GetPhiBoundariesOfSMGap(4,phimin4,phimax4);
287 Int_t n0 = fGrid0->GetNEntries();
288 Int_t n1 = fGrid1->GetNEntries();
289 Int_t n2 = fGrid2->GetNEntries();
290 Int_t n3 = fGrid3->GetNEntries();
292 if(phi >= phimin0 && phi <= phimax0){
293 Int_t id0 = fGrid0->GetIndex(phi,eta)-1;
294 AliJetUnitArray *uArray0 = (AliJetUnitArray*)fUnitArray->At(id0+fNumUnits+nElements);
295 TRefArray *reference0 = uArray0->GetUnitTrackRef();
296 if (reference0->GetEntries() == 0) {
297 new(reference0) TRefArray(TProcessID::GetProcessWithUID(track));
300 reference0->Add(track);
301 uArray0->SetUnitTrackID(it);
303 Float_t uEnergy0 = uArray0->GetUnitEnergy();
305 Bool_t ref0 = kFALSE;
311 uArray0->SetUnitEnergy(uEnergy0+pt);
313 if(uArray0->GetUnitEnergy()<ptMin)
314 uArray0->SetUnitCutFlag(kPtSmaller);
316 uArray0->SetUnitCutFlag(kPtHigher);
317 if(ok0) nTracksEmcalDZCut++;
319 if(sflag[goodTrack] == 1) {
320 uArray0->SetUnitSignalFlag(kGood);
321 uArray0->SetUnitSignalFlagC(kFALSE,kGood);
322 } else uArray0->SetUnitSignalFlagC(kFALSE,kBad);
324 if(uArray0->GetUnitEnergy()>0 && ref0){
326 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray0));
329 fRefArray->Add(uArray0);
333 if(phi >= phimin1 && phi <= phimax1){
334 Int_t id1 = fGrid1->GetIndex(phi,eta)-1+n0;
335 AliJetUnitArray *uArray1 = (AliJetUnitArray*)fUnitArray->At(id1+fNumUnits+nElements);
336 TRefArray *reference1 = uArray1->GetUnitTrackRef();
337 if (reference1->GetEntries() == 0) {
338 new(reference1) TRefArray(TProcessID::GetProcessWithUID(track));
341 reference1->Add(track);
343 Float_t uEnergy1 = uArray1->GetUnitEnergy();
345 Bool_t ref1 = kFALSE;
351 uArray1->SetUnitEnergy(uEnergy1+pt);
353 if(uArray1->GetUnitEnergy()<ptMin)
354 uArray1->SetUnitCutFlag(kPtSmaller);
356 uArray1->SetUnitCutFlag(kPtHigher);
357 if(ok1) nTracksEmcalDZCut++;
359 if(sflag[goodTrack] == 1) {
360 uArray1->SetUnitSignalFlag(kGood);
361 uArray1->SetUnitSignalFlagC(kFALSE,kGood);
362 } else uArray1->SetUnitSignalFlagC(kFALSE,kBad);
364 if(uArray1->GetUnitEnergy()>0 && ref1){
366 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray1));
369 fRefArray->Add(uArray1);
373 if(phi >= phimin2 && phi <= phimax2){
374 Int_t id2 = fGrid2->GetIndex(phi,eta)-1+n0+n1;
375 AliJetUnitArray *uArray2 = (AliJetUnitArray*)fUnitArray->At(id2+fNumUnits+nElements);
376 TRefArray *reference2 = uArray2->GetUnitTrackRef();
377 if (reference2->GetEntries() == 0) {
378 new(reference2) TRefArray(TProcessID::GetProcessWithUID(track));
381 reference2->Add(track);
383 Float_t uEnergy2 = uArray2->GetUnitEnergy();
385 Bool_t ref2 = kFALSE;
391 uArray2->SetUnitEnergy(uEnergy2+pt);
393 if(uArray2->GetUnitEnergy()<ptMin)
394 uArray2->SetUnitCutFlag(kPtSmaller);
396 uArray2->SetUnitCutFlag(kPtHigher);
397 if(ok2) nTracksEmcalDZCut++;
399 if(sflag[goodTrack] == 1) {
400 uArray2->SetUnitSignalFlag(kGood);
401 uArray2->SetUnitSignalFlagC(kFALSE,kGood);
402 } else uArray2->SetUnitSignalFlagC(kFALSE,kBad);
404 if(uArray2->GetUnitEnergy()>0 && ref2){
406 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray2));
409 fRefArray->Add(uArray2);
413 if(phi >= phimin3 && phi <= phimax3){
414 Int_t id3 = fGrid3->GetIndex(phi,eta)-1+n0+n1+n2;
415 AliJetUnitArray *uArray3 = (AliJetUnitArray*)fUnitArray->At(id3+fNumUnits+nElements);
416 TRefArray *reference3 = uArray3->GetUnitTrackRef();
417 if (reference3->GetEntries() == 0) {
418 new(reference3) TRefArray(TProcessID::GetProcessWithUID(track));
421 reference3->Add(track);
423 Float_t uEnergy3 = uArray3->GetUnitEnergy();
425 Bool_t ref3 = kFALSE;
431 uArray3->SetUnitEnergy(uEnergy3+pt);
433 if(uArray3->GetUnitEnergy()<ptMin)
434 uArray3->SetUnitCutFlag(kPtSmaller);
436 uArray3->SetUnitCutFlag(kPtHigher);
437 if(ok3) nTracksEmcalDZCut++;
439 if(sflag[goodTrack] == 1) {
440 uArray3->SetUnitSignalFlag(kGood);
441 uArray3->SetUnitSignalFlagC(kFALSE,kGood);
442 } else uArray3->SetUnitSignalFlagC(kFALSE,kBad);
444 if(uArray3->GetUnitEnergy()>0 && ref3){
446 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray3));
449 fRefArray->Add(uArray3);
453 if(phi >= phimin4 && phi <= phimax4){
454 Int_t id4 = fGrid4->GetIndex(phi,eta)-1+n0+n1+n2+n3;
455 AliJetUnitArray *uArray4 = (AliJetUnitArray*)fUnitArray->At(id4+fNumUnits+nElements);
456 TRefArray *reference4 = uArray4->GetUnitTrackRef();
457 if (reference4->GetEntries() == 0) {
458 new(reference4) TRefArray(TProcessID::GetProcessWithUID(track));
461 reference4->Add(track);
463 Float_t uEnergy4 = uArray4->GetUnitEnergy();
465 Bool_t ref4 = kFALSE;
471 uArray4->SetUnitEnergy(uEnergy4+pt);
473 if(uArray4->GetUnitEnergy()<ptMin)
474 uArray4->SetUnitCutFlag(kPtSmaller);
476 uArray4->SetUnitCutFlag(kPtHigher);
477 if(ok4) nTracksEmcalDZCut++;
479 if(sflag[goodTrack] == 1) {
480 uArray4->SetUnitSignalFlag(kGood);
481 uArray4->SetUnitSignalFlagC(kFALSE,kGood);
482 } else uArray4->SetUnitSignalFlagC(kFALSE,kBad);
484 if(uArray4->GetUnitEnergy()>0 && ref4){
486 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray4));
489 fRefArray->Add(uArray4);
495 fGeom->GetAbsCellIdFromEtaPhi(eta,phi,towerID);
496 if(towerID==-1) continue;
498 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(towerID);
499 TRefArray *reference = uArray->GetUnitTrackRef();
500 if (reference->GetEntries() == 0) {
501 new(reference) TRefArray(TProcessID::GetProcessWithUID(track));
504 reference->Add(track);
506 Float_t unitEnergy = uArray->GetUnitEnergy();
515 // Do Hadron Correction
516 // For the moment I apply MIP correction if p >= 0.5 GeV/c
517 // and only for tracks inside EMCal acceptance
518 // End of if(fGrid==1) -> hadron correction for all tracks
519 if (fApplyMIPCorrection != 0 && p3.Mag() >= 0.5)
521 ((AliJetHadronCorrectionv1*)fHadCorr)->SetGeometry("EMCAL_COMPLETE",1.);
523 Double_t etaOut = 0.;
524 Double_t phiOut = 0.;
525 ((AliJetHadronCorrectionv1*)fHadCorr)->TrackPositionEMCal(track,etaOut,phiOut);
527 // One has to make sure that the track hits the calorimeter
528 // We can then correct on average for these particles
529 if((etaOut >= fEtaMin && etaOut <= fEtaMax) &&
530 (phiOut >= fPhiMin && phiOut <= fPhiMax))// &&
532 Double_t hCEnergy = (Double_t)fHadCorr->GetEnergy(p3.Mag(), (Double_t)eta,0);
533 unitEnergy -= hCEnergy*TMath::Sin(2.0*TMath::ATan(TMath::Exp(-eta)));
535 } //end Hadron Correction loop
537 if((unitEnergy + pt) > 0.) uArray->SetUnitEnergy(unitEnergy + pt);
538 else uArray->SetUnitEnergy(0.);
541 if(uArray->GetUnitEnergy()<ptMin){
542 uArray->SetUnitCutFlag(kPtSmaller);
545 uArray->SetUnitCutFlag(kPtHigher);
546 if(ok) nTracksEmcalCut++;
550 if(sflag[goodTrack] == 1) {
551 uArray->SetUnitSignalFlag(kGood);
552 uArray->SetUnitSignalFlagC(kFALSE,kGood);
553 } else uArray->SetUnitSignalFlagC(kFALSE,kBad);
555 if(uArray->GetUnitEnergy()>0 && ref){
557 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
560 fRefArray->Add(uArray);
562 } // end loop on EMCal acceptance cut + tracks quality
564 // Outside EMCal acceptance
565 // Int_t idTPC = GetIndexFromEtaPhi(etaT[i],phiT[i]);
566 Int_t idTPC = fTPCGrid->GetIndex(phi,eta);
568 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(fNumUnits-1+idTPC);
569 TRefArray *reference = uArray->GetUnitTrackRef();
570 if (reference->GetEntries() == 0) {
571 new(reference) TRefArray(TProcessID::GetProcessWithUID(track));
573 reference->Add(track);
575 Float_t unitEnergy2 = uArray->GetUnitEnergy(); // check if fNumUnits or fNumUnits-1
576 Bool_t okout = kFALSE;
577 Bool_t refout = kFALSE;
583 // Fill energy outside emcal acceptance
584 uArray->SetUnitEnergy(unitEnergy2 + pt);
587 if(uArray->GetUnitEnergy()<ptMin){
588 uArray->SetUnitCutFlag(kPtSmaller);
591 uArray->SetUnitCutFlag(kPtHigher);
592 if(okout) nTracksTpcCut++;
595 if(sflag[goodTrack] == 1) {
596 uArray->SetUnitSignalFlag(kGood);
597 uArray->SetUnitSignalFlagC(kFALSE,kGood);
598 } else uArray->SetUnitSignalFlagC(kFALSE,kBad);
600 if(uArray->GetUnitEnergy()>0 && refout){
602 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
605 fRefArray->Add(uArray);
610 // Do hadron correction
611 // In principle, the best thing to do as particles which
612 // eta,phi (at vertex) correspond to EMCal dead zone can deposit energy
613 // in the calorimeter as well as particles with eta,phi (at vertex)
614 // outside the EMCal acceptance
615 cout << "Hadronic correction for all tracks goes here" << endl;
617 // For the moment I apply MIP correction if p >= 0.5 GeV/c
618 if (fApplyMIPCorrection != 0 && p3.Mag() >= 0.5)
620 ((AliJetHadronCorrectionv1*)fHadCorr)->SetGeometry("EMCAL_COMPLETE",1.);
622 Double_t etaOut = 0.;
623 Double_t phiOut = 0.;
625 Float_t unitEnergy = 0.;
626 ((AliJetHadronCorrectionv1*)fHadCorr)->TrackPositionEMCal(track,etaOut,phiOut);
628 // One has to make sure that the track hits the calorimeter
629 // We can then correct on average for these particles
630 if((etaOut >= fEtaMin && etaOut <= fEtaMax) &&
631 (phiOut >= fPhiMin && phiOut <= fPhiMax))// &&
635 fGeom->GetAbsCellIdFromEtaPhi(etaOut,phiOut,towerID2);
636 if(towerID2==-1) continue;
638 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(towerID2);
639 unitEnergy = uArray->GetUnitEnergy();
648 Double_t hCEnergy = (Double_t)fHadCorr->GetEnergy(p3.Mag(), (Double_t)eta,0);
649 unitEnergy -= hCEnergy*TMath::Sin(2.0*TMath::ATan(TMath::Exp(-eta)));
651 // The energy in this unitarray can be negative
652 // If the Digits are then filled, make sure the energy sum
653 // is positive. If not, put the value to zero in AliJetAODUnitArrayEMCalDigit.
654 // If they are not filled, put this value to zero
655 uArray->SetUnitEnergy(unitEnergy);
657 if(uArray->GetUnitEnergy()!=0. && ref){
659 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
662 fRefArray->Add(uArray);
665 else cout << "Track out of EMCal acceptance" << endl;
667 } //end Hadron Correction loop
674 } // end loop on entries (tpc tracks)
679 cout << "End of Tracks %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << endl;
680 cout << "goodTracks: " << goodTrack << endl;
683 fSignalFlag.Set(goodTrack,sflag);
684 fCutFlag.Set(goodTrack,cflag);
689 // } // end loop on entries (tpc tracks)
692 fNTracks = nTracksTpcOnly;
693 fNTracksCut = nTracksTpcOnlyCut;
695 cout << "fNTracks : " << fNTracks << endl;
696 cout << "fNTracksCut : " << fNTracksCut << endl;
700 fNTracks = nTracksEmcal+nTracksEmcalDZ+nTracksTpc;
701 fNTracksCut = nTracksEmcalCut+nTracksEmcalDZCut+nTracksTpcCut;
703 cout << "fNTracks : " << fNTracks << endl;
704 cout << "fNTracksCut : " << fNTracksCut << endl;