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 // --- AliRoot header files ---
31 #include "AliJetUnitArray.h"
32 #include "AliJetESDFillUnitArrayTracks.h"
34 // --- ROOT system ---
39 // --- AliRoot header files ---
42 class AliJetESDReader;
44 ClassImp(AliJetESDFillUnitArrayTracks)
46 //_____________________________________________________________________________
47 AliJetESDFillUnitArrayTracks::AliJetESDFillUnitArrayTracks()
48 : AliJetFillUnitArray(),
51 fApplyMIPCorrection(kTRUE),
62 //_____________________________________________________________________________
63 AliJetESDFillUnitArrayTracks::AliJetESDFillUnitArrayTracks(AliESDEvent* esd)
64 : AliJetFillUnitArray(),
67 fApplyMIPCorrection(kTRUE),
78 //_____________________________________________________________________________
79 AliJetESDFillUnitArrayTracks::AliJetESDFillUnitArrayTracks(const AliJetESDFillUnitArrayTracks &det)
80 : AliJetFillUnitArray(det),
81 fNumUnits(det.fNumUnits),
82 fHadCorr(det.fHadCorr),
83 fApplyMIPCorrection(det.fApplyMIPCorrection),
94 //------------------------------------------------------------------
95 AliJetESDFillUnitArrayTracks& AliJetESDFillUnitArrayTracks::operator=(const AliJetESDFillUnitArrayTracks& other)
99 fNumUnits = other.fNumUnits;
100 fHadCorr = other.fHadCorr;
101 fApplyMIPCorrection = other.fApplyMIPCorrection;
103 fGrid0 = other.fGrid0;
104 fGrid1 = other.fGrid1;
105 fGrid2 = other.fGrid2;
106 fGrid3 = other.fGrid3;
107 fGrid4 = other.fGrid4;
112 //____________________________________________________________________________
113 void AliJetESDFillUnitArrayTracks::InitParameters()
115 fHadCorr = 0; // For hadron correction
116 fNumUnits = fGeom->GetNCells(); // Number of towers in EMCAL
117 cout << "In AliJetESDFillUnitArrayTracks:InitParameters(), Ncells : " << fNumUnits << endl;
119 fTPCGrid->GetAccParam(fNphi,fNeta,fPhiMin,
120 fPhiMax,fEtaMin,fEtaMax);
121 fTPCGrid->GetBinParam(fPhiBinInTPCAcc,fEtaBinInTPCAcc,
122 fPhiBinInEMCalAcc,fEtaBinInEMCalAcc,fNbinPhi);
124 fIndex = fTPCGrid->GetIndexObject();
127 for(Int_t i=0; i<fNphi+1; i++)
128 for(Int_t j=0; j<fNeta+1; j++) {cout << "fIndex[" << i << "," << j << "] : " <<
129 (*fIndex)(i,j) << endl; }
131 if(fDebug>1) printf("\n Parameters initiated ! \n");
134 //_____________________________________________________________________________
135 AliJetESDFillUnitArrayTracks::~AliJetESDFillUnitArrayTracks()
140 //_____________________________________________________________________________
141 void AliJetESDFillUnitArrayTracks::Exec(Option_t* const /*option*/)
145 // Fill the unit array with the charged particle information in ESD
148 fDebug = fReaderHeader->GetDebug();
153 // get number of tracks in event (for the loop)
156 //(not used ?) Int_t nt2 = 0;
162 if(fDebug>1) cout << "fESD in Fill array : " << fESD << endl;
164 nt = fESD->GetNumberOfTracks();
165 if(fDebug>1) cout << "Number of Tracks in ESD : " << nt << endl;
167 // temporary storage of signal and pt cut flag
168 Int_t* sflag = new Int_t[nt];
169 Int_t* cflag = new Int_t[nt];
171 // get cuts set by user
172 Float_t ptMin = fReaderHeader->GetPtCut();
173 Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
174 Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
175 fOpt = fReaderHeader->GetDetector();
176 fDZ = fReaderHeader->GetDZ();
178 Int_t nTracksEmcal = 0;
179 Int_t nTracksEmcalDZ = 0;
180 Int_t nTracksTpc = 0;
181 Int_t nTracksTpcOnly = 0;
182 Int_t nTracksEmcalCut = 0;
183 Int_t nTracksEmcalDZCut = 0;
184 Int_t nTracksTpcCut = 0;
185 Int_t nTracksTpcOnlyCut = 0;
187 //(not used ?) Int_t nElements = fTPCGrid->GetNEntries();
188 //(not used ?) Int_t nElements2 = fEMCalGrid->GetNEntries();
189 fGrid = fTPCGrid->GetGridType();
191 //(not used ?) Int_t nEntries = fUnitArray->GetEntries();
193 //(not used ?) Int_t nRefEnt = fRefArray->GetEntries();
197 for (Int_t it = 0; it < nmax; it++) {
199 track = fESD->GetTrack(it);
200 UInt_t status = track->GetStatus();
203 track->GetPxPyPz(mom);
205 p3.SetXYZ(mom[0],mom[1],mom[2]);
209 mass = track->GetMass();
211 if ((status & AliESDtrack::kTPCrefit) == 0) continue; // quality check
212 if ((status & AliESDtrack::kITSrefit) == 0) continue; // quality check
213 if (((AliJetESDReaderHeader*) fReaderHeader)->ReadSignalOnly()
214 && TMath::Abs(track->GetLabel()) > 10000) continue; // quality check
215 if (((AliJetESDReaderHeader*) fReaderHeader)->ReadBkgdOnly()
216 && TMath::Abs(track->GetLabel()) < 10000) continue; // quality check
218 phi = ( (p3.Phi()) < 0) ? (p3.Phi()) + 2. * TMath::Pi() : (p3.Phi());
220 if ( (eta > etaMax) || (eta < etaMin)) continue; // checking eta cut
223 if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
225 if (pt > ptMin) cflag[goodTrack]=1; // pt cut
229 // Only TPC filled from its grid in its total acceptance
231 Int_t idTPC = fTPCGrid->GetIndex(phi,eta);
235 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(idTPC-1);
236 TRefArray *reference = uArray->GetUnitTrackRef();
237 reference->Add(track);
239 Float_t unitEnergy = 0.;
240 unitEnergy = uArray->GetUnitEnergy();
241 // nTracksTpcOnly is necessary to count the number of candidate cells
242 // but it doesn't give the real multiplicity -> it will be extracted
243 // in the jet finder through the reference to tracks
250 // Fill energy in TPC acceptance
251 uArray->SetUnitEnergy(unitEnergy + pt);
252 uArray->SetUnitMass(mass);
255 if(uArray->GetUnitEnergy()<ptMin){
256 uArray->SetUnitCutFlag(kPtSmaller);
259 uArray->SetUnitCutFlag(kPtHigher);
260 if(ok) nTracksTpcOnlyCut++;
264 if(sflag[goodTrack] == 1) {
265 uArray->SetUnitSignalFlag(kGood);
266 uArray->SetUnitSignalFlagC(kFALSE,kGood);
267 } else uArray->SetUnitSignalFlagC(kFALSE,kBad);
269 if(uArray->GetUnitEnergy()>0 && ref){
272 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
274 fRefArray->Add(uArray);
280 Int_t nElements = fTPCGrid->GetNEntries();
281 // Fill track information in EMCAL acceptance
282 if((eta >= fEtaMin && eta <= fEtaMax) &&
283 (phi >= fPhiMin && phi <= fPhiMax))// &&
286 // Include dead-zones
289 Double_t phimin0 = 0., phimin1 = 0., phimin2 = 0., phimin3 = 0., phimin4 = 0.;
290 Double_t phimax0 = 0., phimax1 = 0., phimax2 = 0., phimax3 = 0., phimax4 = 0.;
291 fGeom->GetPhiBoundariesOfSMGap(0,phimin0,phimax0);
292 fGeom->GetPhiBoundariesOfSMGap(1,phimin1,phimax1);
293 fGeom->GetPhiBoundariesOfSMGap(2,phimin2,phimax2);
294 fGeom->GetPhiBoundariesOfSMGap(3,phimin3,phimax3);
295 fGeom->GetPhiBoundariesOfSMGap(4,phimin4,phimax4);
296 Int_t n0 = fGrid0->GetNEntries();
297 Int_t n1 = fGrid1->GetNEntries();
298 Int_t n2 = fGrid2->GetNEntries();
299 Int_t n3 = fGrid3->GetNEntries();
300 //(not used ?) Int_t n4 = fGrid4->GetNEntries();
302 if(phi >= phimin0 && phi <= phimax0){
303 Int_t id0 = fGrid0->GetIndex(phi,eta)-1;
304 AliJetUnitArray *uArray0 = (AliJetUnitArray*)fUnitArray->At(id0+fNumUnits+nElements);
305 TRefArray *reference0 = uArray0->GetUnitTrackRef();
306 reference0->Add(track);
308 Float_t uEnergy0 = uArray0->GetUnitEnergy();
310 Bool_t ref0 = kFALSE;
316 uArray0->SetUnitEnergy(uEnergy0+pt);
318 if(uArray0->GetUnitEnergy()<ptMin)
319 uArray0->SetUnitCutFlag(kPtSmaller);
321 uArray0->SetUnitCutFlag(kPtHigher);
322 if(ok0) nTracksEmcalDZCut++;
324 if(sflag[goodTrack] == 1) {
325 uArray0->SetUnitSignalFlag(kGood);
326 uArray0->SetUnitSignalFlagC(kFALSE,kGood);
327 } else uArray0->SetUnitSignalFlagC(kFALSE,kBad);
329 if(uArray0->GetUnitEnergy()>0 && ref0){
331 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray0));
334 fRefArray->Add(uArray0);
338 if(phi >= phimin1 && phi <= phimax1){
339 Int_t id1 = fGrid1->GetIndex(phi,eta)-1+n0;
340 AliJetUnitArray *uArray1 = (AliJetUnitArray*)fUnitArray->At(id1+fNumUnits+nElements);
341 TRefArray *reference1 = uArray1->GetUnitTrackRef();
342 reference1->Add(track);
344 Float_t uEnergy1 = uArray1->GetUnitEnergy();
346 Bool_t ref1 = kFALSE;
352 uArray1->SetUnitEnergy(uEnergy1+pt);
354 if(uArray1->GetUnitEnergy()<ptMin)
355 uArray1->SetUnitCutFlag(kPtSmaller);
357 uArray1->SetUnitCutFlag(kPtHigher);
358 if(ok1) nTracksEmcalDZCut++;
360 if(sflag[goodTrack] == 1) {
361 uArray1->SetUnitSignalFlag(kGood);
362 uArray1->SetUnitSignalFlagC(kFALSE,kGood);
363 } else uArray1->SetUnitSignalFlagC(kFALSE,kBad);
365 if(uArray1->GetUnitEnergy()>0 && ref1){
367 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray1));
370 fRefArray->Add(uArray1);
374 if(phi >= phimin2 && phi <= phimax2){
375 Int_t id2 = fGrid2->GetIndex(phi,eta)-1+n0+n1;
376 AliJetUnitArray *uArray2 = (AliJetUnitArray*)fUnitArray->At(id2+fNumUnits+nElements);
377 TRefArray *reference2 = uArray2->GetUnitTrackRef();
378 reference2->Add(track);
380 Float_t uEnergy2 = uArray2->GetUnitEnergy();
382 Bool_t ref2 = kFALSE;
388 uArray2->SetUnitEnergy(uEnergy2+pt);
390 if(uArray2->GetUnitEnergy()<ptMin)
391 uArray2->SetUnitCutFlag(kPtSmaller);
393 uArray2->SetUnitCutFlag(kPtHigher);
394 if(ok2) nTracksEmcalDZCut++;
396 if(sflag[goodTrack] == 1) {
397 uArray2->SetUnitSignalFlag(kGood);
398 uArray2->SetUnitSignalFlagC(kFALSE,kGood);
399 } else uArray2->SetUnitSignalFlagC(kFALSE,kBad);
401 if(uArray2->GetUnitEnergy()>0 && ref2){
403 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray2));
406 fRefArray->Add(uArray2);
410 if(phi >= phimin3 && phi <= phimax3){
411 Int_t id3 = fGrid3->GetIndex(phi,eta)-1+n0+n1+n2;
412 AliJetUnitArray *uArray3 = (AliJetUnitArray*)fUnitArray->At(id3+fNumUnits+nElements);
413 TRefArray *reference3 = uArray3->GetUnitTrackRef();
414 reference3->Add(track);
416 Float_t uEnergy3 = uArray3->GetUnitEnergy();
418 Bool_t ref3 = kFALSE;
424 uArray3->SetUnitEnergy(uEnergy3+pt);
426 if(uArray3->GetUnitEnergy()<ptMin)
427 uArray3->SetUnitCutFlag(kPtSmaller);
429 uArray3->SetUnitCutFlag(kPtHigher);
430 if(ok3) nTracksEmcalDZCut++;
432 if(sflag[goodTrack] == 1) {
433 uArray3->SetUnitSignalFlag(kGood);
434 uArray3->SetUnitSignalFlagC(kFALSE,kGood);
435 } else uArray3->SetUnitSignalFlagC(kFALSE,kBad);
437 if(uArray3->GetUnitEnergy()>0 && ref3){
439 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray3));
442 fRefArray->Add(uArray3);
446 if(phi >= phimin4 && phi <= phimax4){
447 Int_t id4 = fGrid4->GetIndex(phi,eta)-1+n0+n1+n2+n3;
448 AliJetUnitArray *uArray4 = (AliJetUnitArray*)fUnitArray->At(id4+fNumUnits+nElements);
449 TRefArray *reference4 = uArray4->GetUnitTrackRef();
450 reference4->Add(track);
452 Float_t uEnergy4 = uArray4->GetUnitEnergy();
454 Bool_t ref4 = kFALSE;
460 uArray4->SetUnitEnergy(uEnergy4+pt);
462 if(uArray4->GetUnitEnergy()<ptMin)
463 uArray4->SetUnitCutFlag(kPtSmaller);
465 uArray4->SetUnitCutFlag(kPtHigher);
466 if(ok4) nTracksEmcalDZCut++;
468 if(sflag[goodTrack] == 1) {
469 uArray4->SetUnitSignalFlag(kGood);
470 uArray4->SetUnitSignalFlagC(kFALSE,kGood);
471 } else uArray4->SetUnitSignalFlagC(kFALSE,kBad);
473 if(uArray4->GetUnitEnergy()>0 && ref4){
475 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray4));
478 fRefArray->Add(uArray4);
484 fGeom->GetAbsCellIdFromEtaPhi(eta,phi,towerID);
485 if(towerID==-1) continue;
487 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(towerID);
488 TRefArray *reference = uArray->GetUnitTrackRef();
489 reference->Add(track);
491 Float_t unitEnergy = uArray->GetUnitEnergy();
500 // Do Hadron Correction
501 // This is under construction !!!!!!!!!!!!!!!!!!!!!!!
502 // Parametrization to be added
503 if (fApplyMIPCorrection != 0)
505 // Float_t hCEnergy = fHadCorr->GetEnergy(p3.Mag(), (Double_t)eta,0);
506 // unitEnergy -= hCEnergy*TMath::Sin(2.0*TMath::ATan(TMath::Exp(-eta)));
507 } //end Hadron Correction loop
509 uArray->SetUnitEnergy(unitEnergy + pt);
512 if(uArray->GetUnitEnergy()<ptMin){
513 uArray->SetUnitCutFlag(kPtSmaller);
516 uArray->SetUnitCutFlag(kPtHigher);
517 if(ok) nTracksEmcalCut++;
521 if(sflag[goodTrack] == 1) {
522 uArray->SetUnitSignalFlag(kGood);
523 uArray->SetUnitSignalFlagC(kFALSE,kGood);
524 } else uArray->SetUnitSignalFlagC(kFALSE,kBad);
527 if(uArray->GetUnitEnergy()>0 && ref){
529 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
532 fRefArray->Add(uArray);
534 } // end loop on EMCal acceptance cut + tracks quality
536 // Outside EMCal acceptance
538 Int_t idTPC = fTPCGrid->GetIndex(phi,eta);
540 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(fNumUnits-1+idTPC);
541 TRefArray *reference = uArray->GetUnitTrackRef();
542 reference->Add(track);
544 Float_t unitEnergy2 = uArray->GetUnitEnergy();
545 Bool_t okout = kFALSE;
546 Bool_t refout = kFALSE;
552 // Fill energy outside emcal acceptance
553 uArray->SetUnitEnergy(unitEnergy2 + pt);
556 if(uArray->GetUnitEnergy()<ptMin){
557 uArray->SetUnitCutFlag(kPtSmaller);
560 uArray->SetUnitCutFlag(kPtHigher);
561 if(okout) nTracksTpcCut++;
564 if(sflag[goodTrack] == 1) {
565 uArray->SetUnitSignalFlag(kGood);
566 uArray->SetUnitSignalFlagC(kFALSE,kGood);
567 } else uArray->SetUnitSignalFlagC(kFALSE,kBad);
570 if(uArray->GetUnitEnergy()>0 && refout){
572 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
575 fRefArray->Add(uArray);
582 } // end loop on entries (tpc tracks)
586 cout << "End of Tracks %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << endl;
587 cout << "goodTracks: " << goodTrack << endl;
594 fNTracks = nTracksTpcOnly;
595 fNTracksCut = nTracksTpcOnlyCut;
597 cout << "fNTracks : " << fNTracks << endl;
598 cout << "fNTracksCut : " << fNTracksCut << endl;
602 fNTracks = nTracksEmcal+nTracksEmcalDZ+nTracksTpc;
603 fNTracksCut = nTracksEmcalCut+nTracksEmcalDZCut+nTracksTpcCut;
605 cout << "fNTracks : " << fNTracks << endl;
606 cout << "fNTracksCut : " << fNTracksCut << endl;