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