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), fVZEROcentralityBin(-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), fVZEROcentralityBin(-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), fVZEROcentralityBin(-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   fVZEROcentralityBin = event.fVZEROcentralityBin;
 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), fVZEROcentralityBin(-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), fVZEROcentralityBin(-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), fVZEROcentralityBin(-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), fVZEROcentralityBin(-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), fVZEROcentralityBin(-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), fVZEROcentralityBin(-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         AliExternalTrackParam copy(*tpcTrack);
 632         if(hybrid)
 633           copy.PropagateToDCA(vertexSPD,esd->GetMagneticField(),100.,dca,cov);
 634         else
 635           copy.PropagateToDCA(vertexTPC,esd->GetMagneticField(),100.,dca,cov);
 636
 637 //        dca3D = TMath::Sqrt(TMath::Power(dca[0],2)+TMath::Power(dca[1],2));   FIXME unused variable
 638
 639       }
 640
 641       //make new AliFLowTrack
 642       AliFlowTrack* pTrack = new AliFlowTrack(pParticle);
 643
 644       pTrack->SetSource(AliFlowTrack::kFromESD);
 645
 646       //marking the particles used for diff. flow:
 647       if(poiOK && poiCFManager)
 648       {
 649         pTrack->SetForPOISelection(kTRUE);
 650         IncrementNumberOfPOIs(1);
 651       }
 652
 653       if(useTPC)
 654       {
 655         pTrack->SetForRPSelection(kTRUE);
 656         IncrementNumberOfPOIs(0);
 657       }
 658
 659       AddTrack(pTrack);
 660
 661     }//end of while (itrkN < iNumberOfInputTracks)
 662
 663 }
 664
 665 //-----------------------------------------------------------------------
 666 AliFlowEvent::AliFlowEvent( const AliESDEvent* anInput,
 667                             const TH2F* anInputFMDhist,
 668                             const AliCFManager* poiCFManager ):
 669   AliFlowEventSimple(20), fApplyRecentering(-1), fCachedRun(-1), fVZEROcentralityBin(-1)
 670 {
 671     // constructor
 672     for(Int_t i(0); i < 9; i++) {
 673         for(Int_t j(0); j < 2; j++) {
 674             for(Int_t k(0); k < 2; k++) {
 675                 fMeanQ[i][j][k] = 0.;
 676                 fWidthQ[i][j][k] = 0.;
 677                 fMeanQv3[i][j][k] = 0.;
 678                 fWidthQv3[i][j][k] = 0.;
 679             }
 680         }
 681     }
 682
 683   //Select the particles of interest from the ESD
 684   Int_t iNumberOfInputTracks = anInput->GetNumberOfTracks() ;
 685
 686   //loop over tracks
 687   for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)
 688     {
 689       AliESDtrack* pParticle = anInput->GetTrack(itrkN);   //get input particle
 690
 691       //check if pParticle passes the cuts
 692       Bool_t poiOK = kTRUE;
 693       if (poiCFManager)
 694         {
 695           poiOK = ( poiCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&
 696                     poiCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));
 697         }
 698       if (!poiOK) continue;
 699
 700       //make new AliFLowTrack
 701       AliFlowTrack* pTrack = new AliFlowTrack(pParticle);
 702
 703       //marking the particles used for the particle of interest (POI) selection:
 704       if(poiOK && poiCFManager)
 705         {
 706     IncrementNumberOfPOIs(1);
 707           pTrack->SetForPOISelection(kTRUE);
 708           pTrack->SetSource(AliFlowTrack::kFromESD);
 709         }
 710
 711       AddTrack(pTrack);
 712     }//end of while (itrkN < iNumberOfInputTracks)
 713
 714   //Select the reference particles from the FMD hits
 715   //loop over FMD histogram
 716   Int_t iBinsEta = anInputFMDhist->GetNbinsX();
 717   Int_t iBinsPhi = anInputFMDhist->GetNbinsY();
 718
 719   for (Int_t iEta = 1; iEta <= iBinsEta; iEta++){
 720     Double_t etaFMD = anInputFMDhist->GetXaxis()->GetBinCenter(iEta);
 721     for (Int_t iPhi = 1; iPhi <= iBinsPhi; iPhi++){
 722       Double_t phiFMD = anInputFMDhist->GetYaxis()->GetBinCenter(iPhi);
 723       Double_t weightFMD = anInputFMDhist->GetBinContent(iEta,iPhi);
 724
 725       if (weightFMD > 0.0) { //do not add empty bins
 726         //make new AliFLowTrackSimple
 727         AliFlowTrack* pTrack = new AliFlowTrack();
 728         pTrack->SetPt(0.0);
 729         pTrack->SetEta(etaFMD);
 730         pTrack->SetPhi(phiFMD);
 731         pTrack->SetWeight(weightFMD);
 732         //marking the particles used for the reference particle (RP) selection:
 733         pTrack->TagRP();
 734         IncrementNumberOfPOIs(0);
 735         pTrack->SetSource(AliFlowTrack::kFromFMD);
 736
 737         //Add the track to the flowevent
 738         AddTrack(pTrack);
 739
 740       }
 741     }
 742   }
 743
 744 }
 745
 746 //-----------------------------------------------------------------------
 747 void AliFlowEvent::FindDaughters(Bool_t keepDaughtersInRPselection)
 748 {
 749   //each flow track holds it's esd track index as well as its daughters esd index.
 750   //fill the array of daughters for every track with the pointers to flow tracks
 751   //to associate the mothers with daughters directly
 752   for (Int_t iTrack=0; iTrack<fMothersCollection->GetEntriesFast(); iTrack++)
 753   {
 754     AliFlowTrack* mother = static_cast<AliFlowTrack*>(fMothersCollection->At(iTrack));
 755     if (!mother) continue;
 756     if (mother->GetNDaughters()<1) continue;
 757     for (Int_t iDaughterCandidate=0; iDaughterCandidate<fNumberOfTracks; iDaughterCandidate++)
 758     {
 759       AliFlowTrack* daughterCandidate = static_cast<AliFlowTrack*>(fTrackCollection->At(iDaughterCandidate));
 760       Int_t esdIndexDaughterCandidate = daughterCandidate->GetID();
 761       for (Int_t iDaughter=0; iDaughter<mother->GetNDaughters(); iDaughter++)
 762       {
 763         Int_t esdIndexDaughter = mother->GetIDDaughter(iDaughter);
 764         if (esdIndexDaughter==esdIndexDaughterCandidate)
 765         {
 766           mother->SetDaughter(iDaughter,daughterCandidate);
 767           daughterCandidate->SetForRPSelection(keepDaughtersInRPselection);
 768         }
 769       }
 770     }
 771   }
 772 }
 773
 774 //-----------------------------------------------------------------------
 775 void AliFlowEvent::Fill( AliFlowTrackCuts* rpCuts,
 776                          AliFlowTrackCuts* poiCuts )
 777 {
 778   //Fills the event from a vevent: AliESDEvent,AliAODEvent,AliMCEvent
 779   //the input data needs to be attached to the cuts
 780   //we have two cases, if we're cutting the same collection of tracks
 781   //(same param type) then we can have tracks that are both rp and poi
 782   //in the other case we want to have two exclusive sets of rps and pois
 783   //e.g. one tracklets, the other PMD or global - USER IS RESPOSIBLE
 784   //FOR MAKING SURE THEY DONT OVERLAP OR ELSE THE SAME PARTICLE WILL BE
 785   //TAKEN TWICE
 786
 787   ClearFast();
 788
 789   if (!rpCuts || !poiCuts) return;
 790   AliFlowTrackCuts::trackParameterType sourceRP = rpCuts->GetParamType();
 791   AliFlowTrackCuts::trackParameterType sourcePOI = poiCuts->GetParamType();
 792   AliFlowTrack* pTrack=NULL;
 793
 794   if (sourceRP==sourcePOI)
 795   {
 796     //loop over tracks
 797     Int_t numberOfInputObjects = rpCuts->GetNumberOfInputObjects();
 798     for (Int_t i=0; i<numberOfInputObjects; i++)
 799     {
 800       //get input object (particle)
 801       TObject* particle = rpCuts->GetInputObject(i);
 802
 803       Bool_t rp = rpCuts->IsSelected(particle,i);
 804       Bool_t poi = poiCuts->IsSelected(particle,i);
 805
 806       if (!(rp||poi)) continue;
 807
 808       //make new AliFlowTrack
 809       if (rp)
 810       {
 811         pTrack = rpCuts->FillFlowTrack(fTrackCollection,fNumberOfTracks);
 812         if (!pTrack) continue;
 813         pTrack->Tag(0); IncrementNumberOfPOIs(0);
 814         if (poi) {pTrack->Tag(1); IncrementNumberOfPOIs(1);}
 815         if (pTrack->GetNDaughters()>0) fMothersCollection->Add(pTrack);
 816       }
 817       else if (poi)
 818       {
 819         pTrack = poiCuts->FillFlowTrack(fTrackCollection,fNumberOfTracks);
 820         if (!pTrack) continue;
 821         pTrack->Tag(1); IncrementNumberOfPOIs(1);
 822         if (pTrack->GetNDaughters()>0) fMothersCollection->Add(pTrack);
 823       }
 824       fNumberOfTracks++;
 825     }//end of while (i < numberOfTracks)
 826   }
 827   else if (sourceRP!=sourcePOI)
 828   {
 829     //here we have two different sources of particles, so we fill
 830     //them independently
 831     //POI
 832     for (Int_t i=0; i<poiCuts->GetNumberOfInputObjects(); i++)
 833     {
 834       TObject* particle = poiCuts->GetInputObject(i);
 835       Bool_t poi = poiCuts->IsSelected(particle,i);
 836       if (!poi) continue;
 837       pTrack = poiCuts->FillFlowTrack(fTrackCollection,fNumberOfTracks);
 838       if (!pTrack) continue;
 839       pTrack->Tag(1);
 840       IncrementNumberOfPOIs(1);
 841       fNumberOfTracks++;
 842       if (pTrack->GetNDaughters()>0) fMothersCollection->Add(pTrack);
 843     }
 844     //RP
 845     Int_t numberOfInputObjects = rpCuts->GetNumberOfInputObjects();
 846     for (Int_t i=0; i<numberOfInputObjects; i++)
 847       {
 848       TObject* particle = rpCuts->GetInputObject(i);
 849       Bool_t rp = rpCuts->IsSelected(particle,i);
 850       if (!rp) continue;
 851       pTrack = rpCuts->FillFlowTrack(fTrackCollection,fNumberOfTracks);
 852       if (!pTrack) continue;
 853       pTrack->Tag(0);
 854       IncrementNumberOfPOIs(0);
 855       fNumberOfTracks++;
 856       if (pTrack->GetNDaughters()>0) fMothersCollection->Add(pTrack);
 857     }
 858   }
 859 }
 860
 861 //-----------------------------------------------------------------------
 862 void AliFlowEvent::InsertTrack(AliFlowTrack *track) {
 863   // adds a flow track at the end of the container
 864   AliFlowTrack *pTrack = ReuseTrack( fNumberOfTracks++ );
 865   *pTrack = *track;
 866   if (track->GetNDaughters()>0)
 867   {
 868     fMothersCollection->Add(pTrack);
 869   }
 870   return;
 871 }
 872
 873 //-----------------------------------------------------------------------
 874 AliFlowTrack* AliFlowEvent::ReuseTrack(Int_t i)
 875 {
 876   //try to reuse an existing track, if empty, make new one
 877   AliFlowTrack* pTrack = static_cast<AliFlowTrack*>(fTrackCollection->At(i));
 878   if (pTrack)
 879   {
 880     pTrack->Clear();
 881   }
 882   else
 883   {
 884     pTrack = new AliFlowTrack();
 885     fTrackCollection->AddAtAndExpand(pTrack,i);
 886   }
 887   return pTrack;
 888 }
 889
 890 //-----------------------------------------------------------------------
 891 AliFlowEvent::AliFlowEvent( AliFlowTrackCuts* rpCuts,
 892                             AliFlowTrackCuts* poiCuts ):
 893   AliFlowEventSimple(20), fApplyRecentering(kFALSE), fCachedRun(-1), fVZEROcentralityBin(-1)
 894 {
 895     // constructor
 896     for(Int_t i(0); i < 9; i++) {
 897         for(Int_t j(0); j < 2; j++) {
 898             for(Int_t k(0); k < 2; k++) {
 899                 fMeanQ[i][j][k] = 0.;
 900                 fWidthQ[i][j][k] = 0.;
 901                 fMeanQv3[i][j][k] = 0.;
 902                 fWidthQv3[i][j][k] = 0.;
 903             }
 904         }
 905     }
 906   //Fills the event from a vevent: AliESDEvent,AliAODEvent,AliMCEvent
 907   //the input data needs to be attached to the cuts
 908   Fill(rpCuts,poiCuts);
 909 }
 910
 911 //-------------------------------------------------------------------//
 912 //---- Including PMD tracks as RP --------------------------//
 913
 914 AliFlowEvent::AliFlowEvent( const AliESDEvent* anInput,
 915                             const AliESDPmdTrack *pmdtracks,
 916                             const AliCFManager* poiCFManager ):
 917   AliFlowEventSimple(20), fApplyRecentering(kFALSE), fCachedRun(-1), fVZEROcentralityBin(-1)
 918 {
 919     // constructor
 920     for(Int_t i(0); i < 9; i++) {
 921         for(Int_t j(0); j < 2; j++) {
 922             for(Int_t k(0); k < 2; k++) {
 923                 fMeanQ[i][j][k] = 0.;
 924                 fWidthQ[i][j][k] = 0.;
 925                 fMeanQv3[i][j][k] = 0.;
 926                 fWidthQv3[i][j][k] = 0.;
 927             }
 928         }
 929     }
 930   Float_t GetPmdEta(Float_t xPos, Float_t yPos, Float_t zPos);
 931   Float_t GetPmdPhi(Float_t xPos, Float_t yPos);
 932   //Select the particles of interest from the ESD
 933   Int_t iNumberOfInputTracks = anInput->GetNumberOfTracks() ;
 934
 935   //loop over tracks
 936   for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)
 937     {
 938       AliESDtrack* pParticle = anInput->GetTrack(itrkN);   //get input particle
 939       //check if pParticle passes the cuts
 940       Bool_t poiOK = kTRUE;
 941       if (poiCFManager)
 942         {
 943           poiOK = ( poiCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&
 944                     poiCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));
 945         }
 946       if (!poiOK) continue;
 947
 948       //make new AliFLowTrack
 949       AliFlowTrack* pTrack = new AliFlowTrack(pParticle);
 950
 951       //marking the particles used for the particle of interest (POI) selection:
 952       if(poiOK && poiCFManager)
 953         {
 954           IncrementNumberOfPOIs(1);
 955           pTrack->SetForPOISelection(kTRUE);
 956           pTrack->SetSource(AliFlowTrack::kFromESD);
 957         }
 958
 959       AddTrack(pTrack);
 960     }//end of while (itrkN < iNumberOfInputTracks)
 961
 962   //Select the reference particles from the PMD tracks
 963   Int_t npmdcl = anInput->GetNumberOfPmdTracks();
 964   printf("======There are %d PMD tracks in this event\n-------",npmdcl);
 965   //loop over clusters
 966   for(Int_t iclust=0; iclust < npmdcl; iclust++){
 967     //AliESDPmdTrack *pmdtr = anInput->GetPmdTrack(iclust);
 968     pmdtracks = anInput->GetPmdTrack(iclust);
 969     Int_t   det   = pmdtracks->GetDetector();
 970     //Int_t   smn   = pmdtracks->GetSmn();
 971     Float_t clsX  = pmdtracks->GetClusterX();
 972     Float_t clsY  = pmdtracks->GetClusterY();
 973     Float_t clsZ  = pmdtracks->GetClusterZ();
 974     Float_t ncell = pmdtracks->GetClusterCells();
 975     Float_t adc   = pmdtracks->GetClusterADC();
 976     //Float_t pid   = pmdtracks->GetClusterPID();
 977     Float_t etacls = GetPmdEta(clsX,clsY,clsZ);
 978     Float_t phicls = GetPmdPhi(clsX,clsY);
 979     //make new AliFLowTrackSimple
 980     AliFlowTrack* pTrack = new AliFlowTrack();
 981     //if(det == 0){ //selecting preshower plane only
 982     if(det == 0 && adc > 270 && ncell > 1){ //selecting preshower plane only
 983       //pTrack->SetPt(adc);//cluster adc
 984       pTrack->SetPt(0.0);
 985       pTrack->SetEta(etacls);
 986       pTrack->SetPhi(phicls);
 987       //marking the particles used for the reference particle (RP) selection:
 988       IncrementNumberOfPOIs(0);
 989       pTrack->SetForRPSelection(kTRUE);
 990       pTrack->SetSource(AliFlowTrack::kFromPMD);
 991       //Add the track to the flowevent
 992       AddTrack(pTrack);
 993     }//if det
 994   }
 995 }
 996 //----------------------------------------------------------------------------//
 997 Float_t GetPmdEta(Float_t xPos, Float_t yPos, Float_t zPos)
 998 {
 999   Float_t rpxpy, theta, eta;
1000   rpxpy  = TMath::Sqrt(xPos*xPos + yPos*yPos);
1001   theta  = TMath::ATan2(rpxpy,zPos);
1002   eta    = -TMath::Log(TMath::Tan(0.5*theta));
1003   return eta;
1004 }
1005 //--------------------------------------------------------------------------//
1006 Float_t GetPmdPhi(Float_t xPos, Float_t yPos)
1007 {
1008   Float_t pybypx, phi = 0., phi1;
1009   if(xPos==0)
1010     {
1011       if(yPos>0) phi = 90.;
1012       if(yPos<0) phi = 270.;
1013     }
1014   if(xPos != 0)
1015     {
1016       pybypx = yPos/xPos;
1017       if(pybypx < 0) pybypx = - pybypx;
1018       phi1 = TMath::ATan(pybypx)*180./3.14159;
1019
1020       if(xPos > 0 && yPos > 0) phi = phi1;        // 1st Quadrant
1021       if(xPos < 0 && yPos > 0) phi = 180 - phi1;  // 2nd Quadrant
1022       if(xPos < 0 && yPos < 0) phi = 180 + phi1;  // 3rd Quadrant
1023       if(xPos > 0 && yPos < 0) phi = 360 - phi1;  // 4th Quadrant
1024
1025     }
1026   phi = phi*3.14159/180.;
1027   return   phi;
1028 }
1029 //---------------------------------------------------------------//
1030
1031 void AliFlowEvent::Get2Qsub(AliFlowVector* Qarray, Int_t n, TList *weightsList, Bool_t usePhiWeights, Bool_t usePtWeights, Bool_t useEtaWeights)
1032 {
1033   // get q vectors for the subevents. if no recentering is necessary, get the guy from the flow event simple
1034   AliFlowEventSimple::Get2Qsub(Qarray, n, weightsList, usePhiWeights, usePtWeights, useEtaWeights);
1035   // else get the recentering from the cached info
1036   if (fApplyRecentering == 2010)        // 10h style recentering
1037   {
1038     // first retrieve the q-vectors from the AliFlowEventSimple:: routine
1039     AliFlowVector vA = Qarray[0];
1040     AliFlowVector vB = Qarray[1];
1041     // extract the information form the current flow vectors
1042     Double_t Qxc(vA.X());       // IMPORTANT: user is responsible for the sign of eta
1043     Double_t Qyc(vA.Y());       // vzeroC has negative pseudorapidity and is taken as subevent A
1044     Double_t Qxa(vB.X());       // vzeroA has positive pseudorapidity and is taken as subevent B
1045     Double_t Qya(vB.Y());
1046     // init some values for the corrections
1047
1048     // values for vector a (VZEROA)
1049     Double_t Qxamean(0);
1050     Double_t Qxarms(1);
1051     Double_t Qyamean(0);
1052     Double_t Qyarms(1);
1053     // values for vector b (VZEROC)
1054     Double_t Qxcmean(0);
1055     Double_t Qxcrms(1);
1056     Double_t Qycmean(0);
1057     Double_t Qycrms(1);
1058
1059     if( n == 2) {       // second order symmetry
1060         Qxamean = fMeanQ[fVZEROcentralityBin][1][0];
1061         Qxarms  = fWidthQ[fVZEROcentralityBin][1][0];
1062         Qyamean = fMeanQ[fVZEROcentralityBin][1][1];
1063         Qyarms  = fWidthQ[fVZEROcentralityBin][1][1];
1064
1065         Qxcmean = fMeanQ[fVZEROcentralityBin][0][0];
1066         Qxcrms  = fWidthQ[fVZEROcentralityBin][0][0];
1067         Qycmean = fMeanQ[fVZEROcentralityBin][0][1];
1068         Qycrms  = fWidthQ[fVZEROcentralityBin][0][1];
1069     } else if (n == 3) {        // third order symmetry
1070         Qxamean = fMeanQv3[fVZEROcentralityBin][1][0];
1071         Qxarms  = fWidthQv3[fVZEROcentralityBin][1][0];
1072         Qyamean = fMeanQv3[fVZEROcentralityBin][1][1];
1073         Qyarms  = fWidthQv3[fVZEROcentralityBin][1][1];
1074
1075         Qxcmean = fMeanQv3[fVZEROcentralityBin][0][0];
1076         Qxcrms  = fWidthQv3[fVZEROcentralityBin][0][0];
1077         Qycmean = fMeanQv3[fVZEROcentralityBin][0][1];
1078         Qycrms  = fWidthQv3[fVZEROcentralityBin][0][1];
1079     }
1080     // do the correction
1081     Double_t QxaCor = (Qxa - Qxamean)/Qxarms;
1082     Double_t QyaCor = (Qya - Qyamean)/Qyarms;
1083     Double_t QxcCor = (Qxc - Qxcmean)/Qxcrms;
1084     Double_t QycCor = (Qyc - Qycmean)/Qycrms;
1085     // update the vector
1086     vA.Set(QxcCor, QycCor);
1087     vB.Set(QxaCor, QyaCor);
1088   } else if (fApplyRecentering == 2011) { // 11h style recentering
1089     // in this case, the q-vectors are repaced by the ones from
1090     // the event header
1091
1092     // first retrieve the q-vectors from the AliFlowEventSimple:: routine
1093     AliFlowVector vA = Qarray[0];
1094     AliFlowVector vB = Qarray[1];
1095
1096     Double_t QxaCor = 0.;
1097     Double_t QyaCor = 0.;
1098     Double_t QxcCor = 0.;
1099     Double_t QycCor = 0.;
1100
1101     // copy the new q-vectors from the cache
1102     if(n == 2) {
1103        QxaCor = fMeanQ[0][1][0];
1104        QyaCor = fMeanQ[0][1][1];
1105        QxcCor = fMeanQ[0][0][0];
1106        QycCor = fMeanQ[0][0][1];
1107     } else if (n == 3) {
1108        QxaCor = fMeanQv3[0][1][0];
1109        QyaCor = fMeanQv3[0][1][1];
1110        QxcCor = fMeanQv3[0][0][0];
1111        QycCor = fMeanQv3[0][0][1];
1112     }
1113     // set the new q-vectors (which in this case means REPLACING)
1114     vA.Set(QxcCor, QycCor);
1115     vB.Set(QxaCor, QyaCor);
1116   }
1117 }
1118 //_____________________________________________________________________________
1119 void AliFlowEvent::SetVZEROCalibrationForTrackCuts(AliFlowTrackCuts* cuts) {
1120     // open calibration info, copied from AliAnalyisTaskVnV0.cxx
1121     if(!cuts->GetEvent()) return; // coverity. we need to know the event to get the runnumber and centrlaity
1122     // get the vzero centrality percentile (cc dependent calibration)
1123     Float_t v0Centr(cuts->GetEvent()->GetCentrality()->GetCentralityPercentile("V0M"));
1124     if(v0Centr < 5) fVZEROcentralityBin = 0;
1125     else if(v0Centr < 10) fVZEROcentralityBin = 1;
1126     else if(v0Centr < 20) fVZEROcentralityBin = 2;
1127     else if(v0Centr < 30) fVZEROcentralityBin = 3;
1128     else if(v0Centr < 40) fVZEROcentralityBin = 4;
1129     else if(v0Centr < 50) fVZEROcentralityBin = 5;
1130     else if(v0Centr < 60) fVZEROcentralityBin = 6;
1131     else if(v0Centr < 70) fVZEROcentralityBin = 7;
1132     else fVZEROcentralityBin = 8;
1133
1134     // if this event is from the same run as the previous event
1135     // we can use the cached calibration values, no need to re-open the
1136     // aodb file
1137     Int_t run(cuts->GetEvent()->GetRunNumber());
1138 //    printf ( " > run number is %i \n", run);
1139     if(fCachedRun == run) {
1140         // the runnumber did not change, no need to open the database again
1141         // in case of 11h style recentering, update the q-sub vectors
1142         if(fApplyRecentering == 2011) SetVZEROCalibrationForTrackCuts2011(cuts);
1143         return;
1144     }
1145     // set the chached run number
1146     fCachedRun = run;
1147
1148     TString oadbfilename = "$ALICE_ROOT/OADB/PWGCF/VZERO/VZEROcalibEP.root";
1149     TFile *foadb = TFile::Open(oadbfilename.Data());
1150
1151     if(!foadb){
1152         printf("OADB file %s cannot be opened\n",oadbfilename.Data());
1153         return;
1154     }
1155
1156     AliOADBContainer *cont = (AliOADBContainer*) foadb->Get("hMultV0BefCorr");
1157     if(!cont){
1158         printf("OADB object hMultV0BefCorr is not available in the file\n");
1159         return;
1160     }
1161     if(!(cont->GetObject(run))){
1162         // if the multiplicity correction cannot be found for the specified run,
1163         // loop over the 11h runs to see if it's 11h data
1164         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};
1165         for(Int_t r(0); r < 176; r++) {
1166             if(run == runs11h[r]) {
1167                 printf(" > run has been identified as 11h < \n");
1168                 if(cuts->GetVZEROgainEqualizationPerRing()) {
1169                     // enable or disable rings through the weights, weight 1. is enabled, 0. is disabled
1170                     // start with the vzero c rings (segments 0 through 31)
1171                     (cuts->GetUseVZERORing(0)) ? cuts->SetVZEROCpol(0, 1.) : cuts->SetVZEROCpol(0, 0.);
1172                     (cuts->GetUseVZERORing(1)) ? cuts->SetVZEROCpol(1, 1.) : cuts->SetVZEROCpol(1, 0.);
1173                     (cuts->GetUseVZERORing(2)) ? cuts->SetVZEROCpol(2, 1.) : cuts->SetVZEROCpol(2, 0.);
1174                     (cuts->GetUseVZERORing(3)) ? cuts->SetVZEROCpol(3, 1.) : cuts->SetVZEROCpol(3, 0.);
1175                     // same for vzero a
1176                     (cuts->GetUseVZERORing(4)) ? cuts->SetVZEROApol(0, 1.) : cuts->SetVZEROApol(0, 0.);
1177                     (cuts->GetUseVZERORing(5)) ? cuts->SetVZEROApol(1, 1.) : cuts->SetVZEROApol(1, 0.);
1178                     (cuts->GetUseVZERORing(6)) ? cuts->SetVZEROApol(2, 1.) : cuts->SetVZEROApol(2, 0.);
1179                     (cuts->GetUseVZERORing(7)) ? cuts->SetVZEROApol(3, 1.) : cuts->SetVZEROApol(3, 0.);
1180                 } else {
1181                     // else enable all rings
1182                     for(Int_t i(0); i < 4; i++) cuts->SetVZEROCpol(i, 1.);
1183                     for(Int_t i(0); i < 4; i++) cuts->SetVZEROApol(i, 1.);
1184                 }
1185                 // pass a NULL pointer to the track cuts object
1186                 // the NULL pointer will identify 11h runs
1187                 cuts->SetVZEROgainEqualisation(NULL);
1188                 // this will identify the recentering style that is required. flight might be changed if recenetering is disabled
1189                 fApplyRecentering = 2011;
1190                 SetVZEROCalibrationForTrackCuts2011(cuts);
1191                 return; // the rest of the steps are not necessary
1192             }
1193         }
1194         // the run has not been identified as lhc11h data, so we assume a template calibration
1195         printf("OADB object hMultVZEROBefCorr is not available for run %i (used run 137366)\n",run);
1196         run = 137366;
1197     }
1198     printf(" > run has been identified as 10h < \n");
1199     // step 1) get the proper multiplicity weights from the vzero signal
1200     TProfile* fMultVZERO = ((TH2F *) cont->GetObject(run))->ProfileX();
1201
1202     TF1 *fpol0 = new TF1("fpol0","pol0");
1203     if(cuts->GetVZEROgainEqualizationPerRing()) {
1204         // do the calibration per ring
1205         // start with the vzero c rings (segments 0 through 31)
1206         fMultVZERO->Fit(fpol0, "", "", 0, 8);
1207         (cuts->GetUseVZERORing(0)) ? cuts->SetVZEROCpol(0, fpol0->GetParameter(0)) : cuts->SetVZEROCpol(0, 0.);
1208         fMultVZERO->Fit(fpol0, "", "", 8, 16);
1209         (cuts->GetUseVZERORing(1)) ? cuts->SetVZEROCpol(1, fpol0->GetParameter(0)) : cuts->SetVZEROCpol(1, 0.);
1210         fMultVZERO->Fit(fpol0, "", "", 16, 24);
1211         (cuts->GetUseVZERORing(2)) ? cuts->SetVZEROCpol(2, fpol0->GetParameter(0)) : cuts->SetVZEROCpol(2, 0.);
1212         fMultVZERO->Fit(fpol0, "", "", 24, 32);
1213         (cuts->GetUseVZERORing(3)) ? cuts->SetVZEROCpol(3, fpol0->GetParameter(0)) : cuts->SetVZEROCpol(3, 0.);
1214         // same thing for vero A
1215         fMultVZERO->Fit(fpol0, "", "", 32, 40);
1216         (cuts->GetUseVZERORing(4)) ? cuts->SetVZEROApol(0, fpol0->GetParameter(0)) : cuts->SetVZEROApol(0, 0.);
1217         fMultVZERO->Fit(fpol0, "", "", 40, 48);
1218         (cuts->GetUseVZERORing(5)) ? cuts->SetVZEROApol(1, fpol0->GetParameter(0)) : cuts->SetVZEROApol(1, 0.);
1219         fMultVZERO->Fit(fpol0, "", "", 48, 56);
1220         (cuts->GetUseVZERORing(6)) ? cuts->SetVZEROApol(2, fpol0->GetParameter(0)) : cuts->SetVZEROApol(2, 0.);
1221         fMultVZERO->Fit(fpol0, "", "", 56, 64);
1222         (cuts->GetUseVZERORing(7)) ? cuts->SetVZEROApol(3, fpol0->GetParameter(0)) : cuts->SetVZEROApol(3, 0.);
1223     } else {
1224         // do the calibration in one go. the calibration will still be
1225         // stored per ring, but each ring has the same weight now
1226        fMultVZERO->Fit(fpol0,"","",0,31);
1227        for(Int_t i(0); i < 4; i++) cuts->SetVZEROCpol(i, fpol0->GetParameter(0));
1228        fMultVZERO->Fit(fpol0,"","",32,64);
1229        for(Int_t i(0); i < 4; i++) cuts->SetVZEROApol(i, fpol0->GetParameter(0));
1230     }
1231     // the parameters to weigh the vzero track cuts have been extracted now,
1232     // so we can pass them to the current track cuts obect
1233     cuts->SetVZEROgainEqualisation(fMultVZERO);       // passed as a TH1
1234
1235     // step 2) reweight the q-vectors that will be  called by flow methods which use
1236     // subevents
1237     // underlying assumption is that subevent a uses VZEROA
1238     // and subevent b uses VZEROC
1239     for(Int_t iside=0;iside<2;iside++){
1240         for(Int_t icoord=0;icoord<2;icoord++){
1241             for(Int_t i=0;i  < 9;i++){
1242                 char namecont[100];
1243                 if(iside==0 && icoord==0)
1244                   snprintf(namecont,100,"hQxc2_%i",i);
1245                 else if(iside==1 && icoord==0)
1246                   snprintf(namecont,100,"hQxa2_%i",i);
1247                 else if(iside==0 && icoord==1)
1248                   snprintf(namecont,100,"hQyc2_%i",i);
1249                 else if(iside==1 && icoord==1)
1250                   snprintf(namecont,100,"hQya2_%i",i);
1251
1252                 cont = (AliOADBContainer*) foadb->Get(namecont);
1253                 if(!cont){
1254                     printf("OADB object %s is not available in the file\n",namecont);
1255                     return;
1256                 }
1257
1258                 if(!(cont->GetObject(run))){
1259                     printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
1260                     run = 137366;
1261                 }
1262
1263                 // after grabbing all the info, set the CORRECTION TERMS to
1264                 // the 2nd and 3rd order qsub-vectors
1265                 // we do this here for all centralities, so that subsequent events
1266                 // can grab the correction from these cached values
1267                 fMeanQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
1268                 fWidthQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
1269
1270                 //for v3
1271                 if(iside==0 && icoord==0)
1272                   snprintf(namecont,100,"hQxc3_%i",i);
1273                 else if(iside==1 && icoord==0)
1274                   snprintf(namecont,100,"hQxa3_%i",i);
1275                 else if(iside==0 && icoord==1)
1276                   snprintf(namecont,100,"hQyc3_%i",i);
1277                 else if(iside==1 && icoord==1)
1278                   snprintf(namecont,100,"hQya3_%i",i);
1279
1280                 cont = (AliOADBContainer*) foadb->Get(namecont);
1281                 if(!cont){
1282                     printf("OADB object %s is not available in the file\n",namecont);
1283                     return;
1284                 }
1285
1286                 if(!(cont->GetObject(run))){
1287                     printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
1288                     run = 137366;
1289                 }
1290                 fMeanQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
1291                 fWidthQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
1292
1293             }
1294         }
1295     }
1296     // set the recentering style (might be switched back to -1 if recentering is disabeled)
1297     fApplyRecentering = 2010;
1298 }
1299 //_____________________________________________________________________________
1300 void AliFlowEvent::SetVZEROCalibrationForTrackCuts2011(AliFlowTrackCuts* cuts)
1301 {
1302     // load the vzero q-sub vectors
1303     if(!cuts->GetEvent() || !cuts->GetEvent()->GetEventplane()) return;       // coverity
1304     Double_t qxEPa = 0, qyEPa = 0;
1305     Double_t qxEPc = 0, qyEPc = 0;
1306     Double_t qxEPa3 = 0, qyEPa3 = 0;
1307     Double_t qxEPc3 = 0, qyEPc3 = 0;
1308
1309     // get the q-vectors from the header
1310     cuts->GetEvent()->GetEventplane()->CalculateVZEROEventPlane(cuts->GetEvent(), 8, 2, qxEPa, qyEPa);
1311     cuts->GetEvent()->GetEventplane()->CalculateVZEROEventPlane(cuts->GetEvent(), 9, 2, qxEPc, qyEPc);
1312     cuts->GetEvent()->GetEventplane()->CalculateVZEROEventPlane(cuts->GetEvent(), 8, 3, qxEPa3, qyEPa3);
1313     cuts->GetEvent()->GetEventplane()->CalculateVZEROEventPlane(cuts->GetEvent(), 9, 3, qxEPc3, qyEPc3);
1314
1315     // store the values temporarily. this may seem
1316     // inelegant, but we don't want to include
1317     // aliflowtrackcuts or alivevnet in get2qsub
1318
1319     // qx and qy for vzero a, second harmonc
1320     fMeanQ[0][1][0] = qxEPa;
1321     fMeanQ[0][1][1] = qyEPa;
1322     // qx and qx for vzero c, second harmonic
1323     fMeanQ[0][0][0] = qxEPc;
1324     fMeanQ[0][0][1] = qyEPc;
1325     // qx and qy for vzero a, third harmonic
1326     fMeanQv3[0][1][0] = qxEPa3;
1327     fMeanQv3[0][1][1] = qyEPa3;
1328     // qx and qy for vzero c, third harmonic
1329     fMeanQv3[0][0][0] = qxEPc3;
1330     fMeanQv3[0][0][1] = qyEPc3;
1331 }
1332 //_____________________________________________________________________________
1333
1334 void AliFlowEvent::ClearFast()
1335 {
1336   //clear the event without releasing any memory
1337   AliFlowEventSimple::ClearFast();
1338   fApplyRecentering=0;
1339   fCachedRun=0;
1340 }