Adding rho's depenence on trigger track (M. Verweij)
[u/mrichter/AliRoot.git] / JETAN / AliJetAODFillUnitArrayTracks.cxx
CommitLineData
be6e5811 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
16
17//======================================================================
18// ***July 2006
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//======================================================================
24
25
26// --- ROOT system ---
27#include <TVector3.h>
28#include <TProcessID.h>
29
30#include "AliJetUnitArray.h"
e36a3f22 31#include "AliJetHadronCorrectionv1.h"
be6e5811 32#include "AliJetAODFillUnitArrayTracks.h"
33
34// --- ROOT system ---
35class TSystem;
36class TLorentzVector;
37class TGeoManager;
38
39// --- AliRoot header files ---
40class AliJetFinder;
41
3a7af7bd 42using std::cout;
43using std::endl;
be6e5811 44ClassImp(AliJetAODFillUnitArrayTracks)
45
46//_____________________________________________________________________________
47AliJetAODFillUnitArrayTracks::AliJetAODFillUnitArrayTracks()
48 : AliJetFillUnitArray(),
49 fNumUnits(0),
50 fHadCorr(0),
51 fApplyMIPCorrection(kTRUE),
52 fAOD(0x0),
53 fGrid0(0x0),
54 fGrid1(0x0),
55 fGrid2(0x0),
56 fGrid3(0x0),
57 fGrid4(0x0)
58{
59 // constructor
60}
61
62//_____________________________________________________________________________
63AliJetAODFillUnitArrayTracks::AliJetAODFillUnitArrayTracks(AliAODEvent* aod)
64 : AliJetFillUnitArray(),
65 fNumUnits(0),
66 fHadCorr(0),
67 fApplyMIPCorrection(kTRUE),
68 fAOD(aod),
69 fGrid0(0x0),
70 fGrid1(0x0),
71 fGrid2(0x0),
72 fGrid3(0x0),
73 fGrid4(0x0)
74{
75 // constructor
76}
77
78//_____________________________________________________________________________
79AliJetAODFillUnitArrayTracks::AliJetAODFillUnitArrayTracks(const AliJetAODFillUnitArrayTracks &det)
80 : AliJetFillUnitArray(det),
81 fNumUnits(det.fNumUnits),
82 fHadCorr(det.fHadCorr),
83 fApplyMIPCorrection(det.fApplyMIPCorrection),
84 fAOD(det.fAOD),
85 fGrid0(det.fGrid0),
86 fGrid1(det.fGrid1),
87 fGrid2(det.fGrid2),
88 fGrid3(det.fGrid3),
89 fGrid4(det.fGrid4)
90{
91 // Copy constructor
92}
93
94//_____________________________________________________________________________
95AliJetAODFillUnitArrayTracks& AliJetAODFillUnitArrayTracks::operator=(const AliJetAODFillUnitArrayTracks& other)
96{
97 // Assignment
98
99 fNumUnits = other.fNumUnits;
100 fHadCorr = other.fHadCorr;
101 fApplyMIPCorrection = other.fApplyMIPCorrection;
102 fAOD = other.fAOD;
103 fGrid0 = other.fGrid0;
104 fGrid1 = other.fGrid1;
105 fGrid2 = other.fGrid2;
106 fGrid3 = other.fGrid3;
107 fGrid4 = other.fGrid4;
108
109 return (*this);
110}
111
112//____________________________________________________________________________
113void AliJetAODFillUnitArrayTracks::InitParameters()
114{
e36a3f22 115 // fHadCorr = 0; // For hadron correction
f45e6f72 116 fNumUnits = fGeom->GetEMCGeometry()->GetNCells(); // Number of towers in EMCAL
be6e5811 117
118 fTPCGrid->GetAccParam(fNphi,fNeta,fPhiMin,
119 fPhiMax,fEtaMin,fEtaMax);
120 fTPCGrid->GetBinParam(fPhiBinInTPCAcc,fEtaBinInTPCAcc,
121 fPhiBinInEMCalAcc,fEtaBinInEMCalAcc,fNbinPhi);
122
123 fIndex = fTPCGrid->GetIndexObject();
124
125 if(fDebug>20){
126 for(Int_t i=0; i<fNphi+1; i++)
127 for(Int_t j=0; j<fNeta+1; j++) {cout << "fIndex[" << i << "," << j << "] : " <<
128 (*fIndex)(i,j) << endl; }
129 }
130 if(fDebug>1) printf("\n Parameters initiated ! \n");
131}
132
133//_____________________________________________________________________________
134AliJetAODFillUnitArrayTracks::~AliJetAODFillUnitArrayTracks()
135{
136 // destructor
137}
138
139//_____________________________________________________________________________
140void AliJetAODFillUnitArrayTracks::Exec(Option_t* const /*option*/)
141//void AliJetAODFillUnitArrayTracks::Exec(Option_t* option)
142{
143 //
144 // Main method.
145 // Fill the unit array with the charged particle information in AOD
146 //
147
148 fDebug = fReaderHeader->GetDebug();
149
150 // Set parameters
151 InitParameters();
e36a3f22 152 fRef->Clear();
be6e5811 153
154 // get number of tracks in event (for the loop)
155 Int_t goodTrack = 0;
156 Int_t nt = 0;
157//(not used ?) Int_t nt2 = 0;
158 Int_t nmax = 0;
159 Float_t pt,eta,phi;
160 // Int_t sflag = 0;
161 TVector3 p3;
162
163 nt = fAOD->GetNumberOfTracks();
e36a3f22 164 if(fDebug>1) cout << "Number of Tracks in AOD : " << nt << endl;
165
be6e5811 166 // temporary storage of signal and pt cut flag
167 Int_t* sflag = new Int_t[nt];
168 Int_t* cflag = new Int_t[nt];
169
170 // get cuts set by user
171 Float_t ptMin = fReaderHeader->GetPtCut();
172 Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
173 Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
174 fOpt = fReaderHeader->GetDetector();
175 fDZ = fReaderHeader->GetDZ();
176 UInt_t filterMask = ((AliJetAODReaderHeader*)fReaderHeader)->GetTestFilterMask();
177
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;
186
be6e5811 187 fGrid = fTPCGrid->GetGridType();
188
be6e5811 189 //loop over tracks
190 nmax = nt;
191 for (Int_t it = 0; it < nmax; it++) {
192 AliAODTrack *track;
193 track = fAOD->GetTrack(it);
194 UInt_t status = track->GetStatus();
195
196 Double_t mom[3];
197 track->GetPxPyPz(mom);
198
199 p3.SetXYZ(mom[0],mom[1],mom[2]);
200 pt = p3.Pt();
201
be6e5811 202 eta = p3.Eta();
203 phi = ( (p3.Phi()) < 0) ? (p3.Phi()) + 2. * TMath::Pi() : (p3.Phi());
204 if (status == 0) continue;
205 if((filterMask>0)&&!(track->TestFilterBit(filterMask)))continue;
206
207 if ( (eta > etaMax) || (eta < etaMin)) continue; // checking eta cut
208
e36a3f22 209 if (goodTrack == 0) new(fRef) TRefArray(TProcessID::GetProcessWithUID(track));
210
be6e5811 211 sflag[goodTrack]=0;
212 if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
213 cflag[goodTrack]=0;
214 if (pt > ptMin) cflag[goodTrack]=1; // pt cut
e36a3f22 215 fRef->Add(track);
be6e5811 216
217 if(fGrid==0)
218 {
219 // Only TPC filled from its grid in its total acceptance
220
221 Int_t idTPC = fTPCGrid->GetIndex(phi,eta);
222 Bool_t ok = kFALSE;
223 Bool_t ref = kFALSE;
224
225 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(idTPC-1);
226 TRefArray *reference = uArray->GetUnitTrackRef();
a9c14d46 227 if (reference->GetEntries() == 0) {
228 new(reference) TRefArray(TProcessID::GetProcessWithUID(track));
229 }
230
be6e5811 231 reference->Add(track);
232
233 Float_t unitEnergy = 0.;
234 unitEnergy = uArray->GetUnitEnergy();
235 // nTracksTpcOnly is necessary to count the number of candidate cells
236 // but it doesn't give the real multiplicity -> it will be extracted
237 // in the jet finder through the reference to tracks
238 if(unitEnergy==0.){
239 nTracksTpcOnly++;
240 ok = kTRUE;
241 ref = kTRUE;
242 }
243
244 // Fill energy in TPC acceptance
245 uArray->SetUnitEnergy(unitEnergy + pt);
246
247 // Pt cut flag
248 if(uArray->GetUnitEnergy()<ptMin){
249 uArray->SetUnitCutFlag(kPtSmaller);
250 }
251 else {
252 uArray->SetUnitCutFlag(kPtHigher);
253 if(ok) nTracksTpcOnlyCut++;
254 }
255
256 // Signal flag
257 if(sflag[goodTrack] == 1) {
258 uArray->SetUnitSignalFlag(kGood);
259 uArray->SetUnitSignalFlagC(kFALSE,kGood);
260 } else uArray->SetUnitSignalFlagC(kFALSE,kBad);
261
262 if(uArray->GetUnitEnergy()>0 && ref){
263 if(!fProcId) {
264 fProcId = kTRUE;
265 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
266 }
267 fRefArray->Add(uArray);
268 }
269 }
270
271 if(fGrid==1)
272 {
273
274 Int_t nElements = fTPCGrid->GetNEntries();
275 // Fill track information in EMCAL acceptance
276 if((eta >= fEtaMin && eta <= fEtaMax) &&
277 (phi >= fPhiMin && phi <= fPhiMax))// &&
278 {
279 // Include dead-zones
280 if(fDZ)
281 {
282 Double_t phimin0 = 0., phimin1 = 0., phimin2 = 0., phimin3 = 0., phimin4 = 0.;
283 Double_t phimax0 = 0., phimax1 = 0., phimax2 = 0., phimax3 = 0., phimax4 = 0.;
f45e6f72 284 fGeom->GetEMCGeometry()->GetPhiBoundariesOfSMGap(0,phimin0,phimax0);
285 fGeom->GetEMCGeometry()->GetPhiBoundariesOfSMGap(1,phimin1,phimax1);
286 fGeom->GetEMCGeometry()->GetPhiBoundariesOfSMGap(2,phimin2,phimax2);
287 fGeom->GetEMCGeometry()->GetPhiBoundariesOfSMGap(3,phimin3,phimax3);
288 fGeom->GetEMCGeometry()->GetPhiBoundariesOfSMGap(4,phimin4,phimax4);
be6e5811 289 Int_t n0 = fGrid0->GetNEntries();
290 Int_t n1 = fGrid1->GetNEntries();
291 Int_t n2 = fGrid2->GetNEntries();
292 Int_t n3 = fGrid3->GetNEntries();
be6e5811 293
be6e5811 294 if(phi >= phimin0 && phi <= phimax0){
295 Int_t id0 = fGrid0->GetIndex(phi,eta)-1;
296 AliJetUnitArray *uArray0 = (AliJetUnitArray*)fUnitArray->At(id0+fNumUnits+nElements);
297 TRefArray *reference0 = uArray0->GetUnitTrackRef();
a9c14d46 298 if (reference0->GetEntries() == 0) {
299 new(reference0) TRefArray(TProcessID::GetProcessWithUID(track));
300 }
301
be6e5811 302 reference0->Add(track);
303 uArray0->SetUnitTrackID(it);
304
305 Float_t uEnergy0 = uArray0->GetUnitEnergy();
306 Bool_t ok0 = kFALSE;
307 Bool_t ref0 = kFALSE;
308 if(uEnergy0==0.){
309 nTracksEmcalDZ++;
310 ok0 = kTRUE;
311 ref0 = kTRUE;
312 }
313 uArray0->SetUnitEnergy(uEnergy0+pt);
314
315 if(uArray0->GetUnitEnergy()<ptMin)
316 uArray0->SetUnitCutFlag(kPtSmaller);
317 else {
318 uArray0->SetUnitCutFlag(kPtHigher);
319 if(ok0) nTracksEmcalDZCut++;
320 }
321 if(sflag[goodTrack] == 1) {
322 uArray0->SetUnitSignalFlag(kGood);
323 uArray0->SetUnitSignalFlagC(kFALSE,kGood);
324 } else uArray0->SetUnitSignalFlagC(kFALSE,kBad);
325
326 if(uArray0->GetUnitEnergy()>0 && ref0){
327 if(!fProcId){
328 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray0));
329 fProcId = kTRUE;
330 }
331 fRefArray->Add(uArray0);
332 }
333 }
334
335 if(phi >= phimin1 && phi <= phimax1){
336 Int_t id1 = fGrid1->GetIndex(phi,eta)-1+n0;
337 AliJetUnitArray *uArray1 = (AliJetUnitArray*)fUnitArray->At(id1+fNumUnits+nElements);
338 TRefArray *reference1 = uArray1->GetUnitTrackRef();
a9c14d46 339 if (reference1->GetEntries() == 0) {
340 new(reference1) TRefArray(TProcessID::GetProcessWithUID(track));
341 }
342
be6e5811 343 reference1->Add(track);
344
345 Float_t uEnergy1 = uArray1->GetUnitEnergy();
346 Bool_t ok1 = kFALSE;
347 Bool_t ref1 = kFALSE;
348 if(uEnergy1==0.){
349 nTracksEmcalDZ++;
350 ok1 = kTRUE;
351 ref1 = kTRUE;
352 }
353 uArray1->SetUnitEnergy(uEnergy1+pt);
354
355 if(uArray1->GetUnitEnergy()<ptMin)
356 uArray1->SetUnitCutFlag(kPtSmaller);
357 else {
358 uArray1->SetUnitCutFlag(kPtHigher);
359 if(ok1) nTracksEmcalDZCut++;
360 }
361 if(sflag[goodTrack] == 1) {
362 uArray1->SetUnitSignalFlag(kGood);
363 uArray1->SetUnitSignalFlagC(kFALSE,kGood);
364 } else uArray1->SetUnitSignalFlagC(kFALSE,kBad);
365
366 if(uArray1->GetUnitEnergy()>0 && ref1){
367 if(!fProcId){
368 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray1));
369 fProcId = kTRUE;
370 }
371 fRefArray->Add(uArray1);
372 }
373 }
374
375 if(phi >= phimin2 && phi <= phimax2){
376 Int_t id2 = fGrid2->GetIndex(phi,eta)-1+n0+n1;
377 AliJetUnitArray *uArray2 = (AliJetUnitArray*)fUnitArray->At(id2+fNumUnits+nElements);
378 TRefArray *reference2 = uArray2->GetUnitTrackRef();
a9c14d46 379 if (reference2->GetEntries() == 0) {
380 new(reference2) TRefArray(TProcessID::GetProcessWithUID(track));
381 }
382
be6e5811 383 reference2->Add(track);
384
385 Float_t uEnergy2 = uArray2->GetUnitEnergy();
386 Bool_t ok2 = kFALSE;
387 Bool_t ref2 = kFALSE;
388 if(uEnergy2==0.){
389 nTracksEmcalDZ++;
390 ok2 = kTRUE;
391 ref2 = kTRUE;
392 }
393 uArray2->SetUnitEnergy(uEnergy2+pt);
394
395 if(uArray2->GetUnitEnergy()<ptMin)
396 uArray2->SetUnitCutFlag(kPtSmaller);
397 else {
398 uArray2->SetUnitCutFlag(kPtHigher);
399 if(ok2) nTracksEmcalDZCut++;
400 }
401 if(sflag[goodTrack] == 1) {
402 uArray2->SetUnitSignalFlag(kGood);
403 uArray2->SetUnitSignalFlagC(kFALSE,kGood);
404 } else uArray2->SetUnitSignalFlagC(kFALSE,kBad);
405
406 if(uArray2->GetUnitEnergy()>0 && ref2){
407 if(!fProcId){
408 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray2));
409 fProcId = kTRUE;
410 }
411 fRefArray->Add(uArray2);
412 }
413 }
414
415 if(phi >= phimin3 && phi <= phimax3){
416 Int_t id3 = fGrid3->GetIndex(phi,eta)-1+n0+n1+n2;
417 AliJetUnitArray *uArray3 = (AliJetUnitArray*)fUnitArray->At(id3+fNumUnits+nElements);
418 TRefArray *reference3 = uArray3->GetUnitTrackRef();
a9c14d46 419 if (reference3->GetEntries() == 0) {
420 new(reference3) TRefArray(TProcessID::GetProcessWithUID(track));
421 }
422
be6e5811 423 reference3->Add(track);
424
425 Float_t uEnergy3 = uArray3->GetUnitEnergy();
426 Bool_t ok3 = kFALSE;
427 Bool_t ref3 = kFALSE;
428 if(uEnergy3==0.){
429 nTracksEmcalDZ++;
430 ok3 = kTRUE;
431 ref3 = kTRUE;
432 }
433 uArray3->SetUnitEnergy(uEnergy3+pt);
434
435 if(uArray3->GetUnitEnergy()<ptMin)
436 uArray3->SetUnitCutFlag(kPtSmaller);
437 else {
438 uArray3->SetUnitCutFlag(kPtHigher);
439 if(ok3) nTracksEmcalDZCut++;
440 }
441 if(sflag[goodTrack] == 1) {
442 uArray3->SetUnitSignalFlag(kGood);
443 uArray3->SetUnitSignalFlagC(kFALSE,kGood);
444 } else uArray3->SetUnitSignalFlagC(kFALSE,kBad);
445
446 if(uArray3->GetUnitEnergy()>0 && ref3){
447 if(!fProcId){
448 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray3));
449 fProcId = kTRUE;
450 }
451 fRefArray->Add(uArray3);
452 }
453 }
454
455 if(phi >= phimin4 && phi <= phimax4){
456 Int_t id4 = fGrid4->GetIndex(phi,eta)-1+n0+n1+n2+n3;
457 AliJetUnitArray *uArray4 = (AliJetUnitArray*)fUnitArray->At(id4+fNumUnits+nElements);
458 TRefArray *reference4 = uArray4->GetUnitTrackRef();
a9c14d46 459 if (reference4->GetEntries() == 0) {
460 new(reference4) TRefArray(TProcessID::GetProcessWithUID(track));
461 }
462
be6e5811 463 reference4->Add(track);
464
465 Float_t uEnergy4 = uArray4->GetUnitEnergy();
466 Bool_t ok4 = kFALSE;
467 Bool_t ref4 = kFALSE;
468 if(uEnergy4==0.){
469 nTracksEmcalDZ++;
470 ok4 = kTRUE;
471 ref4 = kTRUE;
472 }
473 uArray4->SetUnitEnergy(uEnergy4+pt);
474
475 if(uArray4->GetUnitEnergy()<ptMin)
476 uArray4->SetUnitCutFlag(kPtSmaller);
477 else {
478 uArray4->SetUnitCutFlag(kPtHigher);
479 if(ok4) nTracksEmcalDZCut++;
480 }
481 if(sflag[goodTrack] == 1) {
482 uArray4->SetUnitSignalFlag(kGood);
483 uArray4->SetUnitSignalFlagC(kFALSE,kGood);
484 } else uArray4->SetUnitSignalFlagC(kFALSE,kBad);
485
486 if(uArray4->GetUnitEnergy()>0 && ref4){
487 if(!fProcId){
488 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray4));
489 fProcId = kTRUE;
490 }
491 fRefArray->Add(uArray4);
492 }
493 }
494 } // end fDZ
495
496 Int_t towerID = 0;
497 fGeom->GetAbsCellIdFromEtaPhi(eta,phi,towerID);
498 if(towerID==-1) continue;
499
500 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(towerID);
501 TRefArray *reference = uArray->GetUnitTrackRef();
a9c14d46 502 if (reference->GetEntries() == 0) {
503 new(reference) TRefArray(TProcessID::GetProcessWithUID(track));
504 }
505
be6e5811 506 reference->Add(track);
507
508 Float_t unitEnergy = uArray->GetUnitEnergy();
509 Bool_t ok = kFALSE;
510 Bool_t ref = kFALSE;
511 if(unitEnergy==0.){
512 nTracksEmcal++;
513 ok=kTRUE;
514 ref=kTRUE;
515 }
516
517 // Do Hadron Correction
e36a3f22 518 // For the moment I apply MIP correction if p >= 0.5 GeV/c
519 // and only for tracks inside EMCal acceptance
520 // End of if(fGrid==1) -> hadron correction for all tracks
aa862876 521// Temporarily commented - will be updated removing the AliAODpid dependence
522/*
e36a3f22 523 if (fApplyMIPCorrection != 0 && p3.Mag() >= 0.5)
aa862876 524 {
f445ed00 525 ((AliJetHadronCorrectionv1*)fHadCorr)->SetGeometry("EMCAL_COMPLETEV1",1.);
e36a3f22 526
527 Double_t etaOut = 0.;
528 Double_t phiOut = 0.;
529 ((AliJetHadronCorrectionv1*)fHadCorr)->TrackPositionEMCal(track,etaOut,phiOut);
530
531 // One has to make sure that the track hits the calorimeter
532 // We can then correct on average for these particles
533 if((etaOut >= fEtaMin && etaOut <= fEtaMax) &&
534 (phiOut >= fPhiMin && phiOut <= fPhiMax))// &&
535 {
536 Double_t hCEnergy = (Double_t)fHadCorr->GetEnergy(p3.Mag(), (Double_t)eta,0);
537 unitEnergy -= hCEnergy*TMath::Sin(2.0*TMath::ATan(TMath::Exp(-eta)));
538 }
be6e5811 539 } //end Hadron Correction loop
aa862876 540*/
be6e5811 541
e36a3f22 542 if((unitEnergy + pt) > 0.) uArray->SetUnitEnergy(unitEnergy + pt);
543 else uArray->SetUnitEnergy(0.);
be6e5811 544
545 // Put a pt cut flag
546 if(uArray->GetUnitEnergy()<ptMin){
547 uArray->SetUnitCutFlag(kPtSmaller);
548 }
549 else {
550 uArray->SetUnitCutFlag(kPtHigher);
551 if(ok) nTracksEmcalCut++;
552 }
553
554 // Signal flag
555 if(sflag[goodTrack] == 1) {
556 uArray->SetUnitSignalFlag(kGood);
557 uArray->SetUnitSignalFlagC(kFALSE,kGood);
558 } else uArray->SetUnitSignalFlagC(kFALSE,kBad);
559
560 if(uArray->GetUnitEnergy()>0 && ref){
561 if(!fProcId){
562 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
563 fProcId = kTRUE;
564 }
565 fRefArray->Add(uArray);
566 }
567 } // end loop on EMCal acceptance cut + tracks quality
568 else{
569 // Outside EMCal acceptance
570 // Int_t idTPC = GetIndexFromEtaPhi(etaT[i],phiT[i]);
571 Int_t idTPC = fTPCGrid->GetIndex(phi,eta);
572
573 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(fNumUnits-1+idTPC);
574 TRefArray *reference = uArray->GetUnitTrackRef();
a9c14d46 575 if (reference->GetEntries() == 0) {
576 new(reference) TRefArray(TProcessID::GetProcessWithUID(track));
577 }
be6e5811 578 reference->Add(track);
579
580 Float_t unitEnergy2 = uArray->GetUnitEnergy(); // check if fNumUnits or fNumUnits-1
581 Bool_t okout = kFALSE;
582 Bool_t refout = kFALSE;
583 if(unitEnergy2==0.){
584 nTracksTpc++;
585 okout=kTRUE;
586 refout=kTRUE;
587 }
588 // Fill energy outside emcal acceptance
589 uArray->SetUnitEnergy(unitEnergy2 + pt);
590
591 // Pt cut flag
592 if(uArray->GetUnitEnergy()<ptMin){
593 uArray->SetUnitCutFlag(kPtSmaller);
594 }
595 else {
596 uArray->SetUnitCutFlag(kPtHigher);
597 if(okout) nTracksTpcCut++;
598 }
599 // Signal flag
600 if(sflag[goodTrack] == 1) {
601 uArray->SetUnitSignalFlag(kGood);
602 uArray->SetUnitSignalFlagC(kFALSE,kGood);
603 } else uArray->SetUnitSignalFlagC(kFALSE,kBad);
604
605 if(uArray->GetUnitEnergy()>0 && refout){
606 if(!fProcId){
607 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
608 fProcId = kTRUE;
609 }
610 fRefArray->Add(uArray);
611 }
612 }
e36a3f22 613
614 /*
615 // Do hadron correction
616 // In principle, the best thing to do as particles which
617 // eta,phi (at vertex) correspond to EMCal dead zone can deposit energy
618 // in the calorimeter as well as particles with eta,phi (at vertex)
619 // outside the EMCal acceptance
620 cout << "Hadronic correction for all tracks goes here" << endl;
621
622 // For the moment I apply MIP correction if p >= 0.5 GeV/c
623 if (fApplyMIPCorrection != 0 && p3.Mag() >= 0.5)
624 {
f445ed00 625 ((AliJetHadronCorrectionv1*)fHadCorr)->SetGeometry("EMCAL_COMPLETEV1",1.);
e36a3f22 626
627 Double_t etaOut = 0.;
628 Double_t phiOut = 0.;
629 Int_t towerID2 = 0;
630 Float_t unitEnergy = 0.;
631 ((AliJetHadronCorrectionv1*)fHadCorr)->TrackPositionEMCal(track,etaOut,phiOut);
632
633 // One has to make sure that the track hits the calorimeter
634 // We can then correct on average for these particles
635 if((etaOut >= fEtaMin && etaOut <= fEtaMax) &&
636 (phiOut >= fPhiMin && phiOut <= fPhiMax))// &&
637 {
638 Int_t towerID = 0;
639 Int_t towerID2 = 0;
640 fGeom->GetAbsCellIdFromEtaPhi(etaOut,phiOut,towerID2);
641 if(towerID2==-1) continue;
642
643 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(towerID2);
644 unitEnergy = uArray->GetUnitEnergy();
645 Bool_t ok = kFALSE;
646 Bool_t ref = kFALSE;
647 if(unitEnergy==0.){
648 nTracksEmcal++;
649 ok=kTRUE;
650 ref=kTRUE;
651 }
652
653 Double_t hCEnergy = (Double_t)fHadCorr->GetEnergy(p3.Mag(), (Double_t)eta,0);
654 unitEnergy -= hCEnergy*TMath::Sin(2.0*TMath::ATan(TMath::Exp(-eta)));
655
656 // The energy in this unitarray can be negative
657 // If the Digits are then filled, make sure the energy sum
658 // is positive. If not, put the value to zero in AliJetAODUnitArrayEMCalDigit.
659 // If they are not filled, put this value to zero
660 uArray->SetUnitEnergy(unitEnergy);
661
662 if(uArray->GetUnitEnergy()!=0. && ref){
663 if(!fProcId){
664 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
665 fProcId = kTRUE;
666 }
667 fRefArray->Add(uArray);
668 }
669 }
670 else cout << "Track out of EMCal acceptance" << endl;
671
672 } //end Hadron Correction loop
673 */
674
be6e5811 675 } // end fGrid==1
676
677 goodTrack++;
678
679 } // end loop on entries (tpc tracks)
680
681
682 if(fDebug>0)
683 {
684 cout << "End of Tracks %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << endl;
685 cout << "goodTracks: " << goodTrack << endl;
686 }
687
e36a3f22 688 fSignalFlag.Set(goodTrack,sflag);
689 fCutFlag.Set(goodTrack,cflag);
690
be6e5811 691 delete[] sflag;
692 delete[] cflag;
693
694 // } // end loop on entries (tpc tracks)
695
696 if(fGrid==0) {
697 fNTracks = nTracksTpcOnly;
698 fNTracksCut = nTracksTpcOnlyCut;
699 if(fDebug>10){
700 cout << "fNTracks : " << fNTracks << endl;
701 cout << "fNTracksCut : " << fNTracksCut << endl;
702 }
703 }
704 if(fGrid==1) {
705 fNTracks = nTracksEmcal+nTracksEmcalDZ+nTracksTpc;
706 fNTracksCut = nTracksEmcalCut+nTracksEmcalDZCut+nTracksTpcCut;
707 if(fDebug>10){
708 cout << "fNTracks : " << fNTracks << endl;
709 cout << "fNTracksCut : " << fNTracksCut << endl;
710 }
711 }
712
713}
714