]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PWG2/FLOW/AliFlowCommon/AliStarEventReader.cxx
changed behavior cut on integer values
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowCommon / AliStarEventReader.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 // Title:   Class for ccessing data from STAR NTuples
17 //          produces data encapsulated in AliStarEvent and AliStarTrack classes
18 //
19 // Origin:  Jim Thomas,        jhthomas@lbl.gov
20 //          Mikolaj Krzewicki, mikolaj.krzewicki@cern.ch
21
22 #include <Riostream.h>
23
24 #include <TSystem.h>
25 #include <TFile.h>
26 #include <TList.h>
27 #include <TLeaf.h>
28 #include <TMath.h>
29 #include <TNtuple.h>
30
31 #include "AliStarEventReader.h"
32 #include "AliStarEvent.h"
33 #include "AliStarTrack.h"
34
35 ClassImp(AliStarEventReader)
36
37 //______________________________________________________________________________
38 AliStarEventReader::AliStarEventReader():
39   TObject(),
40   fFileList(NULL),
41   fEventHeader(NULL),
42   fTracks(NULL),
43   fEvent(NULL)
44 {
45   //ctor
46 }
47 //______________________________________________________________________________
48 AliStarEventReader::AliStarEventReader( const char* inputFileDirectory ):
49   TObject(),
50   fFileList(NULL),
51   fEventHeader( new TNtuple("event","event",
52                             "runId:eventNumber:VtxX:VtxY:VtxZ:BField:refMult:centralityId:numberOfPrimaryTracks:numberOfParticles:h0:h1:h2:h3:h4" )),
53   fTracks( new TNtuple("tracks","tracks",
54                        "ID:Charge:Eta:Phi:Pt:Dca:nHits:nHitsFit:nHitsPoss:nHitsDedx:dEdx:nSigElect:nSigPi:nSigK:nSigProton" )),
55   fEvent(new AliStarEvent(1024))
56 {
57   //ctor
58   MakeFileList ( inputFileDirectory ) ;
59 }
60
61 //______________________________________________________________________________
62 AliStarEventReader::~AliStarEventReader()
63 {
64   //dtor
65   delete fEventHeader;
66   fEventHeader = NULL ;
67   delete fTracks;
68   fTracks = NULL ;
69   delete fFileList;
70   fFileList = NULL ;
71   delete fEvent;
72 }
73
74 //______________________________________________________________________________
75 Bool_t AliStarEventReader::GetNextEvent( )
76 {
77   //gets next event
78   static TFile*    nextFile    = NULL ;
79   static TNtuple*  ntData      = NULL ;
80   static Int_t     DoOnce      = 0 ;
81   static Long64_t  NextEntry   = 0 ;
82   static Long64_t  entries     = 0 ;
83   static Long64_t  FileCounter = 0 ;
84
85   if ( DoOnce == 0 )
86   {
87     DoOnce = 1     ;
88     nextFile = (TFile*)  fFileList->First() ;
89     if ( nextFile == 0 ) return false      ;
90     ntData        = (TNtuple*) ( nextFile->Get("NTtracks") ) ;
91     entries       = ntData->GetEntriesFast() ;
92     Int_t columns = ntData->GetNvar() ;
93     if ( columns != 15 )
94     {
95       cout << "Error in reading Ntuple file: no. columns != 15" << endl ;
96       return false ;
97     }
98     FileCounter++ ;
99     cout << "Start New File " << FileCounter << endl ;
100   }
101
102   while ( nextFile )
103   {
104     while ( NextEntry < entries )
105     {
106       Float_t* header = NULL;
107       Int_t numberOfParticles =  0 ;                   // Number of particle tracks in the next event
108       Long64_t HeaderEntry    =  0 ;                   // Store position of Header and Set Flag in case of EOF or error
109       Long64_t SkipEvent      =  0 ;                   // Flag in case of wrong number of tracks in this event
110
111       ////////////////////////////////////////////
112       fEvent->Reset();           //reset the event
113       ////////////////////////////////////////////
114
115       // Search for the first "Event" record
116
117       delete fEventHeader  ;
118       fEventHeader  = NULL ;               // Note that a first 'event'  ntuple was created in constructor
119       fEventHeader  = new TNtuple("event","event",
120                                   "runId:eventNumber:VtxX:VtxY:VtxZ:BField:refMult:centralityId:numberOfPrimaryTracks:numberOfParticles:h0:h1:h2:h3:h4" ) ;
121
122       for ( Long64_t j = NextEntry ; j < entries ; j++ )
123       {
124         Long64_t BytesRead = ntData->GetEntry(j) ;
125         if ( BytesRead < 60 )
126         {
127           cout << "Warning: error in file or EOF " <<  endl ;
128           HeaderEntry = -1 ;
129           break ;
130         }
131         header = ntData->GetArgs() ;
132         if ( (int) header[10] == -1 && (int) header[11] == -1 && (int) header[12] == -1 &&
133              (int) header[13] == -1 && (int) header[14] == -1 )
134         {
135           fEventHeader->Fill(header) ;
136
137           //////////////////////////////////////////////////
138           fEvent->SetParams(header);  //set the event params
139           //////////////////////////////////////////////////
140
141           numberOfParticles = (int) header[9]  ;   // # of particles passing track cuts, thus in ntuple
142           HeaderEntry = j ;
143           break ;
144         }
145         cout << "Warning: no header entries found in this file" << endl ;
146         HeaderEntry = -1 ;
147       }
148
149       if ( HeaderEntry == -1 ) break ;                 // Break out of main loop if I/O error
150
151       // Get subsequent "track" data
152
153       delete fTracks ;
154       fTracks = NULL ;                 // Note that a first 'tracks' ntuple was created in constructor
155       fTracks = new TNtuple("tracks","tracks",
156                             "ID:Charge:Eta:Phi:Pt:Dca:nHits:nHitsFit:nHitsPoss:nHitsDedx:dEdx:nSigElect:nSigPi:nSigK:nSigProton" ) ;
157
158       for ( Long64_t j = HeaderEntry + 1 ; j < HeaderEntry + 1 + numberOfParticles  ; j++ )
159       {
160         Long64_t BytesRead = ntData->GetEntry(j) ;
161         if ( BytesRead < 60 )
162         {
163           cout << "Warning: error in file sequence or EOF" << endl ;
164           NextEntry = -1 ;
165           break ;
166         }
167         header = ntData->GetArgs() ;
168
169         if ( TMath::IsNaN(header[10]) == 1 )
170         {
171           cout << "IsNan ... dEdx will be zeroed out" << endl ;
172           header[10] = 0 ;
173           header[11] = 999 ;
174           header[12] = 999 ;
175           header[13] = 999 ;
176           header[14] = 999 ;
177           cout << header[0]  << " " << header[1]  << " " << header[2]  << " " << header[3]  << " "
178                << header[4]  << " " << header[5]  << " " << header[6]  << " " << header[7]  << " "
179                << header[8]  << " " << header[9]  << " " << header[10] << " " << header[11] << " "
180                << header[12] << " " << header[13] << " " << header[14] << endl ;  // JT test
181         }
182
183         if ( (int) header[10] == -1 && (int) header[11] == -1 && (int) header[12] == -1 &&
184              (int) header[13] == -1 && (int) header[14] == -1 )
185         {
186           cout << "Warning: Header in the wrong place, skipping event" << endl ;
187           SkipEvent = 1 ;
188           NextEntry = j ;          // Skip event and freeze NextEntry counter
189           break ;
190         }
191
192         fTracks->Fill(header) ;
193
194         ////////////////////////////////////////////////////////////////////
195         fEvent->AddTrack( new AliStarTrack(header) );    //add the new track
196         ////////////////////////////////////////////////////////////////////
197
198         NextEntry = j+1 ;
199       }
200       if ( NextEntry == -1 ) break ;      // Bad record in file, go to next file in fFileList
201       if ( SkipEvent ==  1 ) continue ;   // Bad event, go to next event in this file
202       return true ;                       // Success: Event read OK, note unusual location for a successful return
203     }
204
205     NextEntry = 0 ; // this entry goes before nextFile
206     nextFile = (TFile*) fFileList->After(nextFile) ;
207     if ( nextFile == 0 ) break ;
208     if (ntData) delete ntData;
209     ntData        = (TNtuple*) ( nextFile->Get("NTtracks") ) ;
210     entries       = ntData->GetEntriesFast() ;
211     Int_t columns = ntData->GetNvar() ;
212     if ( columns != 15 )
213     {
214       cout << "Error in reading Ntuple file: no. columns != 15" << endl ;
215       break ;
216     }
217     FileCounter++ ;
218     cout << "Start New File " << FileCounter << endl ;
219   }
220
221   return false ;  // Failure: Error or EOF
222 }
223
224 //______________________________________________________________________________
225 Bool_t AliStarEventReader::AcceptEvent( AliStarEvent* event )
226 {
227   // Cut parameters for each event
228
229   const Float_t VertexXMin  =  -1.0 ;  // cm
230   const Float_t VertexXMax  =   1.0 ;  // cm
231   const Float_t VertexYMin  =  -1.0 ;  // cm
232   const Float_t VertexYMax  =   1.0 ;  // cm
233   const Float_t VertexZMin  = -30.0 ;  // cm
234   const Float_t VertexZMax  =  30.0 ;  // cm
235   const Int_t   MultMin     =    10 ;  // Note: this is a cut on refMult which is not your ordinary multiplicity.
236   const Int_t   MultMax     =  1000 ;  // Refmult corresponds to Zhangbu's ascii table of cuts for centrality bins.
237   const Int_t   BlackEvent  =  3000 ;  // Maximum number of primary tracks allowed in one event (more is an error)
238
239   Int_t refMult               = event->GetRefMult();  // Reference Multiplicity for Centrality determination
240   Int_t NumberOfPrimaryTracks = event->GetNumberOfPrimaryTracks();  // Number of primary tracks in full event (not just those in main TPC)
241
242   if ( refMult < MultMin          || refMult > MultMax )                   return false ;
243   if ( NumberOfPrimaryTracks <= 0 || NumberOfPrimaryTracks > BlackEvent )  return false ;
244
245   // Cut on Vertex location
246   Float_t vertex[3] = { event->GetVtxX(),
247                         event->GetVtxY(),
248                         event->GetVtxZ() };
249
250   if ( vertex[0] < VertexXMin || vertex[0] > VertexXMax )    return false ;  // Skip events that fall outside Vtx cuts
251   if ( vertex[1] < VertexYMin || vertex[1] > VertexYMax )    return false ;
252   if ( vertex[2] < VertexZMin || vertex[2] > VertexZMax )    return false ;
253
254   return true ;
255 }
256
257 //______________________________________________________________________________
258 Bool_t AliStarEventReader::AcceptTrack( AliStarTrack* track )
259 {
260   // Cut Parameters for individual tracks
261
262   const Float_t dcaCut      =   3.0  ;      // cm
263   const Float_t PtMin       =   0.15 ;      // GeV
264   const Float_t PtMax       =   2.0  ;      // GeV
265   const Float_t EtaMin      =  -1.1  ;
266   const Float_t EtaMax      =   1.1  ;
267   const Float_t FitRatio    =  0.52  ;      // Number of hits over number of hits possible
268   const Int_t   nHitMin     =    15  ;      // 15 is typical but sometimes goes as high as 25
269   const Int_t   nHitMax     =   100  ;      // 45 pad rows in the TPC and so anything bigger than 45+Silicon is infinite
270   const Int_t   nHitPossMin =     5  ;      // Don't bother to fit tracks if # possible hits is too low, also protect / 0
271
272   if ( track->GetDCA() >  dcaCut ) return false ;                 // magnitude of 3D DCA for global tracks
273   if ( track->GetNHitsPoss() <  nHitPossMin ||
274        track->GetNHitsPoss() >  nHitMax ) return false ;   // Minimum number of Possible hits, see above.
275   if ( track->GetEta() <  EtaMin  || track->GetEta() >  EtaMax  ) return false ;
276   if ( track->GetPt()  <  PtMin   || track->GetEta() >  PtMax   ) return false ;
277   if ( track->GetNHitsFit() < nHitMin || track->GetNHitsFit() > nHitMax ) return false ;
278   if ( ((1.0*track->GetNHitsFit()) /
279         (1.0*track->GetNHitsPoss())) < FitRatio )  return false ;
280
281   return true ;
282 }
283
284 //______________________________________________________________________________
285 Int_t AliStarEventReader::ParticleID( AliStarTrack* track )
286 {
287   // Note: This is a very simple PID selection scheme.  More elaborate methods (with multiple cuts) may be required.
288   // When you *are* using dEdx information, you must chose a finite number of good Dedx hits ... but the limit should
289   // be about 2/3 of nHitsMin.  This is because some clusters do not form good dEdx hits due to track
290   // merging, etc., and so nHitsDedx is always less than nHitsFit.  A rule of thumb says ~2/3 ratio.
291
292   Int_t ID = 0 ;
293
294   const Int_t   nHitDedxMin =    15  ;       // 10 to 20 is typical.  nHitDedxMin is often chosen to be about 2/3 of nHitMin.
295   const Float_t nSigmaPID   =    2.0 ;       // Number of Sigma cut to apply to PID bands
296
297   // Test on Number of dE/dx hits required, return 0 if not enough hits
298   if ( track->GetNHitsDedx() <  nHitDedxMin ) return ID;
299
300   // Begin PID
301
302   if ( TMath::Abs( track->GetNSigElect() ) >= nSigmaPID )
303   {
304     if ( TMath::Abs( track->GetNSigK()  ) <= nSigmaPID )
305     {
306       ID = 321  ;
307     }
308     if ( TMath::Abs( track->GetNSigProton()  ) <= nSigmaPID )
309     {
310       ID = 2212 ;
311     }
312     if ( TMath::Abs( track->GetNSigPi()  ) <= nSigmaPID )
313     {
314       ID = 211  ;
315     }
316   }
317
318   // Pion is the default in case of ambiguity because it is most abundent. Don't re-arrange order, above.
319
320   return ID ;
321 }
322
323 //______________________________________________________________________________
324 Int_t AliStarEventReader::Centrality( Int_t referenceMultiplicity )
325
326 {
327   // Note Carefully:  Centrality is based on refMult.  This is the 'reference' multiplicity that is measured
328   // independpently from the TPC.  Selecting the centrality bins according to the refMult is something that
329   // is calibrated for each year and each run.  You can get the basic information off the web:
330   // For Example .... http://www.star.bnl.gov/protected/common/common2004/trigger2004/200gev/200gevFaq.html
331   // An index pointing to FAQs, Trigger and Centrality data, for all years, is available at:
332   // http://www.star.bnl.gov/public/all
333   //
334   // Note: Add 0.5 to the (int) that is returned by centrality() when using it as an argument for a histogram
335   // that expects (float) or (double) as input parameters.  This will place the data point in the center of
336   // the bin, avoids ambiguities, and is best for plotting scatter plots and contour plots.
337   // For example histogram2D[1] -> Fill ( (float)CentralityID + 0.5 , SumData )   ;
338   //
339   // The refMult quoted in the Centrality bins array is the lower limit on refMult
340
341
342   Int_t   CentralityBins  [] = { 14 , 31 , 57 , 96 , 150 , 222 , 319 , 441 , 520 , 1000 } ;  // Run4 200 GeV
343   Int_t   MiddleBinID     [] = {  0 ,  1 ,  2 ,  3 ,   4 ,   5 ,   6 ,   7 ,   8 ,    9 } ;  // ID Number
344   //Info  MiddleBinPercent[] = { 85., 75., 65., 55.,  45.,  35.,  25.,  15., 7.5 ,  2.5 } ;  // Percent
345   Int_t   myCentrality  ;
346
347   if      ( referenceMultiplicity < CentralityBins[0] )
348   {
349     myCentrality = MiddleBinID[0] ;
350   }
351   else if ( referenceMultiplicity < CentralityBins[1] )
352   {
353     myCentrality = MiddleBinID[1] ;
354   }
355   else if ( referenceMultiplicity < CentralityBins[2] )
356   {
357     myCentrality = MiddleBinID[2] ;
358   }
359   else if ( referenceMultiplicity < CentralityBins[3] )
360   {
361     myCentrality = MiddleBinID[3] ;
362   }
363   else if ( referenceMultiplicity < CentralityBins[4] )
364   {
365     myCentrality = MiddleBinID[4] ;
366   }
367   else if ( referenceMultiplicity < CentralityBins[5] )
368   {
369     myCentrality = MiddleBinID[5] ;
370   }
371   else if ( referenceMultiplicity < CentralityBins[6] )
372   {
373     myCentrality = MiddleBinID[6] ;
374   }
375   else if ( referenceMultiplicity < CentralityBins[7] )
376   {
377     myCentrality = MiddleBinID[7] ;
378   }
379   else if ( referenceMultiplicity < CentralityBins[8] )
380   {
381     myCentrality = MiddleBinID[8] ;
382   }
383   else
384   {
385     myCentrality = MiddleBinID[9] ;
386   }
387
388   return myCentrality ;
389 }
390
391 //______________________________________________________________________________
392 void AliStarEventReader::PrintEventHeader ( )
393 {
394   // TNtuple* event: names are documented in the next 2 lines
395   // event  = new TNtuple("event","event",
396   //   "runId:eventNumber:VtxX:VtxY:VtxZ:BField:refMult:centralityId:numberOfPrimaryTracks:numberOfParticles:h0:h1:h2:h3:h4" ) ;
397   //
398   Float_t  *eventhd ;
399   eventhd = fEventHeader->GetArgs() ;
400
401   Int_t   runId                  = (int)   eventhd[0]  ;
402   Int_t   eventNumber            = (int)   eventhd[1]  ;
403   Float_t   primaryVertexPosition[3] = { (float) eventhd[2],  // (cm)
404                                          (float) eventhd[3],  // (cm)
405                                          (float) eventhd[4] };  // (cm)
406   Float_t magneticField          = (float) eventhd[5]  ;  // kilogauss
407   Int_t   refMult                = (int)   eventhd[6]  ;  // Raw Mult into 0.5 unit: a relative number, not total Mult.
408   Int_t   centralityId           = (int)   eventhd[7]  ;  // STAR centrality bin # based on refMult
409   Int_t   numberOfPrimaryTracks  = (int)   eventhd[8]  ;  // # of primaries, including FTPC tracks which are not in ntuple
410   Int_t   numberOfParticles      = (int)   eventhd[9]  ;  // # of particles passing track cuts, thus in ntuple
411
412   printf("\n  runId eventNo    VtxX    VtxY    VtxZ  MagFld  refMult  CentBin  #PrimeTrak  #Tracks \n") ;
413   printf("%7d  %6d %7.4f %7.4f %7.3f  %6.3f   %6d   %6d      %6d   %6d \n\n",
414          runId, eventNumber, primaryVertexPosition[0], primaryVertexPosition[1], primaryVertexPosition[2],
415          magneticField, refMult, centralityId, numberOfPrimaryTracks, numberOfParticles ) ;
416
417   //Int_t newCentralityID ;
418   //newCentralityID = Centrality( refMult) ;              // Should be the same as "centralityID" from tape
419   //cout << "Test: should be the same " << newCentralityID << " " << centralityId << endl ;  // JT test
420 }
421
422 //______________________________________________________________________________
423 void AliStarEventReader::PrintTrack ( Int_t counter )
424 {
425   // TNtuple* track: names are documented in the next 2 lines
426   // tracks = new TNtuple("tracks","tracks",
427   //   "ID:Charge:Eta:Phi:Pt:Dca:nHits:nHitsFit:nHitsPoss:nHitsDedx:dEdx:nSigElect:nSigPi:nSigK:nSigProton" ) ;
428   //
429   if ( counter == 0 )
430   {
431     printf(
432       "    id charge     eta     phi      pt     dca  nHits  nFit nPoss ndEdx   dEdx nSigElec nSigPi  nSigK nSigPr\n") ;
433   }
434   Float_t* data = fTracks -> GetArgs()       ;  // Extract data from the track
435   Int_t   id             = (int)   data[0]   ;  // id - a unique integer for each track in this event
436   Int_t   charge         = (int)   data[1]   ;  // +1 or -1
437   Float_t eta            = (float) data[2]   ;  // Pseudo-rapidity at the vertex
438   Float_t phi            = (float) data[3]   ;  // Phi emission angle at the vertexcd
439   Float_t pt             = (float) data[4]   ;  // Pt of the track at the vertex
440   Float_t dca            = (float) data[5]   ;  // Magnitude of 3D DCA to vertex
441   Int_t   nHits          = (int)   data[6]   ;  // Number of clusters available to the Kalman Filter
442   Int_t   nHitsFit       = (int)   data[7]   ;  // Number of clusters used in the fit (after cuts)
443   Int_t   nHitsPoss      = (int)   data[8]   ;  // Number of possible cluster on track (theoretical max)
444   Int_t   nHitsDedx      = (int)   data[9]   ;  // Number of clusters used in the fit (after dEdx cuts)
445   Float_t dEdx           = 1.e6*(float)data[10]  ;  // Measured dEdx (Note: GeV/cm so convert to keV/cm!!)
446   Float_t nSigmaElectron = (float) data[11]  ;  // Number of sigma from electron Bethe-Bloch curve
447   Float_t nSigmaPion     = (float) data[12]  ;  // Number of sigma from pion Bethe-Bloch curve
448   Float_t nSigmaKaon     = (float) data[13]  ;  // Number of sigma from kaon Bethe-Bloch curve
449   Float_t nSigmaProton   = (float) data[14]  ;  // Number of sigma from proton Bethe-Bloch curve
450
451   // Alternative way to access the data
452   nHitsPoss      = (int) ( fTracks->GetLeaf("nHitsPoss")->GetValue() ) ;  // Note alternative method to retrieve data
453   // Using the definition of the original NTuple
454   // TrackTuple      = new TNtuple("NTtracks","NTtracks",
455   // "ID:Charge:Eta:Phi:Pt:Dca:nHits:nHitsFit:nHitsPoss:nHitsDedx:dEdx:nSigElect:nSigPi:nSigK:nSigProton" )
456
457   printf("%6d %4d   %7.3f %7.3f %7.3f %7.4f %6d %5d %5d %5d %6.2f   %6.2f %6.2f %6.2f %6.2f \n",
458          id, charge, eta, phi, pt, dca, nHits, nHitsFit, nHitsPoss, nHitsDedx, dEdx,
459          nSigmaElectron, nSigmaPion, nSigmaKaon, nSigmaProton ) ;
460 }
461
462 //______________________________________________________________________________
463 Bool_t AliStarEventReader::MakeFileList ( const char* inputFileDirectory )
464 {
465   //get the files to process
466   Int_t  Count        = 0 ;
467   static Int_t DoOnce = 0 ;
468   fFileList =  new TList() ;
469   void*   directory = gSystem->OpenDirectory(inputFileDirectory) ;
470   const char* entry = gSystem->GetDirEntry(directory) ;
471
472   if ( entry == 0 )
473   {
474     cout << endl <<  "Error: \"" << inputFileDirectory << "\" does not exist" << endl << endl ;
475     return false ;
476   }
477   else cout << endl ;
478
479   while(entry != 0)
480   {
481     int len = strlen(entry);
482     if( len >= 5 && strcmp( &entry[len - 5], ".root" ) == 0 )
483     {
484       TString fileName ;
485       fileName = inputFileDirectory ;
486       if( !fileName.EndsWith("/") ) fileName += "/" ;
487       fileName += entry;
488       fFileList->Add ( TFile::Open(fileName) ) ;
489       if ( DoOnce == 0 )
490       {
491         cout << "Add: " << fileName << endl ;
492         DoOnce = 1 ;
493       }
494       Count ++ ;
495     }
496     entry = gSystem->GetDirEntry(directory) ;
497   }
498
499   cout << "Add: " << Count-1 << " more file(s) from this directory for a total of " << Count << " files." << endl ;
500   cout << "Finished creating file list ... preparing to open first file." << endl << endl ;
501   return true ;
502 }
503