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