PWG/FLOW/Tasks/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   mods:     Redmer A. Bertens (rbertens@cern.ch)
  21 *****************************************************************/
  22
  23 #include "Riostream.h"
  24 #include "TFile.h"
  25 #include "TList.h"
  26 #include "TH1.h"
  27 #include "TH2F.h"
  28 #include "TProfile.h"
  29 #include "AliMCEvent.h"
  30 #include "AliMCParticle.h"
  31 #include "AliCFManager.h"
  32 #include "AliESDtrack.h"
  33 #include "AliESDPmdTrack.h"
  34 #include "AliESDEvent.h"
  35 #include "AliAODEvent.h"
  36 #include "AliOADBContainer.h"
  37 #include "AliGenCocktailEventHeader.h"
  38 #include "AliGenEposEventHeader.h"
  39 #include "AliGenHijingEventHeader.h"
  40 #include "AliGenGeVSimEventHeader.h"
  41 #include "AliCollisionGeometry.h"
  42 #include "AliMultiplicity.h"
  43 #include "AliFlowTrackCuts.h"
  44 #include "AliFlowEventSimple.h"
  45 #include "AliFlowTrack.h"
  46 #include "AliFlowVector.h"
  47 #include "AliFlowEvent.h"
  48 #include "AliLog.h"
  49
  50 using std::endl;
  51 using std::cout;
  52 ClassImp(AliFlowEvent)
  53
  54 //-----------------------------------------------------------------------
  55 AliFlowEvent::AliFlowEvent():
  56   AliFlowEventSimple(), fApplyRecentering(-1), fCachedRun(-1), fCurrentCentrality(-1)
  57 {
  58     // constructor
  59     for(Int_t i(0); i < 9; i++) {
  60         for(Int_t j(0); j < 2; j++) {
  61             for(Int_t k(0); k < 2; k++) {
  62                 fMeanQ[i][j][k] = 0.;
  63                 fWidthQ[i][j][k] = 0.;
  64                 fMeanQv3[i][j][k] = 0.;
  65                 fWidthQv3[i][j][k] = 0.;
  66             }
  67         }
  68     }
  69   //ctor
  70   cout << "AliFlowEvent: Default constructor to be used only by root for io" << endl;
  71 }
  72
  73 //-----------------------------------------------------------------------
  74 AliFlowEvent::AliFlowEvent(Int_t n):
  75   AliFlowEventSimple(n), fApplyRecentering(-1), fCachedRun(-1), fCurrentCentrality(-1)
  76 {
  77     // constructor
  78     for(Int_t i(0); i < 9; i++) {
  79         for(Int_t j(0); j < 2; j++) {
  80             for(Int_t k(0); k < 2; k++) {
  81                 fMeanQ[i][j][k] = 0.;
  82                 fWidthQ[i][j][k] = 0.;
  83                 fMeanQv3[i][j][k] = 0.;
  84                 fWidthQv3[i][j][k] = 0.;
  85             }
  86         }
  87     }
  88 }
  89
  90 //-----------------------------------------------------------------------
  91 AliFlowEvent::AliFlowEvent(const AliFlowEvent& event):
  92   AliFlowEventSimple(event), fApplyRecentering(event.fApplyRecentering), fCachedRun(-1), fCurrentCentrality(-1)
  93 {
  94     // copy constructor
  95     for(Int_t i(0); i < 9; i++) {
  96         for(Int_t j(0); j < 2; j++) {
  97             for(Int_t k(0); k < 2; k++) {
  98                 fMeanQ[i][j][k] = 0.;
  99                 fWidthQ[i][j][k] = 0.;
 100                 fMeanQv3[i][j][k] = 0.;
 101                 fWidthQv3[i][j][k] = 0.;
 102             }
 103         }
 104     }
 105 }
 106
 107 //-----------------------------------------------------------------------
 108 AliFlowEvent& AliFlowEvent::operator=(const AliFlowEvent& event)
 109 {
 110   //assignment operator
 111   if (&event==this) return *this;       // check self-assignment
 112
 113   fApplyRecentering = event.fApplyRecentering;
 114   fCachedRun = event.fCachedRun;
 115   fCurrentCentrality = event.fCurrentCentrality;
 116   for(Int_t i(0); i < 9; i++) {
 117       for(Int_t j(0); j < 2; j++) {
 118           for(Int_t k(0); k < 2; k++) {
 119               fMeanQ[i][j][k] = event.fMeanQ[i][j][k];
 120               fWidthQ[i][j][k] = event.fWidthQ[i][j][k];
 121               fMeanQv3[i][j][k] = event.fMeanQv3[i][j][k];
 122               fWidthQv3[i][j][k] = event.fWidthQv3[i][j][k];
 123           }
 124       }
 125   }
 126   AliFlowEventSimple::operator=(event);
 127   return *this;
 128 }
 129
 130 //-----------------------------------------------------------------------
 131 AliFlowTrack* AliFlowEvent::GetTrack(Int_t i)
 132 {
 133   //get track i from collection
 134   if (i>=fNumberOfTracks) return NULL;
 135   AliFlowTrack* pTrack = static_cast<AliFlowTrack*>(fTrackCollection->At(i)) ;
 136   return pTrack;
 137 }
 138
 139 //-----------------------------------------------------------------------
 140 void AliFlowEvent::SetMCReactionPlaneAngle(const AliMCEvent* mcEvent)
 141 {
 142   //sets the event plane angle from the proper header in the MC
 143
 144   //COCKTAIL with HIJING
 145   if (!strcmp(mcEvent-> GenEventHeader()->GetName(),"Cocktail Header"))   //returns 0 if matches
 146   {
 147     AliGenCocktailEventHeader *headerC = dynamic_cast<AliGenCocktailEventHeader *> (mcEvent-> GenEventHeader());
 148     if (headerC)
 149     {
 150       TList *lhd = headerC->GetHeaders();
 151       if (lhd)
 152       {
 153         AliGenHijingEventHeader *hdh = dynamic_cast<AliGenHijingEventHeader *> (lhd->At(0));
 154         if (hdh) AliFlowEventSimple::SetMCReactionPlaneAngle( hdh->ReactionPlaneAngle() );
 155       }
 156     }
 157   }
 158   //THERMINATOR
 159   else if (!strcmp(mcEvent-> GenEventHeader()->GetName(),"Therminator"))   //returns 0 if matches
 160   {
 161     AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
 162     if (headerH) AliFlowEventSimple::SetMCReactionPlaneAngle( headerH->ReactionPlaneAngle() );
 163   }
 164   //GEVSIM
 165   else if (!strcmp(mcEvent-> GenEventHeader()->GetName(),"GeVSim header"))   //returns 0 if matches
 166   {
 167     AliGenGeVSimEventHeader* headerG = dynamic_cast<AliGenGeVSimEventHeader*>(mcEvent->GenEventHeader());
 168     if (headerG) AliFlowEventSimple::SetMCReactionPlaneAngle( headerG->GetEventPlane() );
 169   }
 170   //HIJING
 171   else if (!strcmp(mcEvent-> GenEventHeader()->GetName(),"Hijing"))   //returns 0 if matches
 172   {
 173     AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
 174     if (headerH) AliFlowEventSimple::SetMCReactionPlaneAngle( headerH->ReactionPlaneAngle() );
 175   }
 176   //AMPT
 177   else if (!strcmp(mcEvent-> GenEventHeader()->GetName(),"Ampt"))   //returns 0 if matches
 178   {
 179     AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
 180     if (headerH) AliFlowEventSimple::SetMCReactionPlaneAngle( headerH->ReactionPlaneAngle() );
 181   }
 182   //EPOS
 183   else if (!strcmp(mcEvent->GenEventHeader()->GetName(),"EPOS"))
 184   {
 185     AliGenEposEventHeader* headerE = dynamic_cast<AliGenEposEventHeader*>(mcEvent->GenEventHeader());
 186     if (headerE) AliFlowEventSimple::SetMCReactionPlaneAngle( headerE->ReactionPlaneAngle() );
 187   }
 188   //Hydjet
 189   else
 190   {
 191     AliCollisionGeometry* header = dynamic_cast<AliCollisionGeometry*>(mcEvent->GenEventHeader());
 192     if (header) AliFlowEventSimple::SetMCReactionPlaneAngle( header->ReactionPlaneAngle() );
 193   }
 194 }
 195
 196 //-----------------------------------------------------------------------
 197 AliFlowEvent::AliFlowEvent( const AliMCEvent* anInput,
 198                             const AliCFManager* rpCFManager,
 199                             const AliCFManager* poiCFManager):
 200   AliFlowEventSimple(20), fApplyRecentering(-1), fCachedRun(-1), fCurrentCentrality(-1)
 201 {
 202     // constructor
 203     for(Int_t i(0); i < 9; i++) {
 204         for(Int_t j(0); j < 2; j++) {
 205             for(Int_t k(0); k < 2; k++) {
 206                 fMeanQ[i][j][k] = 0.;
 207                 fWidthQ[i][j][k] = 0.;
 208                 fMeanQv3[i][j][k] = 0.;
 209                 fWidthQv3[i][j][k] = 0.;
 210             }
 211         }
 212     }
 213   //Fills the event from the MC kinematic information
 214
 215   Int_t iNumberOfInputTracks = anInput->GetNumberOfTracks() ;
 216
 217   //loop over tracks
 218   for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)
 219   {
 220     //get input particle
 221     AliMCParticle* pParticle = dynamic_cast<AliMCParticle*>(anInput->GetTrack(itrkN));
 222     if (!pParticle) continue;
 223
 224     //check if pParticle passes the cuts
 225     Bool_t rpOK = kTRUE;
 226     Bool_t poiOK = kTRUE;
 227     if (rpCFManager && poiCFManager)
 228     {
 229       rpOK = rpCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,pParticle);
 230       poiOK = poiCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,pParticle);
 231     }
 232     if (!(rpOK||poiOK)) continue;
 233
 234     AliFlowTrack* pTrack = new AliFlowTrack(pParticle);
 235     pTrack->SetSource(AliFlowTrack::kFromMC);
 236
 237     if (rpOK && rpCFManager)
 238     {
 239       pTrack->SetForRPSelection(kTRUE);
 240       IncrementNumberOfPOIs(0);
 241     }
 242     if (poiOK && poiCFManager)
 243     {
 244       pTrack->SetForPOISelection(kTRUE);
 245       IncrementNumberOfPOIs(1);
 246     }
 247
 248     AddTrack(pTrack) ;
 249   }//for all tracks
 250   SetMCReactionPlaneAngle(anInput);
 251 }
 252
 253 //-----------------------------------------------------------------------
 254 AliFlowEvent::AliFlowEvent( const AliESDEvent* anInput,
 255                             const AliCFManager* rpCFManager,
 256                             const AliCFManager* poiCFManager ):
 257   AliFlowEventSimple(20),  fApplyRecentering(-1), fCachedRun(-1), fCurrentCentrality(-1)
 258 {
 259     // constructor
 260     for(Int_t i(0); i < 9; i++) {
 261         for(Int_t j(0); j < 2; j++) {
 262             for(Int_t k(0); k < 2; k++) {
 263                 fMeanQ[i][j][k] = 0.;
 264                 fWidthQ[i][j][k] = 0.;
 265                 fMeanQv3[i][j][k] = 0.;
 266                 fWidthQv3[i][j][k] = 0.;
 267             }
 268         }
 269     }
 270
 271   //Fills the event from the ESD
 272
 273   Int_t iNumberOfInputTracks = anInput->GetNumberOfTracks() ;
 274
 275   //loop over tracks
 276   for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)
 277   {
 278     AliESDtrack* pParticle = anInput->GetTrack(itrkN);   //get input particle
 279
 280     //check if pParticle passes the cuts
 281     Bool_t rpOK = kTRUE;
 282     Bool_t poiOK = kTRUE;
 283     if (rpCFManager && poiCFManager)
 284     {
 285       rpOK = ( rpCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&
 286                rpCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));
 287       poiOK = ( poiCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&
 288                 poiCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));
 289     }
 290     if (!(rpOK || poiOK)) continue;
 291
 292     //make new AliFLowTrack
 293     AliFlowTrack* pTrack = new AliFlowTrack(pParticle);
 294     pTrack->SetSource(AliFlowTrack::kFromESD);
 295
 296     //marking the particles used for int. flow:
 297     if(rpOK && rpCFManager)
 298     {
 299       pTrack->SetForRPSelection(kTRUE);
 300       IncrementNumberOfPOIs(0);
 301     }
 302     //marking the particles used for diff. flow:
 303     if(poiOK && poiCFManager)
 304     {
 305       pTrack->SetForPOISelection(kTRUE);
 306       IncrementNumberOfPOIs(1);
 307     }
 308
 309     AddTrack(pTrack);
 310   }//end of while (itrkN < iNumberOfInputTracks)
 311 }
 312
 313 //-----------------------------------------------------------------------
 314 AliFlowEvent::AliFlowEvent( const AliAODEvent* anInput,
 315                             const AliCFManager* rpCFManager,
 316                             const AliCFManager* poiCFManager):
 317   AliFlowEventSimple(20), fApplyRecentering(-1), fCachedRun(-1), fCurrentCentrality(-1)
 318 {
 319     // constructor
 320     for(Int_t i(0); i < 9; i++) {
 321         for(Int_t j(0); j < 2; j++) {
 322             for(Int_t k(0); k < 2; k++) {
 323                 fMeanQ[i][j][k] = 0.;
 324                 fWidthQ[i][j][k] = 0.;
 325                 fMeanQv3[i][j][k] = 0.;
 326                 fWidthQv3[i][j][k] = 0.;
 327             }
 328         }
 329     }
 330   //Fills the event from the AOD
 331   Int_t iNumberOfInputTracks = anInput->GetNumberOfTracks() ;
 332
 333   //loop over tracks
 334   for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)
 335   {
 336     AliAODTrack* pParticle = anInput->GetTrack(itrkN);   //get input particle
 337
 338     //check if pParticle passes the cuts
 339     Bool_t rpOK = kTRUE;
 340     Bool_t poiOK = kTRUE;
 341     if (rpCFManager && poiCFManager)
 342     {
 343       rpOK = ( rpCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&
 344                rpCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));
 345       poiOK = ( poiCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&
 346                 poiCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));
 347     }
 348     if (!(rpOK || poiOK)) continue;
 349
 350     //make new AliFlowTrack
 351     AliFlowTrack* pTrack = new AliFlowTrack(pParticle);
 352     pTrack->SetSource(AliFlowTrack::kFromAOD);
 353
 354     if (rpOK /* && rpCFManager */ ) // to be fixed - with CF managers uncommented only empty events (NULL in header files)
 355     {
 356       pTrack->SetForRPSelection(kTRUE);
 357       IncrementNumberOfPOIs(0);
 358     }
 359     if (poiOK /* && poiCFManager*/ )
 360     {
 361       pTrack->SetForPOISelection(kTRUE);
 362       IncrementNumberOfPOIs(1);
 363     }
 364     AddTrack(pTrack);
 365   }
 366
 367   //  if (iSelParticlesRP >= fMinMult && iSelParticlesRP <= fMaxMult)
 368   //  {
 369   //    if ( (++fCount % 100) == 0)
 370   //    {
 371   //      if (!fMCReactionPlaneAngle == 0) cout<<" MC Reaction Plane Angle = "<<  fMCReactionPlaneAngle << endl;
 372   //      else cout<<" MC Reaction Plane Angle = unknown "<< endl;
 373   //      cout<<" iGoodTracks = "<<iGoodTracks<<endl;
 374   //      cout<<" # of RP selected tracks = "<<iSelParticlesRP<<endl;
 375   //      cout<<" # of POI selected tracks = "<<iSelParticlesPOI<<endl;
 376   //      cout << "# " << fCount << " events processed" << endl;
 377   //    }
 378   //    return pEvent;
 379   //  }
 380   //  else
 381   //  {
 382   //    cout<<"Not enough tracks in the FlowEventSimple"<<endl;
 383   //    return 0;
 384   //  }
 385   //}
 386   //else
 387   //{
 388   //  cout<<"Event does not pass multiplicity cuts"<<endl;
 389   //  return 0;
 390   //}
 391
 392 }
 393
 394 //-----------------------------------------------------------------------
 395 AliFlowEvent::AliFlowEvent( const AliESDEvent* anInput,
 396                             const AliMCEvent* anInputMc,
 397                             KineSource anOption,
 398                             const AliCFManager* rpCFManager,
 399                             const AliCFManager* poiCFManager ):
 400   AliFlowEventSimple(20), fApplyRecentering(-1), fCachedRun(-1), fCurrentCentrality(-1)
 401 {
 402     // constructor
 403     for(Int_t i(0); i < 9; i++) {
 404         for(Int_t j(0); j < 2; j++) {
 405             for(Int_t k(0); k < 2; k++) {
 406                 fMeanQ[i][j][k] = 0.;
 407                 fWidthQ[i][j][k] = 0.;
 408                 fMeanQv3[i][j][k] = 0.;
 409                 fWidthQv3[i][j][k] = 0.;
 410             }
 411         }
 412     }
 413   //fills the event with tracks from the ESD and kinematics from the MC info via the track label
 414   if (anOption==kNoKine)
 415   {
 416     AliFatal("WRONG OPTION IN AliFlowEventMaker::FillTracks(AliESDEvent* anInput, AliMCEvent* anInputMc, KineSource anOption)");
 417     exit(1);
 418   }
 419
 420   Int_t iNumberOfInputTracks = anInput->GetNumberOfTracks() ;
 421
 422   Int_t iNumberOfInputTracksMC = anInputMc->GetNumberOfTracks() ;
 423   if (iNumberOfInputTracksMC==-1)
 424   {
 425     AliError("Skipping Event -- No MC information available for this event");
 426     return;
 427   }
 428
 429   //loop over ESD tracks
 430   for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)
 431   {
 432     AliESDtrack* pParticle = anInput->GetTrack(itrkN);   //get input particle
 433     //get Label
 434     Int_t iLabel = pParticle->GetLabel();
 435     //match to mc particle
 436     AliMCParticle* pMcParticle = (AliMCParticle*) anInputMc->GetTrack(TMath::Abs(iLabel));
 437
 438     //check
 439     if (TMath::Abs(pParticle->GetLabel())!=pMcParticle->Label())
 440       AliWarning(Form("pParticle->GetLabel()!=pMcParticle->Label(), %i, %i", pParticle->GetLabel(), pMcParticle->Label()));
 441
 442     //check if pParticle passes the cuts
 443     Bool_t rpOK = kFALSE;
 444     Bool_t poiOK = kFALSE;
 445     if (rpCFManager && poiCFManager)
 446     {
 447       if(anOption == kESDkine)
 448       {
 449         if (rpCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,pMcParticle,"mcGenCuts1") &&
 450             rpCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle))
 451           rpOK=kTRUE;
 452         if (poiCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,pMcParticle,"mcGenCuts2") &&
 453             poiCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle))
 454           poiOK=kTRUE;
 455       }
 456       else if (anOption == kMCkine)
 457       {
 458         if (rpCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,pMcParticle))
 459           rpOK=kTRUE;
 460         if (poiCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,pMcParticle))
 461           poiOK=kTRUE;
 462       }
 463     }
 464
 465     if (!(rpOK || poiOK)) continue;
 466
 467     //make new AliFlowTrack
 468     AliFlowTrack* pTrack = NULL;
 469     if(anOption == kESDkine)   //take the PID from the MC & the kinematics from the ESD
 470     {
 471       pTrack = new AliFlowTrack(pParticle);
 472     }
 473     else if (anOption == kMCkine)   //take the PID and kinematics from the MC
 474     {
 475       pTrack = new AliFlowTrack(pMcParticle);
 476     }
 477
 478     if (rpOK && rpCFManager)
 479     {
 480       IncrementNumberOfPOIs(0);
 481       pTrack->SetForRPSelection();
 482     }
 483     if (poiOK && poiCFManager)
 484     {
 485       IncrementNumberOfPOIs(1);
 486       pTrack->SetForPOISelection();
 487     }
 488
 489     AddTrack(pTrack);
 490   }
 491   SetMCReactionPlaneAngle(anInputMc);
 492 }
 493
 494 //-----------------------------------------------------------------------
 495 AliFlowEvent::AliFlowEvent( const AliESDEvent* anInput,
 496                             const AliMultiplicity* anInputTracklets,
 497                             const AliCFManager* poiCFManager ):
 498   AliFlowEventSimple(20), fApplyRecentering(-1), fCachedRun(-1), fCurrentCentrality(-1)
 499 {
 500     // constructor
 501     for(Int_t i(0); i < 9; i++) {
 502         for(Int_t j(0); j < 2; j++) {
 503             for(Int_t k(0); k < 2; k++) {
 504                 fMeanQ[i][j][k] = 0.;
 505                 fWidthQ[i][j][k] = 0.;
 506                 fMeanQv3[i][j][k] = 0.;
 507                 fWidthQv3[i][j][k] = 0.;
 508             }
 509         }
 510     }
 511
 512   //Select the particles of interest from the ESD
 513   Int_t iNumberOfInputTracks = anInput->GetNumberOfTracks() ;
 514
 515   //loop over tracks
 516   for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)
 517     {
 518       AliESDtrack* pParticle = anInput->GetTrack(itrkN);   //get input particle
 519
 520       //check if pParticle passes the cuts
 521       Bool_t poiOK = kTRUE;
 522       if (poiCFManager)
 523         {
 524           poiOK = ( poiCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&
 525                     poiCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));
 526         }
 527       if (!poiOK) continue;
 528
 529       //make new AliFLowTrack
 530       AliFlowTrack* pTrack = new AliFlowTrack(pParticle);
 531
 532       //marking the particles used for the particle of interest (POI) selection:
 533       if(poiOK && poiCFManager)
 534         {
 535           IncrementNumberOfPOIs(1);
 536           pTrack->SetForPOISelection(kTRUE);
 537           pTrack->SetSource(AliFlowTrack::kFromESD);
 538         }
 539
 540       AddTrack(pTrack);
 541     }//end of while (itrkN < iNumberOfInputTracks)
 542
 543   //Select the reference particles from the SPD tracklets
 544   anInputTracklets = anInput->GetMultiplicity();
 545   Int_t multSPD = anInputTracklets->GetNumberOfTracklets();
 546
 547   //loop over tracklets
 548   for (Int_t itracklet=0; itracklet<multSPD; ++itracklet) {
 549     Float_t thetaTr= anInputTracklets->GetTheta(itracklet);
 550     Float_t phiTr= anInputTracklets->GetPhi(itracklet);
 551     // calculate eta
 552     Float_t etaTr = -TMath::Log(TMath::Tan(thetaTr/2.));
 553
 554     //make new AliFLowTrackSimple
 555     AliFlowTrack* pTrack = new AliFlowTrack();
 556     pTrack->SetPt(0.0);
 557     pTrack->SetEta(etaTr);
 558     pTrack->SetPhi(phiTr);
 559     //marking the particles used for the reference particle (RP) selection:
 560     IncrementNumberOfPOIs(0);
 561     pTrack->SetForRPSelection(kTRUE);
 562     pTrack->SetSource(AliFlowTrack::kFromTracklet);
 563
 564     //Add the track to the flowevent
 565     AddTrack(pTrack);
 566   }
 567
 568 }
 569
 570 //-----------------------------------------------------------------------
 571 AliFlowEvent::AliFlowEvent( const AliESDEvent* esd,
 572                             const AliCFManager* poiCFManager,
 573                             Bool_t hybrid):
 574   AliFlowEventSimple(20), fApplyRecentering(-1), fCachedRun(-1), fCurrentCentrality(-1)
 575 {
 576     // constructor
 577     for(Int_t i(0); i < 9; i++) {
 578         for(Int_t j(0); j < 2; j++) {
 579             for(Int_t k(0); k < 2; k++) {
 580                 fMeanQ[i][j][k] = 0.;
 581                 fWidthQ[i][j][k] = 0.;
 582                 fMeanQv3[i][j][k] = 0.;
 583                 fWidthQv3[i][j][k] = 0.;
 584             }
 585         }
 586     }
 587
 588   //Select the particles of interest from the ESD
 589   Int_t iNumberOfInputTracks = esd->GetNumberOfTracks() ;
 590
 591   //Double_t gPt = 0.0, gP = 0.0;
 592   Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};  //The impact parameters and their covariance.
 593 //  Double_t dca3D = 0.0;       FIXME unused variable
 594
 595   AliESDtrack trackTPC;
 596
 597   //loop over tracks
 598   for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)
 599     {
 600
 601       if (!esd->GetTrack(itrkN)) continue;
 602
 603       Bool_t useTPC = kFALSE;
 604
 605       AliESDtrack* pParticle = esd->GetTrack(itrkN);   //get input particle
 606
 607       //check if pParticle passes the cuts
 608       Bool_t poiOK = kTRUE;
 609
 610       if (poiCFManager)
 611       {
 612         poiOK = ( poiCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&
 613                   poiCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));
 614       }
 615
 616       if (!(poiOK)) continue;
 617
 618       AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)pParticle->GetTPCInnerParam();
 619
 620       if (tpcTrack)
 621       {
 622
 623 //      gPt = tpcTrack->Pt();
 624 //      gP = tpcTrack->P();
 625
 626         useTPC = kTRUE;
 627
 628         const AliESDVertex *vertexSPD = esd->GetPrimaryVertexSPD();
 629         const AliESDVertex *vertexTPC = esd->GetPrimaryVertexTPC();
 630
 631         if(hybrid)
 632           tpcTrack->PropagateToDCA(vertexSPD,esd->GetMagneticField(),100.,dca,cov);
 633         else
 634           tpcTrack->PropagateToDCA(vertexTPC,esd->GetMagneticField(),100.,dca,cov);
 635
 636 //        dca3D = TMath::Sqrt(TMath::Power(dca[0],2)+TMath::Power(dca[1],2));   FIXME unused variable
 637
 638       }
 639
 640       //make new AliFLowTrack
 641       AliFlowTrack* pTrack = new AliFlowTrack(pParticle);
 642
 643       pTrack->SetSource(AliFlowTrack::kFromESD);
 644
 645       //marking the particles used for diff. flow:
 646       if(poiOK && poiCFManager)
 647       {
 648         pTrack->SetForPOISelection(kTRUE);
 649         IncrementNumberOfPOIs(1);
 650       }
 651
 652       if(useTPC)
 653       {
 654         pTrack->SetForRPSelection(kTRUE);
 655         IncrementNumberOfPOIs(0);
 656       }
 657
 658       AddTrack(pTrack);
 659
 660     }//end of while (itrkN < iNumberOfInputTracks)
 661
 662 }
 663
 664 //-----------------------------------------------------------------------
 665 AliFlowEvent::AliFlowEvent( const AliESDEvent* anInput,
 666                             const TH2F* anInputFMDhist,
 667                             const AliCFManager* poiCFManager ):
 668   AliFlowEventSimple(20), fApplyRecentering(-1), fCachedRun(-1), fCurrentCentrality(-1)
 669 {
 670     // constructor
 671     for(Int_t i(0); i < 9; i++) {
 672         for(Int_t j(0); j < 2; j++) {
 673             for(Int_t k(0); k < 2; k++) {
 674                 fMeanQ[i][j][k] = 0.;
 675                 fWidthQ[i][j][k] = 0.;
 676                 fMeanQv3[i][j][k] = 0.;
 677                 fWidthQv3[i][j][k] = 0.;
 678             }
 679         }
 680     }
 681
 682   //Select the particles of interest from the ESD
 683   Int_t iNumberOfInputTracks = anInput->GetNumberOfTracks() ;
 684
 685   //loop over tracks
 686   for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)
 687     {
 688       AliESDtrack* pParticle = anInput->GetTrack(itrkN);   //get input particle
 689
 690       //check if pParticle passes the cuts
 691       Bool_t poiOK = kTRUE;
 692       if (poiCFManager)
 693         {
 694           poiOK = ( poiCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&
 695                     poiCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));
 696         }
 697       if (!poiOK) continue;
 698
 699       //make new AliFLowTrack
 700       AliFlowTrack* pTrack = new AliFlowTrack(pParticle);
 701
 702       //marking the particles used for the particle of interest (POI) selection:
 703       if(poiOK && poiCFManager)
 704         {
 705           IncrementNumberOfPOIs(1);
 706           pTrack->SetForPOISelection(kTRUE);
 707           pTrack->SetSource(AliFlowTrack::kFromESD);
 708         }
 709
 710       AddTrack(pTrack);
 711     }//end of while (itrkN < iNumberOfInputTracks)
 712
 713   //Select the reference particles from the FMD hits
 714   //loop over FMD histogram
 715   Int_t iBinsEta = anInputFMDhist->GetNbinsX();
 716   Int_t iBinsPhi = anInputFMDhist->GetNbinsY();
 717
 718   for (Int_t iEta = 1; iEta <= iBinsEta; iEta++){
 719     Double_t etaFMD = anInputFMDhist->GetXaxis()->GetBinCenter(iEta);
 720     for (Int_t iPhi = 1; iPhi <= iBinsPhi; iPhi++){
 721       Double_t phiFMD = anInputFMDhist->GetYaxis()->GetBinCenter(iPhi);
 722       Double_t weightFMD = anInputFMDhist->GetBinContent(iEta,iPhi);
 723
 724       if (weightFMD > 0.0) { //do not add empty bins
 725         //make new AliFLowTrackSimple
 726         AliFlowTrack* pTrack = new AliFlowTrack();
 727         pTrack->SetPt(0.0);
 728         pTrack->SetEta(etaFMD);
 729         pTrack->SetPhi(phiFMD);
 730         pTrack->SetWeight(weightFMD);
 731         //marking the particles used for the reference particle (RP) selection:
 732         pTrack->TagRP();
 733         IncrementNumberOfPOIs(0);
 734         pTrack->SetSource(AliFlowTrack::kFromFMD);
 735
 736         //Add the track to the flowevent
 737         AddTrack(pTrack);
 738
 739       }
 740     }
 741   }
 742
 743 }
 744
 745 //-----------------------------------------------------------------------
 746 void AliFlowEvent::FindDaughters(Bool_t keepDaughtersInRPselection)
 747 {
 748   //each flow track holds it's esd track index as well as its daughters esd index.
 749   //fill the array of daughters for every track with the pointers to flow tracks
 750   //to associate the mothers with daughters directly
 751   for (Int_t iTrack=0; iTrack<fMothersCollection->GetEntriesFast(); iTrack++)
 752   {
 753     AliFlowTrack* mother = static_cast<AliFlowTrack*>(fMothersCollection->At(iTrack));
 754     if (!mother) continue;
 755     if (mother->GetNDaughters()<1) continue;
 756     for (Int_t iDaughterCandidate=0; iDaughterCandidate<fNumberOfTracks; iDaughterCandidate++)
 757     {
 758       AliFlowTrack* daughterCandidate = static_cast<AliFlowTrack*>(fTrackCollection->At(iDaughterCandidate));
 759       Int_t esdIndexDaughterCandidate = daughterCandidate->GetID();
 760       for (Int_t iDaughter=0; iDaughter<mother->GetNDaughters(); iDaughter++)
 761       {
 762         Int_t esdIndexDaughter = mother->GetIDDaughter(iDaughter);
 763         if (esdIndexDaughter==esdIndexDaughterCandidate)
 764         {
 765           mother->SetDaughter(iDaughter,daughterCandidate);
 766           daughterCandidate->SetForRPSelection(keepDaughtersInRPselection);
 767         }
 768       }
 769     }
 770   }
 771 }
 772
 773 //-----------------------------------------------------------------------
 774 void AliFlowEvent::Fill( AliFlowTrackCuts* rpCuts,
 775                          AliFlowTrackCuts* poiCuts )
 776 {
 777   //Fills the event from a vevent: AliESDEvent,AliAODEvent,AliMCEvent
 778   //the input data needs to be attached to the cuts
 779   //we have two cases, if we're cutting the same collection of tracks
 780   //(same param type) then we can have tracks that are both rp and poi
 781   //in the other case we want to have two exclusive sets of rps and pois
 782   //e.g. one tracklets, the other PMD or global - USER IS RESPOSIBLE
 783   //FOR MAKING SURE THEY DONT OVERLAP OR ELSE THE SAME PARTICLE WILL BE
 784   //TAKEN TWICE
 785
 786   ClearFast();
 787
 788   if (!rpCuts || !poiCuts) return;
 789   AliFlowTrackCuts::trackParameterType sourceRP = rpCuts->GetParamType();
 790   AliFlowTrackCuts::trackParameterType sourcePOI = poiCuts->GetParamType();
 791   AliFlowTrack* pTrack=NULL;
 792
 793   if (sourceRP==sourcePOI)
 794   {
 795     //loop over tracks
 796     Int_t numberOfInputObjects = rpCuts->GetNumberOfInputObjects();
 797     for (Int_t i=0; i<numberOfInputObjects; i++)
 798     {
 799       //get input object (particle)
 800       TObject* particle = rpCuts->GetInputObject(i);
 801
 802       Bool_t rp = rpCuts->IsSelected(particle,i);
 803       Int_t poiClass = poiCuts->IsSelected(particle,i);
 804
 805       if (!(rp||poiClass>0)) continue;
 806
 807       //make new AliFlowTrack
 808       if (rp)
 809       {
 810         pTrack = rpCuts->FillFlowTrack(fTrackCollection,fNumberOfTracks);
 811         if (!pTrack) continue;
 812         pTrack->TagRP(); IncrementNumberOfPOIs(0);
 813         if (poiClass>0) {pTrack->Tag(poiClass); IncrementNumberOfPOIs(poiClass);}
 814         if (pTrack->GetNDaughters()>0) fMothersCollection->Add(pTrack);
 815       }
 816       else if (poiClass>0)
 817       {
 818         pTrack = poiCuts->FillFlowTrack(fTrackCollection,fNumberOfTracks);
 819         if (!pTrack) continue;
 820         pTrack->Tag(poiClass); IncrementNumberOfPOIs(poiClass);
 821         if (pTrack->GetNDaughters()>0) fMothersCollection->Add(pTrack);
 822       }
 823       fNumberOfTracks++;
 824     }//end of while (i < numberOfTracks)
 825   }
 826   else if (sourceRP!=sourcePOI)
 827   {
 828     //here we have two different sources of particles, so we fill
 829     //them independently
 830     //POI
 831     for (Int_t i=0; i<poiCuts->GetNumberOfInputObjects(); i++)
 832     {
 833       TObject* particle = poiCuts->GetInputObject(i);
 834       Int_t poiClass = poiCuts->IsSelected(particle,i);
 835       if (poiClass<=0) continue;
 836       pTrack = poiCuts->FillFlowTrack(fTrackCollection,fNumberOfTracks);
 837       if (!pTrack) continue;
 838       pTrack->Tag(poiClass);
 839       IncrementNumberOfPOIs(poiClass);
 840       fNumberOfTracks++;
 841       if (pTrack->GetNDaughters()>0) fMothersCollection->Add(pTrack);
 842     }
 843     //RP
 844     Int_t numberOfInputObjects = rpCuts->GetNumberOfInputObjects();
 845     for (Int_t i=0; i<numberOfInputObjects; i++)
 846       {
 847       TObject* particle = rpCuts->GetInputObject(i);
 848       Int_t rp = rpCuts->IsSelected(particle,i);
 849       if (rp<1) continue;
 850       pTrack = rpCuts->FillFlowTrack(fTrackCollection,fNumberOfTracks);
 851       if (!pTrack) continue;
 852       pTrack->TagRP();
 853       IncrementNumberOfPOIs(0);
 854       fNumberOfTracks++;
 855       if (pTrack->GetNDaughters()>0) fMothersCollection->Add(pTrack);
 856     }
 857   }
 858 }
 859
 860 //-----------------------------------------------------------------------
 861 void AliFlowEvent::InsertTrack(AliFlowTrack *track) {
 862   // adds a flow track at the end of the container
 863   AliFlowTrack *pTrack = ReuseTrack( fNumberOfTracks++ );
 864   *pTrack = *track;
 865   if (track->GetNDaughters()>0)
 866   {
 867     fMothersCollection->Add(track);
 868   }
 869   //pTrack->SetPt( track->Pt() );
 870   //pTrack->SetPhi( track->Phi() );
 871   //pTrack->SetEta( track->Eta() );
 872   //pTrack->SetWeight( track->Weight() );
 873   //pTrack->SetCharge( track->Charge() );
 874   //pTrack->SetMass( track->Mass() );
 875   //pTrack->SetForRPSelection( track->InRPSelection() );
 876   //pTrack->SetForPOISelection( track->InPOISelection() );
 877   //if(track->InSubevent(0)) pTrack->SetForSubevent(0);
 878   //if(track->InSubevent(1)) pTrack->SetForSubevent(1);
 879   //pTrack->SetID( track->GetID() );
 880   return;
 881 }
 882
 883 //-----------------------------------------------------------------------
 884 AliFlowTrack* AliFlowEvent::ReuseTrack(Int_t i)
 885 {
 886   //try to reuse an existing track, if empty, make new one
 887   AliFlowTrack* pTrack = static_cast<AliFlowTrack*>(fTrackCollection->At(i));
 888   if (pTrack)
 889   {
 890     pTrack->Clear();
 891   }
 892   else
 893   {
 894     pTrack = new AliFlowTrack();
 895     fTrackCollection->AddAtAndExpand(pTrack,i);
 896   }
 897   return pTrack;
 898 }
 899
 900 //-----------------------------------------------------------------------
 901 AliFlowEvent::AliFlowEvent( AliFlowTrackCuts* rpCuts,
 902                             AliFlowTrackCuts* poiCuts ):
 903   AliFlowEventSimple(20), fApplyRecentering(kFALSE), fCachedRun(-1), fCurrentCentrality(-1)
 904 {
 905     // constructor
 906     for(Int_t i(0); i < 9; i++) {
 907         for(Int_t j(0); j < 2; j++) {
 908             for(Int_t k(0); k < 2; k++) {
 909                 fMeanQ[i][j][k] = 0.;
 910                 fWidthQ[i][j][k] = 0.;
 911                 fMeanQv3[i][j][k] = 0.;
 912                 fWidthQv3[i][j][k] = 0.;
 913             }
 914         }
 915     }
 916   //Fills the event from a vevent: AliESDEvent,AliAODEvent,AliMCEvent
 917   //the input data needs to be attached to the cuts
 918   //we have two cases, if we're cutting the same collection of tracks
 919   //(same param type) then we can have tracks that are both rp and poi
 920   //in the other case we want to have two exclusive sets of rps and pois
 921   //e.g. one tracklets, the other PMD or global - USER IS RESPOSIBLE
 922   //FOR MAKING SURE THEY DONT OVERLAP OR ELSE THE SAME PARTICLE WILL BE
 923   //TAKEN TWICE
 924
 925   if (!rpCuts || !poiCuts) return;
 926   AliFlowTrackCuts::trackParameterType sourceRP = rpCuts->GetParamType();
 927   AliFlowTrackCuts::trackParameterType sourcePOI = poiCuts->GetParamType();
 928
 929   if (sourceRP==sourcePOI)
 930   {
 931     //loop over tracks
 932     Int_t numberOfInputObjects = rpCuts->GetNumberOfInputObjects();
 933     for (Int_t i=0; i<numberOfInputObjects; i++)
 934     {
 935       //get input object (particle)
 936       TObject* particle = rpCuts->GetInputObject(i);
 937
 938       Bool_t rp = rpCuts->IsSelected(particle,i);
 939       Int_t poiClass = poiCuts->IsSelected(particle,i);
 940
 941       if (!(rp||poiClass>0)) continue;
 942
 943       //make new AliFLowTrack
 944       AliFlowTrack* pTrack = NULL;
 945       if (rp)
 946       {
 947         pTrack = rpCuts->FillFlowTrack(fTrackCollection,fNumberOfTracks);
 948         if (!pTrack) continue;
 949         pTrack->TagRP(); IncrementNumberOfPOIs(0);
 950         if (poiClass>0) {pTrack->Tag(poiClass); IncrementNumberOfPOIs(poiClass);}
 951       }
 952       else
 953       if (poiClass>0)
 954       {
 955         pTrack = poiCuts->FillFlowTrack(fTrackCollection,fNumberOfTracks);
 956         if (!pTrack) continue;
 957         pTrack->Tag(poiClass); IncrementNumberOfPOIs(poiClass);
 958       }
 959       TrackAdded();
 960     }//end of while (i < numberOfTracks)
 961   }
 962   else if (sourceRP!=sourcePOI)
 963   {
 964     //here we have two different sources of particles, so we fill
 965     //them independently
 966     AliFlowTrack* pTrack = NULL;
 967     //RP
 968     Int_t numberOfInputObjects = rpCuts->GetNumberOfInputObjects();
 969     for (Int_t i=0; i<numberOfInputObjects; i++)
 970     {
 971       TObject* particle = rpCuts->GetInputObject(i);
 972       Bool_t rp = rpCuts->IsSelected(particle,i);
 973       if (!rp) continue;
 974       pTrack = rpCuts->FillFlowTrack(fTrackCollection,fNumberOfTracks);
 975       if (!pTrack) continue;
 976       pTrack->TagRP(); IncrementNumberOfPOIs(0);
 977       TrackAdded();
 978     }
 979     //POI
 980     numberOfInputObjects = poiCuts->GetNumberOfInputObjects();
 981     for (Int_t i=0; i<numberOfInputObjects; i++)
 982     {
 983       TObject* particle = poiCuts->GetInputObject(i);
 984       Int_t poiClass = poiCuts->IsSelected(particle,i);
 985       if (poiClass<=0) continue;
 986       pTrack = poiCuts->FillFlowTrack(fTrackCollection,fNumberOfTracks);
 987       if (!pTrack) continue;
 988       pTrack->Tag(poiClass); IncrementNumberOfPOIs(poiClass);
 989       TrackAdded();
 990     }
 991   }
 992 }
 993
 994 //-------------------------------------------------------------------//
 995 //---- Including PMD tracks as RP --------------------------//
 996
 997 AliFlowEvent::AliFlowEvent( const AliESDEvent* anInput,
 998                             const AliESDPmdTrack *pmdtracks,
 999                             const AliCFManager* poiCFManager ):
1000   AliFlowEventSimple(20), fApplyRecentering(kFALSE), fCachedRun(-1), fCurrentCentrality(-1)
1001 {
1002     // constructor
1003     for(Int_t i(0); i < 9; i++) {
1004         for(Int_t j(0); j < 2; j++) {
1005             for(Int_t k(0); k < 2; k++) {
1006                 fMeanQ[i][j][k] = 0.;
1007                 fWidthQ[i][j][k] = 0.;
1008                 fMeanQv3[i][j][k] = 0.;
1009                 fWidthQv3[i][j][k] = 0.;
1010             }
1011         }
1012     }
1013   Float_t GetPmdEta(Float_t xPos, Float_t yPos, Float_t zPos);
1014   Float_t GetPmdPhi(Float_t xPos, Float_t yPos);
1015   //Select the particles of interest from the ESD
1016   Int_t iNumberOfInputTracks = anInput->GetNumberOfTracks() ;
1017
1018   //loop over tracks
1019   for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)
1020     {
1021       AliESDtrack* pParticle = anInput->GetTrack(itrkN);   //get input particle
1022       //check if pParticle passes the cuts
1023       Bool_t poiOK = kTRUE;
1024       if (poiCFManager)
1025         {
1026           poiOK = ( poiCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&
1027                     poiCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));
1028         }
1029       if (!poiOK) continue;
1030
1031       //make new AliFLowTrack
1032       AliFlowTrack* pTrack = new AliFlowTrack(pParticle);
1033
1034       //marking the particles used for the particle of interest (POI) selection:
1035       if(poiOK && poiCFManager)
1036         {
1037           IncrementNumberOfPOIs(1);
1038           pTrack->SetForPOISelection(kTRUE);
1039           pTrack->SetSource(AliFlowTrack::kFromESD);
1040         }
1041
1042       AddTrack(pTrack);
1043     }//end of while (itrkN < iNumberOfInputTracks)
1044
1045   //Select the reference particles from the PMD tracks
1046   Int_t npmdcl = anInput->GetNumberOfPmdTracks();
1047   printf("======There are %d PMD tracks in this event\n-------",npmdcl);
1048   //loop over clusters
1049   for(Int_t iclust=0; iclust < npmdcl; iclust++){
1050     //AliESDPmdTrack *pmdtr = anInput->GetPmdTrack(iclust);
1051     pmdtracks = anInput->GetPmdTrack(iclust);
1052     Int_t   det   = pmdtracks->GetDetector();
1053     //Int_t   smn   = pmdtracks->GetSmn();
1054     Float_t clsX  = pmdtracks->GetClusterX();
1055     Float_t clsY  = pmdtracks->GetClusterY();
1056     Float_t clsZ  = pmdtracks->GetClusterZ();
1057     Float_t ncell = pmdtracks->GetClusterCells();
1058     Float_t adc   = pmdtracks->GetClusterADC();
1059     //Float_t pid   = pmdtracks->GetClusterPID();
1060     Float_t etacls = GetPmdEta(clsX,clsY,clsZ);
1061     Float_t phicls = GetPmdPhi(clsX,clsY);
1062     //make new AliFLowTrackSimple
1063     AliFlowTrack* pTrack = new AliFlowTrack();
1064     //if(det == 0){ //selecting preshower plane only
1065     if(det == 0 && adc > 270 && ncell > 1){ //selecting preshower plane only
1066       //pTrack->SetPt(adc);//cluster adc
1067       pTrack->SetPt(0.0);
1068       pTrack->SetEta(etacls);
1069       pTrack->SetPhi(phicls);
1070       //marking the particles used for the reference particle (RP) selection:
1071       IncrementNumberOfPOIs(0);
1072       pTrack->SetForRPSelection(kTRUE);
1073       pTrack->SetSource(AliFlowTrack::kFromPMD);
1074       //Add the track to the flowevent
1075       AddTrack(pTrack);
1076     }//if det
1077   }
1078 }
1079 //----------------------------------------------------------------------------//
1080 Float_t GetPmdEta(Float_t xPos, Float_t yPos, Float_t zPos)
1081 {
1082   Float_t rpxpy, theta, eta;
1083   rpxpy  = TMath::Sqrt(xPos*xPos + yPos*yPos);
1084   theta  = TMath::ATan2(rpxpy,zPos);
1085   eta    = -TMath::Log(TMath::Tan(0.5*theta));
1086   return eta;
1087 }
1088 //--------------------------------------------------------------------------//
1089 Float_t GetPmdPhi(Float_t xPos, Float_t yPos)
1090 {
1091   Float_t pybypx, phi = 0., phi1;
1092   if(xPos==0)
1093     {
1094       if(yPos>0) phi = 90.;
1095       if(yPos<0) phi = 270.;
1096     }
1097   if(xPos != 0)
1098     {
1099       pybypx = yPos/xPos;
1100       if(pybypx < 0) pybypx = - pybypx;
1101       phi1 = TMath::ATan(pybypx)*180./3.14159;
1102
1103       if(xPos > 0 && yPos > 0) phi = phi1;        // 1st Quadrant
1104       if(xPos < 0 && yPos > 0) phi = 180 - phi1;  // 2nd Quadrant
1105       if(xPos < 0 && yPos < 0) phi = 180 + phi1;  // 3rd Quadrant
1106       if(xPos > 0 && yPos < 0) phi = 360 - phi1;  // 4th Quadrant
1107
1108     }
1109   phi = phi*3.14159/180.;
1110   return   phi;
1111 }
1112 //---------------------------------------------------------------//
1113
1114 void AliFlowEvent::Get2Qsub(AliFlowVector* Qarray, Int_t n, TList *weightsList, Bool_t usePhiWeights, Bool_t usePtWeights, Bool_t useEtaWeights)
1115 {
1116   // get q vectors for the subevents. if no recentering is necessary, get the guy from the flow event simple
1117   AliFlowEventSimple::Get2Qsub(Qarray, n, weightsList, usePhiWeights, usePtWeights, useEtaWeights);
1118   // else get the recentering from the cached info
1119   if (fApplyRecentering == 2010)        // 10h style recentering
1120   {
1121     // first retrieve the q-vectors from the AliFlowEventSimple:: routine
1122     AliFlowVector vA = Qarray[0];
1123     AliFlowVector vB = Qarray[1];
1124     // extract the information form the current flow vectors
1125     Double_t Qxc(vA.X());       // IMPORTANT: user is responsible for the sign of eta
1126     Double_t Qyc(vA.Y());       // vzeroC has negative pseudorapidity and is taken as subevent A
1127     Double_t Qxa(vB.X());       // vzeroA has positive pseudorapidity and is taken as subevent B
1128     Double_t Qya(vB.Y());
1129     // init some values for the corrections
1130
1131     // values for vector a (VZEROA)
1132     Double_t Qxamean(0);
1133     Double_t Qxarms(1);
1134     Double_t Qyamean(0);
1135     Double_t Qyarms(1);
1136     // values for vector b (VZEROC)
1137     Double_t Qxcmean(0);
1138     Double_t Qxcrms(1);
1139     Double_t Qycmean(0);
1140     Double_t Qycrms(1);
1141
1142     if( n == 2) {       // second order symmetry
1143         Qxamean = fMeanQ[fCurrentCentrality][1][0];
1144         Qxarms  = fWidthQ[fCurrentCentrality][1][0];
1145         Qyamean = fMeanQ[fCurrentCentrality][1][1];
1146         Qyarms  = fWidthQ[fCurrentCentrality][1][1];
1147
1148         Qxcmean = fMeanQ[fCurrentCentrality][0][0];
1149         Qxcrms  = fWidthQ[fCurrentCentrality][0][0];
1150         Qycmean = fMeanQ[fCurrentCentrality][0][1];
1151         Qycrms  = fWidthQ[fCurrentCentrality][0][1];
1152     } else if (n == 3) {        // third order symmetry
1153         Qxamean = fMeanQv3[fCurrentCentrality][1][0];
1154         Qxarms  = fWidthQv3[fCurrentCentrality][1][0];
1155         Qyamean = fMeanQv3[fCurrentCentrality][1][1];
1156         Qyarms  = fWidthQv3[fCurrentCentrality][1][1];
1157
1158         Qxcmean = fMeanQv3[fCurrentCentrality][0][0];
1159         Qxcrms  = fWidthQv3[fCurrentCentrality][0][0];
1160         Qycmean = fMeanQv3[fCurrentCentrality][0][1];
1161         Qycrms  = fWidthQv3[fCurrentCentrality][0][1];
1162     }
1163     // do the correction
1164     Double_t QxaCor = (Qxa - Qxamean)/Qxarms;
1165     Double_t QyaCor = (Qya - Qyamean)/Qyarms;
1166     Double_t QxcCor = (Qxc - Qxcmean)/Qxcrms;
1167     Double_t QycCor = (Qyc - Qycmean)/Qycrms;
1168     // update the vector
1169     vA.Set(QxcCor, QycCor);
1170     vB.Set(QxaCor, QyaCor);
1171   } else if (fApplyRecentering == 2011) { // 11h style recentering
1172     // in this case, the q-vectors are repaced by the ones from
1173     // the event header
1174
1175     // first retrieve the q-vectors from the AliFlowEventSimple:: routine
1176     AliFlowVector vA = Qarray[0];
1177     AliFlowVector vB = Qarray[1];
1178
1179     Double_t QxaCor = 0.;
1180     Double_t QyaCor = 0.;
1181     Double_t QxcCor = 0.;
1182     Double_t QycCor = 0.;
1183
1184     // copy the new q-vectors from the cache
1185     if(n == 2) {
1186        QxaCor = fMeanQ[0][1][0];
1187        QyaCor = fMeanQ[0][1][1];
1188        QxcCor = fMeanQ[0][0][0];
1189        QycCor = fMeanQ[0][0][1];
1190     } else if (n == 3) {
1191        QxaCor = fMeanQv3[0][1][0];
1192        QyaCor = fMeanQv3[0][1][1];
1193        QxcCor = fMeanQv3[0][0][0];
1194        QycCor = fMeanQv3[0][0][1];
1195     }
1196     // set the new q-vectors (which in this case means REPLACING)
1197     vA.Set(QxcCor, QycCor);
1198     vB.Set(QxaCor, QyaCor);
1199   }
1200 }
1201 //_____________________________________________________________________________
1202 void AliFlowEvent::SetVZEROCalibrationForTrackCuts(AliFlowTrackCuts* cuts) {
1203     // open calibration info, copied from AliAnalyisTaskVnV0.cxx
1204     if(!cuts->GetEvent()) return; // coverity. we need to know the event to get the runnumber and centrlaity
1205     // get the vzero centrality percentile (cc dependent calibration)
1206     Float_t v0Centr(cuts->GetEvent()->GetCentrality()->GetCentralityPercentile("V0M"));
1207     if(v0Centr < 5) fCurrentCentrality = 0;
1208     else if(v0Centr < 10) fCurrentCentrality = 1;
1209     else if(v0Centr < 20) fCurrentCentrality = 2;
1210     else if(v0Centr < 30) fCurrentCentrality = 3;
1211     else if(v0Centr < 40) fCurrentCentrality = 4;
1212     else if(v0Centr < 50) fCurrentCentrality = 5;
1213     else if(v0Centr < 60) fCurrentCentrality = 6;
1214     else if(v0Centr < 70) fCurrentCentrality = 7;
1215     else fCurrentCentrality = 8;
1216
1217     // if this event is from the same run as the previous event
1218     // we can use the cached calibration values, no need to re-open the
1219     // aodb file
1220     Int_t run(cuts->GetEvent()->GetRunNumber());
1221 //    printf ( " > run number is %i \n", run);
1222     if(fCachedRun == run) {
1223         // the runnumber did not change, no need to open the database again
1224         // in case of 11h style recentering, update the q-sub vectors
1225         if(fApplyRecentering == 2011) SetVZEROCalibrationForTrackCuts2011(cuts);
1226         return;
1227     }
1228     // set the chached run number
1229     fCachedRun = run;
1230
1231     TString oadbfilename = "$ALICE_ROOT/OADB/PWGCF/VZERO/VZEROcalibEP.root";
1232     TFile *foadb = TFile::Open(oadbfilename.Data());
1233
1234     if(!foadb){
1235         printf("OADB file %s cannot be opened\n",oadbfilename.Data());
1236         return;
1237     }
1238
1239     AliOADBContainer *cont = (AliOADBContainer*) foadb->Get("hMultV0BefCorr");
1240     if(!cont){
1241         printf("OADB object hMultV0BefCorr is not available in the file\n");
1242         return;
1243     }
1244     if(!(cont->GetObject(run))){
1245         // if the multiplicity correction cannot be found for the specified run,
1246         // loop over the 11h runs to see if it's 11h data
1247         Int_t runs11h[] = {170593, 170572, 170556, 170552, 170546, 170390, 170389, 170388, 170387, 170315, 170313, 170312, 170311, 170309, 170308, 170306, 170270, 170269, 170268, 170267, 170264, 170230, 170228, 170208, 170207, 170205, 170204, 170203, 170195, 170193, 170163, 170162, 170159, 170155, 170152, 170091, 170089, 170088, 170085, 170084, 170083, 170081, 170040, 170038, 170036, 170027, 169981, 169975, 169969, 169965, 169961, 169956, 169926, 169924, 169923, 169922, 169919, 169918, 169914, 169859, 169858, 169855, 169846, 169838, 169837, 169835, 169683, 169628, 169591, 169590, 169588, 169587, 169586, 169584, 169557, 169555, 169554, 169553, 169550, 169515, 169512, 169506, 169504, 169498, 169475, 169420, 169419, 169418, 169417, 169415, 169411, 169238, 169236, 169167, 169160, 169156, 169148, 169145, 169144, 169143, 169138, 169099, 169094, 169091, 169045, 169044, 169040, 169035, 168992, 168988, 168984, 168826, 168777, 168514, 168512, 168511, 168467, 168464, 168461, 168460, 168458, 168362, 168361, 168356, 168342, 168341, 168325, 168322, 168318, 168311, 168310, 168213, 168212, 168208, 168207, 168206, 168205, 168204, 168203, 168181, 168177, 168175, 168173, 168172, 168171, 168115, 168108, 168107, 168105, 168104, 168103, 168076, 168069, 168068, 168066, 167988, 167987, 167986, 167985, 167921, 167920, 167915, 167909, 167903, 167902, 167818, 167814, 167813, 167808, 167807, 167806, 167713, 167712, 167711, 167706, 167693};
1248         for(Int_t r(0); r < 176; r++) {
1249             if(run == runs11h[r]) {
1250                 printf(" > run has been identified as 11h < \n");
1251                 if(cuts->GetVZEROgainEqualizationPerRing()) {
1252                     // enable or disable rings through the weights, weight 1. is enabled, 0. is disabled
1253                     // start with the vzero c rings (segments 0 through 31)
1254                     (cuts->GetUseVZERORing(0)) ? cuts->SetVZEROCpol(0, 1.) : cuts->SetVZEROCpol(0, 0.);
1255                     (cuts->GetUseVZERORing(1)) ? cuts->SetVZEROCpol(1, 1.) : cuts->SetVZEROCpol(1, 0.);
1256                     (cuts->GetUseVZERORing(2)) ? cuts->SetVZEROCpol(2, 1.) : cuts->SetVZEROCpol(2, 0.);
1257                     (cuts->GetUseVZERORing(3)) ? cuts->SetVZEROCpol(3, 1.) : cuts->SetVZEROCpol(3, 0.);
1258                     // same for vzero a
1259                     (cuts->GetUseVZERORing(4)) ? cuts->SetVZEROApol(0, 1.) : cuts->SetVZEROApol(0, 0.);
1260                     (cuts->GetUseVZERORing(5)) ? cuts->SetVZEROApol(1, 1.) : cuts->SetVZEROApol(1, 0.);
1261                     (cuts->GetUseVZERORing(6)) ? cuts->SetVZEROApol(2, 1.) : cuts->SetVZEROApol(2, 0.);
1262                     (cuts->GetUseVZERORing(7)) ? cuts->SetVZEROApol(3, 1.) : cuts->SetVZEROApol(3, 0.);
1263                 } else {
1264                     // else enable all rings
1265                     for(Int_t i(0); i < 4; i++) cuts->SetVZEROCpol(i, 1.);
1266                     for(Int_t i(0); i < 4; i++) cuts->SetVZEROApol(i, 1.);
1267                 }
1268                 // pass a NULL pointer to the track cuts object
1269                 // the NULL pointer will identify 11h runs
1270                 cuts->SetVZEROgainEqualisation(NULL);
1271                 // this will identify the recentering style that is required. flight might be changed if recenetering is disabled
1272                 fApplyRecentering = 2011;
1273                 SetVZEROCalibrationForTrackCuts2011(cuts);
1274                 return; // the rest of the steps are not necessary
1275             }
1276         }
1277         // the run has not been identified as lhc11h data, so we assume a template calibration
1278         printf("OADB object hMultVZEROBefCorr is not available for run %i (used run 137366)\n",run);
1279         run = 137366;
1280     }
1281     printf(" > run has been identified as 10h < \n");
1282     // step 1) get the proper multiplicity weights from the vzero signal
1283     TProfile* fMultVZERO = ((TH2F *) cont->GetObject(run))->ProfileX();
1284
1285     TF1 *fpol0 = new TF1("fpol0","pol0");
1286     if(cuts->GetVZEROgainEqualizationPerRing()) {
1287         // do the calibration per ring
1288         // start with the vzero c rings (segments 0 through 31)
1289         fMultVZERO->Fit(fpol0, "", "", 0, 8);
1290         (cuts->GetUseVZERORing(0)) ? cuts->SetVZEROCpol(0, fpol0->GetParameter(0)) : cuts->SetVZEROCpol(0, 0.);
1291         fMultVZERO->Fit(fpol0, "", "", 8, 16);
1292         (cuts->GetUseVZERORing(1)) ? cuts->SetVZEROCpol(1, fpol0->GetParameter(0)) : cuts->SetVZEROCpol(1, 0.);
1293         fMultVZERO->Fit(fpol0, "", "", 16, 24);
1294         (cuts->GetUseVZERORing(2)) ? cuts->SetVZEROCpol(2, fpol0->GetParameter(0)) : cuts->SetVZEROCpol(2, 0.);
1295         fMultVZERO->Fit(fpol0, "", "", 24, 32);
1296         (cuts->GetUseVZERORing(3)) ? cuts->SetVZEROCpol(3, fpol0->GetParameter(0)) : cuts->SetVZEROCpol(3, 0.);
1297         // same thing for vero A
1298         fMultVZERO->Fit(fpol0, "", "", 32, 40);
1299         (cuts->GetUseVZERORing(4)) ? cuts->SetVZEROApol(0, fpol0->GetParameter(0)) : cuts->SetVZEROApol(0, 0.);
1300         fMultVZERO->Fit(fpol0, "", "", 40, 48);
1301         (cuts->GetUseVZERORing(5)) ? cuts->SetVZEROApol(1, fpol0->GetParameter(0)) : cuts->SetVZEROApol(1, 0.);
1302         fMultVZERO->Fit(fpol0, "", "", 48, 56);
1303         (cuts->GetUseVZERORing(6)) ? cuts->SetVZEROApol(2, fpol0->GetParameter(0)) : cuts->SetVZEROApol(2, 0.);
1304         fMultVZERO->Fit(fpol0, "", "", 56, 64);
1305         (cuts->GetUseVZERORing(7)) ? cuts->SetVZEROApol(3, fpol0->GetParameter(0)) : cuts->SetVZEROApol(3, 0.);
1306     } else {
1307         // do the calibration in one go. the calibration will still be
1308         // stored per ring, but each ring has the same weight now
1309        fMultVZERO->Fit(fpol0,"","",0,31);
1310        for(Int_t i(0); i < 4; i++) cuts->SetVZEROCpol(i, fpol0->GetParameter(0));
1311        fMultVZERO->Fit(fpol0,"","",32,64);
1312        for(Int_t i(0); i < 4; i++) cuts->SetVZEROApol(i, fpol0->GetParameter(0));
1313     }
1314     // the parameters to weigh the vzero track cuts have been extracted now,
1315     // so we can pass them to the current track cuts obect
1316     cuts->SetVZEROgainEqualisation(fMultVZERO);       // passed as a TH1
1317
1318     // step 2) reweight the q-vectors that will be  called by flow methods which use
1319     // subevents
1320     // underlying assumption is that subevent a uses VZEROA
1321     // and subevent b uses VZEROC
1322     for(Int_t iside=0;iside<2;iside++){
1323         for(Int_t icoord=0;icoord<2;icoord++){
1324             for(Int_t i=0;i  < 9;i++){
1325                 char namecont[100];
1326                 if(iside==0 && icoord==0)
1327                   snprintf(namecont,100,"hQxc2_%i",i);
1328                 else if(iside==1 && icoord==0)
1329                   snprintf(namecont,100,"hQxa2_%i",i);
1330                 else if(iside==0 && icoord==1)
1331                   snprintf(namecont,100,"hQyc2_%i",i);
1332                 else if(iside==1 && icoord==1)
1333                   snprintf(namecont,100,"hQya2_%i",i);
1334
1335                 cont = (AliOADBContainer*) foadb->Get(namecont);
1336                 if(!cont){
1337                     printf("OADB object %s is not available in the file\n",namecont);
1338                     return;
1339                 }
1340
1341                 if(!(cont->GetObject(run))){
1342                     printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
1343                     run = 137366;
1344                 }
1345
1346                 // after grabbing all the info, set the CORRECTION TERMS to
1347                 // the 2nd and 3rd order qsub-vectors
1348                 // we do this here for all centralities, so that subsequent events
1349                 // can grab the correction from these cached values
1350                 fMeanQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
1351                 fWidthQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
1352
1353                 //for v3
1354                 if(iside==0 && icoord==0)
1355                   snprintf(namecont,100,"hQxc3_%i",i);
1356                 else if(iside==1 && icoord==0)
1357                   snprintf(namecont,100,"hQxa3_%i",i);
1358                 else if(iside==0 && icoord==1)
1359                   snprintf(namecont,100,"hQyc3_%i",i);
1360                 else if(iside==1 && icoord==1)
1361                   snprintf(namecont,100,"hQya3_%i",i);
1362
1363                 cont = (AliOADBContainer*) foadb->Get(namecont);
1364                 if(!cont){
1365                     printf("OADB object %s is not available in the file\n",namecont);
1366                     return;
1367                 }
1368
1369                 if(!(cont->GetObject(run))){
1370                     printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
1371                     run = 137366;
1372                 }
1373                 fMeanQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
1374                 fWidthQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
1375
1376             }
1377         }
1378     }
1379     // set the recentering style (might be switched back to -1 if recentering is disabeled)
1380     fApplyRecentering = 2010;
1381 }
1382 //_____________________________________________________________________________
1383 void AliFlowEvent::SetVZEROCalibrationForTrackCuts2011(AliFlowTrackCuts* cuts)
1384 {
1385     // load the vzero q-sub vectors
1386     if(!cuts->GetEvent() || !cuts->GetEvent()->GetEventplane()) return;       // coverity
1387     Double_t qxEPa = 0, qyEPa = 0;
1388     Double_t qxEPc = 0, qyEPc = 0;
1389     Double_t qxEPa3 = 0, qyEPa3 = 0;
1390     Double_t qxEPc3 = 0, qyEPc3 = 0;
1391
1392     // get the q-vectors from the header
1393     cuts->GetEvent()->GetEventplane()->CalculateVZEROEventPlane(cuts->GetEvent(), 8, 2, qxEPa, qyEPa);
1394     cuts->GetEvent()->GetEventplane()->CalculateVZEROEventPlane(cuts->GetEvent(), 9, 2, qxEPc, qyEPc);
1395     cuts->GetEvent()->GetEventplane()->CalculateVZEROEventPlane(cuts->GetEvent(), 8, 3, qxEPa3, qyEPa3);
1396     cuts->GetEvent()->GetEventplane()->CalculateVZEROEventPlane(cuts->GetEvent(), 9, 3, qxEPc3, qyEPc3);
1397
1398     // store the values temporarily. this may seem
1399     // inelegant, but we don't want to include
1400     // aliflowtrackcuts or alivevnet in get2qsub
1401
1402     // qx and qy for vzero a, second harmonc
1403     fMeanQ[0][1][0] = qxEPa;
1404     fMeanQ[0][1][1] = qyEPa;
1405     // qx and qx for vzero c, second harmonic
1406     fMeanQ[0][0][0] = qxEPc;
1407     fMeanQ[0][0][1] = qyEPc;
1408     // qx and qy for vzero a, third harmonic
1409     fMeanQv3[0][1][0] = qxEPa3;
1410     fMeanQv3[0][1][1] = qyEPa3;
1411     // qx and qy for vzero c, third harmonic
1412     fMeanQv3[0][0][0] = qxEPc3;
1413     fMeanQv3[0][0][1] = qyEPc3;
1414 }
1415 //_____________________________________________________________________________
1416