]> git.uio.no Git - u/mrichter/AliRoot.git/blame - JETAN/AliJetFillUnitArrayTracks.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Coding violations corrected.
[u/mrichter/AliRoot.git] / JETAN / AliJetFillUnitArrayTracks.cxx
CommitLineData
4a01bb2c 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
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 **************************************************************************/
15
ee7de0dd 16
8838ab7a 17//======================================================================
18// ***July 2006
4a01bb2c 19// Fill Unit Array class
8838ab7a 20// Class used by AliJetESDReader to fill a UnitArray from the information extracted
21// from the particle tracks
4a01bb2c 22// Author: magali.estienne@ires.in2p3.fr
8838ab7a 23//======================================================================
4a01bb2c 24
25
26// --- Standard library ---
27#include <Riostream.h>
28
29// --- ROOT system ---
30#include <TSystem.h>
31#include <TLorentzVector.h>
ee7de0dd 32#include <TRefArray.h>
4a01bb2c 33#include <TVector3.h>
4a01bb2c 34#include "TTask.h"
35#include <TGeoManager.h>
36#include <TMatrixD.h>
37#include <TArrayD.h>
ee7de0dd 38#include <TMath.h>
39#include <TClonesArray.h>
8838ab7a 40#include <TProcessID.h>
4a01bb2c 41
42// --- AliRoot header files ---
43#include "AliJetFinder.h"
44#include "AliJetReaderHeader.h"
45#include "AliJetReader.h"
46#include "AliJetESDReader.h"
47#include "AliJetESDReaderHeader.h"
ee7de0dd 48#include "AliESDEvent.h"
b92e2ccf 49#include "AliJetDummyGeo.h"
4a01bb2c 50#include "AliJetUnitArray.h"
51#include "AliJetFillUnitArrayTracks.h"
7aec0427 52#include "AliJetHadronCorrectionv1.h"
4a01bb2c 53#include "AliJetGrid.h"
54
55ClassImp(AliJetFillUnitArrayTracks)
56
4a01bb2c 57//_____________________________________________________________________________
58AliJetFillUnitArrayTracks::AliJetFillUnitArrayTracks()
ee7de0dd 59 : TTask("AliJetFillUnitArrayTracks","Fill Unit Array with tpc/its and emcal information"),
60 fNumUnits(0),
61 fEtaMinCal(0),
62 fEtaMaxCal(0),
63 fPhiMinCal(0),
64 fPhiMaxCal(0),
65 fHadCorr(0),
66 fHCorrection(0),
67 fNTracks(0),
68 fNTracksCut(0),
69 fOpt(0),
70 fDZ(0),
71 fDebug(0),
72 fReaderHeader(0x0),
73 fMomentumArray(0x0),
74 fUnitArray(0x0),
75 fRefArray(0x0),
8838ab7a 76 fProcId(kFALSE),
ee7de0dd 77 fTPCGrid(0x0),
78 fEMCalGrid(0x0),
79 fGeom(0x0),
80 fESD(0x0),
81 fGrid0(0x0),
82 fGrid1(0x0),
83 fGrid2(0x0),
84 fGrid3(0x0),
85 fGrid4(0x0),
86 fNphi(0),
87 fNeta(0),
88 fPhi2(0x0),
89 fEta2(0x0),
90 fPhi(0x0),
91 fEta(0x0),
92 fIndex(0x0),
93 fParams(0x0),
94 fGrid(0),
95 fPhiMin(0),
96 fPhiMax(0),
97 fEtaMin(0),
98 fEtaMax(0),
99 fEtaBinInTPCAcc(0),
100 fPhiBinInTPCAcc(0),
101 fEtaBinInEMCalAcc(0),
102 fPhiBinInEMCalAcc(0),
103 fNbinPhi(0)
4a01bb2c 104{
105 // constructor
4a01bb2c 106}
107
ee7de0dd 108//_____________________________________________________________________________
8838ab7a 109AliJetFillUnitArrayTracks::AliJetFillUnitArrayTracks(AliESDEvent* esd)
ee7de0dd 110 : TTask("AliJetFillUnitArrayTracks","Fill Unit Array with tpc/its and emcal information"),
111 fNumUnits(0),
112 fEtaMinCal(0),
113 fEtaMaxCal(0),
114 fPhiMinCal(0),
115 fPhiMaxCal(0),
116 fHadCorr(0),
117 fHCorrection(0),
118 fNTracks(0),
119 fNTracksCut(0),
120 fOpt(0),
121 fDZ(0),
122 fDebug(0),
123 fReaderHeader(0x0),
124 fMomentumArray(0x0),
125 fUnitArray(0x0),
126 fRefArray(0x0),
8838ab7a 127 fProcId(kFALSE),
ee7de0dd 128 fTPCGrid(0x0),
129 fEMCalGrid(0x0),
130 fGeom(0x0),
8838ab7a 131 fESD(esd),
ee7de0dd 132 fGrid0(0x0),
133 fGrid1(0x0),
134 fGrid2(0x0),
135 fGrid3(0x0),
8838ab7a 136 fGrid4(0x0),// first: loop over tracks in ESD
ee7de0dd 137 fNphi(0),
138 fNeta(0),
139 fPhi2(0x0),
140 fEta2(0x0),
141 fPhi(0x0),
142 fEta(0x0),
143 fIndex(0x0),
144 fParams(0x0),
145 fGrid(0),
146 fPhiMin(0),
147 fPhiMax(0),
148 fEtaMin(0),
149 fEtaMax(0),
150 fEtaBinInTPCAcc(0),
151 fPhiBinInTPCAcc(0),
152 fEtaBinInEMCalAcc(0),
153 fPhiBinInEMCalAcc(0),
154 fNbinPhi(0)
4a01bb2c 155{
8838ab7a 156 // constructor 2
157}
158
159//_____________________________________________________________________________
160AliJetFillUnitArrayTracks::AliJetFillUnitArrayTracks(const AliJetFillUnitArrayTracks &det)
161 : TTask(det),//"AliJetFillUnitArrayTracks","Fill Unit Array with tpc/its and emcal information"),
162 fNumUnits(det.fNumUnits),
163 fEtaMinCal(det.fEtaMinCal),
164 fEtaMaxCal(det.fEtaMaxCal),
165 fPhiMinCal(det.fPhiMinCal),
166 fPhiMaxCal(det.fPhiMaxCal),
167 fHadCorr(det.fHadCorr),
168 fHCorrection(det.fHCorrection),
169 fNTracks(det.fNTracks),
170 fNTracksCut(det.fNTracksCut),
171 fOpt(det.fOpt),
172 fDZ(det.fDZ),
173 fDebug(det.fDebug),
174 fReaderHeader(det.fReaderHeader),
175 fMomentumArray(det.fMomentumArray),
176 fUnitArray(det.fUnitArray),
177 fRefArray(det.fRefArray),
178 fProcId(det.fProcId),
179 fTPCGrid(det.fTPCGrid),
180 fEMCalGrid(det.fEMCalGrid),
181 fGeom(det.fGeom),
182 fESD(det.fESD),
183 fGrid0(det.fGrid0),
184 fGrid1(det.fGrid1),
185 fGrid2(det.fGrid2),
186 fGrid3(det.fGrid3),
187 fGrid4(det.fGrid4),
188 fNphi(det.fNphi),
189 fNeta(det.fNeta),
190 fPhi2(det.fPhi2),
191 fEta2(det.fEta2),
192 fPhi(det.fPhi),
193 fEta(det.fEta),
194 fIndex(det.fIndex),
195 fParams(det.fParams),
196 fGrid(det.fGrid),
197 fPhiMin(det.fPhiMin),
198 fPhiMax(det.fPhiMax),
199 fEtaMin(det.fEtaMin),
200 fEtaMax(det.fEtaMax),
201 fEtaBinInTPCAcc(det.fEtaBinInTPCAcc),
202 fPhiBinInTPCAcc(det.fPhiBinInTPCAcc),
203 fEtaBinInEMCalAcc(det.fEtaBinInEMCalAcc),
204 fPhiBinInEMCalAcc(det.fPhiBinInEMCalAcc),
205 fNbinPhi(det.fNbinPhi)
206{
207 // Copy constructor
208}
209
210//------------------------------------------------------------------
211AliJetFillUnitArrayTracks& AliJetFillUnitArrayTracks::operator=(const AliJetFillUnitArrayTracks& other)
212{
213 // Assignment
214
215 fNumUnits = other.fNumUnits;
216 fEtaMinCal = other.fEtaMinCal;
217 fEtaMaxCal = other.fEtaMaxCal;
218 fPhiMinCal = other.fPhiMinCal;
219 fPhiMaxCal = other.fPhiMaxCal;
220 fHadCorr = other.fHadCorr;
221 fHCorrection = other.fHCorrection;
222 fNTracks = other.fNTracks;
223 fNTracksCut = other.fNTracksCut;
224 fOpt = other.fOpt;
225 fDZ = other.fDZ;
226 fDebug = other.fDebug;
227 fReaderHeader = other.fReaderHeader;
228 fMomentumArray = other.fMomentumArray;
229 fUnitArray = other.fUnitArray;
230 fRefArray = other.fRefArray;
231 fProcId = other.fProcId;
232 fTPCGrid = other.fTPCGrid;
233 fEMCalGrid = other.fEMCalGrid;
234 fGeom = other.fGeom;
235 fESD = other.fESD;
236 fGrid0 = other.fGrid0;
237 fGrid1 = other.fGrid1;
238 fGrid2 = other.fGrid2;
239 fGrid3 = other.fGrid3;
240 fGrid4 = other.fGrid4;
241 fNphi = other.fNphi;
242 fNeta = other.fNeta;
243 fPhi2 = other.fPhi2;
244 fEta2 = other.fEta2;
245 fPhi = other.fPhi;
246 fEta = other.fEta;
247 fIndex = other.fIndex;
248 fParams = other.fParams;
249 fGrid = other.fGrid;
250 fPhiMin = other.fPhiMin;
251 fPhiMax = other.fPhiMax;
252 fEtaMin = other.fEtaMin;
253 fEtaMax = other.fEtaMax;
254 fEtaBinInTPCAcc = other.fEtaBinInTPCAcc;
255 fPhiBinInTPCAcc = other.fPhiBinInTPCAcc;
256 fEtaBinInEMCalAcc = other.fEtaBinInEMCalAcc;
257 fPhiBinInEMCalAcc = other.fPhiBinInEMCalAcc;
258 fNbinPhi = other.fNbinPhi;
259
260 return (*this);
4a01bb2c 261}
262
263//____________________________________________________________________________
264void AliJetFillUnitArrayTracks::InitParameters()
265{
8838ab7a 266 fHadCorr = 0; // For hadron correction
267 fNumUnits = fGeom->GetNCells(); // Number of towers in EMCAL
268 //fDebug = fReaderHeader->GetDebug();
4a01bb2c 269
270 fEtaMinCal = fGeom->GetArm1EtaMin();
271 fEtaMaxCal = fGeom->GetArm1EtaMax();
272 fPhiMinCal = fGeom->GetArm1PhiMin();
8838ab7a 273 fPhiMaxCal = fGeom->GetArm1PhiMax()-10.;
4a01bb2c 274
4a01bb2c 275 fTPCGrid->GetAccParam(fNphi,fNeta,fPhiMin,
276 fPhiMax,fEtaMin,fEtaMax);
277 fTPCGrid->GetBinParam(fPhiBinInTPCAcc,fEtaBinInTPCAcc,
278 fPhiBinInEMCalAcc,fEtaBinInEMCalAcc,fNbinPhi);
279
280 fEta = fTPCGrid->GetArrayEta();
281 fPhi = fTPCGrid->GetArrayPhi();
282 fIndex = fTPCGrid->GetIndexObject();
283
284 if(fDebug>20){
285 for(Int_t i=0; i<fNphi+1; i++) cout << "phi[" << i << "] : " << (*fPhi)[i] << endl;
286 for(Int_t i=0; i<fNeta+1; i++) cout << "eta[" << i << "] : " << (*fEta)[i] << endl;
287
288 for(Int_t i=0; i<fNphi+1; i++)
289 for(Int_t j=0; j<fNeta+1; j++) {cout << "fIndex[" << i << "," << j << "] : " <<
290 (*fIndex)(i,j) << endl; }
291 }
292 if(fDebug>1) printf("\n Parameters initiated ! \n");
293}
294
295//_____________________________________________________________________________
296AliJetFillUnitArrayTracks::~AliJetFillUnitArrayTracks()
297{
298 // destructor
299}
300
301//_____________________________________________________________________________
ee7de0dd 302void AliJetFillUnitArrayTracks::Exec(Option_t* /*option*/)
4a01bb2c 303{
304 //
305 // Main method.
ee7de0dd 306 //
4a01bb2c 307
308 fDebug = fReaderHeader->GetDebug();
8838ab7a 309
4a01bb2c 310 // Set parameters
311 InitParameters();
312
ee7de0dd 313 // get number of tracks in event (for the loop)
314 Int_t goodTrack = 0;
315 Int_t nt = 0;
8838ab7a 316 Int_t nmax = 0;
317 Float_t pt,eta,phi;
318 // Int_t sflag = 0;
ee7de0dd 319 TVector3 p3;
8838ab7a 320
321 if(fDebug>1) cout << "fESD in Fill array : " << fESD << endl;
322
ee7de0dd 323 nt = fESD->GetNumberOfTracks();
324 if(fDebug>1) cout << "Number of Tracks in ESD : " << nt << endl;
8838ab7a 325
326 // temporary storage of signal and pt cut flag
ee7de0dd 327 Int_t* sflag = new Int_t[nt];
328 Int_t* cflag = new Int_t[nt];
329
330 // get cuts set by user
331 Float_t ptMin = fReaderHeader->GetPtCut();
332 Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
333 Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
8838ab7a 334 Float_t phiMin = fReaderHeader->GetFiducialPhiMin();
335 Float_t phiMax = fReaderHeader->GetFiducialPhiMax();
ee7de0dd 336 fOpt = fReaderHeader->GetDetector();
337 fDZ = fReaderHeader->GetDZ();
4a01bb2c 338
ee7de0dd 339 Int_t nTracksEmcal = 0;
340 Int_t nTracksEmcalDZ = 0;
341 Int_t nTracksTpc = 0;
342 Int_t nTracksTpcOnly = 0;
343 Int_t nTracksEmcalCut = 0;
344 Int_t nTracksEmcalDZCut = 0;
345 Int_t nTracksTpcCut = 0;
346 Int_t nTracksTpcOnlyCut = 0;
4a01bb2c 347
ee7de0dd 348 fGrid = fTPCGrid->GetGridType();
4a01bb2c 349
8838ab7a 350 // Loop over tracks
351 nmax = nt;
352 for (Int_t it = 0; it < nmax; it++) {
353 AliESDtrack *track;
354 track = fESD->GetTrack(it);
ee7de0dd 355 UInt_t status = track->GetStatus();
356
357 Double_t mom[3];
358 track->GetPxPyPz(mom);
8838ab7a 359
ee7de0dd 360 p3.SetXYZ(mom[0],mom[1],mom[2]);
8838ab7a 361 vector<Float_t> vtmp(3);
362 vtmp[0] = mom[0]; vtmp[1] = mom[1]; vtmp[2] = mom[2];
ee7de0dd 363 pt = p3.Pt();
8838ab7a 364
ee7de0dd 365 Float_t mass = 0.;
366 mass = track->GetMass();
8838ab7a 367
368 if ((status & AliESDtrack::kTPCrefit) == 0) continue; // quality check
369 if ((status & AliESDtrack::kITSrefit) == 0) continue; // quality check
ee7de0dd 370 if (((AliJetESDReaderHeader*) fReaderHeader)->ReadSignalOnly()
371 && TMath::Abs(track->GetLabel()) > 10000) continue; // quality check
372 if (((AliJetESDReaderHeader*) fReaderHeader)->ReadBkgdOnly()
373 && TMath::Abs(track->GetLabel()) < 10000) continue; // quality check
374 eta = p3.Eta();
375 phi = ( (p3.Phi()) < 0) ? (p3.Phi()) + 2. * TMath::Pi() : (p3.Phi());
8838ab7a 376
ee7de0dd 377 if ( (eta > etaMax) || (eta < etaMin)) continue; // checking eta cut
8838ab7a 378 if ( (phi > phiMax) || (phi < phiMin)) continue; // checking phi cut
379
380 sflag[goodTrack]=0;
381 if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
382 cflag[goodTrack]=0;
383 if (pt > ptMin) cflag[goodTrack]=1; // pt cut
384
385 if(fDebug>20) cout << "In Pythia: Track:" << it << ", eta: " << eta << ", phi: " << phi << ", pt: " << pt << endl;
ee7de0dd 386
ee7de0dd 387 if(fGrid==0)
388 {
389 // Only TPC filled from its grid in its total acceptance
390
391 Int_t idTPC = fTPCGrid->GetIndex(phi,eta);
392 Bool_t ok = kFALSE;
8838ab7a 393 Bool_t ref = kFALSE;
394
ee7de0dd 395 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(idTPC-1);
8838ab7a 396
ee7de0dd 397 uArray->SetUnitTrackID(it);
8838ab7a 398
ee7de0dd 399 Float_t unitEnergy = 0.;
400 unitEnergy = uArray->GetUnitEnergy();
8838ab7a 401 // nTracksTpcOnly is necessary to count the number of candidate cells
402 // but it doesn't give the real multiplicity -> it will be extracted
403 // in the jet finder and for that, it is necessary to fill the Px,Py,Pz
404 // information for each tracks stored in a given unitcell
ee7de0dd 405 if(unitEnergy==0.){
406 nTracksTpcOnly++;
407 ok = kTRUE;
8838ab7a 408 ref = kTRUE;
ee7de0dd 409 }
8838ab7a 410
ee7de0dd 411 // Fill energy in TPC acceptance
412 uArray->SetUnitEnergy(unitEnergy + pt);
8838ab7a 413 uArray->SetUnitPxPyPz(kFALSE,vtmp);
414 uArray->SetUnitMass(mass);
ee7de0dd 415
416 // Pt cut flag
417 if(uArray->GetUnitEnergy()<ptMin){
418 uArray->SetUnitCutFlag(kPtSmaller);
419 }
420 else {
421 uArray->SetUnitCutFlag(kPtHigher);
422 if(ok) nTracksTpcOnlyCut++;
4a01bb2c 423 }
424
8838ab7a 425 // Signal flag
426 if(sflag[goodTrack] == 1) {
427 uArray->SetUnitSignalFlag(kGood);
428 uArray->SetUnitSignalFlagC(kFALSE,kGood);
429 } else uArray->SetUnitSignalFlagC(kFALSE,kBad);
430
431 if(uArray->GetUnitEnergy()>0 && ref){
432 if(!fProcId) {
433 fProcId = kTRUE;
434 // delete fRefArray;
435 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
436 }
ee7de0dd 437 fRefArray->Add(uArray);
438 }
8838ab7a 439
ee7de0dd 440 }
441
442 if(fGrid==1)
443 {
444 Int_t nElements = fTPCGrid->GetNEntries();
ee7de0dd 445 // Fill track information in EMCAL acceptance
446 if((eta >= fEtaMin && eta <= fEtaMax) &&
447 (phi >= fPhiMin && phi <= fPhiMax))// &&
448 {
ee7de0dd 449 // Include dead-zones
450 if(fDZ)
451 {
452 Double_t phimin0 = 0., phimin1 = 0., phimin2 = 0., phimin3 = 0., phimin4 = 0.;
453 Double_t phimax0 = 0., phimax1 = 0., phimax2 = 0., phimax3 = 0., phimax4 = 0.;
454 fGeom->GetPhiBoundariesOfSMGap(0,phimin0,phimax0);
455 fGeom->GetPhiBoundariesOfSMGap(1,phimin1,phimax1);
456 fGeom->GetPhiBoundariesOfSMGap(2,phimin2,phimax2);
457 fGeom->GetPhiBoundariesOfSMGap(3,phimin3,phimax3);
458 fGeom->GetPhiBoundariesOfSMGap(4,phimin4,phimax4);
459 Int_t n0 = fGrid0->GetNEntries();
460 Int_t n1 = fGrid1->GetNEntries();
461 Int_t n2 = fGrid2->GetNEntries();
462 Int_t n3 = fGrid3->GetNEntries();
8838ab7a 463
ee7de0dd 464 if(phi >= phimin0 && phi <= phimax0){
465 Int_t id0 = fGrid0->GetIndex(phi,eta)-1;
466 AliJetUnitArray *uArray0 = (AliJetUnitArray*)fUnitArray->At(id0+fNumUnits+nElements);
8838ab7a 467
ee7de0dd 468 uArray0->SetUnitTrackID(it);
8838ab7a 469
ee7de0dd 470 Float_t uEnergy0 = uArray0->GetUnitEnergy();
471 Bool_t ok0 = kFALSE;
8838ab7a 472 Bool_t ref0 = kFALSE;
ee7de0dd 473 if(uEnergy0==0.){
474 nTracksEmcalDZ++;
475 ok0 = kTRUE;
8838ab7a 476 ref0 = kTRUE;
ee7de0dd 477 }
478 uArray0->SetUnitEnergy(uEnergy0+pt);
8838ab7a 479 uArray0->SetUnitPxPyPz(kFALSE,vtmp);
ee7de0dd 480 if(uArray0->GetUnitEnergy()<ptMin)
481 uArray0->SetUnitCutFlag(kPtSmaller);
482 else {
483 uArray0->SetUnitCutFlag(kPtHigher);
484 if(ok0) nTracksEmcalDZCut++;
485 }
8838ab7a 486 if(sflag[goodTrack] == 1) {
487 uArray0->SetUnitSignalFlag(kGood);
488 uArray0->SetUnitSignalFlagC(kFALSE,kGood);
489 } else uArray0->SetUnitSignalFlagC(kFALSE,kBad);
490
491 if(uArray0->GetUnitEnergy()>0 && ref0){
492 if(!fProcId){
493 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray0));
494 fProcId = kTRUE;
495 }
ee7de0dd 496 fRefArray->Add(uArray0);
8838ab7a 497 }
4a01bb2c 498 }
ee7de0dd 499 if(phi >= phimin1 && phi <= phimax1){
500 Int_t id1 = fGrid1->GetIndex(phi,eta)-1+n0;
501 AliJetUnitArray *uArray1 = (AliJetUnitArray*)fUnitArray->At(id1+fNumUnits+nElements);
502 uArray1->SetUnitTrackID(it);
8838ab7a 503
ee7de0dd 504 Float_t uEnergy1 = uArray1->GetUnitEnergy();
505 Bool_t ok1 = kFALSE;
8838ab7a 506 Bool_t ref1 = kFALSE;
ee7de0dd 507 if(uEnergy1==0.){
508 nTracksEmcalDZ++;
509 ok1 = kTRUE;
8838ab7a 510 ref1 = kTRUE;
ee7de0dd 511 }
512 uArray1->SetUnitEnergy(uEnergy1+pt);
8838ab7a 513 uArray1->SetUnitPxPyPz(kFALSE,vtmp);
ee7de0dd 514 if(uArray1->GetUnitEnergy()<ptMin)
515 uArray1->SetUnitCutFlag(kPtSmaller);
4a01bb2c 516 else {
ee7de0dd 517 uArray1->SetUnitCutFlag(kPtHigher);
518 if(ok1) nTracksEmcalDZCut++;
4a01bb2c 519 }
8838ab7a 520 if(sflag[goodTrack] == 1) {
521 uArray1->SetUnitSignalFlag(kGood);
522 uArray1->SetUnitSignalFlagC(kFALSE,kGood);
523 } else uArray1->SetUnitSignalFlagC(kFALSE,kBad);
524
525 if(uArray1->GetUnitEnergy()>0 && ref1){
526 if(!fProcId){
527 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray1));
528 fProcId = kTRUE;
529 }
ee7de0dd 530 fRefArray->Add(uArray1);
8838ab7a 531 }
ee7de0dd 532 }
533 if(phi >= phimin2 && phi <= phimax2){
534 Int_t id2 = fGrid2->GetIndex(phi,eta)-1+n0+n1;
535 AliJetUnitArray *uArray2 = (AliJetUnitArray*)fUnitArray->At(id2+fNumUnits+nElements);
536 uArray2->SetUnitTrackID(it);
8838ab7a 537
ee7de0dd 538 Float_t uEnergy2 = uArray2->GetUnitEnergy();
539 Bool_t ok2 = kFALSE;
8838ab7a 540 Bool_t ref2 = kFALSE;
ee7de0dd 541 if(uEnergy2==0.){
542 nTracksEmcalDZ++;
543 ok2 = kTRUE;
8838ab7a 544 ref2 = kTRUE;
ee7de0dd 545 }
546 uArray2->SetUnitEnergy(uEnergy2+pt);
8838ab7a 547 uArray2->SetUnitPxPyPz(kFALSE,vtmp);
ee7de0dd 548 if(uArray2->GetUnitEnergy()<ptMin)
549 uArray2->SetUnitCutFlag(kPtSmaller);
550 else {
551 uArray2->SetUnitCutFlag(kPtHigher);
552 if(ok2) nTracksEmcalDZCut++;
4a01bb2c 553 }
8838ab7a 554 if(sflag[goodTrack] == 1) {
555 uArray2->SetUnitSignalFlag(kGood);
556 uArray2->SetUnitSignalFlagC(kFALSE,kGood);
557 } else uArray2->SetUnitSignalFlagC(kFALSE,kBad);
558
559 if(uArray2->GetUnitEnergy()>0 && ref2){
560 if(!fProcId){
561 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray2));
562 fProcId = kTRUE;
563 }
ee7de0dd 564 fRefArray->Add(uArray2);
8838ab7a 565 }
4a01bb2c 566 }
ee7de0dd 567 if(phi >= phimin3 && phi <= phimax3){
568 Int_t id3 = fGrid3->GetIndex(phi,eta)-1+n0+n1+n2;
569 AliJetUnitArray *uArray3 = (AliJetUnitArray*)fUnitArray->At(id3+fNumUnits+nElements);
570 uArray3->SetUnitTrackID(it);
8838ab7a 571
ee7de0dd 572 Float_t uEnergy3 = uArray3->GetUnitEnergy();
573 Bool_t ok3 = kFALSE;
8838ab7a 574 Bool_t ref3 = kFALSE;
ee7de0dd 575 if(uEnergy3==0.){
576 nTracksEmcalDZ++;
577 ok3 = kTRUE;
8838ab7a 578 ref3 = kTRUE;
ee7de0dd 579 }
580 uArray3->SetUnitEnergy(uEnergy3+pt);
8838ab7a 581 uArray3->SetUnitPxPyPz(kFALSE,vtmp);
ee7de0dd 582 if(uArray3->GetUnitEnergy()<ptMin)
583 uArray3->SetUnitCutFlag(kPtSmaller);
584 else {
585 uArray3->SetUnitCutFlag(kPtHigher);
586 if(ok3) nTracksEmcalDZCut++;
587 }
8838ab7a 588 if(sflag[goodTrack] == 1) {
589 uArray3->SetUnitSignalFlag(kGood);
590 uArray3->SetUnitSignalFlagC(kFALSE,kGood);
591 } else uArray3->SetUnitSignalFlagC(kFALSE,kBad);
592
593 if(uArray3->GetUnitEnergy()>0 && ref3){
594 if(!fProcId){
595 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray3));
596 fProcId = kTRUE;
597 }
ee7de0dd 598 fRefArray->Add(uArray3);
8838ab7a 599 }
ee7de0dd 600 }
601 if(phi >= phimin4 && phi <= phimax4){
602 Int_t id4 = fGrid4->GetIndex(phi,eta)-1+n0+n1+n2+n3;
603 AliJetUnitArray *uArray4 = (AliJetUnitArray*)fUnitArray->At(id4+fNumUnits+nElements);
604 uArray4->SetUnitTrackID(it);
8838ab7a 605
ee7de0dd 606 Float_t uEnergy4 = uArray4->GetUnitEnergy();
607 Bool_t ok4 = kFALSE;
8838ab7a 608 Bool_t ref4 = kFALSE;
ee7de0dd 609 if(uEnergy4==0.){
610 nTracksEmcalDZ++;
611 ok4 = kTRUE;
8838ab7a 612 ref4 = kTRUE;
ee7de0dd 613 }
614 uArray4->SetUnitEnergy(uEnergy4+pt);
8838ab7a 615 uArray4->SetUnitPxPyPz(kFALSE,vtmp);
ee7de0dd 616 if(uArray4->GetUnitEnergy()<ptMin)
617 uArray4->SetUnitCutFlag(kPtSmaller);
618 else {
619 uArray4->SetUnitCutFlag(kPtHigher);
620 if(ok4) nTracksEmcalDZCut++;
621 }
8838ab7a 622 if(sflag[goodTrack] == 1) {
623 uArray4->SetUnitSignalFlag(kGood);
624 uArray4->SetUnitSignalFlagC(kFALSE,kGood);
625 } else uArray4->SetUnitSignalFlagC(kFALSE,kBad);
626
627 if(uArray4->GetUnitEnergy()>0 && ref4){
628 if(!fProcId){
629 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray4));
630 fProcId = kTRUE;
631 }
ee7de0dd 632 fRefArray->Add(uArray4);
8838ab7a 633 }
ee7de0dd 634 }
635 } // end fDZ
636
637 Int_t towerID = 0;
638 fGeom->GetAbsCellIdFromEtaPhi(eta,phi,towerID);
ee7de0dd 639 if(towerID==-1) continue;
8838ab7a 640
ee7de0dd 641 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(towerID);
642 uArray->SetUnitTrackID(it);
8838ab7a 643
ee7de0dd 644 Float_t unitEnergy = uArray->GetUnitEnergy();
645 Bool_t ok = kFALSE;
8838ab7a 646 Bool_t ref = kFALSE;
ee7de0dd 647 if(unitEnergy==0.){
648 nTracksEmcal++;
649 ok=kTRUE;
8838ab7a 650 ref=kTRUE;
4a01bb2c 651 }
652
ee7de0dd 653 // Do Hadron Correction
654 // Parametrization to be added
655 if (fHCorrection != 0)
656 {
ee7de0dd 657 Float_t hCEnergy = fHadCorr->GetEnergy(p3.Mag(), (Double_t)eta);
658 unitEnergy -= hCEnergy*TMath::Sin(2.0*TMath::ATan(TMath::Exp(-eta)));
659
660 } //end Hadron Correction loop
8838ab7a 661
ee7de0dd 662 uArray->SetUnitEnergy(unitEnergy + pt);
8838ab7a 663 uArray->SetUnitPxPyPz(kFALSE,vtmp);
ee7de0dd 664 // Put a pt cut flag
8838ab7a 665 if(uArray->GetUnitEnergy()<ptMin){
ee7de0dd 666 uArray->SetUnitCutFlag(kPtSmaller);
8838ab7a 667 }
ee7de0dd 668 else {
669 uArray->SetUnitCutFlag(kPtHigher);
670 if(ok) nTracksEmcalCut++;
671 }
ee7de0dd 672
8838ab7a 673 // Signal flag
674 if(sflag[goodTrack] == 1) {
675 uArray->SetUnitSignalFlag(kGood);
676 uArray->SetUnitSignalFlagC(kFALSE,kGood);
677 } else uArray->SetUnitSignalFlagC(kFALSE,kBad);
ee7de0dd 678
8838ab7a 679 if(uArray->GetUnitEnergy()>0 && ref){
680 if(!fProcId){
681 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
682 fProcId = kTRUE;
683 }
684 fRefArray->Add(uArray);
685 }
ee7de0dd 686 } // end loop on EMCal acceptance cut + tracks quality
687 else{
688 // Outside EMCal acceptance
689
ee7de0dd 690 Int_t idTPC = fTPCGrid->GetIndex(phi,eta);
8838ab7a 691
ee7de0dd 692 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(fNumUnits-1+idTPC);
693 uArray->SetUnitTrackID(it);
694
695 Float_t unitEnergy2 = uArray->GetUnitEnergy(); // check if fNumUnits or fNumUnits-1
8838ab7a 696 Bool_t okout = kFALSE;
697 Bool_t refout = kFALSE;
ee7de0dd 698 if(unitEnergy2==0.){
699 nTracksTpc++;
8838ab7a 700 okout=kTRUE;
701 refout=kTRUE;
ee7de0dd 702 }
703 // Fill energy outside emcal acceptance
704 uArray->SetUnitEnergy(unitEnergy2 + pt);
8838ab7a 705 uArray->SetUnitPxPyPz(kFALSE,vtmp);
ee7de0dd 706
707 // Pt cut flag
708 if(uArray->GetUnitEnergy()<ptMin){
709 uArray->SetUnitCutFlag(kPtSmaller);
710 }
711 else {
712 uArray->SetUnitCutFlag(kPtHigher);
8838ab7a 713 if(okout) nTracksTpcCut++;
ee7de0dd 714 }
8838ab7a 715 // Signal flag
716 if(sflag[goodTrack] == 1) {
717 uArray->SetUnitSignalFlag(kGood);
718 uArray->SetUnitSignalFlagC(kFALSE,kGood);
719 } else uArray->SetUnitSignalFlagC(kFALSE,kBad);
ee7de0dd 720
8838ab7a 721 if(uArray->GetUnitEnergy()>0 && refout){
722 if(!fProcId){
723 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
724 fProcId = kTRUE;
725 }
726 fRefArray->Add(uArray);
727 }
ee7de0dd 728 }
729 } // end fGrid==1
8838ab7a 730
731 goodTrack++;
732
ee7de0dd 733 } // end loop on entries (tpc tracks)
4a01bb2c 734
8838ab7a 735 if(fDebug>0)
736 {
737 cout << "End of Tracks %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << endl;
738 cout << "goodTracks: " << goodTrack << endl;
739 }
ee7de0dd 740
8838ab7a 741 delete[] sflag;
742 delete[] cflag;
ee7de0dd 743
ee7de0dd 744 if(fGrid==0) {
745 fNTracks = nTracksTpcOnly;
746 fNTracksCut = nTracksTpcOnlyCut;
747 if(fDebug>10){
748 cout << "fNTracks : " << fNTracks << endl;
749 cout << "fNTracksCut : " << fNTracksCut << endl;
750 }
751 }
752 if(fGrid==1) {
753 fNTracks = nTracksEmcal+nTracksEmcalDZ+nTracksTpc;
754 fNTracksCut = nTracksEmcalCut+nTracksEmcalDZCut+nTracksTpcCut;
755 if(fDebug>10){
756 cout << "fNTracks : " << fNTracks << endl;
757 cout << "fNTracksCut : " << fNTracksCut << endl;
758 }
759 }
8838ab7a 760
4a01bb2c 761}
762
763//__________________________________________________________
764void AliJetFillUnitArrayTracks::GetEtaPhiFromIndex(Int_t index, Float_t &eta, Float_t &phi)
765{
8838ab7a 766 for(Int_t j=0; j<fNphi+1; j++)
767 {
768 for(Int_t i=0; i<fNeta+1; i++)
769 {
770
771 // TPC grid only
772 //-------------------------------------
773 if(fGrid==0) {
774 if(j*(fNeta+1)+i == index) {
775 eta = fEta2->At(i);
776 phi = fPhi2->At(j);
777 }
778 }
779
780 // TPC-EMCAL grid
781 //-------------------------------------
782 Int_t ii = 0;
783 if(i==0) ii = 0;
784 if(i>0 && i<(fEtaBinInTPCAcc-fEtaBinInEMCalAcc)/2) ii = i;
785 if(i>=(fEtaBinInTPCAcc+fEtaBinInEMCalAcc)/2 && i<fNeta+1) ii = i-fEtaBinInEMCalAcc;
786
787 if(fGrid==1) {
788 if(j<(fNbinPhi+1) && j*(fNeta+1)+i == index) {
789 eta = fEta2->At(i);
790 phi = fPhi2->At(j);
791 }
792
793 if((j>=(fNbinPhi+1) && j<(fNbinPhi+1+fPhiBinInEMCalAcc)) &&
794 ((fNbinPhi+1)*(fNeta+1) + (j-fNbinPhi-1)*(fEtaBinInTPCAcc-fEtaBinInEMCalAcc) + ii)== index ) {
795 if(ii==0) {Int_t ind = 0; eta = fEta2->At(ind);}
796 else eta = fEta2->At(i);
797 phi = fPhi2->At(j);
798 }
799
800 if(j>=(fNbinPhi+1+fPhiBinInEMCalAcc) && ((fNbinPhi+1)*(fNeta+1)+fPhiBinInEMCalAcc*((fEtaBinInTPCAcc-fEtaBinInEMCalAcc))+(j-(fNbinPhi+1+fPhiBinInEMCalAcc))*(fNeta+1)+i == index)) {
801 eta = fEta2->At(i);
802 phi = fPhi2->At(j);
803 }
804 }
4a01bb2c 805 }
4a01bb2c 806 }
4a01bb2c 807}
808
809
810
811
812
813
814
815
816
817
818
819