1671804c0bb3d68428dfd18c5d26a08b51ede100
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowKineMaker.cxx
1 //////////////////////////////////////////////////////////////////////
2 //
3 // $Id$
4 //
5 // Author: Emanuele Simili
6 //
7 //////////////////////////////////////////////////////////////////////
8 //_____________________________________________________________
9 //
10 // Description: 
11 //        AliFlowKineMaker provides the method to create AliFlowEvent(s)
12 // creates AliFlowEvent from the KineTree . 
13 // TParticle(s) is translated into AliFlowTrack(s), with exact momentum, 
14 // P.Id., etc. Very basic track cuts are applyed (like primaries). 
15 // The present class can be used in a simple AliRoot macro or in a 
16 // more complex enviroment such as AliSelector or AliTask.
17 //
18 //////////////////////////////////////////////////////////////////////
19
20 #ifndef ALIFLOWKINEMAKER_CXX
21 #define ALIFLOWKINEMAKER_CXX
22
23 // ROOT things
24 #include <TROOT.h>
25 #include <TFile.h>
26 #include <TString.h>
27 #include <TMath.h>
28 #include <TTree.h>
29 #include "TClonesArray.h"
30 #include "TParticle.h"
31 #include "TParticlePDG.h"
32 //#include "TDatabasePDG.h"
33
34 // AliRoot things (...not used here, but in the macro)
35 //#include "AliRun.h"
36 //#include "AliRunLoader.h"
37 //#include "AliStack.h"
38
39 // Flow things
40 #include "AliFlowEvent.h"
41 #include "AliFlowTrack.h"
42 #include "AliFlowV0.h"
43 #include "AliFlowConstants.h"
44 #include "AliFlowKineMaker.h"
45
46 // ANSI things
47 #include <stdlib.h>
48 using namespace std; //required for resolving the 'cout' symbol
49
50 ClassImp(AliFlowKineMaker) 
51 //-----------------------------------------------------------------------
52 AliFlowKineMaker::AliFlowKineMaker():
53   fKTree(0x0), fParticle(0x0), fParticlePDG(0x0),
54   fFlowEvent(0x0), fFlowTrack(0x0), fFlowV0(0x0)
55 {
56  // default constructor 
57  // resets counters , sets defaults
58  
59  fNewAli = kFALSE ;
60  // flags
61  fLoopParts = kTRUE ;
62  fCounter = 0 ; 
63  // loop variable
64  fRunID = 0 ;       
65  fNumberOfParticles = 0 ;
66  fEventNumber = 0 ; 
67  fPartNumber = 0 ; 
68  fEventNumber = 0 ;
69  fMagField = 0. ;
70  // counters
71  fGoodTracks = 0  ;  
72  fGoodV0s    = 0  ;     
73  fGoodTracksEta = 0  ;
74  fPosiTracks = 0  ;  
75  fNegaTracks = 0  ;  
76  fUnconstrained = 0  ;
77  fCutParts   = 0  ;     
78  for(Int_t bb=0;bb<5;bb++) { fBayesianAll[bb] = 0 ; } ;
79  fSumAll = 0 ; 
80  fCharge = 0 ;
81
82  // trak cuts
83  fPrimary = kTRUE ;
84  fAbsEta = 2.1 ;                        
85  fElow = 0.001 ; fEup = 1000. ;   
86  fLabel[0] = 0 ; fLabel[1] = -1 ;
87
88 //  // TGeant3::AddParticlesToPdgDataBase() ---  Stolen From TGeant3.cxx ----(
89 //  TDatabasePDG *pdgDB = TDatabasePDG::Instance();
90 //  const Int_t kion=10000000;
91 //  const Int_t kspe=50000000;
92 //  const Double_t kAu2Gev=0.9314943228;
93 //  const Double_t khSlash = 1.0545726663e-27;
94 //  const Double_t kErg2Gev = 1/1.6021773349e-3;
95 //  const Double_t khShGev = khSlash*kErg2Gev;
96 //  const Double_t kYear2Sec = 3600*24*365.25;
97 //  // Ions
98 //  pdgDB->AddParticle("Deuteron","Deuteron",2*kAu2Gev+8.071e-3,kTRUE,0,3,"Ion",kion+10020);
99 //  pdgDB->AddParticle("Triton","Triton",3*kAu2Gev+14.931e-3,kFALSE,khShGev/(12.33*kYear2Sec),3,"Ion",kion+10030);
100 //  pdgDB->AddParticle("Alpha","Alpha",4*kAu2Gev+2.424e-3,kTRUE,khShGev/(12.33*kYear2Sec),6,"Ion",kion+20040);
101 //  pdgDB->AddParticle("HE3","HE3",3*kAu2Gev+14.931e-3,kFALSE,0,6,"Ion",kion+20030);
102 //  // Special particles
103 //  pdgDB->AddParticle("Cherenkov","Cherenkov",0,kFALSE,0,0,"Special",kspe+50);
104 //  pdgDB->AddParticle("FeedbackPhoton","FeedbackPhoton",0,kFALSE,0,0,"Special",kspe+51);
105 //  // ----------------------------------------------------------------------)
106 }
107 //-----------------------------------------------------------------------
108 AliFlowKineMaker::~AliFlowKineMaker()
109 {
110  // default destructor (no actions) 
111 }
112 //-----------------------------------------------------------------------
113
114 //-----------------------------------------------------------------------
115 AliFlowEvent* AliFlowKineMaker::FillFlowEvent(TTree* fKTree)
116 {
117  // From the MC KineTree (input) fills the AliFlowEvent (output) . 
118  // It loops on the stored TParticles and calls the methods to fill the 
119  // arrays in the AliFlowEvent (charged -> tracks , neutral -> v0s) . 
120
121  fFlowEvent = new AliFlowEvent() ; if(!fFlowEvent) { return 0 ; }
122  //cout << " -evt- " << fFlowEvent << endl ;
123
124  fRunID = -1 ;
125  fEventNumber = -1 ;
126  fNumberOfParticles = fKTree->GetEntries() ; 
127  //
128  cout << " *evt n. " << fEventNumber << " (run " << fRunID << ")  -  tracks/v0s : " << fNumberOfParticles << endl ;
129
130  // Event id 
131  fFlowEvent->SetRunID(fRunID) ;                
132  fFlowEvent->SetEventID(fEventNumber) ;         
133  fFlowEvent->SetOrigMult((UInt_t)fNumberOfParticles) ;
134
135  // Run information (fixed - ???)
136  fMagField = 4 ; // (?)
137  fFlowEvent->SetMagneticField(fMagField) ;      
138  fFlowEvent->SetCenterOfMassEnergy(AliFlowConstants::fgCenterOfMassEnergy) ;  
139  fFlowEvent->SetBeamMassNumberEast(AliFlowConstants::fgBeamMassNumberEast) ;  
140  fFlowEvent->SetBeamMassNumberWest(AliFlowConstants::fgBeamMassNumberWest) ;  
141
142  // Trigger information (now is: ULon64_t - some trigger mask)
143  fFlowEvent->SetL0TriggerWord(-1); 
144  
145  // Get primary vertex position
146  fVertex[0] = 0. ; fVertex[1] = 0. ; fVertex[2] = 0. ;         // fVertex = // ?! how to get primary vertex !?
147  fFlowEvent->SetVertexPos(fVertex[0],fVertex[1],fVertex[2]) ; 
148
149  // Zero Degree Calorimeter information
150  Int_t zdcp = (Int_t)(TMath::Sqrt(TMath::Sqrt(fNumberOfParticles))) ;
151  Float_t zdce[3] ; zdce[0] = -1 ; zdce[1] = -1 ; zdce[2] = -1 ;
152  fFlowEvent->SetZDCpart(zdcp);                          
153  fFlowEvent->SetZDCenergy(zdce[0],zdce[1],zdce[2]);     
154
155  fKTree->SetBranchAddress("Particles",&fParticle) ;
156
157  // Track (& V0) loop
158  if(fLoopParts)
159  {
160   Int_t badPart = 0 ;   
161   for(fPartNumber=0;fPartNumber<fNumberOfParticles;fPartNumber++) 
162   {
163    fKTree->GetEntry(fPartNumber) ;
164    if(CheckTrack(fParticle))
165    {
166     // fParticlePDG = fParticle->GetPDG() ; 
167     fCharge = (Int_t)((fParticle->GetPDG()->Charge())/3) ; // cout << fCharge << endl ;  
168     // fCharge = (Int_t)(TMath::Sign(1,(fParticle->GetPdgCode()))) ;             
169   
170     if(TMath::Abs(fCharge) > 0)
171     {
172      FillFlowTrack(fParticle) ;   
173      fGoodTracks++ ;                                   
174     }
175     else if(fCharge == 0)
176     {
177      FillFlowV0(fParticle) ;
178      fGoodV0s++ ;                        
179     }
180    }
181    else { badPart++ ; continue ; }
182   }
183   fCutParts += badPart ;
184  } // cout << " -particle number- :  " << fPartNumber << endl ;
185
186  // Evt setting stuff
187  fFlowEvent->SetCentrality();   
188                                 
189  fCounter++ ; 
190  return fFlowEvent ;
191 }
192 //----------------------------------------------------------------------
193 AliFlowTrack* AliFlowKineMaker::FillFlowTrack(TParticle* fParticle)
194 {
195  // From a charged TParticle (input) fills the AliFlowTrack (output) .
196
197  TString name = "" ; name += fPartNumber ;
198  Int_t idx = fFlowEvent->TrackCollection()->GetEntries() ;
199  fFlowTrack = (AliFlowTrack*)(fFlowEvent->TrackCollection()->New(idx)) ;
200  fFlowTrack->SetName(name.Data()) ;
201  
202  // cout << " -tr- " << name.Data() << "(" << idx << ")"  << endl ;
203
204  // TParticle label (link: KineTree-ESD)
205  Int_t label = TMath::Abs(fPartNumber);
206  fFlowTrack->SetLabel(label) ;                  
207
208  // signed DCA from ESDtrack
209  Float_t x = fParticle->Vx() ; 
210  Float_t y = fParticle->Vy() ; 
211  Float_t z = fParticle->Vz() ; 
212  Float_t xy = TMath::Sqrt(x*x + y*y) ; 
213  fFlowTrack->SetDcaSigned(xy,z) ;                   
214
215  // error on the DCA = 0
216  fFlowTrack->SetDcaError(0.,0.,0.) ;                
217
218  // UnConstrained (global) first
219  Double_t gD[3] ;                               
220  gD[0] = fParticle->Px() ; gD[1] = fParticle->Py() ; gD[2] = fParticle->Pz() ;                  
221  // -
222  Float_t phiGl = (Float_t)Phi(gD) ;  
223  if(phiGl<0) { phiGl += 2*TMath::Pi() ; }
224  fFlowTrack->SetPhiGlobal(phiGl) ;              
225  Float_t ptGl = (Float_t)Pt(gD) ;  if(ptGl<=0) { cout << " !!! ptGlobal = " << ptGl << endl ; }
226  fFlowTrack->SetPtGlobal(ptGl) ;                
227  Float_t etaGl = (Float_t)Eta(gD) ; 
228  fFlowTrack->SetEtaGlobal(etaGl) ;              
229
230  // Constrained (same, if primary)
231  if((fParticle->IsPrimary()) && (Norm(gD)!=0.))  
232  {                                                      
233   fFlowTrack->SetPhi(phiGl) ;                           
234   fFlowTrack->SetPt(ptGl) ;                             
235   fFlowTrack->SetEta(etaGl) ;                           
236
237   // number of constrainable tracks with |eta| < AliFlowConstants::fgEtaGood (0.9)
238   if(TMath::Abs(etaGl) < AliFlowConstants::fgEtaGood)  { fGoodTracksEta++ ; }
239  }
240  else  // in case Constriction impossible for track, fill the UnConstrained (global)
241  {
242   fUnconstrained++ ;    
243   fFlowTrack->SetPhi(0.) ;
244   fFlowTrack->SetPt(0.) ;   
245   fFlowTrack->SetEta(0.) ; 
246  }           
247
248  // positive - negative tracks
249  
250  //Int_t fCharge = (Int_t)(fParticle->GetPDG()->Charge()) ; 
251  fFlowTrack->SetCharge(fCharge) ;               
252  if(fCharge>0)          { fPosiTracks++ ; }
253  else if(fCharge<0)     { fNegaTracks++ ; }
254  else                   { return 0 ; }
255
256  // Track parametrization (p at, hits, clusters, dE/dx) 
257  Double_t pVecAt[3] ; 
258  for(Int_t gg=0;gg<3;gg++) { pVecAt[gg] = gD[gg] ; }
259  Float_t pAt = (Float_t)Norm(pVecAt) ;
260  
261  Int_t nClus[9] ; Int_t nHits[9] ; 
262  nClus[0] = 1 ;   nHits[0] = 0 ;                            // ITS - pixel
263  nClus[1] = 1 ;   nHits[1] = 0 ;  
264  nClus[2] = 1 ;   nHits[2] = 0 ;                            // ITS - drift
265  nClus[3] = 1 ;   nHits[3] = 0 ;  
266  nClus[4] = 1 ;   nHits[4] = 0 ;                            // ITS - strips
267  nClus[5] = 1 ;   nHits[5] = 0 ;  
268  nClus[6] = 160 ; nHits[6] = 0 ;                            // TPC
269  nClus[7] = 130 ; nHits[7] = 0 ;                            // TRD
270  nClus[8] = 1 ;   nHits[8] = 0 ;                            // TOF
271
272  Int_t pdgcode = 0 ;
273  Float_t dEdx[4] ; for(Int_t de=0;de<4;de++) { dEdx[de] = -1. ; }
274  Float_t detResFun[6] ; for(Int_t de=0;de<6;de++) { detResFun[de] = 0. ; }
275  Float_t zFirst = fVertex[2] ; 
276  Float_t zLast  = 1000. ;
277  Float_t rFirst = TMath::Sqrt((fVertex[0]*fVertex[0]) + (fVertex[1]*fVertex[1])) ; 
278  Float_t rLast  = 1000. ;
279
280 // // Geometrical acceptance (calculated assuming straight tracks) 
281 //
282 //  Float_t rDet[9][2] ; Float_t zDet[9] ;      
283 //  rDet[0][0] = 3.9 ;   rDet[0][1] = 3.9 ;   zDet[0] = 14.1/2. ;     // ITS - pixel
284 //  rDet[1][0] = 7.6 ;   rDet[1][1] = 7.6 ;   zDet[1] = 14.1/2. ;    
285 //  rDet[2][0] = 15.0 ;  rDet[2][1] = 15.0 ;  zDet[2] = 22.2/2. ;     // ITS - drift
286 //  rDet[3][0] = 23.9 ;  rDet[3][1] = 23.9 ;  zDet[3] = 29.7/2. ;    
287 //  rDet[4][0] = 37.8 ;  rDet[4][1] = 38.4 ;  zDet[4] = 43.1/2. ;     // ITS - strips
288 //  rDet[5][0] = 42.8 ;  rDet[5][1] = 43.4 ;  zDet[5] = 48.9/2. ;    
289 //  rDet[6][0] = 84.5 ;  rDet[6][1] = 246.6 ; zDet[6] = 500./2. ;     // TPC
290 //  rDet[7][0] = 290. ;  rDet[7][1] = 370. ;  zDet[7] = 700./2. ;     // TRD
291 //  rDet[8][0] = 370. ;  rDet[8][1] = 399. ;  zDet[8] = 745./2. ;     // TOF
292 //
293 //  Float_t Atheta = fParticle->Pz()/fParticle->Pt() ; 
294 //  Float_t z0 = fParticle->Vz() ; 
295 //  Float_t r0 = TMath::Sqrt((fParticle->Vx()*fParticle->Vx())+(fParticle->Vy()*fParticle->Vy())) ;
296 //  if((fParticle->Vx()*fParticle->Px()+fParticle->Vy()*fParticle->Py())>0)      { r0 *= 1. ; }   // sign given basing on track direction in respect to position
297 //  else if((fParticle->Vx()*fParticle->Px()+fParticle->Vy()*fParticle->Py())<0) { r0 *= -1.; }
298 //  else                                                                         { r0  = 0. ; }
299 //
300 //  // rFirst = rDet[0][0] ; rLast  = rDet[0][0] ;
301 //  zFirst = z0 + Atheta * (rDet[0][0] - r0) ; 
302 //  zLast  = z0 + Atheta * (rDet[4][1] - r0) ;
303 //  Float_t Pin  = 0. ; Float_t Pout = 0. ; Float_t Rout = 0. ; 
304 //  for(int dd=0;dd<9;dd++)
305 //  {
306 //   Pin  = z0 + Atheta * (rDet[dd][0] - r0) ; 
307 //   if(Pin<zDet[dd]) 
308 //   {
309 //    Pout = z0 + Atheta * (rDet[dd][1] - r0) ; 
310 //    if(TMath::Abs(Pout<zDet[dd]))                     // track gets in and out inside acceptance -> full hits
311 //    { 
312 //     nHits[dd] = nClus[dd] ; 
313 //     Rout = rDet[dd][1] ; 
314 //     rLast = TMath::Abs(Rout) ;               
315 //    } 
316 //    else                                              // track goes out from one side -> SOME hits (...)
317 //    {
318 //     Rout = r0 + ((TMath::Sign(zDet[dd],eta)-z0)/Atheta) ; 
319 //     rLast = TMath::Abs(Rout) ;       
320 //     Float_t proportion = TMath::Abs((rLast-rDet[dd][0])/(rDet[dd][1]-rDet[dd][0])) ; proportion *= nClus[dd] ;       
321 //     nHits[dd] = (Int_t)proportion ; 
322 //    }                         
323 //   }
324 //   else                                               // track does not get in -> zero hits
325 //   {
326 //    nHits[0] = 0 ; rFirst = 0. ; //rLast = 0. ; 
327 //   }  
328 //  }
329 //
330 //  if(nHits[7])                { dEdx[0] = 1. ; }  // implement bethe-block for TPC
331 //  if(nHits[5] || nHits[6])    { dEdx[1] = 1. ; }  // implement bethe-block for ITS
332 //  if(nHits[8])                { dEdx[2] = 1. ; }  // implement transition-radiation for TRD
333 //  if(nHits[9])                { dEdx[3] = 1. ; }  // implement time of flight for TOF
334 //
335
336  // P.id. (basing on the pdg code from MC -> an exact P.Id (probability=1) is given for [e , mu , pi , K , p , d], others are 0) 
337  pdgcode = fParticle->GetPdgCode() ;                         // cout << FlowDebug << "PDG code = " << pdgcode << endl ; 
338  if(TMath::Abs(pdgcode) == 11)               { detResFun[0] = 1. ; fBayesianAll[0]++ ; }
339  else if(TMath::Abs(pdgcode) == 13)          { detResFun[1] = 1. ; fBayesianAll[1]++ ; }
340  else if(TMath::Abs(pdgcode) == 211)         { detResFun[2] = 1. ; fBayesianAll[2]++ ; } 
341  else if(TMath::Abs(pdgcode) == 321)         { detResFun[3] = 1. ; fBayesianAll[3]++ ; } 
342  else if(TMath::Abs(pdgcode) == 2212)        { detResFun[4] = 1. ; fBayesianAll[4]++ ; }
343  else if(TMath::Abs(pdgcode) == 10010020)    { detResFun[5] = 1. ; fBayesianAll[5]++ ; }
344  else           { for(Int_t de=0;de<6;de++)  { detResFun[de] = 0.2 ; } }
345  fSumAll++ ;
346
347  // Fill the (fake) track parameters (fit , P.Id. ...)
348  fFlowTrack->SetMostLikelihoodPID(pdgcode); 
349
350  fFlowTrack->SetElectronPositronProb(detResFun[0]);               
351  fFlowTrack->SetMuonPlusMinusProb(detResFun[1]);
352  fFlowTrack->SetPionPlusMinusProb(detResFun[2]);
353  fFlowTrack->SetKaonPlusMinusProb(detResFun[3]);
354  fFlowTrack->SetProtonPbarProb(detResFun[4]);
355  fFlowTrack->SetDeuteriumAntiDeuteriumProb(detResFun[5]);    // *!* implement P.Id. for Deuterium
356
357  fFlowTrack->SetZFirstPoint(zFirst) ;  fFlowTrack->SetZLastPoint(zLast) ;            
358  fFlowTrack->SetTrackLength(TMath::Sqrt((zLast-zFirst)*(zLast-zFirst)+(rLast-rFirst)*(rLast-rFirst))) ;
359  fFlowTrack->SetChi2(0.) ;
360
361  // Fill the (fake) detector information (nHits, dE/dx, det, resp. func.) 
362  fFlowTrack->SetFitPtsTPC(nHits[6]) ;                // cout << FlowDebug << "nHits TPC = " << nHits[6] << endl ; 
363  fFlowTrack->SetMaxPtsTPC(nClus[6]) ;                // cout << FlowDebug << "nClus = " << nClus[6] << endl ;  
364  fFlowTrack->SetChi2TPC(nHits[6]/nClus[6]) ;         // cout << FlowDebug << "Chi2 = " << nHits[6]/nClus[6] << endl ;        
365  fFlowTrack->SetDedxTPC(dEdx[0]) ;                   // cout << FlowDebug << "Dedx = " << dEdx << endl ; 
366  fFlowTrack->SetPatTPC(pAt) ;                        // cout << FlowDebug << "p = " << pAt << endl ;                            
367  fFlowTrack->SetRespFunTPC(detResFun) ;              // cout << FlowDebug << "response function = " << detResFun << endl ;  
368  // -
369  Int_t nITShits = 0 ; for(int dd=0;dd<6;dd++) { nITShits += nHits[dd] ; } 
370  Int_t nITSclus = 6 ; 
371  fFlowTrack->SetFitPtsITS(nITShits) ;                // cout << FlowDebug << "nHits ITS = " << nITShits << endl ;                                   
372  fFlowTrack->SetMaxPtsITS(nITSclus) ;                // cout << FlowDebug << "nClus = " << nITSclus << endl ; 
373  fFlowTrack->SetChi2ITS(nITShits/nITSclus) ;         // cout << FlowDebug << "Chi2 = " << nITShits/nITSclus << endl ; 
374  fFlowTrack->SetDedxITS(dEdx[1]) ;                   // cout << FlowDebug << "Dedx = " << dEdx << endl ; 
375  fFlowTrack->SetPatITS(pAt) ;                        // cout << FlowDebug << "p = " << pAt << endl ;                                    
376  fFlowTrack->SetRespFunITS(detResFun) ;              // cout << FlowDebug << "response function = " << detResFun << endl ;  
377  // -
378  fFlowTrack->SetNhitsTRD(nHits[7]) ;                 // cout << FlowDebug << "nHits TOF = " << nHits[7] << endl ;                                                   
379  fFlowTrack->SetMaxPtsTRD(nClus[7]) ;                // cout << FlowDebug << "nClus = " << nClus[7] << endl ;  
380  fFlowTrack->SetChi2TRD(nHits[7]/nClus[7]) ;         // cout << FlowDebug << "Chi2 = " << nHits[7]/nClus[7] << endl ; 
381  fFlowTrack->SetSigTRD(dEdx[2]) ;                    // cout << FlowDebug << "Dedx = " << dEdx << endl ; 
382  fFlowTrack->SetPatTRD(pAt) ;                        // cout << FlowDebug << "p = " << pAt << endl ;                                    
383  fFlowTrack->SetRespFunTRD(detResFun) ;              // cout << FlowDebug << "response function = " << detResFun << endl ;  
384  // -
385  fFlowTrack->SetNhitsTOF(nHits[8]) ;                 // cout << FlowDebug << "nHits TOF = " << nHits[8] << endl ;                                                    
386  fFlowTrack->SetMaxPtsTOF(nClus[8]) ;                // cout << FlowDebug << "nClus = " << nClus[8] << endl ; 
387  fFlowTrack->SetChi2TOF(nHits[8]/nClus[8]) ;         // cout << FlowDebug << "Chi2 = " << nHits[8]/nClus[8] << endl ; 
388  fFlowTrack->SetTofTOF(dEdx[3]) ;                    // cout << FlowDebug << "Dedx = " << 1. << endl ; 
389  fFlowTrack->SetPatTOF(pAt) ;                        // cout << FlowDebug << "p = " << pAt << endl ;                                    
390  fFlowTrack->SetRespFunTOF(detResFun) ;              // cout << FlowDebug << "response function = " << detResFun << endl ;  
391  
392  return fFlowTrack ; 
393 }
394 //-----------------------------------------------------------------------
395 AliFlowV0* AliFlowKineMaker::FillFlowV0(TParticle* fParticle)
396 {
397  // From a neutral TParticle (input) fills the AliFlowV0 (output) .
398
399  TString name = "" ; name += fPartNumber ;
400  Int_t idx = fFlowEvent->V0Collection()->GetEntries() ;
401  fFlowV0 = (AliFlowV0*)(fFlowEvent->V0Collection()->New(idx)) ;
402  fFlowV0->SetName(name.Data()) ;
403  
404  // cout << " -v0- " << name.Data() << "(" << idx << ")"  << endl ;
405
406  // TParticle label (link: KineTree-ESD)
407  Int_t label = TMath::Abs(fPartNumber);
408  fFlowV0->SetLabel(label) ;
409
410  // reconstructed position of the V0 
411  Double_t xyz[3] ;              
412  xyz[0] = fParticle->Vx() ; 
413  xyz[1] = fParticle->Vy() ; 
414  xyz[2] = fParticle->Vz() ; 
415  fFlowV0->SetCrossPoint(xyz[0],xyz[1],xyz[2]) ;
416
417  // V0's impact parameter & error (chi2 , DCA , sigma , pointing angle)  
418  fFlowV0->SetDca((Float_t)Norm(xyz)) ;
419  fFlowV0->SetSigma(0.) ;
420  fFlowV0->SetCosPointingAngle(1.) ;  
421  fFlowV0->SetDaughtersDca(0.) ;
422  fFlowV0->SetChi2(0.) ;
423
424  // reconstructed momentum of the V0
425  Double_t pxyz[3] ;
426  pxyz[0] = fParticle->Px() ; pxyz[1] = fParticle->Py() ; pxyz[2] = fParticle->Pz() ;                    
427  Float_t phi = (Float_t)Phi(pxyz) ; if(phi<0) { phi += 2*TMath::Pi() ; }
428  fFlowV0->SetPhi(phi) ;         
429  Float_t pt = (Float_t)Pt(pxyz) ; 
430  fFlowV0->SetPt(pt) ;           
431  Float_t eta = (Float_t)Eta(pxyz) ; 
432  fFlowV0->SetEta(eta) ;         
433
434  // P.id. 
435  Int_t pdgCode = fParticle->GetPdgCode() ;
436  fFlowV0->SetMostLikelihoodPID(pdgCode);        
437
438  // mass 
439  fFlowV0->SetVmass((Float_t)fParticle->GetMass()) ; 
440  
441  // daughters (should be taken on puprose from the KineTree, and wrote into the flow event)
442  Int_t nDaughters = fParticle->GetNDaughters() ;
443  Int_t d1 = fParticle->GetDaughter(nDaughters-1) ;
444  Int_t d2 = fParticle->GetDaughter(nDaughters-2) ;
445 //
446 // AliFlowTrack* pos = (AliFlowTrack*)fFlowEvent->TrackCollection()->At(pN) ; 
447 // AliFlowTrack* neg = (AliFlowTrack*)fFlowEvent->TrackCollection()->At(nN) ; 
448 // fFlowV0->SetDaughters(pos,neg) ;
449 //
450  // d1 + d2 ; // warning: statement with no effect :)
451
452  return fFlowV0 ;
453 }
454 //-----------------------------------------------------------------------
455 Bool_t AliFlowKineMaker::CheckTrack(TParticle* fParticle) const
456 {
457  // applies track cuts (pE , eta , label)
458
459  Float_t eta   = (Float_t)fParticle->Eta() ;
460  Float_t pE    = (Float_t)fParticle->P() ;
461  Int_t   label = -1 ;   // check how to assign this !
462  Bool_t  prim  = fParticle->IsPrimary() ;
463  
464  if(fAbsEta && (eta > fAbsEta))                                          { return kFALSE ; }
465  if((fElow < fEup) && ((pE<fElow) || (pE>fEup)))                         { return kFALSE ; }
466  if((fLabel[0] < fLabel[1]) && ((label<fLabel[0]) || (label>fLabel[1]))) { return kFALSE ; }
467  if(fPrimary && !prim)                                                   { return kFALSE ; }
468
469  return kTRUE ;
470 }
471 //-----------------------------------------------------------------------
472 Bool_t AliFlowKineMaker::CheckEvent(TTree* fKTree) const
473 {
474  // applies event cuts (dummy)
475  
476  if(!fKTree) { return kFALSE ; }
477
478  return kTRUE ;
479 }
480 //-----------------------------------------------------------------------
481 void AliFlowKineMaker::PrintCutList()
482
483  // Prints the list of Cuts 
484
485  cout << " * ESD cuts list * " << endl ; 
486  if(fLabel[0]<fLabel[1])
487  { 
488   cout << "  Track Label [ " << fLabel[0] << " , " << fLabel[1] << " ] " << endl ; 
489  }
490  if(fAbsEta)
491  { 
492   cout << "  |eta| < " << fAbsEta << endl ; 
493  }
494  if(fElow<fEup)
495  { 
496   cout << "  Track Energy (P_total) [ " << fElow << " , " << fEup << " ] " << endl ; 
497  }
498  if(fPrimary)                                                    
499  { 
500   cout << "  Primary Particles " << endl ;
501  }
502  cout << " *               * " << endl ;
503 }
504 //-----------------------------------------------------------------------
505
506 //-----------------------------------------------------------------------
507 //*** USEFULL METHODS for a 3-array of double (~ TVector3) ***
508 //-----------------------------------------------------------------------
509 Double_t AliFlowKineMaker::Norm(Double_t nu[3])
510
511  // returns the norm of a double[3] 
512
513  Double_t norm2 = nu[0]*nu[0] + nu[1]*nu[1] + nu[2]*nu[2] ;
514  return TMath::Sqrt(norm2) ; 
515 }
516 //-----------------------------------------------------------------------
517 Double_t AliFlowKineMaker::Phi(Double_t nu[3])
518 {
519  // returns the azimuthal angle of a double[3] 
520
521  if(nu[0]==0 && nu[1]==0) { return 0. ; }
522  else                     { return TMath::ATan2(nu[1],nu[0]) ; }
523 }
524 //-----------------------------------------------------------------------
525 Double_t AliFlowKineMaker::Pt(Double_t nu[3])
526 {
527  // returns the transvers momentum of a double[3] 
528
529  Double_t trans = nu[0]*nu[0] + nu[1]*nu[1] ;
530  return TMath::Sqrt(trans) ; 
531 }
532 //-----------------------------------------------------------------------
533 Double_t AliFlowKineMaker::Eta(Double_t nu[3])
534 {
535  // returns the PseudoRapidity of a double[3] 
536  // if transvers momentum = 0 --> returns +/- 1.000
537
538  Double_t m = Norm(nu) ;
539  if(nu[0]!=0 || nu[1]!=0) { return 0.5*TMath::Log((m+nu[2])/(m-nu[2])) ; }
540  else                     { return TMath::Sign((Double_t)1000.,nu[2]) ; }
541 }
542 //-----------------------------------------------------------------------
543
544 #endif
545
546 //////////////////////////////////////////////////////////////////////
547 // - one way to open the alice KineTree -
548 //////////////////////////////////////////////////////////////////////
549 //
550 //   TString fileName = "galice.root" ;
551 //   rl = AliRunLoader::Open(fileName.Data(),"MyEvent","read");
552 //   rl->LoadgAlice();
553 //   gAlice = rl->GetAliRun();
554 //   rl->LoadHeader();
555 //   rl->LoadKinematics();
556 //   fNumberOfEvents = rl->GetNumberOfEvents() ;
557 //
558 //   Int_t exitStatus = rl->GetEvent(getEv) ; if(exitStatus!=0) { return kFALSE ; }
559 //
560 //   TTree* pKTree = (TTree*)rl->TreeK();         // Particles' TTree (KineTree)
561 //   AliStack* pStack = gAlice->Stack();          // Particles' Stack - "Label()" to get the number in the stack 
562 //
563 // // else if(rl)      // opens files one by one (unload and reload)
564 // // {
565 // //  rl->UnloadgAlice() ;
566 // //  rl->UnloadHeader() ;
567 // //  rl->UnloadKinematics() ;
568 // //  delete rl ; rl = 0 ;
569 // // }
570 //
571 //   fNumberOfParticles = pKTree->GetEntries() ;
572 //   nPart = pStack->GetNtrack() ;
573 //   cout << " Event : " << evtN << " :  particles : " << fNumberOfParticles << "  (stack: " << nPart << ") . " << endl ; }
574 //
575 //////////////////////////////////////////////////////////////////////
576