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 "AliJetAODFillUnitArrayTracks.h"
33 // --- ROOT system ---
38 // --- AliRoot header files ---
41 ClassImp(AliJetAODFillUnitArrayTracks)
43 //_____________________________________________________________________________
44 AliJetAODFillUnitArrayTracks::AliJetAODFillUnitArrayTracks()
45 : AliJetFillUnitArray(),
48 fApplyMIPCorrection(kTRUE),
59 //_____________________________________________________________________________
60 AliJetAODFillUnitArrayTracks::AliJetAODFillUnitArrayTracks(AliAODEvent* aod)
61 : AliJetFillUnitArray(),
64 fApplyMIPCorrection(kTRUE),
75 //_____________________________________________________________________________
76 AliJetAODFillUnitArrayTracks::AliJetAODFillUnitArrayTracks(const AliJetAODFillUnitArrayTracks &det)
77 : AliJetFillUnitArray(det),
78 fNumUnits(det.fNumUnits),
79 fHadCorr(det.fHadCorr),
80 fApplyMIPCorrection(det.fApplyMIPCorrection),
91 //_____________________________________________________________________________
92 AliJetAODFillUnitArrayTracks& AliJetAODFillUnitArrayTracks::operator=(const AliJetAODFillUnitArrayTracks& other)
96 fNumUnits = other.fNumUnits;
97 fHadCorr = other.fHadCorr;
98 fApplyMIPCorrection = other.fApplyMIPCorrection;
100 fGrid0 = other.fGrid0;
101 fGrid1 = other.fGrid1;
102 fGrid2 = other.fGrid2;
103 fGrid3 = other.fGrid3;
104 fGrid4 = other.fGrid4;
109 //____________________________________________________________________________
110 void AliJetAODFillUnitArrayTracks::InitParameters()
112 fHadCorr = 0; // For hadron correction
113 fNumUnits = fGeom->GetNCells(); // Number of towers in EMCAL
115 fTPCGrid->GetAccParam(fNphi,fNeta,fPhiMin,
116 fPhiMax,fEtaMin,fEtaMax);
117 fTPCGrid->GetBinParam(fPhiBinInTPCAcc,fEtaBinInTPCAcc,
118 fPhiBinInEMCalAcc,fEtaBinInEMCalAcc,fNbinPhi);
120 fIndex = fTPCGrid->GetIndexObject();
123 for(Int_t i=0; i<fNphi+1; i++)
124 for(Int_t j=0; j<fNeta+1; j++) {cout << "fIndex[" << i << "," << j << "] : " <<
125 (*fIndex)(i,j) << endl; }
127 if(fDebug>1) printf("\n Parameters initiated ! \n");
130 //_____________________________________________________________________________
131 AliJetAODFillUnitArrayTracks::~AliJetAODFillUnitArrayTracks()
136 //_____________________________________________________________________________
137 void AliJetAODFillUnitArrayTracks::Exec(Option_t* const /*option*/)
138 //void AliJetAODFillUnitArrayTracks::Exec(Option_t* option)
142 // Fill the unit array with the charged particle information in AOD
145 fDebug = fReaderHeader->GetDebug();
150 // get number of tracks in event (for the loop)
153 //(not used ?) Int_t nt2 = 0;
159 nt = fAOD->GetNumberOfTracks();
160 if(fDebug>1) cout << "Number of Tracks in ESD : " << nt << endl;
162 // temporary storage of signal and pt cut flag
163 Int_t* sflag = new Int_t[nt];
164 Int_t* cflag = new Int_t[nt];
166 // get cuts set by user
167 Float_t ptMin = fReaderHeader->GetPtCut();
168 Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
169 Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
170 fOpt = fReaderHeader->GetDetector();
171 fDZ = fReaderHeader->GetDZ();
172 UInt_t filterMask = ((AliJetAODReaderHeader*)fReaderHeader)->GetTestFilterMask();
174 Int_t nTracksEmcal = 0;
175 Int_t nTracksEmcalDZ = 0;
176 Int_t nTracksTpc = 0;
177 Int_t nTracksTpcOnly = 0;
178 Int_t nTracksEmcalCut = 0;
179 Int_t nTracksEmcalDZCut = 0;
180 Int_t nTracksTpcCut = 0;
181 Int_t nTracksTpcOnlyCut = 0;
183 //(not used ?) Int_t nElements = fTPCGrid->GetNEntries();
184 //(not used ?) Int_t nElements2 = fEMCalGrid->GetNEntries();
185 fGrid = fTPCGrid->GetGridType();
187 //(not used ?) Int_t nEntries = fUnitArray->GetEntries();
189 //(not used ?) Int_t nRefEnt = fRefArray->GetEntries();
193 for (Int_t it = 0; it < nmax; it++) {
195 track = fAOD->GetTrack(it);
196 UInt_t status = track->GetStatus();
199 track->GetPxPyPz(mom);
201 p3.SetXYZ(mom[0],mom[1],mom[2]);
204 //(not used ?) Float_t mass = 0.;
207 phi = ( (p3.Phi()) < 0) ? (p3.Phi()) + 2. * TMath::Pi() : (p3.Phi());
208 if (status == 0) continue;
209 if((filterMask>0)&&!(track->TestFilterBit(filterMask)))continue;
211 if ( (eta > etaMax) || (eta < etaMin)) continue; // checking eta cut
214 if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
216 if (pt > ptMin) cflag[goodTrack]=1; // pt cut
220 // Only TPC filled from its grid in its total acceptance
222 Int_t idTPC = fTPCGrid->GetIndex(phi,eta);
226 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(idTPC-1);
227 TRefArray *reference = uArray->GetUnitTrackRef();
228 reference->Add(track);
230 Float_t unitEnergy = 0.;
231 unitEnergy = uArray->GetUnitEnergy();
232 // nTracksTpcOnly is necessary to count the number of candidate cells
233 // but it doesn't give the real multiplicity -> it will be extracted
234 // in the jet finder through the reference to tracks
241 // Fill energy in TPC acceptance
242 uArray->SetUnitEnergy(unitEnergy + pt);
245 if(uArray->GetUnitEnergy()<ptMin){
246 uArray->SetUnitCutFlag(kPtSmaller);
249 uArray->SetUnitCutFlag(kPtHigher);
250 if(ok) nTracksTpcOnlyCut++;
254 if(sflag[goodTrack] == 1) {
255 uArray->SetUnitSignalFlag(kGood);
256 uArray->SetUnitSignalFlagC(kFALSE,kGood);
257 } else uArray->SetUnitSignalFlagC(kFALSE,kBad);
259 if(uArray->GetUnitEnergy()>0 && ref){
262 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
264 fRefArray->Add(uArray);
271 Int_t nElements = fTPCGrid->GetNEntries();
272 // Fill track information in EMCAL acceptance
273 if((eta >= fEtaMin && eta <= fEtaMax) &&
274 (phi >= fPhiMin && phi <= fPhiMax))// &&
276 // Include dead-zones
279 Double_t phimin0 = 0., phimin1 = 0., phimin2 = 0., phimin3 = 0., phimin4 = 0.;
280 Double_t phimax0 = 0., phimax1 = 0., phimax2 = 0., phimax3 = 0., phimax4 = 0.;
281 fGeom->GetPhiBoundariesOfSMGap(0,phimin0,phimax0);
282 fGeom->GetPhiBoundariesOfSMGap(1,phimin1,phimax1);
283 fGeom->GetPhiBoundariesOfSMGap(2,phimin2,phimax2);
284 fGeom->GetPhiBoundariesOfSMGap(3,phimin3,phimax3);
285 fGeom->GetPhiBoundariesOfSMGap(4,phimin4,phimax4);
286 Int_t n0 = fGrid0->GetNEntries();
287 Int_t n1 = fGrid1->GetNEntries();
288 Int_t n2 = fGrid2->GetNEntries();
289 Int_t n3 = fGrid3->GetNEntries();
290 //(not used ?) Int_t n4 = fGrid4->GetNEntries();
292 // cout << "phimin0: " << phimin0 << ", phimax0: " << phimax0 << endl;
293 if(phi >= phimin0 && phi <= phimax0){
294 Int_t id0 = fGrid0->GetIndex(phi,eta)-1;
295 AliJetUnitArray *uArray0 = (AliJetUnitArray*)fUnitArray->At(id0+fNumUnits+nElements);
296 TRefArray *reference0 = uArray0->GetUnitTrackRef();
297 reference0->Add(track);
298 uArray0->SetUnitTrackID(it);
300 Float_t uEnergy0 = uArray0->GetUnitEnergy();
302 Bool_t ref0 = kFALSE;
308 uArray0->SetUnitEnergy(uEnergy0+pt);
310 if(uArray0->GetUnitEnergy()<ptMin)
311 uArray0->SetUnitCutFlag(kPtSmaller);
313 uArray0->SetUnitCutFlag(kPtHigher);
314 if(ok0) nTracksEmcalDZCut++;
316 if(sflag[goodTrack] == 1) {
317 uArray0->SetUnitSignalFlag(kGood);
318 uArray0->SetUnitSignalFlagC(kFALSE,kGood);
319 } else uArray0->SetUnitSignalFlagC(kFALSE,kBad);
321 if(uArray0->GetUnitEnergy()>0 && ref0){
323 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray0));
326 fRefArray->Add(uArray0);
330 if(phi >= phimin1 && phi <= phimax1){
331 Int_t id1 = fGrid1->GetIndex(phi,eta)-1+n0;
332 AliJetUnitArray *uArray1 = (AliJetUnitArray*)fUnitArray->At(id1+fNumUnits+nElements);
333 TRefArray *reference1 = uArray1->GetUnitTrackRef();
334 reference1->Add(track);
336 Float_t uEnergy1 = uArray1->GetUnitEnergy();
338 Bool_t ref1 = kFALSE;
344 uArray1->SetUnitEnergy(uEnergy1+pt);
346 if(uArray1->GetUnitEnergy()<ptMin)
347 uArray1->SetUnitCutFlag(kPtSmaller);
349 uArray1->SetUnitCutFlag(kPtHigher);
350 if(ok1) nTracksEmcalDZCut++;
352 if(sflag[goodTrack] == 1) {
353 uArray1->SetUnitSignalFlag(kGood);
354 uArray1->SetUnitSignalFlagC(kFALSE,kGood);
355 } else uArray1->SetUnitSignalFlagC(kFALSE,kBad);
357 if(uArray1->GetUnitEnergy()>0 && ref1){
359 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray1));
362 fRefArray->Add(uArray1);
366 if(phi >= phimin2 && phi <= phimax2){
367 Int_t id2 = fGrid2->GetIndex(phi,eta)-1+n0+n1;
368 AliJetUnitArray *uArray2 = (AliJetUnitArray*)fUnitArray->At(id2+fNumUnits+nElements);
369 TRefArray *reference2 = uArray2->GetUnitTrackRef();
370 reference2->Add(track);
372 Float_t uEnergy2 = uArray2->GetUnitEnergy();
374 Bool_t ref2 = kFALSE;
380 uArray2->SetUnitEnergy(uEnergy2+pt);
382 if(uArray2->GetUnitEnergy()<ptMin)
383 uArray2->SetUnitCutFlag(kPtSmaller);
385 uArray2->SetUnitCutFlag(kPtHigher);
386 if(ok2) nTracksEmcalDZCut++;
388 if(sflag[goodTrack] == 1) {
389 uArray2->SetUnitSignalFlag(kGood);
390 uArray2->SetUnitSignalFlagC(kFALSE,kGood);
391 } else uArray2->SetUnitSignalFlagC(kFALSE,kBad);
393 if(uArray2->GetUnitEnergy()>0 && ref2){
395 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray2));
398 fRefArray->Add(uArray2);
402 if(phi >= phimin3 && phi <= phimax3){
403 Int_t id3 = fGrid3->GetIndex(phi,eta)-1+n0+n1+n2;
404 AliJetUnitArray *uArray3 = (AliJetUnitArray*)fUnitArray->At(id3+fNumUnits+nElements);
405 TRefArray *reference3 = uArray3->GetUnitTrackRef();
406 reference3->Add(track);
408 Float_t uEnergy3 = uArray3->GetUnitEnergy();
410 Bool_t ref3 = kFALSE;
416 uArray3->SetUnitEnergy(uEnergy3+pt);
418 if(uArray3->GetUnitEnergy()<ptMin)
419 uArray3->SetUnitCutFlag(kPtSmaller);
421 uArray3->SetUnitCutFlag(kPtHigher);
422 if(ok3) nTracksEmcalDZCut++;
424 if(sflag[goodTrack] == 1) {
425 uArray3->SetUnitSignalFlag(kGood);
426 uArray3->SetUnitSignalFlagC(kFALSE,kGood);
427 } else uArray3->SetUnitSignalFlagC(kFALSE,kBad);
429 if(uArray3->GetUnitEnergy()>0 && ref3){
431 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray3));
434 fRefArray->Add(uArray3);
438 if(phi >= phimin4 && phi <= phimax4){
439 Int_t id4 = fGrid4->GetIndex(phi,eta)-1+n0+n1+n2+n3;
440 AliJetUnitArray *uArray4 = (AliJetUnitArray*)fUnitArray->At(id4+fNumUnits+nElements);
441 TRefArray *reference4 = uArray4->GetUnitTrackRef();
442 reference4->Add(track);
444 Float_t uEnergy4 = uArray4->GetUnitEnergy();
446 Bool_t ref4 = kFALSE;
452 uArray4->SetUnitEnergy(uEnergy4+pt);
454 if(uArray4->GetUnitEnergy()<ptMin)
455 uArray4->SetUnitCutFlag(kPtSmaller);
457 uArray4->SetUnitCutFlag(kPtHigher);
458 if(ok4) nTracksEmcalDZCut++;
460 if(sflag[goodTrack] == 1) {
461 uArray4->SetUnitSignalFlag(kGood);
462 uArray4->SetUnitSignalFlagC(kFALSE,kGood);
463 } else uArray4->SetUnitSignalFlagC(kFALSE,kBad);
465 if(uArray4->GetUnitEnergy()>0 && ref4){
467 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray4));
470 fRefArray->Add(uArray4);
476 fGeom->GetAbsCellIdFromEtaPhi(eta,phi,towerID);
477 if(towerID==-1) continue;
479 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(towerID);
480 TRefArray *reference = uArray->GetUnitTrackRef();
481 reference->Add(track);
483 Float_t unitEnergy = uArray->GetUnitEnergy();
492 // Do Hadron Correction
493 // This is under construction !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
494 // Parametrization to be added
495 if (fApplyMIPCorrection != 0)
497 // Float_t hCEnergy = fHadCorr->GetEnergy(p3.Mag(), (Double_t)eta,0);
498 // unitEnergy -= hCEnergy*TMath::Sin(2.0*TMath::ATan(TMath::Exp(-eta)));
499 } //end Hadron Correction loop
501 uArray->SetUnitEnergy(unitEnergy + pt);
504 if(uArray->GetUnitEnergy()<ptMin){
505 uArray->SetUnitCutFlag(kPtSmaller);
508 uArray->SetUnitCutFlag(kPtHigher);
509 if(ok) nTracksEmcalCut++;
513 if(sflag[goodTrack] == 1) {
514 uArray->SetUnitSignalFlag(kGood);
515 uArray->SetUnitSignalFlagC(kFALSE,kGood);
516 } else uArray->SetUnitSignalFlagC(kFALSE,kBad);
518 if(uArray->GetUnitEnergy()>0 && ref){
520 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
523 fRefArray->Add(uArray);
525 } // end loop on EMCal acceptance cut + tracks quality
527 // Outside EMCal acceptance
528 // Int_t idTPC = GetIndexFromEtaPhi(etaT[i],phiT[i]);
529 Int_t idTPC = fTPCGrid->GetIndex(phi,eta);
531 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(fNumUnits-1+idTPC);
532 TRefArray *reference = uArray->GetUnitTrackRef();
533 reference->Add(track);
535 Float_t unitEnergy2 = uArray->GetUnitEnergy(); // check if fNumUnits or fNumUnits-1
536 Bool_t okout = kFALSE;
537 Bool_t refout = kFALSE;
543 // Fill energy outside emcal acceptance
544 uArray->SetUnitEnergy(unitEnergy2 + pt);
547 if(uArray->GetUnitEnergy()<ptMin){
548 uArray->SetUnitCutFlag(kPtSmaller);
551 uArray->SetUnitCutFlag(kPtHigher);
552 if(okout) nTracksTpcCut++;
555 if(sflag[goodTrack] == 1) {
556 uArray->SetUnitSignalFlag(kGood);
557 uArray->SetUnitSignalFlagC(kFALSE,kGood);
558 } else uArray->SetUnitSignalFlagC(kFALSE,kBad);
560 if(uArray->GetUnitEnergy()>0 && refout){
562 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
565 fRefArray->Add(uArray);
572 } // end loop on entries (tpc tracks)
577 cout << "End of Tracks %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << endl;
578 cout << "goodTracks: " << goodTrack << endl;
584 // } // end loop on entries (tpc tracks)
587 fNTracks = nTracksTpcOnly;
588 fNTracksCut = nTracksTpcOnlyCut;
590 cout << "fNTracks : " << fNTracks << endl;
591 cout << "fNTracksCut : " << fNTracksCut << endl;
595 fNTracks = nTracksEmcal+nTracksEmcalDZ+nTracksTpc;
596 fNTracksCut = nTracksEmcalCut+nTracksEmcalDZCut+nTracksTpcCut;
598 cout << "fNTracks : " << fNTracks << endl;
599 cout << "fNTracksCut : " << fNTracksCut << endl;