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   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)
  60 {
  61  // default constructor
  62  // resets counters , sets defaults
  63
  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
  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 //-----------------------------------------------------------------------
  90 AliFlowKineMaker::~AliFlowKineMaker()
  91 {
  92  // default destructor (no actions)
  93 }
  94 //-----------------------------------------------------------------------
  95
  96 //-----------------------------------------------------------------------
  97 AliFlowEvent* 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
 103  fFlowEvent = new AliFlowEvent(10000) ; if(!fFlowEvent) { return 0 ; }
 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 //----------------------------------------------------------------------
 175 AliFlowTrack* 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) ;
 195  fFlowTrack->SetDcaSigned(xy,z) ;     // cout << "DCA :   xy = " << xy << "   z = " << z << endl ;
 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)
 220   if(TMath::Abs(etaGl) < AliFlowConstants::fgEtaMid)  { fGoodTracksEta++ ; }
 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] ;
 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
 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 //-----------------------------------------------------------------------
 377 AliFlowV0* 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)
 424  //Int_t nDaughters = fParticle->GetNDaughters() ;
 425  //Int_t d1 = fParticle->GetDaughter(nDaughters-1) ;
 426  //Int_t d2 = fParticle->GetDaughter(nDaughters-2) ;
 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 //-----------------------------------------------------------------------
 437 Bool_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 //-----------------------------------------------------------------------
 454 Bool_t AliFlowKineMaker::CheckEvent(TTree* fKTree) const
 455 {
 456  // applies event cuts (dummy)
 457
 458  if(!fKTree) { return kFALSE ; }
 459
 460  return kTRUE ;
 461 }
 462 //-----------------------------------------------------------------------
 463 void 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 //-----------------------------------------------------------------------
 491 Double_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 //-----------------------------------------------------------------------
 499 Double_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 //-----------------------------------------------------------------------
 507 Double_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 //-----------------------------------------------------------------------
 515 Double_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