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 //----------------------------------------------------------------------
18 // Fill Unit Array class
19 // Class used by AliJetESDReader to fill a UnitArray from the information
20 // extracted from the particle tracks
21 // Author: magali.estienne@ires.in2p3.fr
22 //----------------------------------------------------------------------
25 // --- Standard library ---
26 #include <Riostream.h>
28 // --- ROOT system ---
30 #include <TLorentzVector.h>
31 #include <TRefArray.h>
34 #include <TGeoManager.h>
38 #include <TClonesArray.h>
40 // --- AliRoot header files ---
41 #include "AliJetFinder.h"
42 #include "AliJetReaderHeader.h"
43 #include "AliJetReader.h"
44 #include "AliJetESDReader.h"
45 #include "AliJetESDReaderHeader.h"
47 #include "AliESDEvent.h"
48 #include "AliJetDummyGeo.h"
49 #include "AliJetUnitArray.h"
50 #include "AliJetFillUnitArrayTracks.h"
51 #include "AliJetHadronCorrectionv1.h"
52 #include "AliJetGrid.h"
54 ClassImp(AliJetFillUnitArrayTracks)
56 //_____________________________________________________________________________
57 AliJetFillUnitArrayTracks::AliJetFillUnitArrayTracks()
58 : TTask("AliJetFillUnitArrayTracks","Fill Unit Array with tpc/its and emcal information"),
100 fPhiBinInEMCalAcc(0),
106 //_____________________________________________________________________________
107 AliJetFillUnitArrayTracks::AliJetFillUnitArrayTracks(AliESDEvent* /*esd*/)
108 : TTask("AliJetFillUnitArrayTracks","Fill Unit Array with tpc/its and emcal information"),
149 fEtaBinInEMCalAcc(0),
150 fPhiBinInEMCalAcc(0),
156 //____________________________________________________________________________
157 void AliJetFillUnitArrayTracks::InitParameters()
159 // fHCorrection = 0; // For hadron correction
160 fHadCorr = 0; // For hadron correction
161 fNumUnits = fGeom->GetNCells(); // Number of towers in EMCAL
162 fDebug = fReaderHeader->GetDebug();
164 fEtaMinCal = fGeom->GetArm1EtaMin();
165 fEtaMaxCal = fGeom->GetArm1EtaMax();
166 fPhiMinCal = fGeom->GetArm1PhiMin();
167 fPhiMaxCal = fGeom->GetArm1PhiMax()-10.; // A verifier quelle doit etre la derniere valeur !
169 fTPCGrid->GetAccParam(fNphi,fNeta,fPhiMin,
170 fPhiMax,fEtaMin,fEtaMax);
171 fTPCGrid->GetBinParam(fPhiBinInTPCAcc,fEtaBinInTPCAcc,
172 fPhiBinInEMCalAcc,fEtaBinInEMCalAcc,fNbinPhi);
174 fEta = fTPCGrid->GetArrayEta();
175 fPhi = fTPCGrid->GetArrayPhi();
176 fIndex = fTPCGrid->GetIndexObject();
179 for(Int_t i=0; i<fNphi+1; i++) cout << "phi[" << i << "] : " << (*fPhi)[i] << endl;
180 for(Int_t i=0; i<fNeta+1; i++) cout << "eta[" << i << "] : " << (*fEta)[i] << endl;
182 for(Int_t i=0; i<fNphi+1; i++)
183 for(Int_t j=0; j<fNeta+1; j++) {cout << "fIndex[" << i << "," << j << "] : " <<
184 (*fIndex)(i,j) << endl; }
186 if(fDebug>1) printf("\n Parameters initiated ! \n");
189 //_____________________________________________________________________________
190 AliJetFillUnitArrayTracks::~AliJetFillUnitArrayTracks()
195 //_____________________________________________________________________________
196 void AliJetFillUnitArrayTracks::Exec(Option_t* /*option*/)
202 fDebug = fReaderHeader->GetDebug();
207 // get number of tracks in event (for the loop)
212 nt = fESD->GetNumberOfTracks();
213 if(fDebug>1) cout << "Number of Tracks in ESD : " << nt << endl;
215 // temporary storage of signal and pt cut flag
216 Int_t* sflag = new Int_t[nt];
217 Int_t* cflag = new Int_t[nt];
219 // get cuts set by user
220 Float_t ptMin = fReaderHeader->GetPtCut();
221 Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
222 Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
223 fOpt = fReaderHeader->GetDetector();
224 fDZ = fReaderHeader->GetDZ();
226 Int_t nTracksEmcal = 0;
227 Int_t nTracksEmcalDZ = 0;
228 Int_t nTracksTpc = 0;
229 Int_t nTracksTpcOnly = 0;
230 Int_t nTracksEmcalCut = 0;
231 Int_t nTracksEmcalDZCut = 0;
232 Int_t nTracksTpcCut = 0;
233 Int_t nTracksTpcOnlyCut = 0;
235 fGrid = fTPCGrid->GetGridType();
239 for (Int_t it = 0; it < nt; it++) {
240 AliESDtrack *track = fESD->GetTrack(it);
241 UInt_t status = track->GetStatus();
244 track->GetPxPyPz(mom);
245 p3.SetXYZ(mom[0],mom[1],mom[2]);
248 mass = track->GetMass();
250 if (((status & AliESDtrack::kITSrefit) == 0) ||
251 ((status & AliESDtrack::kTPCrefit) == 0)) continue; // quality check
252 if (((AliJetESDReaderHeader*) fReaderHeader)->ReadSignalOnly()
253 && TMath::Abs(track->GetLabel()) > 10000) continue; // quality check
254 if (((AliJetESDReaderHeader*) fReaderHeader)->ReadBkgdOnly()
255 && TMath::Abs(track->GetLabel()) < 10000) continue; // quality check
257 phi = ( (p3.Phi()) < 0) ? (p3.Phi()) + 2. * TMath::Pi() : (p3.Phi());
259 if ( (eta > etaMax) || (eta < etaMin)) continue; // checking eta cut
261 // sflag -> not yet implemented !!!!
265 // Only TPC filled from its grid in its total acceptance
267 Int_t idTPC = fTPCGrid->GetIndex(phi,eta);
270 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(idTPC-1);
271 uArray->SetUnitTrackID(it);
273 Float_t unitEnergy = 0.;
274 unitEnergy = uArray->GetUnitEnergy();
279 // Fill energy in TPC acceptance
280 uArray->SetUnitEnergy(unitEnergy + pt);
281 uArray->SetUnitPxPyPz(mom);
282 uArray->SetUnitMass(mass);
285 if(uArray->GetUnitEnergy()<ptMin){
286 uArray->SetUnitCutFlag(kPtSmaller);
289 uArray->SetUnitCutFlag(kPtHigher);
290 if(ok) nTracksTpcOnlyCut++;
295 uArray->SetUnitDetectorFlag(kAll);
297 if(uArray->GetUnitEnergy()>0){
298 fRefArray->Add(uArray);
302 if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
304 if (pt > ptMin) cflag[goodTrack]=1; // pt cut
311 Int_t nElements = fTPCGrid->GetNEntries();
313 // Fill track information in EMCAL acceptance
314 if((eta >= fEtaMin && eta <= fEtaMax) &&
315 (phi >= fPhiMin && phi <= fPhiMax))// &&
318 // Include dead-zones
321 Double_t phimin0 = 0., phimin1 = 0., phimin2 = 0., phimin3 = 0., phimin4 = 0.;
322 Double_t phimax0 = 0., phimax1 = 0., phimax2 = 0., phimax3 = 0., phimax4 = 0.;
323 fGeom->GetPhiBoundariesOfSMGap(0,phimin0,phimax0);
324 fGeom->GetPhiBoundariesOfSMGap(1,phimin1,phimax1);
325 fGeom->GetPhiBoundariesOfSMGap(2,phimin2,phimax2);
326 fGeom->GetPhiBoundariesOfSMGap(3,phimin3,phimax3);
327 fGeom->GetPhiBoundariesOfSMGap(4,phimin4,phimax4);
328 Int_t n0 = fGrid0->GetNEntries();
329 Int_t n1 = fGrid1->GetNEntries();
330 Int_t n2 = fGrid2->GetNEntries();
331 Int_t n3 = fGrid3->GetNEntries();
333 if(phi >= phimin0 && phi <= phimax0){
334 Int_t id0 = fGrid0->GetIndex(phi,eta)-1;
335 AliJetUnitArray *uArray0 = (AliJetUnitArray*)fUnitArray->At(id0+fNumUnits+nElements);
336 uArray0->SetUnitTrackID(it);
337 Float_t uEnergy0 = uArray0->GetUnitEnergy();
343 uArray0->SetUnitEnergy(uEnergy0+pt);
344 if(uArray0->GetUnitEnergy()<ptMin)
345 uArray0->SetUnitCutFlag(kPtSmaller);
347 uArray0->SetUnitCutFlag(kPtHigher);
348 if(ok0) nTracksEmcalDZCut++;
350 if(uArray0->GetUnitEnergy()>0)
351 fRefArray->Add(uArray0);
353 if(phi >= phimin1 && phi <= phimax1){
354 Int_t id1 = fGrid1->GetIndex(phi,eta)-1+n0;
355 AliJetUnitArray *uArray1 = (AliJetUnitArray*)fUnitArray->At(id1+fNumUnits+nElements);
356 uArray1->SetUnitTrackID(it);
357 Float_t uEnergy1 = uArray1->GetUnitEnergy();
363 uArray1->SetUnitEnergy(uEnergy1+pt);
364 if(uArray1->GetUnitEnergy()<ptMin)
365 uArray1->SetUnitCutFlag(kPtSmaller);
367 uArray1->SetUnitCutFlag(kPtHigher);
368 if(ok1) nTracksEmcalDZCut++;
370 if(uArray1->GetUnitEnergy()>0)
371 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 uArray2->SetUnitTrackID(it);
377 Float_t uEnergy2 = uArray2->GetUnitEnergy();
383 uArray2->SetUnitEnergy(uEnergy2+pt);
384 if(uArray2->GetUnitEnergy()<ptMin)
385 uArray2->SetUnitCutFlag(kPtSmaller);
387 uArray2->SetUnitCutFlag(kPtHigher);
388 if(ok2) nTracksEmcalDZCut++;
390 if(uArray2->GetUnitEnergy()>0)
391 fRefArray->Add(uArray2);
393 if(phi >= phimin3 && phi <= phimax3){
394 Int_t id3 = fGrid3->GetIndex(phi,eta)-1+n0+n1+n2;
395 AliJetUnitArray *uArray3 = (AliJetUnitArray*)fUnitArray->At(id3+fNumUnits+nElements);
396 uArray3->SetUnitTrackID(it);
397 Float_t uEnergy3 = uArray3->GetUnitEnergy();
403 uArray3->SetUnitEnergy(uEnergy3+pt);
404 if(uArray3->GetUnitEnergy()<ptMin)
405 uArray3->SetUnitCutFlag(kPtSmaller);
407 uArray3->SetUnitCutFlag(kPtHigher);
408 if(ok3) nTracksEmcalDZCut++;
410 if(uArray3->GetUnitEnergy()>0)
411 fRefArray->Add(uArray3);
413 if(phi >= phimin4 && phi <= phimax4){
414 Int_t id4 = fGrid4->GetIndex(phi,eta)-1+n0+n1+n2+n3;
415 AliJetUnitArray *uArray4 = (AliJetUnitArray*)fUnitArray->At(id4+fNumUnits+nElements);
416 uArray4->SetUnitTrackID(it);
417 Float_t uEnergy4 = uArray4->GetUnitEnergy();
423 uArray4->SetUnitEnergy(uEnergy4+pt);
424 if(uArray4->GetUnitEnergy()<ptMin)
425 uArray4->SetUnitCutFlag(kPtSmaller);
427 uArray4->SetUnitCutFlag(kPtHigher);
428 if(ok4) nTracksEmcalDZCut++;
430 if(uArray4->GetUnitEnergy()>0)
431 fRefArray->Add(uArray4);
436 fGeom->GetAbsCellIdFromEtaPhi(eta,phi,towerID);
438 if(towerID==-1) continue;
440 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(towerID);
441 uArray->SetUnitTrackID(it);
442 Float_t unitEnergy = uArray->GetUnitEnergy();
449 // Do Hadron Correction
450 // Parametrization to be added
451 if (fHCorrection != 0)
453 // Float_t hCEnergy = fHadCorr->GetEnergy(enT[i], (Double_t)etaT[i]);
454 Float_t hCEnergy = fHadCorr->GetEnergy(p3.Mag(), (Double_t)eta);
455 unitEnergy -= hCEnergy*TMath::Sin(2.0*TMath::ATan(TMath::Exp(-eta)));
457 } //end Hadron Correction loop
459 uArray->SetUnitEnergy(unitEnergy + pt);
462 if(uArray->GetUnitEnergy()<ptMin)
463 uArray->SetUnitCutFlag(kPtSmaller);
465 uArray->SetUnitCutFlag(kPtHigher);
466 if(ok) nTracksEmcalCut++;
470 uArray->SetUnitDetectorFlag(kAll);
472 if(uArray->GetUnitEnergy()>0)
473 fRefArray->Add(uArray);
476 if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
478 if (pt > ptMin) cflag[goodTrack]=1; // pt cut
481 } // end loop on EMCal acceptance cut + tracks quality
483 // Outside EMCal acceptance
485 // Int_t idTPC = GetIndexFromEtaPhi(etaT[i],phiT[i]);
486 Int_t idTPC = fTPCGrid->GetIndex(phi,eta);
488 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(fNumUnits-1+idTPC);
489 uArray->SetUnitTrackID(it);
491 Float_t unitEnergy2 = uArray->GetUnitEnergy(); // check if fNumUnits or fNumUnits-1
497 // Fill energy outside emcal acceptance
498 uArray->SetUnitEnergy(unitEnergy2 + pt);
501 if(uArray->GetUnitEnergy()<ptMin){
502 uArray->SetUnitCutFlag(kPtSmaller);
505 uArray->SetUnitCutFlag(kPtHigher);
506 if(ok2) nTracksTpcCut++;
510 uArray->SetUnitDetectorFlag(kTpc);
511 if(uArray->GetUnitEnergy()>0)
512 fRefArray->Add(uArray);
515 if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
517 if (pt > ptMin) cflag[goodTrack]=1; // pt cut
522 } // end loop on entries (tpc tracks)
524 // // set the signal flags
525 // fSignalFlag.Set(goodTrack,sflag);
526 // fCutFlag.Set(goodTrack,cflag);
531 // } // end loop on entries (tpc tracks)
534 fNTracks = nTracksTpcOnly;
535 fNTracksCut = nTracksTpcOnlyCut;
537 cout << "fNTracks : " << fNTracks << endl;
538 cout << "fNTracksCut : " << fNTracksCut << endl;
542 fNTracks = nTracksEmcal+nTracksEmcalDZ+nTracksTpc;
543 fNTracksCut = nTracksEmcalCut+nTracksEmcalDZCut+nTracksTpcCut;
545 cout << "fNTracks : " << fNTracks << endl;
546 cout << "fNTracksCut : " << fNTracksCut << endl;
552 //__________________________________________________________
553 void AliJetFillUnitArrayTracks::GetEtaPhiFromIndex(Int_t index, Float_t &eta, Float_t &phi)
556 // Obtain the index from eta and phi
558 for(Int_t j=0; j<fNphi+1; j++) {
559 for(Int_t i=0; i<fNeta+1; i++) {
562 //-------------------------------------
564 if(j*(fNeta+1)+i == index) {
571 //-------------------------------------
574 if(i>0 && i<(fEtaBinInTPCAcc-fEtaBinInEMCalAcc)/2) ii = i;
575 if(i>=(fEtaBinInTPCAcc+fEtaBinInEMCalAcc)/2 && i<fNeta+1) ii = i-fEtaBinInEMCalAcc;
578 if(j<(fNbinPhi+1) && j*(fNeta+1)+i == index) {
583 if((j>=(fNbinPhi+1) && j<(fNbinPhi+1+fPhiBinInEMCalAcc)) &&
584 ((fNbinPhi+1)*(fNeta+1) + (j-fNbinPhi-1)*(fEtaBinInTPCAcc-fEtaBinInEMCalAcc) + ii)== index ) {
585 if(ii==0) {Int_t ind = 0; eta = fEta2->At(ind);}
586 else eta = fEta2->At(i);
590 if(j>=(fNbinPhi+1+fPhiBinInEMCalAcc) && ((fNbinPhi+1)*(fNeta+1)+fPhiBinInEMCalAcc*((fEtaBinInTPCAcc-fEtaBinInEMCalAcc))+(j-(fNbinPhi+1+fPhiBinInEMCalAcc))*(fNeta+1)+i == index)) {