Improved version, kinetree included
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowKineMaker.cxx
CommitLineData
4e566f2f 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>
48using namespace std; //required for resolving the 'cout' symbol
49
50ClassImp(AliFlowKineMaker)
51//-----------------------------------------------------------------------
52AliFlowKineMaker::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//-----------------------------------------------------------------------
108AliFlowKineMaker::~AliFlowKineMaker()
109{
110 // default destructor (no actions)
111}
112//-----------------------------------------------------------------------
113
114//-----------------------------------------------------------------------
115AliFlowEvent* 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//----------------------------------------------------------------------
193AliFlowTrack* 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//-----------------------------------------------------------------------
395AliFlowV0* 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//-----------------------------------------------------------------------
455Bool_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//-----------------------------------------------------------------------
472Bool_t AliFlowKineMaker::CheckEvent(TTree* fKTree) const
473{
474 // applies event cuts (dummy)
475
476 if(!fKTree) { return kFALSE ; }
477
478 return kTRUE ;
479}
480//-----------------------------------------------------------------------
481void 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//-----------------------------------------------------------------------
509Double_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//-----------------------------------------------------------------------
517Double_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//-----------------------------------------------------------------------
525Double_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//-----------------------------------------------------------------------
533Double_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