PWG2/FLOW/AliFlowTasks/AliFlowEvent.cxx

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