PWG/FLOW/Tasks/AliFlowEvent.cxx

   1 /**************************************************************************\r
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
   3  *                                                                        *\r
   4  * Author: The ALICE Off-line Project.                                    *\r
   5  * Contributors are mentioned in the code where appropriate.              *\r
   6  *                                                                        *\r
   7  * Permission to use, copy, modify and distribute this software and its   *\r
   8  * documentation strictly for non-commercial purposes is hereby granted   *\r
   9  * without fee, provided that the above copyright notice appears in all   *\r
  10  * copies and that both the copyright notice and this permission notice   *\r
  11  * appear in the supporting documentation. The authors make no claims     *\r
  12  * about the suitability of this software for any purpose. It is          *\r
  13  * provided "as is" without express or implied warranty.                  *\r
  14  **************************************************************************/\r
  15 \r
  16 /*****************************************************************\r
  17   AliFlowEvent: Event container for flow analysis\r
  18 \r
  19   origin:   Mikolaj Krzewicki  (mikolaj.krzewicki@cern.ch)\r
  20 *****************************************************************/\r
  21 \r
  22 #include "Riostream.h"\r
  23 #include "TList.h"\r
  24 #include "TH2F.h"\r
  25 #include "AliMCEvent.h"\r
  26 #include "AliMCParticle.h"\r
  27 #include "AliCFManager.h"\r
  28 #include "AliESDtrack.h"\r
  29 #include "AliESDPmdTrack.h"\r
  30 #include "AliESDEvent.h"\r
  31 #include "AliAODEvent.h"\r
  32 #include "AliGenCocktailEventHeader.h"\r
  33 #include "AliGenEposEventHeader.h"\r
  34 #include "AliGenHijingEventHeader.h"\r
  35 #include "AliGenGeVSimEventHeader.h"\r
  36 #include "AliCollisionGeometry.h"\r
  37 #include "AliMultiplicity.h"\r
  38 #include "AliFlowTrackCuts.h"\r
  39 #include "AliFlowEventSimple.h"\r
  40 #include "AliFlowTrack.h"\r
  41 #include "AliFlowEvent.h"\r
  42 #include "AliLog.h"\r
  43 \r
  44 using std::endl;\r
  45 using std::cout;\r
  46 ClassImp(AliFlowEvent)\r
  47 \r
  48 //-----------------------------------------------------------------------\r
  49 AliFlowEvent::AliFlowEvent():\r
  50   AliFlowEventSimple()\r
  51 {\r
  52   //ctor\r
  53   cout << "AliFlowEvent: Default constructor to be used only by root for io" << endl;\r
  54 }\r
  55 \r
  56 //-----------------------------------------------------------------------\r
  57 AliFlowEvent::AliFlowEvent(Int_t n):\r
  58   AliFlowEventSimple(n)\r
  59 {\r
  60   //ctor\r
  61 }\r
  62 \r
  63 //-----------------------------------------------------------------------\r
  64 AliFlowEvent::AliFlowEvent(const AliFlowEvent& event):\r
  65   AliFlowEventSimple(event)\r
  66 {\r
  67   //cpy ctor\r
  68 }\r
  69 \r
  70 //-----------------------------------------------------------------------\r
  71 AliFlowEvent& AliFlowEvent::operator=(const AliFlowEvent& event)\r
  72 {\r
  73   //assignment operator\r
  74   AliFlowEventSimple::operator=(event);\r
  75   return *this;\r
  76 }\r
  77 \r
  78 //-----------------------------------------------------------------------\r
  79 AliFlowTrack* AliFlowEvent::GetTrack(Int_t i)\r
  80 {\r
  81   //get track i from collection\r
  82   if (i>=fNumberOfTracks) return NULL;\r
  83   AliFlowTrack* pTrack = static_cast<AliFlowTrack*>(fTrackCollection->At(i)) ;\r
  84   return pTrack;\r
  85 }\r
  86 \r
  87 //-----------------------------------------------------------------------\r
  88 void AliFlowEvent::SetMCReactionPlaneAngle(const AliMCEvent* mcEvent)\r
  89 {\r
  90   //sets the event plane angle from the proper header in the MC\r
  91 \r
  92   //COCKTAIL with HIJING\r
  93   if (!strcmp(mcEvent-> GenEventHeader()->GetName(),"Cocktail Header"))   //returns 0 if matches\r
  94   {\r
  95     AliGenCocktailEventHeader *headerC = dynamic_cast<AliGenCocktailEventHeader *> (mcEvent-> GenEventHeader());\r
  96     if (headerC)\r
  97     {\r
  98       TList *lhd = headerC->GetHeaders();\r
  99       if (lhd)\r
 100       {\r
 101         AliGenHijingEventHeader *hdh = dynamic_cast<AliGenHijingEventHeader *> (lhd->At(0));\r
 102         if (hdh) AliFlowEventSimple::SetMCReactionPlaneAngle( hdh->ReactionPlaneAngle() );\r
 103       }\r
 104     }\r
 105   }\r
 106   //THERMINATOR\r
 107   else if (!strcmp(mcEvent-> GenEventHeader()->GetName(),"Therminator"))   //returns 0 if matches\r
 108   {\r
 109     AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());\r
 110     if (headerH) AliFlowEventSimple::SetMCReactionPlaneAngle( headerH->ReactionPlaneAngle() );\r
 111   }\r
 112   //GEVSIM\r
 113   else if (!strcmp(mcEvent-> GenEventHeader()->GetName(),"GeVSim header"))   //returns 0 if matches\r
 114   {\r
 115     AliGenGeVSimEventHeader* headerG = dynamic_cast<AliGenGeVSimEventHeader*>(mcEvent->GenEventHeader());\r
 116     if (headerG) AliFlowEventSimple::SetMCReactionPlaneAngle( headerG->GetEventPlane() );\r
 117   }\r
 118   //HIJING\r
 119   else if (!strcmp(mcEvent-> GenEventHeader()->GetName(),"Hijing"))   //returns 0 if matches\r
 120   {\r
 121     AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());\r
 122     if (headerH) AliFlowEventSimple::SetMCReactionPlaneAngle( headerH->ReactionPlaneAngle() );\r
 123   }\r
 124   //AMPT\r
 125   else if (!strcmp(mcEvent-> GenEventHeader()->GetName(),"Ampt"))   //returns 0 if matches\r
 126   {\r
 127     AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());\r
 128     if (headerH) AliFlowEventSimple::SetMCReactionPlaneAngle( headerH->ReactionPlaneAngle() );\r
 129   }\r
 130   //EPOS\r
 131   else if (!strcmp(mcEvent->GenEventHeader()->GetName(),"EPOS"))\r
 132   {\r
 133     AliGenEposEventHeader* headerE = dynamic_cast<AliGenEposEventHeader*>(mcEvent->GenEventHeader());\r
 134     if (headerE) AliFlowEventSimple::SetMCReactionPlaneAngle( headerE->ReactionPlaneAngle() );\r
 135   }\r
 136   //Hydjet\r
 137   else\r
 138   {\r
 139     AliCollisionGeometry* header = dynamic_cast<AliCollisionGeometry*>(mcEvent->GenEventHeader());\r
 140     if (header) AliFlowEventSimple::SetMCReactionPlaneAngle( header->ReactionPlaneAngle() );\r
 141   }\r
 142 }\r
 143 \r
 144 //-----------------------------------------------------------------------\r
 145 AliFlowEvent::AliFlowEvent( const AliMCEvent* anInput,\r
 146                             const AliCFManager* rpCFManager,\r
 147                             const AliCFManager* poiCFManager):\r
 148   AliFlowEventSimple(20)\r
 149 {\r
 150   //Fills the event from the MC kinematic information\r
 151 \r
 152   Int_t iNumberOfInputTracks = anInput->GetNumberOfTracks() ;\r
 153 \r
 154   //loop over tracks\r
 155   for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)\r
 156   {\r
 157     //get input particle\r
 158     AliMCParticle* pParticle = dynamic_cast<AliMCParticle*>(anInput->GetTrack(itrkN));\r
 159     if (!pParticle) continue;\r
 160 \r
 161     //check if pParticle passes the cuts\r
 162     Bool_t rpOK = kTRUE;\r
 163     Bool_t poiOK = kTRUE;\r
 164     if (rpCFManager && poiCFManager)\r
 165     {\r
 166       rpOK = rpCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,pParticle);\r
 167       poiOK = poiCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,pParticle);\r
 168     }\r
 169     if (!(rpOK||poiOK)) continue;\r
 170 \r
 171     AliFlowTrack* pTrack = new AliFlowTrack(pParticle);\r
 172     pTrack->SetSource(AliFlowTrack::kFromMC);\r
 173 \r
 174     if (rpOK && rpCFManager)\r
 175     {\r
 176       pTrack->SetForRPSelection(kTRUE);\r
 177       fNumberOfRPs++;\r
 178     }\r
 179     if (poiOK && poiCFManager)\r
 180     {\r
 181       pTrack->SetForPOISelection(kTRUE);\r
 182       fNumberOfPOIs++;\r
 183     }\r
 184 \r
 185     AddTrack(pTrack) ;\r
 186   }//for all tracks\r
 187   SetMCReactionPlaneAngle(anInput);\r
 188 }\r
 189 \r
 190 //-----------------------------------------------------------------------\r
 191 AliFlowEvent::AliFlowEvent( const AliESDEvent* anInput,\r
 192                             const AliCFManager* rpCFManager,\r
 193                             const AliCFManager* poiCFManager ):\r
 194   AliFlowEventSimple(20)\r
 195 {\r
 196   //Fills the event from the ESD\r
 197 \r
 198   Int_t iNumberOfInputTracks = anInput->GetNumberOfTracks() ;\r
 199 \r
 200   //loop over tracks\r
 201   for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)\r
 202   {\r
 203     AliESDtrack* pParticle = anInput->GetTrack(itrkN);   //get input particle\r
 204 \r
 205     //check if pParticle passes the cuts\r
 206     Bool_t rpOK = kTRUE;\r
 207     Bool_t poiOK = kTRUE;\r
 208     if (rpCFManager && poiCFManager)\r
 209     {\r
 210       rpOK = ( rpCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&\r
 211                rpCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));\r
 212       poiOK = ( poiCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&\r
 213                 poiCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));\r
 214     }\r
 215     if (!(rpOK || poiOK)) continue;\r
 216 \r
 217     //make new AliFLowTrack\r
 218     AliFlowTrack* pTrack = new AliFlowTrack(pParticle);\r
 219     pTrack->SetSource(AliFlowTrack::kFromESD);\r
 220 \r
 221     //marking the particles used for int. flow:\r
 222     if(rpOK && rpCFManager)\r
 223     {\r
 224       pTrack->SetForRPSelection(kTRUE);\r
 225       fNumberOfRPs++;\r
 226     }\r
 227     //marking the particles used for diff. flow:\r
 228     if(poiOK && poiCFManager)\r
 229     {\r
 230       pTrack->SetForPOISelection(kTRUE);\r
 231       fNumberOfPOIs++;\r
 232     }\r
 233 \r
 234     AddTrack(pTrack);\r
 235   }//end of while (itrkN < iNumberOfInputTracks)\r
 236 }\r
 237 \r
 238 //-----------------------------------------------------------------------\r
 239 AliFlowEvent::AliFlowEvent( const AliAODEvent* anInput,\r
 240                             const AliCFManager* rpCFManager,\r
 241                             const AliCFManager* poiCFManager):\r
 242   AliFlowEventSimple(20)\r
 243 {\r
 244   //Fills the event from the AOD\r
 245   Int_t iNumberOfInputTracks = anInput->GetNumberOfTracks() ;\r
 246 \r
 247   //loop over tracks\r
 248   for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)\r
 249   {\r
 250     AliAODTrack* pParticle = anInput->GetTrack(itrkN);   //get input particle\r
 251 \r
 252     //check if pParticle passes the cuts\r
 253     Bool_t rpOK = kTRUE;\r
 254     Bool_t poiOK = kTRUE;\r
 255     if (rpCFManager && poiCFManager)\r
 256     {\r
 257       rpOK = ( rpCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&\r
 258                rpCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));\r
 259       poiOK = ( poiCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&\r
 260                 poiCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));\r
 261     }\r
 262     if (!(rpOK || poiOK)) continue;\r
 263 \r
 264     //make new AliFlowTrack\r
 265     AliFlowTrack* pTrack = new AliFlowTrack(pParticle);\r
 266     pTrack->SetSource(AliFlowTrack::kFromAOD);\r
 267 \r
 268     if (rpOK /* && rpCFManager */ ) // to be fixed - with CF managers uncommented only empty events (NULL in header files)\r
 269     {\r
 270       pTrack->SetForRPSelection(kTRUE);\r
 271       fNumberOfRPs++;\r
 272     }\r
 273     if (poiOK /* && poiCFManager*/ )\r
 274     {\r
 275       pTrack->SetForPOISelection(kTRUE);\r
 276       fNumberOfPOIs++;\r
 277     }\r
 278     AddTrack(pTrack);\r
 279   }\r
 280 \r
 281   //  if (iSelParticlesRP >= fMinMult && iSelParticlesRP <= fMaxMult)\r
 282   //  {\r
 283   //    if ( (++fCount % 100) == 0)\r
 284   //    {\r
 285   //      if (!fMCReactionPlaneAngle == 0) cout<<" MC Reaction Plane Angle = "<<  fMCReactionPlaneAngle << endl;\r
 286   //      else cout<<" MC Reaction Plane Angle = unknown "<< endl;\r
 287   //      cout<<" iGoodTracks = "<<iGoodTracks<<endl;\r
 288   //      cout<<" # of RP selected tracks = "<<iSelParticlesRP<<endl;\r
 289   //      cout<<" # of POI selected tracks = "<<iSelParticlesPOI<<endl;\r
 290   //      cout << "# " << fCount << " events processed" << endl;\r
 291   //    }\r
 292   //    return pEvent;\r
 293   //  }\r
 294   //  else\r
 295   //  {\r
 296   //    cout<<"Not enough tracks in the FlowEventSimple"<<endl;\r
 297   //    return 0;\r
 298   //  }\r
 299   //}\r
 300   //else\r
 301   //{\r
 302   //  cout<<"Event does not pass multiplicity cuts"<<endl;\r
 303   //  return 0;\r
 304   //}\r
 305 \r
 306 }\r
 307 \r
 308 //-----------------------------------------------------------------------\r
 309 AliFlowEvent::AliFlowEvent( const AliESDEvent* anInput,\r
 310                             const AliMCEvent* anInputMc,\r
 311                             KineSource anOption,\r
 312                             const AliCFManager* rpCFManager,\r
 313                             const AliCFManager* poiCFManager ):\r
 314   AliFlowEventSimple(20)\r
 315 {\r
 316   //fills the event with tracks from the ESD and kinematics from the MC info via the track label\r
 317   if (anOption==kNoKine)\r
 318   {\r
 319     AliFatal("WRONG OPTION IN AliFlowEventMaker::FillTracks(AliESDEvent* anInput, AliMCEvent* anInputMc, KineSource anOption)");\r
 320     exit(1);\r
 321   }\r
 322 \r
 323   Int_t iNumberOfInputTracks = anInput->GetNumberOfTracks() ;\r
 324 \r
 325   Int_t iNumberOfInputTracksMC = anInputMc->GetNumberOfTracks() ;\r
 326   if (iNumberOfInputTracksMC==-1)\r
 327   {\r
 328     AliError("Skipping Event -- No MC information available for this event");\r
 329     return;\r
 330   }\r
 331 \r
 332   //loop over ESD tracks\r
 333   for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)\r
 334   {\r
 335     AliESDtrack* pParticle = anInput->GetTrack(itrkN);   //get input particle\r
 336     //get Label\r
 337     Int_t iLabel = pParticle->GetLabel();\r
 338     //match to mc particle\r
 339     AliMCParticle* pMcParticle = (AliMCParticle*) anInputMc->GetTrack(TMath::Abs(iLabel));\r
 340 \r
 341     //check\r
 342     if (TMath::Abs(pParticle->GetLabel())!=pMcParticle->Label())\r
 343       AliWarning(Form("pParticle->GetLabel()!=pMcParticle->Label(), %i, %i", pParticle->GetLabel(), pMcParticle->Label()));\r
 344 \r
 345     //check if pParticle passes the cuts\r
 346     Bool_t rpOK = kFALSE;\r
 347     Bool_t poiOK = kFALSE;\r
 348     if (rpCFManager && poiCFManager)\r
 349     {\r
 350       if(anOption == kESDkine)\r
 351       {\r
 352         if (rpCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,pMcParticle,"mcGenCuts1") &&\r
 353             rpCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle))\r
 354           rpOK=kTRUE;\r
 355         if (poiCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,pMcParticle,"mcGenCuts2") &&\r
 356             poiCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle))\r
 357           poiOK=kTRUE;\r
 358       }\r
 359       else if (anOption == kMCkine)\r
 360       {\r
 361         if (rpCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,pMcParticle))\r
 362           rpOK=kTRUE;\r
 363         if (poiCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,pMcParticle))\r
 364           poiOK=kTRUE;\r
 365       }\r
 366     }\r
 367 \r
 368     if (!(rpOK || poiOK)) continue;\r
 369 \r
 370     //make new AliFlowTrack\r
 371     AliFlowTrack* pTrack = NULL;\r
 372     if(anOption == kESDkine)   //take the PID from the MC & the kinematics from the ESD\r
 373     {\r
 374       pTrack = new AliFlowTrack(pParticle);\r
 375     }\r
 376     else if (anOption == kMCkine)   //take the PID and kinematics from the MC\r
 377     {\r
 378       pTrack = new AliFlowTrack(pMcParticle);\r
 379     }\r
 380 \r
 381     if (rpOK && rpCFManager)\r
 382     {\r
 383       fNumberOfRPs++;\r
 384       pTrack->SetForRPSelection();\r
 385     }\r
 386     if (poiOK && poiCFManager) \r
 387     { \r
 388       fNumberOfPOIs++;\r
 389       pTrack->SetForPOISelection();\r
 390     }\r
 391 \r
 392     AddTrack(pTrack);\r
 393   }\r
 394   SetMCReactionPlaneAngle(anInputMc);\r
 395 }\r
 396 \r
 397 //-----------------------------------------------------------------------\r
 398 AliFlowEvent::AliFlowEvent( const AliESDEvent* anInput,\r
 399                             const AliMultiplicity* anInputTracklets,\r
 400                             const AliCFManager* poiCFManager ):\r
 401   AliFlowEventSimple(20)\r
 402 {\r
 403 \r
 404   //Select the particles of interest from the ESD\r
 405   Int_t iNumberOfInputTracks = anInput->GetNumberOfTracks() ;\r
 406 \r
 407   //loop over tracks\r
 408   for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)\r
 409     {\r
 410       AliESDtrack* pParticle = anInput->GetTrack(itrkN);   //get input particle\r
 411 \r
 412       //check if pParticle passes the cuts\r
 413       Bool_t poiOK = kTRUE;\r
 414       if (poiCFManager)\r
 415         {\r
 416           poiOK = ( poiCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&\r
 417                     poiCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));\r
 418         }\r
 419       if (!poiOK) continue;\r
 420       \r
 421       //make new AliFLowTrack\r
 422       AliFlowTrack* pTrack = new AliFlowTrack(pParticle);\r
 423           \r
 424       //marking the particles used for the particle of interest (POI) selection:\r
 425       if(poiOK && poiCFManager)\r
 426         {\r
 427           fNumberOfPOIs++;\r
 428           pTrack->SetForPOISelection(kTRUE);\r
 429           pTrack->SetSource(AliFlowTrack::kFromESD);\r
 430         }\r
 431 \r
 432       AddTrack(pTrack);\r
 433     }//end of while (itrkN < iNumberOfInputTracks)\r
 434 \r
 435   //Select the reference particles from the SPD tracklets\r
 436   anInputTracklets = anInput->GetMultiplicity();\r
 437   Int_t multSPD = anInputTracklets->GetNumberOfTracklets();\r
 438   \r
 439   //loop over tracklets\r
 440   for (Int_t itracklet=0; itracklet<multSPD; ++itracklet) {\r
 441     Float_t thetaTr= anInputTracklets->GetTheta(itracklet);\r
 442     Float_t phiTr= anInputTracklets->GetPhi(itracklet);\r
 443     // calculate eta\r
 444     Float_t etaTr = -TMath::Log(TMath::Tan(thetaTr/2.));\r
 445     \r
 446     //make new AliFLowTrackSimple\r
 447     AliFlowTrack* pTrack = new AliFlowTrack();\r
 448     pTrack->SetPt(0.0);\r
 449     pTrack->SetEta(etaTr);\r
 450     pTrack->SetPhi(phiTr);\r
 451     //marking the particles used for the reference particle (RP) selection:\r
 452     fNumberOfRPs++;\r
 453     pTrack->SetForRPSelection(kTRUE);\r
 454     pTrack->SetSource(AliFlowTrack::kFromTracklet);\r
 455 \r
 456     //Add the track to the flowevent\r
 457     AddTrack(pTrack);\r
 458   }\r
 459 \r
 460 }\r
 461 \r
 462 //-----------------------------------------------------------------------\r
 463 AliFlowEvent::AliFlowEvent( const AliESDEvent* esd,\r
 464                             const AliCFManager* poiCFManager,\r
 465                             Bool_t hybrid):\r
 466   AliFlowEventSimple(20)\r
 467 {\r
 468 \r
 469   //Select the particles of interest from the ESD\r
 470   Int_t iNumberOfInputTracks = esd->GetNumberOfTracks() ;\r
 471 \r
 472   //Double_t gPt = 0.0, gP = 0.0;\r
 473   Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};  //The impact parameters and their covariance.\r
 474   Double_t dca3D = 0.0;\r
 475 \r
 476   AliESDtrack trackTPC;\r
 477 \r
 478   //loop over tracks\r
 479   for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)\r
 480     {\r
 481 \r
 482       if (!esd->GetTrack(itrkN)) continue;\r
 483 \r
 484       Bool_t useTPC = kFALSE;\r
 485 \r
 486       AliESDtrack* pParticle = esd->GetTrack(itrkN);   //get input particle\r
 487 \r
 488       //check if pParticle passes the cuts\r
 489       Bool_t poiOK = kTRUE;\r
 490 \r
 491       if (poiCFManager)\r
 492       {\r
 493         poiOK = ( poiCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&\r
 494                   poiCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));\r
 495       }\r
 496 \r
 497       if (!(poiOK)) continue;\r
 498 \r
 499       AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)pParticle->GetTPCInnerParam();\r
 500 \r
 501       if (tpcTrack)\r
 502       {\r
 503 \r
 504 //      gPt = tpcTrack->Pt();\r
 505 //      gP = tpcTrack->P();\r
 506 \r
 507         useTPC = kTRUE;\r
 508 \r
 509         const AliESDVertex *vertexSPD = esd->GetPrimaryVertexSPD();\r
 510         const AliESDVertex *vertexTPC = esd->GetPrimaryVertexTPC();\r
 511 \r
 512         if(hybrid)\r
 513           tpcTrack->PropagateToDCA(vertexSPD,esd->GetMagneticField(),100.,dca,cov);\r
 514         else\r
 515           tpcTrack->PropagateToDCA(vertexTPC,esd->GetMagneticField(),100.,dca,cov);\r
 516 \r
 517         dca3D = TMath::Sqrt(TMath::Power(dca[0],2)+TMath::Power(dca[1],2));\r
 518 \r
 519       }\r
 520 \r
 521       //make new AliFLowTrack\r
 522       AliFlowTrack* pTrack = new AliFlowTrack(pParticle);\r
 523 \r
 524       pTrack->SetSource(AliFlowTrack::kFromESD);\r
 525 \r
 526       //marking the particles used for diff. flow:\r
 527       if(poiOK && poiCFManager)\r
 528       {\r
 529         pTrack->SetForPOISelection(kTRUE);\r
 530         fNumberOfPOIs++;\r
 531       }\r
 532 \r
 533       if(useTPC)\r
 534       {\r
 535         pTrack->SetForRPSelection(kTRUE);\r
 536         fNumberOfRPs++;\r
 537       }\r
 538 \r
 539       AddTrack(pTrack);\r
 540 \r
 541     }//end of while (itrkN < iNumberOfInputTracks)\r
 542 \r
 543 }\r
 544 \r
 545 //-----------------------------------------------------------------------\r
 546 AliFlowEvent::AliFlowEvent( const AliESDEvent* anInput,\r
 547                             const TH2F* anInputFMDhist,\r
 548                             const AliCFManager* poiCFManager ):\r
 549   AliFlowEventSimple(20)\r
 550 {\r
 551 \r
 552   //Select the particles of interest from the ESD\r
 553   Int_t iNumberOfInputTracks = anInput->GetNumberOfTracks() ;\r
 554 \r
 555   //loop over tracks\r
 556   for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)\r
 557     {\r
 558       AliESDtrack* pParticle = anInput->GetTrack(itrkN);   //get input particle\r
 559 \r
 560       //check if pParticle passes the cuts\r
 561       Bool_t poiOK = kTRUE;\r
 562       if (poiCFManager)\r
 563         {\r
 564           poiOK = ( poiCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&\r
 565                     poiCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));\r
 566         }\r
 567       if (!poiOK) continue;\r
 568  \r
 569       //make new AliFLowTrack\r
 570       AliFlowTrack* pTrack = new AliFlowTrack(pParticle);\r
 571           \r
 572       //marking the particles used for the particle of interest (POI) selection:\r
 573       if(poiOK && poiCFManager)\r
 574         {\r
 575           fNumberOfPOIs++; \r
 576           pTrack->SetForPOISelection(kTRUE);\r
 577           pTrack->SetSource(AliFlowTrack::kFromESD);\r
 578         }\r
 579 \r
 580       AddTrack(pTrack);\r
 581     }//end of while (itrkN < iNumberOfInputTracks)\r
 582 \r
 583   //Select the reference particles from the FMD hits\r
 584   //loop over FMD histogram\r
 585   Int_t iBinsEta = anInputFMDhist->GetNbinsX();\r
 586   Int_t iBinsPhi = anInputFMDhist->GetNbinsY();\r
 587   \r
 588   for (Int_t iEta = 1; iEta <= iBinsEta; iEta++){\r
 589     Double_t etaFMD = anInputFMDhist->GetXaxis()->GetBinCenter(iEta);\r
 590     for (Int_t iPhi = 1; iPhi <= iBinsPhi; iPhi++){\r
 591       Double_t phiFMD = anInputFMDhist->GetYaxis()->GetBinCenter(iPhi);\r
 592       Double_t weightFMD = anInputFMDhist->GetBinContent(iEta,iPhi);\r
 593     \r
 594       if (weightFMD > 0.0) { //do not add empty bins\r
 595         //make new AliFLowTrackSimple\r
 596         AliFlowTrack* pTrack = new AliFlowTrack();\r
 597         pTrack->SetPt(0.0);\r
 598         pTrack->SetEta(etaFMD);\r
 599         pTrack->SetPhi(phiFMD);\r
 600         pTrack->SetWeight(weightFMD);\r
 601         //marking the particles used for the reference particle (RP) selection:\r
 602         pTrack->TagRP();\r
 603         fNumberOfRPs++;\r
 604         pTrack->SetSource(AliFlowTrack::kFromFMD);\r
 605 \r
 606         //Add the track to the flowevent\r
 607         AddTrack(pTrack);\r
 608         \r
 609       }\r
 610     }\r
 611   }\r
 612 \r
 613 }\r
 614 \r
 615 //-----------------------------------------------------------------------\r
 616 void AliFlowEvent::Fill( AliFlowTrackCuts* rpCuts,\r
 617                          AliFlowTrackCuts* poiCuts )\r
 618 {\r
 619   //Fills the event from a vevent: AliESDEvent,AliAODEvent,AliMCEvent\r
 620   //the input data needs to be attached to the cuts\r
 621   //we have two cases, if we're cutting the same collection of tracks\r
 622   //(same param type) then we can have tracks that are both rp and poi\r
 623   //in the other case we want to have two exclusive sets of rps and pois\r
 624   //e.g. one tracklets, the other PMD or global - USER IS RESPOSIBLE\r
 625   //FOR MAKING SURE THEY DONT OVERLAP OR ELSE THE SAME PARTICLE WILL BE\r
 626   //TAKEN TWICE\r
 627 \r
 628   ClearFast();\r
 629 \r
 630   if (!rpCuts || !poiCuts) return;\r
 631   AliFlowTrackCuts::trackParameterType sourceRP = rpCuts->GetParamType();\r
 632   AliFlowTrackCuts::trackParameterType sourcePOI = poiCuts->GetParamType();\r
 633   AliFlowTrack* pTrack=NULL;\r
 634 \r
 635   if (sourceRP==sourcePOI)\r
 636   {\r
 637     //loop over tracks\r
 638     Int_t numberOfInputObject = rpCuts->GetNumberOfInputObjects();\r
 639     for (Int_t i=0; i<numberOfInputObject; i++)\r
 640     {\r
 641       //get input object (particle)\r
 642       TObject* particle = rpCuts->GetInputObject(i);\r
 643 \r
 644       Bool_t rp = rpCuts->IsSelected(particle,i);\r
 645       Bool_t poi = poiCuts->IsSelected(particle,i);\r
 646       \r
 647       if (!(rp||poi)) continue;\r
 648 \r
 649       //make new AliFLowTrack\r
 650       if (rp)\r
 651       {\r
 652         pTrack = ReuseTrack(fNumberOfTracks);\r
 653         if (!rpCuts->FillFlowTrack(pTrack)) continue;\r
 654         pTrack->TagRP(); fNumberOfRPs++;\r
 655         if (poi) {pTrack->TagPOI(); fNumberOfPOIs++;}\r
 656       }\r
 657       else if (poi)\r
 658       {\r
 659         pTrack = ReuseTrack(fNumberOfTracks);\r
 660         if (!poiCuts->FillFlowTrack(pTrack)) continue;\r
 661         pTrack->TagPOI(); fNumberOfPOIs++;\r
 662       }\r
 663       fNumberOfTracks++;\r
 664     }//end of while (i < numberOfTracks)\r
 665   }\r
 666   else if (sourceRP!=sourcePOI)\r
 667   {\r
 668     //here we have two different sources of particles, so we fill\r
 669     //them independently\r
 670     //RP\r
 671     Int_t numberOfInputObject = rpCuts->GetNumberOfInputObjects();\r
 672     for (Int_t i=0; i<numberOfInputObject; i++)\r
 673       {\r
 674       TObject* particle = rpCuts->GetInputObject(i);\r
 675       Bool_t rp = rpCuts->IsSelected(particle,i);\r
 676       if (!rp) continue;\r
 677       pTrack = ReuseTrack(fNumberOfTracks);\r
 678       if (!rpCuts->FillFlowTrack(pTrack)) continue;\r
 679       pTrack->TagRP();\r
 680       fNumberOfRPs++;\r
 681       fNumberOfTracks++;\r
 682     }\r
 683     //POI\r
 684     numberOfInputObject = poiCuts->GetNumberOfInputObjects();\r
 685     for (Int_t i=0; i<numberOfInputObject; i++)\r
 686     {\r
 687       TObject* particle = poiCuts->GetInputObject(i);\r
 688       Bool_t poi = poiCuts->IsSelected(particle,i);\r
 689       if (!poi) continue;\r
 690       pTrack = ReuseTrack(fNumberOfTracks);\r
 691       if (!poiCuts->FillFlowTrack(pTrack)) continue;\r
 692       pTrack->TagPOI();\r
 693       fNumberOfPOIs++;\r
 694       fNumberOfTracks++;\r
 695     }\r
 696   }\r
 697 }\r
 698 \r
 699 //-----------------------------------------------------------------------\r
 700 void AliFlowEvent::InsertTrack(AliFlowTrack *thisTrack) {\r
 701   // adds a flow track at the end of the container\r
 702   AliFlowTrack *pTrack = ReuseTrack( fNumberOfTracks++ );\r
 703   pTrack->SetPt( thisTrack->Pt() );\r
 704   pTrack->SetPhi( thisTrack->Phi() );\r
 705   pTrack->SetEta( thisTrack->Eta() );\r
 706   pTrack->SetWeight( thisTrack->Weight() );\r
 707   pTrack->SetCharge( thisTrack->Charge() );\r
 708   pTrack->SetMass( thisTrack->Mass() );\r
 709   pTrack->SetForRPSelection( thisTrack->InRPSelection() );\r
 710   pTrack->SetForPOISelection( thisTrack->InPOISelection() );\r
 711   if(thisTrack->InSubevent(0)) pTrack->SetForSubevent(0);\r
 712   if(thisTrack->InSubevent(1)) pTrack->SetForSubevent(1);\r
 713   pTrack->SetID( thisTrack->GetID() );\r
 714   return;\r
 715 }\r
 716 \r
 717 //-----------------------------------------------------------------------\r
 718 AliFlowTrack* AliFlowEvent::ReuseTrack(Int_t i)\r
 719 {\r
 720   //try to reuse an existing track, if empty, make new one\r
 721   AliFlowTrack* pTrack = static_cast<AliFlowTrack*>(fTrackCollection->At(i));\r
 722   if (pTrack)\r
 723   {\r
 724     pTrack->Clear();\r
 725   }\r
 726   else \r
 727   {\r
 728     pTrack = new AliFlowTrack();\r
 729     fTrackCollection->AddAtAndExpand(pTrack,i);\r
 730   }\r
 731   return pTrack;\r
 732 }\r
 733 \r
 734 //-----------------------------------------------------------------------\r
 735 AliFlowEvent::AliFlowEvent( AliFlowTrackCuts* rpCuts,\r
 736                             AliFlowTrackCuts* poiCuts ):\r
 737   AliFlowEventSimple(20)\r
 738 {\r
 739   //Fills the event from a vevent: AliESDEvent,AliAODEvent,AliMCEvent\r
 740   //the input data needs to be attached to the cuts\r
 741   //we have two cases, if we're cutting the same collection of tracks\r
 742   //(same param type) then we can have tracks that are both rp and poi\r
 743   //in the other case we want to have two exclusive sets of rps and pois\r
 744   //e.g. one tracklets, the other PMD or global - USER IS RESPOSIBLE\r
 745   //FOR MAKING SURE THEY DONT OVERLAP OR ELSE THE SAME PARTICLE WILL BE\r
 746   //TAKEN TWICE\r
 747 \r
 748   if (!rpCuts || !poiCuts) return;\r
 749   AliFlowTrackCuts::trackParameterType sourceRP = rpCuts->GetParamType();\r
 750   AliFlowTrackCuts::trackParameterType sourcePOI = poiCuts->GetParamType();\r
 751 \r
 752   if (sourceRP==sourcePOI)\r
 753   {\r
 754     //loop over tracks\r
 755     Int_t numberOfInputObject = rpCuts->GetNumberOfInputObjects();\r
 756     for (Int_t i=0; i<numberOfInputObject; i++)\r
 757     {\r
 758       //get input object (particle)\r
 759       TObject* particle = rpCuts->GetInputObject(i);\r
 760 \r
 761       Bool_t rp = rpCuts->IsSelected(particle,i);\r
 762       Bool_t poi = poiCuts->IsSelected(particle,i);\r
 763       \r
 764       if (!(rp||poi)) continue;\r
 765 \r
 766       //make new AliFLowTrack\r
 767       AliFlowTrack* pTrack = NULL;\r
 768       if (rp)\r
 769       {\r
 770         pTrack = rpCuts->MakeFlowTrack();\r
 771         if (!pTrack) continue;\r
 772         pTrack->TagRP(); fNumberOfRPs++;\r
 773         if (poi) {pTrack->TagPOI(); fNumberOfPOIs++;}\r
 774       }\r
 775       else\r
 776       if (poi)\r
 777       {\r
 778         pTrack = poiCuts->MakeFlowTrack();\r
 779         if (!pTrack) continue;\r
 780         pTrack->TagPOI(); fNumberOfPOIs++;\r
 781       }\r
 782       AddTrack(pTrack);\r
 783     }//end of while (i < numberOfTracks)\r
 784   }\r
 785   else if (sourceRP!=sourcePOI)\r
 786   {\r
 787     //here we have two different sources of particles, so we fill\r
 788     //them independently\r
 789     AliFlowTrack* pTrack = NULL;\r
 790     //RP\r
 791     Int_t numberOfInputObject = rpCuts->GetNumberOfInputObjects();\r
 792     for (Int_t i=0; i<numberOfInputObject; i++)\r
 793     {\r
 794       TObject* particle = rpCuts->GetInputObject(i);\r
 795       Bool_t rp = rpCuts->IsSelected(particle,i);\r
 796       if (!rp) continue;\r
 797       pTrack = rpCuts->MakeFlowTrack();\r
 798       if (!pTrack) continue;\r
 799       pTrack->TagRP(); fNumberOfRPs++;\r
 800       AddTrack(pTrack);\r
 801     }\r
 802     //POI\r
 803     numberOfInputObject = poiCuts->GetNumberOfInputObjects();\r
 804     for (Int_t i=0; i<numberOfInputObject; i++)\r
 805     {\r
 806       TObject* particle = poiCuts->GetInputObject(i);\r
 807       Bool_t poi = poiCuts->IsSelected(particle,i);\r
 808       if (!poi) continue;\r
 809       pTrack = poiCuts->MakeFlowTrack();\r
 810       if (!pTrack) continue;\r
 811       pTrack->TagPOI(); fNumberOfPOIs++;\r
 812       AddTrack(pTrack);\r
 813     }\r
 814   }\r
 815 }\r
 816 \r
 817 //-------------------------------------------------------------------//\r
 818 //---- Including PMD tracks as RP --------------------------//\r
 819 \r
 820 AliFlowEvent::AliFlowEvent( const AliESDEvent* anInput,\r
 821                             const AliESDPmdTrack *pmdtracks,\r
 822                             const AliCFManager* poiCFManager ):\r
 823   AliFlowEventSimple(20)\r
 824 {\r
 825   Float_t GetPmdEta(Float_t xPos, Float_t yPos, Float_t zPos);\r
 826   Float_t GetPmdPhi(Float_t xPos, Float_t yPos);\r
 827   //Select the particles of interest from the ESD\r
 828   Int_t iNumberOfInputTracks = anInput->GetNumberOfTracks() ;\r
 829   \r
 830   //loop over tracks\r
 831   for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)\r
 832     {\r
 833       AliESDtrack* pParticle = anInput->GetTrack(itrkN);   //get input particle\r
 834       //check if pParticle passes the cuts\r
 835       Bool_t poiOK = kTRUE;\r
 836       if (poiCFManager)\r
 837         {\r
 838           poiOK = ( poiCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&\r
 839                     poiCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));\r
 840         }\r
 841       if (!poiOK) continue;\r
 842       \r
 843       //make new AliFLowTrack\r
 844       AliFlowTrack* pTrack = new AliFlowTrack(pParticle);\r
 845       \r
 846       //marking the particles used for the particle of interest (POI) selection:\r
 847       if(poiOK && poiCFManager)\r
 848         {\r
 849           fNumberOfPOIs++;\r
 850           pTrack->SetForPOISelection(kTRUE);\r
 851           pTrack->SetSource(AliFlowTrack::kFromESD);\r
 852         }\r
 853       \r
 854       AddTrack(pTrack);\r
 855     }//end of while (itrkN < iNumberOfInputTracks)\r
 856   \r
 857   //Select the reference particles from the PMD tracks\r
 858   Int_t npmdcl = anInput->GetNumberOfPmdTracks();\r
 859   printf("======There are %d PMD tracks in this event\n-------",npmdcl);\r
 860   //loop over clusters \r
 861   for(Int_t iclust=0; iclust < npmdcl; iclust++){\r
 862     //AliESDPmdTrack *pmdtr = anInput->GetPmdTrack(iclust);\r
 863     pmdtracks = anInput->GetPmdTrack(iclust);\r
 864     Int_t   det   = pmdtracks->GetDetector();\r
 865     //Int_t   smn   = pmdtracks->GetSmn();\r
 866     Float_t clsX  = pmdtracks->GetClusterX();\r
 867     Float_t clsY  = pmdtracks->GetClusterY();\r
 868     Float_t clsZ  = pmdtracks->GetClusterZ();\r
 869     Float_t ncell = pmdtracks->GetClusterCells();\r
 870     Float_t adc   = pmdtracks->GetClusterADC();\r
 871     //Float_t pid   = pmdtracks->GetClusterPID();\r
 872     Float_t etacls = GetPmdEta(clsX,clsY,clsZ);\r
 873     Float_t phicls = GetPmdPhi(clsX,clsY);\r
 874     //make new AliFLowTrackSimple\r
 875     AliFlowTrack* pTrack = new AliFlowTrack();\r
 876     //if(det == 0){ //selecting preshower plane only\r
 877     if(det == 0 && adc > 270 && ncell > 1){ //selecting preshower plane only\r
 878       //pTrack->SetPt(adc);//cluster adc\r
 879       pTrack->SetPt(0.0);\r
 880       pTrack->SetEta(etacls);\r
 881       pTrack->SetPhi(phicls);\r
 882       //marking the particles used for the reference particle (RP) selection:\r
 883       fNumberOfRPs++;\r
 884       pTrack->SetForRPSelection(kTRUE);\r
 885       pTrack->SetSource(AliFlowTrack::kFromPMD);\r
 886       //Add the track to the flowevent\r
 887       AddTrack(pTrack);\r
 888     }//if det\r
 889   }\r
 890 }\r
 891 //----------------------------------------------------------------------------//\r
 892 Float_t GetPmdEta(Float_t xPos, Float_t yPos, Float_t zPos)\r
 893 {\r
 894   Float_t rpxpy, theta, eta;\r
 895   rpxpy  = TMath::Sqrt(xPos*xPos + yPos*yPos);\r
 896   theta  = TMath::ATan2(rpxpy,zPos);\r
 897   eta    = -TMath::Log(TMath::Tan(0.5*theta));\r
 898   return eta;\r
 899 }\r
 900 //--------------------------------------------------------------------------//\r
 901 Float_t GetPmdPhi(Float_t xPos, Float_t yPos)\r
 902 {\r
 903   Float_t pybypx, phi = 0., phi1;\r
 904   if(xPos==0)\r
 905     {\r
 906       if(yPos>0) phi = 90.;\r
 907       if(yPos<0) phi = 270.;\r
 908     }\r
 909   if(xPos != 0)\r
 910     {\r
 911       pybypx = yPos/xPos;\r
 912       if(pybypx < 0) pybypx = - pybypx;\r
 913       phi1 = TMath::ATan(pybypx)*180./3.14159;\r
 914       \r
 915       if(xPos > 0 && yPos > 0) phi = phi1;        // 1st Quadrant\r
 916       if(xPos < 0 && yPos > 0) phi = 180 - phi1;  // 2nd Quadrant\r
 917       if(xPos < 0 && yPos < 0) phi = 180 + phi1;  // 3rd Quadrant\r
 918       if(xPos > 0 && yPos < 0) phi = 360 - phi1;  // 4th Quadrant\r
 919       \r
 920     }\r
 921   phi = phi*3.14159/180.;\r
 922   return   phi;\r
 923 }\r
 924 //---------------------------------------------------------------//\r
 925 \r
 926 \r