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 "AliJetHadronCorrectionv1.h"
33 #include "AliJetESDFillUnitArrayTracks.h"
35 // --- ROOT system ---
40 // --- AliRoot header files ---
43 class AliJetESDReader;
45 ClassImp(AliJetESDFillUnitArrayTracks)
47 //_____________________________________________________________________________
48 AliJetESDFillUnitArrayTracks::AliJetESDFillUnitArrayTracks()
49 : AliJetFillUnitArray(),
52 fApplyMIPCorrection(kTRUE),
63 //_____________________________________________________________________________
64 AliJetESDFillUnitArrayTracks::AliJetESDFillUnitArrayTracks(AliESDEvent* esd)
65 : AliJetFillUnitArray(),
68 fApplyMIPCorrection(kTRUE),
79 //_____________________________________________________________________________
80 AliJetESDFillUnitArrayTracks::AliJetESDFillUnitArrayTracks(const AliJetESDFillUnitArrayTracks &det)
81 : AliJetFillUnitArray(det),
82 fNumUnits(det.fNumUnits),
83 fHadCorr(det.fHadCorr),
84 fApplyMIPCorrection(det.fApplyMIPCorrection),
95 //------------------------------------------------------------------
96 AliJetESDFillUnitArrayTracks& AliJetESDFillUnitArrayTracks::operator=(const AliJetESDFillUnitArrayTracks& other)
100 fNumUnits = other.fNumUnits;
101 fHadCorr = other.fHadCorr;
102 fApplyMIPCorrection = other.fApplyMIPCorrection;
104 fGrid0 = other.fGrid0;
105 fGrid1 = other.fGrid1;
106 fGrid2 = other.fGrid2;
107 fGrid3 = other.fGrid3;
108 fGrid4 = other.fGrid4;
113 //____________________________________________________________________________
114 void AliJetESDFillUnitArrayTracks::InitParameters()
116 // fHadCorr = 0; // For hadron correction
117 fNumUnits = fGeom->GetNCells(); // Number of towers in EMCAL
118 cout << "In AliJetESDFillUnitArrayTracks:InitParameters(), Ncells : " << fNumUnits << endl;
120 fTPCGrid->GetAccParam(fNphi,fNeta,fPhiMin,
121 fPhiMax,fEtaMin,fEtaMax);
122 fTPCGrid->GetBinParam(fPhiBinInTPCAcc,fEtaBinInTPCAcc,
123 fPhiBinInEMCalAcc,fEtaBinInEMCalAcc,fNbinPhi);
125 fIndex = fTPCGrid->GetIndexObject();
128 for(Int_t i=0; i<fNphi+1; i++)
129 for(Int_t j=0; j<fNeta+1; j++) {cout << "fIndex[" << i << "," << j << "] : " <<
130 (*fIndex)(i,j) << endl; }
132 if(fDebug>1) printf("\n Parameters initiated ! \n");
135 //_____________________________________________________________________________
136 AliJetESDFillUnitArrayTracks::~AliJetESDFillUnitArrayTracks()
141 //_____________________________________________________________________________
142 void AliJetESDFillUnitArrayTracks::Exec(Option_t* const /*option*/)
146 // Fill the unit array with the charged particle information in ESD
149 fDebug = fReaderHeader->GetDebug();
155 // get number of tracks in event (for the loop)
158 //(not used ?) Int_t nt2 = 0;
164 if(fDebug>1) cout << "fESD in Fill array : " << fESD << endl;
166 nt = fESD->GetNumberOfTracks();
167 if(fDebug>1) cout << "Number of Tracks in ESD : " << nt << endl;
169 // temporary storage of signal and pt cut flag
170 Int_t* sflag = new Int_t[nt];
171 Int_t* cflag = new Int_t[nt];
173 // get cuts set by user
174 Float_t ptMin = fReaderHeader->GetPtCut();
175 Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
176 Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
177 fOpt = fReaderHeader->GetDetector();
178 fDZ = fReaderHeader->GetDZ();
180 Int_t nTracksEmcal = 0;
181 Int_t nTracksEmcalDZ = 0;
182 Int_t nTracksTpc = 0;
183 Int_t nTracksTpcOnly = 0;
184 Int_t nTracksEmcalCut = 0;
185 Int_t nTracksEmcalDZCut = 0;
186 Int_t nTracksTpcCut = 0;
187 Int_t nTracksTpcOnlyCut = 0;
189 //(not used ?) Int_t nElements = fTPCGrid->GetNEntries();
190 //(not used ?) Int_t nElements2 = fEMCalGrid->GetNEntries();
191 fGrid = fTPCGrid->GetGridType();
193 //(not used ?) Int_t nEntries = fUnitArray->GetEntries();
195 //(not used ?) Int_t nRefEnt = fRefArray->GetEntries();
199 for (Int_t it = 0; it < nmax; it++) {
201 track = fESD->GetTrack(it);
202 UInt_t status = track->GetStatus();
205 track->GetPxPyPz(mom);
207 p3.SetXYZ(mom[0],mom[1],mom[2]);
211 mass = track->GetMass();
213 if ((status & AliESDtrack::kTPCrefit) == 0) continue; // quality check
214 if ((status & AliESDtrack::kITSrefit) == 0) continue; // quality check
215 if (((AliJetESDReaderHeader*) fReaderHeader)->ReadSignalOnly()
216 && TMath::Abs(track->GetLabel()) > 10000) continue; // quality check
217 if (((AliJetESDReaderHeader*) fReaderHeader)->ReadBkgdOnly()
218 && TMath::Abs(track->GetLabel()) < 10000) continue; // quality check
220 phi = ( (p3.Phi()) < 0) ? (p3.Phi()) + 2. * TMath::Pi() : (p3.Phi());
222 if ( (eta > etaMax) || (eta < etaMin)) continue; // checking eta cut
225 if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
227 if (pt > ptMin) cflag[goodTrack]=1; // pt cut
232 // Only TPC filled from its grid in its total acceptance
234 Int_t idTPC = fTPCGrid->GetIndex(phi,eta);
238 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(idTPC-1);
239 TRefArray *reference = uArray->GetUnitTrackRef();
240 if (reference->GetEntries() == 0) {
241 new(reference) TRefArray(TProcessID::GetProcessWithUID(track));
244 reference->Add(track);
246 Float_t unitEnergy = 0.;
247 unitEnergy = uArray->GetUnitEnergy();
248 // nTracksTpcOnly is necessary to count the number of candidate cells
249 // but it doesn't give the real multiplicity -> it will be extracted
250 // in the jet finder through the reference to tracks
257 // Fill energy in TPC acceptance
258 uArray->SetUnitEnergy(unitEnergy + pt);
259 uArray->SetUnitMass(mass);
262 if(uArray->GetUnitEnergy()<ptMin){
263 uArray->SetUnitCutFlag(kPtSmaller);
266 uArray->SetUnitCutFlag(kPtHigher);
267 if(ok) nTracksTpcOnlyCut++;
271 if(sflag[goodTrack] == 1) {
272 uArray->SetUnitSignalFlag(kGood);
273 uArray->SetUnitSignalFlagC(kFALSE,kGood);
274 } else uArray->SetUnitSignalFlagC(kFALSE,kBad);
276 if(uArray->GetUnitEnergy()>0 && ref){
279 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
281 fRefArray->Add(uArray);
287 Int_t nElements = fTPCGrid->GetNEntries();
288 // Fill track information in EMCAL acceptance
289 if((eta >= fEtaMin && eta <= fEtaMax) &&
290 (phi >= fPhiMin && phi <= fPhiMax))// &&
293 // Include dead-zones
296 Double_t phimin0 = 0., phimin1 = 0., phimin2 = 0., phimin3 = 0., phimin4 = 0.;
297 Double_t phimax0 = 0., phimax1 = 0., phimax2 = 0., phimax3 = 0., phimax4 = 0.;
298 fGeom->GetPhiBoundariesOfSMGap(0,phimin0,phimax0);
299 fGeom->GetPhiBoundariesOfSMGap(1,phimin1,phimax1);
300 fGeom->GetPhiBoundariesOfSMGap(2,phimin2,phimax2);
301 fGeom->GetPhiBoundariesOfSMGap(3,phimin3,phimax3);
302 fGeom->GetPhiBoundariesOfSMGap(4,phimin4,phimax4);
303 Int_t n0 = fGrid0->GetNEntries();
304 Int_t n1 = fGrid1->GetNEntries();
305 Int_t n2 = fGrid2->GetNEntries();
306 Int_t n3 = fGrid3->GetNEntries();
307 //(not used ?) Int_t n4 = fGrid4->GetNEntries();
309 if(phi >= phimin0 && phi <= phimax0){
310 Int_t id0 = fGrid0->GetIndex(phi,eta)-1;
311 AliJetUnitArray *uArray0 = (AliJetUnitArray*)fUnitArray->At(id0+fNumUnits+nElements);
312 TRefArray *reference0 = uArray0->GetUnitTrackRef();
313 if (reference0->GetEntries() == 0) {
314 new(reference0) TRefArray(TProcessID::GetProcessWithUID(track));
317 reference0->Add(track);
319 Float_t uEnergy0 = uArray0->GetUnitEnergy();
321 Bool_t ref0 = kFALSE;
327 uArray0->SetUnitEnergy(uEnergy0+pt);
329 if(uArray0->GetUnitEnergy()<ptMin)
330 uArray0->SetUnitCutFlag(kPtSmaller);
332 uArray0->SetUnitCutFlag(kPtHigher);
333 if(ok0) nTracksEmcalDZCut++;
335 if(sflag[goodTrack] == 1) {
336 uArray0->SetUnitSignalFlag(kGood);
337 uArray0->SetUnitSignalFlagC(kFALSE,kGood);
338 } else uArray0->SetUnitSignalFlagC(kFALSE,kBad);
340 if(uArray0->GetUnitEnergy()>0 && ref0){
342 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray0));
345 fRefArray->Add(uArray0);
349 if(phi >= phimin1 && phi <= phimax1){
350 Int_t id1 = fGrid1->GetIndex(phi,eta)-1+n0;
351 AliJetUnitArray *uArray1 = (AliJetUnitArray*)fUnitArray->At(id1+fNumUnits+nElements);
352 TRefArray *reference1 = uArray1->GetUnitTrackRef();
353 if (reference1->GetEntries() == 0) {
354 new(reference1) TRefArray(TProcessID::GetProcessWithUID(track));
357 reference1->Add(track);
359 Float_t uEnergy1 = uArray1->GetUnitEnergy();
361 Bool_t ref1 = kFALSE;
367 uArray1->SetUnitEnergy(uEnergy1+pt);
369 if(uArray1->GetUnitEnergy()<ptMin)
370 uArray1->SetUnitCutFlag(kPtSmaller);
372 uArray1->SetUnitCutFlag(kPtHigher);
373 if(ok1) nTracksEmcalDZCut++;
375 if(sflag[goodTrack] == 1) {
376 uArray1->SetUnitSignalFlag(kGood);
377 uArray1->SetUnitSignalFlagC(kFALSE,kGood);
378 } else uArray1->SetUnitSignalFlagC(kFALSE,kBad);
380 if(uArray1->GetUnitEnergy()>0 && ref1){
382 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray1));
385 fRefArray->Add(uArray1);
389 if(phi >= phimin2 && phi <= phimax2){
390 Int_t id2 = fGrid2->GetIndex(phi,eta)-1+n0+n1;
391 AliJetUnitArray *uArray2 = (AliJetUnitArray*)fUnitArray->At(id2+fNumUnits+nElements);
392 TRefArray *reference2 = uArray2->GetUnitTrackRef();
393 if (reference2->GetEntries() == 0) {
394 new(reference2) TRefArray(TProcessID::GetProcessWithUID(track));
397 reference2->Add(track);
399 Float_t uEnergy2 = uArray2->GetUnitEnergy();
401 Bool_t ref2 = kFALSE;
407 uArray2->SetUnitEnergy(uEnergy2+pt);
409 if(uArray2->GetUnitEnergy()<ptMin)
410 uArray2->SetUnitCutFlag(kPtSmaller);
412 uArray2->SetUnitCutFlag(kPtHigher);
413 if(ok2) nTracksEmcalDZCut++;
415 if(sflag[goodTrack] == 1) {
416 uArray2->SetUnitSignalFlag(kGood);
417 uArray2->SetUnitSignalFlagC(kFALSE,kGood);
418 } else uArray2->SetUnitSignalFlagC(kFALSE,kBad);
420 if(uArray2->GetUnitEnergy()>0 && ref2){
422 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray2));
425 fRefArray->Add(uArray2);
429 if(phi >= phimin3 && phi <= phimax3){
430 Int_t id3 = fGrid3->GetIndex(phi,eta)-1+n0+n1+n2;
431 AliJetUnitArray *uArray3 = (AliJetUnitArray*)fUnitArray->At(id3+fNumUnits+nElements);
432 TRefArray *reference3 = uArray3->GetUnitTrackRef();
433 if (reference3->GetEntries() == 0) {
434 new(reference3) TRefArray(TProcessID::GetProcessWithUID(track));
437 reference3->Add(track);
439 Float_t uEnergy3 = uArray3->GetUnitEnergy();
441 Bool_t ref3 = kFALSE;
447 uArray3->SetUnitEnergy(uEnergy3+pt);
449 if(uArray3->GetUnitEnergy()<ptMin)
450 uArray3->SetUnitCutFlag(kPtSmaller);
452 uArray3->SetUnitCutFlag(kPtHigher);
453 if(ok3) nTracksEmcalDZCut++;
455 if(sflag[goodTrack] == 1) {
456 uArray3->SetUnitSignalFlag(kGood);
457 uArray3->SetUnitSignalFlagC(kFALSE,kGood);
458 } else uArray3->SetUnitSignalFlagC(kFALSE,kBad);
460 if(uArray3->GetUnitEnergy()>0 && ref3){
462 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray3));
465 fRefArray->Add(uArray3);
469 if(phi >= phimin4 && phi <= phimax4){
470 Int_t id4 = fGrid4->GetIndex(phi,eta)-1+n0+n1+n2+n3;
471 AliJetUnitArray *uArray4 = (AliJetUnitArray*)fUnitArray->At(id4+fNumUnits+nElements);
472 TRefArray *reference4 = uArray4->GetUnitTrackRef();
473 if (reference4->GetEntries() == 0) {
474 new(reference4) TRefArray(TProcessID::GetProcessWithUID(track));
477 reference4->Add(track);
479 Float_t uEnergy4 = uArray4->GetUnitEnergy();
481 Bool_t ref4 = kFALSE;
487 uArray4->SetUnitEnergy(uEnergy4+pt);
489 if(uArray4->GetUnitEnergy()<ptMin)
490 uArray4->SetUnitCutFlag(kPtSmaller);
492 uArray4->SetUnitCutFlag(kPtHigher);
493 if(ok4) nTracksEmcalDZCut++;
495 if(sflag[goodTrack] == 1) {
496 uArray4->SetUnitSignalFlag(kGood);
497 uArray4->SetUnitSignalFlagC(kFALSE,kGood);
498 } else uArray4->SetUnitSignalFlagC(kFALSE,kBad);
500 if(uArray4->GetUnitEnergy()>0 && ref4){
502 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray4));
505 fRefArray->Add(uArray4);
511 fGeom->GetAbsCellIdFromEtaPhi(eta,phi,towerID);
512 if(towerID==-1) continue;
514 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(towerID);
515 TRefArray *reference = uArray->GetUnitTrackRef();
516 if (reference->GetEntries() == 0) {
517 new(reference) TRefArray(TProcessID::GetProcessWithUID(track));
520 reference->Add(track);
522 Float_t unitEnergy = uArray->GetUnitEnergy();
531 // Do Hadron Correction
532 // This is under construction !!!!!!!!!!!!!!!!!!!!!!!
533 // For the moment I apply MIP correction if p >= 0.5 GeV/c
534 if (fApplyMIPCorrection != 0 && p3.Mag() >= 0.5)
536 ((AliJetHadronCorrectionv1*)fHadCorr)->SetGeometry("EMCAL_COMPLETE",1.);
538 // Get track position at the outer part of the reconstruction ~ TRD
539 Double_t phiOut = track->GetOuterParam()->Phi();
540 Double_t etaOut = track->GetOuterParam()->Eta();
542 // If the track in the outer part of the TPC/TDR ? is inside
543 // the calorimeter, it can deposit part of its energy
544 // We can then correct on average for these particles
545 if((etaOut >= fEtaMin && etaOut <= fEtaMax) &&
546 (phiOut >= fPhiMin && phiOut <= fPhiMax))// &&
548 Double_t hCEnergy = (Double_t)fHadCorr->GetEnergy(p3.Mag(), (Double_t)eta,0);
549 unitEnergy -= hCEnergy*TMath::Sin(2.0*TMath::ATan(TMath::Exp(-eta)));
551 } //end Hadron Correction loop
553 cout << "unitEnergy + pt = " << unitEnergy << " + " << pt << " = " << unitEnergy + pt << endl;
555 if((unitEnergy + pt) > 0.) uArray->SetUnitEnergy(unitEnergy + pt);
556 else uArray->SetUnitEnergy(0.);
559 if(uArray->GetUnitEnergy()<ptMin){
560 uArray->SetUnitCutFlag(kPtSmaller);
563 uArray->SetUnitCutFlag(kPtHigher);
564 if(ok) nTracksEmcalCut++;
568 if(sflag[goodTrack] == 1) {
569 uArray->SetUnitSignalFlag(kGood);
570 uArray->SetUnitSignalFlagC(kFALSE,kGood);
571 } else uArray->SetUnitSignalFlagC(kFALSE,kBad);
574 if(uArray->GetUnitEnergy()>0 && ref){
576 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
579 fRefArray->Add(uArray);
581 } // end loop on EMCal acceptance cut + tracks quality
583 // Outside EMCal acceptance
585 Int_t idTPC = fTPCGrid->GetIndex(phi,eta);
587 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(fNumUnits-1+idTPC);
588 TRefArray *reference = uArray->GetUnitTrackRef();
589 if (reference->GetEntries() == 0) {
590 new(reference) TRefArray(TProcessID::GetProcessWithUID(track));
593 reference->Add(track);
595 Float_t unitEnergy2 = uArray->GetUnitEnergy();
596 Bool_t okout = kFALSE;
597 Bool_t refout = kFALSE;
603 // Fill energy outside emcal acceptance
604 uArray->SetUnitEnergy(unitEnergy2 + pt);
607 if(uArray->GetUnitEnergy()<ptMin){
608 uArray->SetUnitCutFlag(kPtSmaller);
611 uArray->SetUnitCutFlag(kPtHigher);
612 if(okout) nTracksTpcCut++;
615 if(sflag[goodTrack] == 1) {
616 uArray->SetUnitSignalFlag(kGood);
617 uArray->SetUnitSignalFlagC(kFALSE,kGood);
618 } else uArray->SetUnitSignalFlagC(kFALSE,kBad);
621 if(uArray->GetUnitEnergy()>0 && refout){
623 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
626 fRefArray->Add(uArray);
633 } // end loop on entries (tpc tracks)
637 cout << "End of Tracks %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << endl;
638 cout << "goodTracks: " << goodTrack << endl;
641 // fSignalFlag.Set(goodTrack,sflag);
642 // fCutFlag.Set(goodTrack,cflag);
648 fNTracks = nTracksTpcOnly;
649 fNTracksCut = nTracksTpcOnlyCut;
651 cout << "fNTracks : " << fNTracks << endl;
652 cout << "fNTracksCut : " << fNTracksCut << endl;
656 fNTracks = nTracksEmcal+nTracksEmcalDZ+nTracksTpc;
657 fNTracksCut = nTracksEmcalCut+nTracksEmcalDZCut+nTracksTpcCut;
659 cout << "fNTracks : " << fNTracks << endl;
660 cout << "fNTracksCut : " << fNTracksCut << endl;