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), fEvent(0x0)
  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), fEvent(0x0)
  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), fEvent(0x0)
  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   fEvent = 0x0;                         // should never be copied
 117   for(Int_t i(0); i < 9; i++) {
 118       for(Int_t j(0); j < 2; j++) {
 119           for(Int_t k(0); k < 2; k++) {
 120               fMeanQ[i][j][k] = event.fMeanQ[i][j][k];
 121               fWidthQ[i][j][k] = event.fWidthQ[i][j][k];
 122               fMeanQv3[i][j][k] = event.fMeanQv3[i][j][k];
 123               fWidthQv3[i][j][k] = event.fWidthQv3[i][j][k];
 124           }
 125       }
 126   }
 127   AliFlowEventSimple::operator=(event);
 128   return *this;
 129 }
 130
 131 //-----------------------------------------------------------------------
 132 AliFlowTrack* AliFlowEvent::GetTrack(Int_t i)
 133 {
 134   //get track i from collection
 135   if (i>=fNumberOfTracks) return NULL;
 136   AliFlowTrack* pTrack = static_cast<AliFlowTrack*>(fTrackCollection->At(i)) ;
 137   return pTrack;
 138 }
 139
 140 //-----------------------------------------------------------------------
 141 void AliFlowEvent::SetMCReactionPlaneAngle(const AliMCEvent* mcEvent)
 142 {
 143   //sets the event plane angle from the proper header in the MC
 144
 145   //COCKTAIL with HIJING
 146   if (!strcmp(mcEvent-> GenEventHeader()->GetName(),"Cocktail Header"))   //returns 0 if matches
 147   {
 148     AliGenCocktailEventHeader *headerC = dynamic_cast<AliGenCocktailEventHeader *> (mcEvent-> GenEventHeader());
 149     if (headerC)
 150     {
 151       TList *lhd = headerC->GetHeaders();
 152       if (lhd)
 153       {
 154         AliGenHijingEventHeader *hdh = dynamic_cast<AliGenHijingEventHeader *> (lhd->At(0));
 155         if (hdh) AliFlowEventSimple::SetMCReactionPlaneAngle( hdh->ReactionPlaneAngle() );
 156       }
 157     }
 158   }
 159   //THERMINATOR
 160   else if (!strcmp(mcEvent-> GenEventHeader()->GetName(),"Therminator"))   //returns 0 if matches
 161   {
 162     AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
 163     if (headerH) AliFlowEventSimple::SetMCReactionPlaneAngle( headerH->ReactionPlaneAngle() );
 164   }
 165   //GEVSIM
 166   else if (!strcmp(mcEvent-> GenEventHeader()->GetName(),"GeVSim header"))   //returns 0 if matches
 167   {
 168     AliGenGeVSimEventHeader* headerG = dynamic_cast<AliGenGeVSimEventHeader*>(mcEvent->GenEventHeader());
 169     if (headerG) AliFlowEventSimple::SetMCReactionPlaneAngle( headerG->GetEventPlane() );
 170   }
 171   //HIJING
 172   else if (!strcmp(mcEvent-> GenEventHeader()->GetName(),"Hijing"))   //returns 0 if matches
 173   {
 174     AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
 175     if (headerH) AliFlowEventSimple::SetMCReactionPlaneAngle( headerH->ReactionPlaneAngle() );
 176   }
 177   //AMPT
 178   else if (!strcmp(mcEvent-> GenEventHeader()->GetName(),"Ampt"))   //returns 0 if matches
 179   {
 180     AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
 181     if (headerH) AliFlowEventSimple::SetMCReactionPlaneAngle( headerH->ReactionPlaneAngle() );
 182   }
 183   //EPOS
 184   else if (!strcmp(mcEvent->GenEventHeader()->GetName(),"EPOS"))
 185   {
 186     AliGenEposEventHeader* headerE = dynamic_cast<AliGenEposEventHeader*>(mcEvent->GenEventHeader());
 187     if (headerE) AliFlowEventSimple::SetMCReactionPlaneAngle( headerE->ReactionPlaneAngle() );
 188   }
 189   //Hydjet
 190   else
 191   {
 192     AliCollisionGeometry* header = dynamic_cast<AliCollisionGeometry*>(mcEvent->GenEventHeader());
 193     if (header) AliFlowEventSimple::SetMCReactionPlaneAngle( header->ReactionPlaneAngle() );
 194   }
 195 }
 196
 197 //-----------------------------------------------------------------------
 198 AliFlowEvent::AliFlowEvent( const AliMCEvent* anInput,
 199                             const AliCFManager* rpCFManager,
 200                             const AliCFManager* poiCFManager):
 201   AliFlowEventSimple(20), fApplyRecentering(-1), fCachedRun(-1), fVZEROcentralityBin(-1), fEvent(0x0)
 202 {
 203     // constructor
 204     for(Int_t i(0); i < 9; i++) {
 205         for(Int_t j(0); j < 2; j++) {
 206             for(Int_t k(0); k < 2; k++) {
 207                 fMeanQ[i][j][k] = 0.;
 208                 fWidthQ[i][j][k] = 0.;
 209                 fMeanQv3[i][j][k] = 0.;
 210                 fWidthQv3[i][j][k] = 0.;
 211             }
 212         }
 213     }
 214   //Fills the event from the MC kinematic information
 215
 216   Int_t iNumberOfInputTracks = anInput->GetNumberOfTracks() ;
 217
 218   //loop over tracks
 219   for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)
 220   {
 221     //get input particle
 222     AliMCParticle* pParticle = dynamic_cast<AliMCParticle*>(anInput->GetTrack(itrkN));
 223     if (!pParticle) continue;
 224
 225     //check if pParticle passes the cuts
 226     Bool_t rpOK = kTRUE;
 227     Bool_t poiOK = kTRUE;
 228     if (rpCFManager && poiCFManager)
 229     {
 230       rpOK = rpCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,pParticle);
 231       poiOK = poiCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,pParticle);
 232     }
 233     if (!(rpOK||poiOK)) continue;
 234
 235     AliFlowTrack* pTrack = new AliFlowTrack(pParticle);
 236     pTrack->SetSource(AliFlowTrack::kFromMC);
 237
 238     if (rpOK && rpCFManager)
 239     {
 240       pTrack->SetForRPSelection(kTRUE);
 241       IncrementNumberOfPOIs(0);
 242     }
 243     if (poiOK && poiCFManager)
 244     {
 245       pTrack->SetForPOISelection(kTRUE);
 246       IncrementNumberOfPOIs(1);
 247     }
 248
 249     AddTrack(pTrack) ;
 250   }//for all tracks
 251   SetMCReactionPlaneAngle(anInput);
 252 }
 253
 254 //-----------------------------------------------------------------------
 255 AliFlowEvent::AliFlowEvent( const AliESDEvent* anInput,
 256                             const AliCFManager* rpCFManager,
 257                             const AliCFManager* poiCFManager ):
 258   AliFlowEventSimple(20),  fApplyRecentering(-1), fCachedRun(-1), fVZEROcentralityBin(-1), fEvent(0x0)
 259 {
 260     // constructor
 261     for(Int_t i(0); i < 9; i++) {
 262         for(Int_t j(0); j < 2; j++) {
 263             for(Int_t k(0); k < 2; k++) {
 264                 fMeanQ[i][j][k] = 0.;
 265                 fWidthQ[i][j][k] = 0.;
 266                 fMeanQv3[i][j][k] = 0.;
 267                 fWidthQv3[i][j][k] = 0.;
 268             }
 269         }
 270     }
 271
 272   //Fills the event from the ESD
 273
 274   Int_t iNumberOfInputTracks = anInput->GetNumberOfTracks() ;
 275
 276   //loop over tracks
 277   for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)
 278   {
 279     AliESDtrack* pParticle = anInput->GetTrack(itrkN);   //get input particle
 280
 281     //check if pParticle passes the cuts
 282     Bool_t rpOK = kTRUE;
 283     Bool_t poiOK = kTRUE;
 284     if (rpCFManager && poiCFManager)
 285     {
 286       rpOK = ( rpCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&
 287                rpCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));
 288       poiOK = ( poiCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&
 289                 poiCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));
 290     }
 291     if (!(rpOK || poiOK)) continue;
 292
 293     //make new AliFLowTrack
 294     AliFlowTrack* pTrack = new AliFlowTrack(pParticle);
 295     pTrack->SetSource(AliFlowTrack::kFromESD);
 296
 297     //marking the particles used for int. flow:
 298     if(rpOK && rpCFManager)
 299     {
 300       pTrack->SetForRPSelection(kTRUE);
 301       IncrementNumberOfPOIs(0);
 302     }
 303     //marking the particles used for diff. flow:
 304     if(poiOK && poiCFManager)
 305     {
 306       pTrack->SetForPOISelection(kTRUE);
 307       IncrementNumberOfPOIs(1);
 308     }
 309
 310     AddTrack(pTrack);
 311   }//end of while (itrkN < iNumberOfInputTracks)
 312 }
 313
 314 //-----------------------------------------------------------------------
 315 AliFlowEvent::AliFlowEvent( const AliAODEvent* anInput,
 316                             const AliCFManager* rpCFManager,
 317                             const AliCFManager* poiCFManager):
 318   AliFlowEventSimple(20), fApplyRecentering(-1), fCachedRun(-1), fVZEROcentralityBin(-1), fEvent(0x0)
 319 {
 320     // constructor
 321     for(Int_t i(0); i < 9; i++) {
 322         for(Int_t j(0); j < 2; j++) {
 323             for(Int_t k(0); k < 2; k++) {
 324                 fMeanQ[i][j][k] = 0.;
 325                 fWidthQ[i][j][k] = 0.;
 326                 fMeanQv3[i][j][k] = 0.;
 327                 fWidthQv3[i][j][k] = 0.;
 328             }
 329         }
 330     }
 331   //Fills the event from the AOD
 332   Int_t iNumberOfInputTracks = anInput->GetNumberOfTracks() ;
 333
 334   //loop over tracks
 335   for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)
 336   {
 337     AliAODTrack* pParticle = dynamic_cast<AliAODTrack*>(anInput->GetTrack(itrkN));
 338     if(!pParticle) AliFatal("Not a standard AOD");   //get input particle
 339
 340     //check if pParticle passes the cuts
 341     Bool_t rpOK = kTRUE;
 342     Bool_t poiOK = kTRUE;
 343     if (rpCFManager && poiCFManager)
 344     {
 345       rpOK = ( rpCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&
 346                rpCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));
 347       poiOK = ( poiCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&
 348                 poiCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));
 349     }
 350     if (!(rpOK || poiOK)) continue;
 351
 352     //make new AliFlowTrack
 353     AliFlowTrack* pTrack = new AliFlowTrack(pParticle);
 354     pTrack->SetSource(AliFlowTrack::kFromAOD);
 355
 356     if (rpOK /* && rpCFManager */ ) // to be fixed - with CF managers uncommented only empty events (NULL in header files)
 357     {
 358       pTrack->SetForRPSelection(kTRUE);
 359       IncrementNumberOfPOIs(0);
 360     }
 361     if (poiOK /* && poiCFManager*/ )
 362     {
 363       pTrack->SetForPOISelection(kTRUE);
 364       IncrementNumberOfPOIs(1);
 365     }
 366     AddTrack(pTrack);
 367   }
 368
 369   //  if (iSelParticlesRP >= fMinMult && iSelParticlesRP <= fMaxMult)
 370   //  {
 371   //    if ( (++fCount % 100) == 0)
 372   //    {
 373   //      if (!fMCReactionPlaneAngle == 0) cout<<" MC Reaction Plane Angle = "<<  fMCReactionPlaneAngle << endl;
 374   //      else cout<<" MC Reaction Plane Angle = unknown "<< endl;
 375   //      cout<<" iGoodTracks = "<<iGoodTracks<<endl;
 376   //      cout<<" # of RP selected tracks = "<<iSelParticlesRP<<endl;
 377   //      cout<<" # of POI selected tracks = "<<iSelParticlesPOI<<endl;
 378   //      cout << "# " << fCount << " events processed" << endl;
 379   //    }
 380   //    return pEvent;
 381   //  }
 382   //  else
 383   //  {
 384   //    cout<<"Not enough tracks in the FlowEventSimple"<<endl;
 385   //    return 0;
 386   //  }
 387   //}
 388   //else
 389   //{
 390   //  cout<<"Event does not pass multiplicity cuts"<<endl;
 391   //  return 0;
 392   //}
 393
 394 }
 395
 396 //-----------------------------------------------------------------------
 397 AliFlowEvent::AliFlowEvent( const AliESDEvent* anInput,
 398                             const AliMCEvent* anInputMc,
 399                             KineSource anOption,
 400                             const AliCFManager* rpCFManager,
 401                             const AliCFManager* poiCFManager ):
 402   AliFlowEventSimple(20), fApplyRecentering(-1), fCachedRun(-1), fVZEROcentralityBin(-1), fEvent(0x0)
 403 {
 404     // constructor
 405     for(Int_t i(0); i < 9; i++) {
 406         for(Int_t j(0); j < 2; j++) {
 407             for(Int_t k(0); k < 2; k++) {
 408                 fMeanQ[i][j][k] = 0.;
 409                 fWidthQ[i][j][k] = 0.;
 410                 fMeanQv3[i][j][k] = 0.;
 411                 fWidthQv3[i][j][k] = 0.;
 412             }
 413         }
 414     }
 415   //fills the event with tracks from the ESD and kinematics from the MC info via the track label
 416   if (anOption==kNoKine)
 417   {
 418     AliFatal("WRONG OPTION IN AliFlowEventMaker::FillTracks(AliESDEvent* anInput, AliMCEvent* anInputMc, KineSource anOption)");
 419     exit(1);
 420   }
 421
 422   Int_t iNumberOfInputTracks = anInput->GetNumberOfTracks() ;
 423
 424   Int_t iNumberOfInputTracksMC = anInputMc->GetNumberOfTracks() ;
 425   if (iNumberOfInputTracksMC==-1)
 426   {
 427     AliError("Skipping Event -- No MC information available for this event");
 428     return;
 429   }
 430
 431   //loop over ESD tracks
 432   for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)
 433   {
 434     AliESDtrack* pParticle = anInput->GetTrack(itrkN);   //get input particle
 435     //get Label
 436     Int_t iLabel = pParticle->GetLabel();
 437     //match to mc particle
 438     AliMCParticle* pMcParticle = (AliMCParticle*) anInputMc->GetTrack(TMath::Abs(iLabel));
 439
 440     //check
 441     if (TMath::Abs(pParticle->GetLabel())!=pMcParticle->Label())
 442       AliWarning(Form("pParticle->GetLabel()!=pMcParticle->Label(), %i, %i", pParticle->GetLabel(), pMcParticle->Label()));
 443
 444     //check if pParticle passes the cuts
 445     Bool_t rpOK = kFALSE;
 446     Bool_t poiOK = kFALSE;
 447     if (rpCFManager && poiCFManager)
 448     {
 449       if(anOption == kESDkine)
 450       {
 451         if (rpCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,pMcParticle,"mcGenCuts1") &&
 452             rpCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle))
 453           rpOK=kTRUE;
 454         if (poiCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,pMcParticle,"mcGenCuts2") &&
 455             poiCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle))
 456           poiOK=kTRUE;
 457       }
 458       else if (anOption == kMCkine)
 459       {
 460         if (rpCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,pMcParticle))
 461           rpOK=kTRUE;
 462         if (poiCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,pMcParticle))
 463           poiOK=kTRUE;
 464       }
 465     }
 466
 467     if (!(rpOK || poiOK)) continue;
 468
 469     //make new AliFlowTrack
 470     AliFlowTrack* pTrack = NULL;
 471     if(anOption == kESDkine)   //take the PID from the MC & the kinematics from the ESD
 472     {
 473       pTrack = new AliFlowTrack(pParticle);
 474     }
 475     else if (anOption == kMCkine)   //take the PID and kinematics from the MC
 476     {
 477       pTrack = new AliFlowTrack(pMcParticle);
 478     }
 479
 480     if (rpOK && rpCFManager)
 481     {
 482       IncrementNumberOfPOIs(0);
 483       pTrack->SetForRPSelection();
 484     }
 485     if (poiOK && poiCFManager)
 486     {
 487       IncrementNumberOfPOIs(1);
 488       pTrack->SetForPOISelection();
 489     }
 490
 491     AddTrack(pTrack);
 492   }
 493   SetMCReactionPlaneAngle(anInputMc);
 494 }
 495
 496 //-----------------------------------------------------------------------
 497 AliFlowEvent::AliFlowEvent( const AliESDEvent* anInput,
 498                             const AliMultiplicity* anInputTracklets,
 499                             const AliCFManager* poiCFManager ):
 500   AliFlowEventSimple(20), fApplyRecentering(-1), fCachedRun(-1), fVZEROcentralityBin(-1), fEvent(0x0)
 501 {
 502     // constructor
 503     for(Int_t i(0); i < 9; i++) {
 504         for(Int_t j(0); j < 2; j++) {
 505             for(Int_t k(0); k < 2; k++) {
 506                 fMeanQ[i][j][k] = 0.;
 507                 fWidthQ[i][j][k] = 0.;
 508                 fMeanQv3[i][j][k] = 0.;
 509                 fWidthQv3[i][j][k] = 0.;
 510             }
 511         }
 512     }
 513
 514   //Select the particles of interest from the ESD
 515   Int_t iNumberOfInputTracks = anInput->GetNumberOfTracks() ;
 516
 517   //loop over tracks
 518   for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)
 519     {
 520       AliESDtrack* pParticle = anInput->GetTrack(itrkN);   //get input particle
 521
 522       //check if pParticle passes the cuts
 523       Bool_t poiOK = kTRUE;
 524       if (poiCFManager)
 525         {
 526           poiOK = ( poiCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&
 527                     poiCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));
 528         }
 529       if (!poiOK) continue;
 530
 531       //make new AliFLowTrack
 532       AliFlowTrack* pTrack = new AliFlowTrack(pParticle);
 533
 534       //marking the particles used for the particle of interest (POI) selection:
 535       if(poiOK && poiCFManager)
 536         {
 537           IncrementNumberOfPOIs(1);
 538           pTrack->SetForPOISelection(kTRUE);
 539           pTrack->SetSource(AliFlowTrack::kFromESD);
 540         }
 541
 542       AddTrack(pTrack);
 543     }//end of while (itrkN < iNumberOfInputTracks)
 544
 545   //Select the reference particles from the SPD tracklets
 546   anInputTracklets = anInput->GetMultiplicity();
 547   Int_t multSPD = anInputTracklets->GetNumberOfTracklets();
 548
 549   //loop over tracklets
 550   for (Int_t itracklet=0; itracklet<multSPD; ++itracklet) {
 551     Float_t thetaTr= anInputTracklets->GetTheta(itracklet);
 552     Float_t phiTr= anInputTracklets->GetPhi(itracklet);
 553     // calculate eta
 554     Float_t etaTr = -TMath::Log(TMath::Tan(thetaTr/2.));
 555
 556     //make new AliFLowTrackSimple
 557     AliFlowTrack* pTrack = new AliFlowTrack();
 558     pTrack->SetPt(0.0);
 559     pTrack->SetEta(etaTr);
 560     pTrack->SetPhi(phiTr);
 561     //marking the particles used for the reference particle (RP) selection:
 562     IncrementNumberOfPOIs(0);
 563     pTrack->SetForRPSelection(kTRUE);
 564     pTrack->SetSource(AliFlowTrack::kFromTracklet);
 565
 566     //Add the track to the flowevent
 567     AddTrack(pTrack);
 568   }
 569
 570 }
 571
 572 //-----------------------------------------------------------------------
 573 AliFlowEvent::AliFlowEvent( const AliESDEvent* esd,
 574                             const AliCFManager* poiCFManager,
 575                             Bool_t hybrid):
 576   AliFlowEventSimple(20), fApplyRecentering(-1), fCachedRun(-1), fVZEROcentralityBin(-1), fEvent(0x0)
 577 {
 578     // constructor
 579     for(Int_t i(0); i < 9; i++) {
 580         for(Int_t j(0); j < 2; j++) {
 581             for(Int_t k(0); k < 2; k++) {
 582                 fMeanQ[i][j][k] = 0.;
 583                 fWidthQ[i][j][k] = 0.;
 584                 fMeanQv3[i][j][k] = 0.;
 585                 fWidthQv3[i][j][k] = 0.;
 586             }
 587         }
 588     }
 589
 590   //Select the particles of interest from the ESD
 591   Int_t iNumberOfInputTracks = esd->GetNumberOfTracks() ;
 592
 593   //Double_t gPt = 0.0, gP = 0.0;
 594   Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};  //The impact parameters and their covariance.
 595 //  Double_t dca3D = 0.0;       FIXME unused variable
 596
 597   AliESDtrack trackTPC;
 598
 599   //loop over tracks
 600   for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)
 601     {
 602
 603       if (!esd->GetTrack(itrkN)) continue;
 604
 605       Bool_t useTPC = kFALSE;
 606
 607       AliESDtrack* pParticle = esd->GetTrack(itrkN);   //get input particle
 608
 609       //check if pParticle passes the cuts
 610       Bool_t poiOK = kTRUE;
 611
 612       if (poiCFManager)
 613       {
 614         poiOK = ( poiCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&
 615                   poiCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));
 616       }
 617
 618       if (!(poiOK)) continue;
 619
 620       AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)pParticle->GetTPCInnerParam();
 621
 622       if (tpcTrack)
 623       {
 624
 625 //      gPt = tpcTrack->Pt();
 626 //      gP = tpcTrack->P();
 627
 628         useTPC = kTRUE;
 629
 630         const AliESDVertex *vertexSPD = esd->GetPrimaryVertexSPD();
 631         const AliESDVertex *vertexTPC = esd->GetPrimaryVertexTPC();
 632
 633         AliExternalTrackParam copy(*tpcTrack);
 634         if(hybrid)
 635           copy.PropagateToDCA(vertexSPD,esd->GetMagneticField(),100.,dca,cov);
 636         else
 637           copy.PropagateToDCA(vertexTPC,esd->GetMagneticField(),100.,dca,cov);
 638
 639 //        dca3D = TMath::Sqrt(TMath::Power(dca[0],2)+TMath::Power(dca[1],2));   FIXME unused variable
 640
 641       }
 642
 643       //make new AliFLowTrack
 644       AliFlowTrack* pTrack = new AliFlowTrack(pParticle);
 645
 646       pTrack->SetSource(AliFlowTrack::kFromESD);
 647
 648       //marking the particles used for diff. flow:
 649       if(poiOK && poiCFManager)
 650       {
 651         pTrack->SetForPOISelection(kTRUE);
 652         IncrementNumberOfPOIs(1);
 653       }
 654
 655       if(useTPC)
 656       {
 657         pTrack->SetForRPSelection(kTRUE);
 658         IncrementNumberOfPOIs(0);
 659       }
 660
 661       AddTrack(pTrack);
 662
 663     }//end of while (itrkN < iNumberOfInputTracks)
 664
 665 }
 666
 667 //-----------------------------------------------------------------------
 668 AliFlowEvent::AliFlowEvent( const AliESDEvent* anInput,
 669                             const TH2F* anInputFMDhist,
 670                             const AliCFManager* poiCFManager ):
 671   AliFlowEventSimple(20), fApplyRecentering(-1), fCachedRun(-1), fVZEROcentralityBin(-1), fEvent(0x0)
 672 {
 673     // constructor
 674     for(Int_t i(0); i < 9; i++) {
 675         for(Int_t j(0); j < 2; j++) {
 676             for(Int_t k(0); k < 2; k++) {
 677                 fMeanQ[i][j][k] = 0.;
 678                 fWidthQ[i][j][k] = 0.;
 679                 fMeanQv3[i][j][k] = 0.;
 680                 fWidthQv3[i][j][k] = 0.;
 681             }
 682         }
 683     }
 684
 685   //Select the particles of interest from the ESD
 686   Int_t iNumberOfInputTracks = anInput->GetNumberOfTracks() ;
 687
 688   //loop over tracks
 689   for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)
 690     {
 691       AliESDtrack* pParticle = anInput->GetTrack(itrkN);   //get input particle
 692
 693       //check if pParticle passes the cuts
 694       Bool_t poiOK = kTRUE;
 695       if (poiCFManager)
 696         {
 697           poiOK = ( poiCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&
 698                     poiCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));
 699         }
 700       if (!poiOK) continue;
 701
 702       //make new AliFLowTrack
 703       AliFlowTrack* pTrack = new AliFlowTrack(pParticle);
 704
 705       //marking the particles used for the particle of interest (POI) selection:
 706       if(poiOK && poiCFManager)
 707         {
 708     IncrementNumberOfPOIs(1);
 709           pTrack->SetForPOISelection(kTRUE);
 710           pTrack->SetSource(AliFlowTrack::kFromESD);
 711         }
 712
 713       AddTrack(pTrack);
 714     }//end of while (itrkN < iNumberOfInputTracks)
 715
 716   //Select the reference particles from the FMD hits
 717   //loop over FMD histogram
 718   Int_t iBinsEta = anInputFMDhist->GetNbinsX();
 719   Int_t iBinsPhi = anInputFMDhist->GetNbinsY();
 720
 721   for (Int_t iEta = 1; iEta <= iBinsEta; iEta++){
 722     Double_t etaFMD = anInputFMDhist->GetXaxis()->GetBinCenter(iEta);
 723     for (Int_t iPhi = 1; iPhi <= iBinsPhi; iPhi++){
 724       Double_t phiFMD = anInputFMDhist->GetYaxis()->GetBinCenter(iPhi);
 725       Double_t weightFMD = anInputFMDhist->GetBinContent(iEta,iPhi);
 726
 727       if (weightFMD > 0.0) { //do not add empty bins
 728         //make new AliFLowTrackSimple
 729         AliFlowTrack* pTrack = new AliFlowTrack();
 730         pTrack->SetPt(0.0);
 731         pTrack->SetEta(etaFMD);
 732         pTrack->SetPhi(phiFMD);
 733         pTrack->SetWeight(weightFMD);
 734         //marking the particles used for the reference particle (RP) selection:
 735         pTrack->TagRP();
 736         IncrementNumberOfPOIs(0);
 737         pTrack->SetSource(AliFlowTrack::kFromFMD);
 738
 739         //Add the track to the flowevent
 740         AddTrack(pTrack);
 741
 742       }
 743     }
 744   }
 745
 746 }
 747
 748 //-----------------------------------------------------------------------
 749 void AliFlowEvent::FindDaughters(Bool_t keepDaughtersInRPselection)
 750 {
 751   //each flow track holds it's esd track index as well as its daughters esd index.
 752   //fill the array of daughters for every track with the pointers to flow tracks
 753   //to associate the mothers with daughters directly
 754   for (Int_t iTrack=0; iTrack<fMothersCollection->GetEntriesFast(); iTrack++)
 755   {
 756     AliFlowTrack* mother = static_cast<AliFlowTrack*>(fMothersCollection->At(iTrack));
 757     if (!mother) continue;
 758     if (mother->GetNDaughters()<1) continue;
 759     for (Int_t iDaughterCandidate=0; iDaughterCandidate<fNumberOfTracks; iDaughterCandidate++)
 760     {
 761       AliFlowTrack* daughterCandidate = static_cast<AliFlowTrack*>(fTrackCollection->At(iDaughterCandidate));
 762       Int_t esdIndexDaughterCandidate = daughterCandidate->GetID();
 763       for (Int_t iDaughter=0; iDaughter<mother->GetNDaughters(); iDaughter++)
 764       {
 765         Int_t esdIndexDaughter = mother->GetIDDaughter(iDaughter);
 766         if (esdIndexDaughter==esdIndexDaughterCandidate)
 767         {
 768           mother->SetDaughter(iDaughter,daughterCandidate);
 769           daughterCandidate->SetForRPSelection(keepDaughtersInRPselection);
 770         }
 771       }
 772     }
 773   }
 774 }
 775
 776 //-----------------------------------------------------------------------
 777 void AliFlowEvent::Fill( AliFlowTrackCuts* rpCuts,
 778                          AliFlowTrackCuts* poiCuts )
 779 {
 780   //Fills the event from a vevent: AliESDEvent,AliAODEvent,AliMCEvent
 781   //the input data needs to be attached to the cuts
 782   //we have two cases, if we're cutting the same collection of tracks
 783   //(same param type) then we can have tracks that are both rp and poi
 784   //in the other case we want to have two exclusive sets of rps and pois
 785   //e.g. one tracklets, the other PMD or global - USER IS RESPOSIBLE
 786   //FOR MAKING SURE THEY DONT OVERLAP OR ELSE THE SAME PARTICLE WILL BE
 787   //TAKEN TWICE
 788
 789   ClearFast();
 790   if (!rpCuts || !poiCuts) return;
 791   AliFlowTrackCuts::trackParameterType sourceRP = rpCuts->GetParamType();
 792   AliFlowTrackCuts::trackParameterType sourcePOI = poiCuts->GetParamType();
 793   AliFlowTrack* pTrack=NULL;
 794
 795   // if the source for rp's or poi's is the VZERO detector, get the calibration
 796   // and set the calibration parameters
 797   if (sourceRP == AliFlowTrackCuts::kVZERO) {
 798       SetVZEROCalibrationForTrackCuts(rpCuts);
 799       if(!rpCuts->GetApplyRecentering()) {
 800           // if the user does not want to recenter, switch the flag
 801           fApplyRecentering = -1;
 802       }
 803       // note: this flag is used in the overloaded implementation of Get2Qsub()
 804       // and tells the function to use as Qsub vectors the recentered Q-vectors
 805       // from the VZERO oadb file or from the event header
 806   }
 807   if (sourcePOI == AliFlowTrackCuts::kVZERO) {
 808       // probably no-one will choose vzero tracks as poi's ...
 809       SetVZEROCalibrationForTrackCuts(poiCuts);
 810   }
 811
 812
 813   if (sourceRP==sourcePOI)
 814   {
 815     //loop over tracks
 816     Int_t numberOfInputObjects = rpCuts->GetNumberOfInputObjects();
 817     for (Int_t i=0; i<numberOfInputObjects; i++)
 818     {
 819       //get input object (particle)
 820       TObject* particle = rpCuts->GetInputObject(i);
 821
 822       Bool_t rp = rpCuts->IsSelected(particle,i);
 823       Bool_t poi = poiCuts->IsSelected(particle,i);
 824
 825       if (!(rp||poi)) continue;
 826
 827       //make new AliFlowTrack
 828       if (rp)
 829       {
 830         pTrack = rpCuts->FillFlowTrack(fTrackCollection,fNumberOfTracks);
 831         if (!pTrack) continue;
 832         pTrack->Tag(0); IncrementNumberOfPOIs(0);
 833         if (poi) {pTrack->Tag(1); IncrementNumberOfPOIs(1);}
 834         if (pTrack->GetNDaughters()>0) fMothersCollection->Add(pTrack);
 835       }
 836       else if (poi)
 837       {
 838         pTrack = poiCuts->FillFlowTrack(fTrackCollection,fNumberOfTracks);
 839         if (!pTrack) continue;
 840         pTrack->Tag(1); IncrementNumberOfPOIs(1);
 841         if (pTrack->GetNDaughters()>0) fMothersCollection->Add(pTrack);
 842       }
 843       fNumberOfTracks++;
 844     }//end of while (i < numberOfTracks)
 845   }
 846   else if (sourceRP!=sourcePOI)
 847   {
 848     //here we have two different sources of particles, so we fill
 849     //them independently
 850     //POI
 851     for (Int_t i=0; i<poiCuts->GetNumberOfInputObjects(); i++)
 852     {
 853       TObject* particle = poiCuts->GetInputObject(i);
 854       Bool_t poi = poiCuts->IsSelected(particle,i);
 855       if (!poi) continue;
 856       pTrack = poiCuts->FillFlowTrack(fTrackCollection,fNumberOfTracks);
 857       if (!pTrack) continue;
 858       pTrack->Tag(1);
 859       IncrementNumberOfPOIs(1);
 860       fNumberOfTracks++;
 861       if (pTrack->GetNDaughters()>0) fMothersCollection->Add(pTrack);
 862     }
 863     //RP
 864     Int_t numberOfInputObjects = rpCuts->GetNumberOfInputObjects();
 865     for (Int_t i=0; i<numberOfInputObjects; i++)
 866       {
 867       TObject* particle = rpCuts->GetInputObject(i);
 868       Bool_t rp = rpCuts->IsSelected(particle,i);
 869       if (!rp) continue;
 870       pTrack = rpCuts->FillFlowTrack(fTrackCollection,fNumberOfTracks);
 871       if (!pTrack) continue;
 872       pTrack->Tag(0);
 873       IncrementNumberOfPOIs(0);
 874       fNumberOfTracks++;
 875       if (pTrack->GetNDaughters()>0) fMothersCollection->Add(pTrack);
 876     }
 877   }
 878 }
 879
 880 //-----------------------------------------------------------------------
 881 void AliFlowEvent::InsertTrack(AliFlowTrack *track) {
 882   // adds a flow track at the end of the container
 883   AliFlowTrack *pTrack = ReuseTrack( fNumberOfTracks++ );
 884   *pTrack = *track;
 885   if (track->GetNDaughters()>0)
 886   {
 887     fMothersCollection->Add(pTrack);
 888   }
 889   return;
 890 }
 891
 892 //-----------------------------------------------------------------------
 893 AliFlowTrack* AliFlowEvent::ReuseTrack(Int_t i)
 894 {
 895   //try to reuse an existing track, if empty, make new one
 896   AliFlowTrack* pTrack = static_cast<AliFlowTrack*>(fTrackCollection->At(i));
 897   if (pTrack)
 898   {
 899     pTrack->Clear();
 900   }
 901   else
 902   {
 903     pTrack = new AliFlowTrack();
 904     fTrackCollection->AddAtAndExpand(pTrack,i);
 905   }
 906   return pTrack;
 907 }
 908
 909 //-----------------------------------------------------------------------
 910 AliFlowEvent::AliFlowEvent( AliFlowTrackCuts* rpCuts,
 911                             AliFlowTrackCuts* poiCuts ):
 912   AliFlowEventSimple(20), fApplyRecentering(kFALSE), fCachedRun(-1), fVZEROcentralityBin(-1), fEvent(0x0)
 913 {
 914     // constructor
 915     for(Int_t i(0); i < 9; i++) {
 916         for(Int_t j(0); j < 2; j++) {
 917             for(Int_t k(0); k < 2; k++) {
 918                 fMeanQ[i][j][k] = 0.;
 919                 fWidthQ[i][j][k] = 0.;
 920                 fMeanQv3[i][j][k] = 0.;
 921                 fWidthQv3[i][j][k] = 0.;
 922             }
 923         }
 924     }
 925   //Fills the event from a vevent: AliESDEvent,AliAODEvent,AliMCEvent
 926   //the input data needs to be attached to the cuts
 927   Fill(rpCuts,poiCuts);
 928 }
 929
 930 //-------------------------------------------------------------------//
 931 //---- Including PMD tracks as RP --------------------------//
 932
 933 AliFlowEvent::AliFlowEvent( const AliESDEvent* anInput,
 934                             const AliESDPmdTrack *pmdtracks,
 935                             const AliCFManager* poiCFManager ):
 936   AliFlowEventSimple(20), fApplyRecentering(kFALSE), fCachedRun(-1), fVZEROcentralityBin(-1), fEvent(0x0)
 937 {
 938     // constructor
 939     for(Int_t i(0); i < 9; i++) {
 940         for(Int_t j(0); j < 2; j++) {
 941             for(Int_t k(0); k < 2; k++) {
 942                 fMeanQ[i][j][k] = 0.;
 943                 fWidthQ[i][j][k] = 0.;
 944                 fMeanQv3[i][j][k] = 0.;
 945                 fWidthQv3[i][j][k] = 0.;
 946             }
 947         }
 948     }
 949   Float_t GetPmdEta(Float_t xPos, Float_t yPos, Float_t zPos);
 950   Float_t GetPmdPhi(Float_t xPos, Float_t yPos);
 951   //Select the particles of interest from the ESD
 952   Int_t iNumberOfInputTracks = anInput->GetNumberOfTracks() ;
 953
 954   //loop over tracks
 955   for (Int_t itrkN=0; itrkN<iNumberOfInputTracks; itrkN++)
 956     {
 957       AliESDtrack* pParticle = anInput->GetTrack(itrkN);   //get input particle
 958       //check if pParticle passes the cuts
 959       Bool_t poiOK = kTRUE;
 960       if (poiCFManager)
 961         {
 962           poiOK = ( poiCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,pParticle) &&
 963                     poiCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,pParticle));
 964         }
 965       if (!poiOK) continue;
 966
 967       //make new AliFLowTrack
 968       AliFlowTrack* pTrack = new AliFlowTrack(pParticle);
 969
 970       //marking the particles used for the particle of interest (POI) selection:
 971       if(poiOK && poiCFManager)
 972         {
 973           IncrementNumberOfPOIs(1);
 974           pTrack->SetForPOISelection(kTRUE);
 975           pTrack->SetSource(AliFlowTrack::kFromESD);
 976         }
 977
 978       AddTrack(pTrack);
 979     }//end of while (itrkN < iNumberOfInputTracks)
 980
 981   //Select the reference particles from the PMD tracks
 982   Int_t npmdcl = anInput->GetNumberOfPmdTracks();
 983   printf("======There are %d PMD tracks in this event\n-------",npmdcl);
 984   //loop over clusters
 985   for(Int_t iclust=0; iclust < npmdcl; iclust++){
 986     //AliESDPmdTrack *pmdtr = anInput->GetPmdTrack(iclust);
 987     pmdtracks = anInput->GetPmdTrack(iclust);
 988     Int_t   det   = pmdtracks->GetDetector();
 989     //Int_t   smn   = pmdtracks->GetSmn();
 990     Float_t clsX  = pmdtracks->GetClusterX();
 991     Float_t clsY  = pmdtracks->GetClusterY();
 992     Float_t clsZ  = pmdtracks->GetClusterZ();
 993     Float_t ncell = pmdtracks->GetClusterCells();
 994     Float_t adc   = pmdtracks->GetClusterADC();
 995     //Float_t pid   = pmdtracks->GetClusterPID();
 996     Float_t etacls = GetPmdEta(clsX,clsY,clsZ);
 997     Float_t phicls = GetPmdPhi(clsX,clsY);
 998     //make new AliFLowTrackSimple
 999     AliFlowTrack* pTrack = new AliFlowTrack();
1000     //if(det == 0){ //selecting preshower plane only
1001     if(det == 0 && adc > 270 && ncell > 1){ //selecting preshower plane only
1002       //pTrack->SetPt(adc);//cluster adc
1003       pTrack->SetPt(0.0);
1004       pTrack->SetEta(etacls);
1005       pTrack->SetPhi(phicls);
1006       //marking the particles used for the reference particle (RP) selection:
1007       IncrementNumberOfPOIs(0);
1008       pTrack->SetForRPSelection(kTRUE);
1009       pTrack->SetSource(AliFlowTrack::kFromPMD);
1010       //Add the track to the flowevent
1011       AddTrack(pTrack);
1012     }//if det
1013   }
1014 }
1015 //----------------------------------------------------------------------------//
1016 Float_t GetPmdEta(Float_t xPos, Float_t yPos, Float_t zPos)
1017 {
1018   Float_t rpxpy, theta, eta;
1019   rpxpy  = TMath::Sqrt(xPos*xPos + yPos*yPos);
1020   theta  = TMath::ATan2(rpxpy,zPos);
1021   eta    = -TMath::Log(TMath::Tan(0.5*theta));
1022   return eta;
1023 }
1024 //--------------------------------------------------------------------------//
1025 Float_t GetPmdPhi(Float_t xPos, Float_t yPos)
1026 {
1027   Float_t pybypx, phi = 0., phi1;
1028   if(xPos==0)
1029     {
1030       if(yPos>0) phi = 90.;
1031       if(yPos<0) phi = 270.;
1032     }
1033   if(xPos != 0)
1034     {
1035       pybypx = yPos/xPos;
1036       if(pybypx < 0) pybypx = - pybypx;
1037       phi1 = TMath::ATan(pybypx)*180./3.14159;
1038
1039       if(xPos > 0 && yPos > 0) phi = phi1;        // 1st Quadrant
1040       if(xPos < 0 && yPos > 0) phi = 180 - phi1;  // 2nd Quadrant
1041       if(xPos < 0 && yPos < 0) phi = 180 + phi1;  // 3rd Quadrant
1042       if(xPos > 0 && yPos < 0) phi = 360 - phi1;  // 4th Quadrant
1043
1044     }
1045   phi = phi*3.14159/180.;
1046   return   phi;
1047 }
1048 //---------------------------------------------------------------//
1049
1050 void AliFlowEvent::Get2Qsub(AliFlowVector* Qarray, Int_t n, TList *weightsList, Bool_t usePhiWeights, Bool_t usePtWeights, Bool_t useEtaWeights)
1051 {
1052   // get q vectors for the subevents. if no recentering is necessary, get the guy from the flow event simple
1053   AliFlowEventSimple::Get2Qsub(Qarray, n, weightsList, usePhiWeights, usePtWeights, useEtaWeights);
1054   AliFlowVector vA = Qarray[0];
1055   AliFlowVector vB = Qarray[1];
1056
1057   // else get the recentering from the cached info
1058   if (fApplyRecentering == 2010)        // 10h style recentering, implemented for n=2 and n=3
1059   {
1060     // first retrieve the q-vectors from the AliFlowEventSimple:: routine
1061     // extract the information form the current flow vectors
1062     Double_t Qxc(vA.X());       // IMPORTANT: user is responsible for the sign of eta
1063     Double_t Qyc(vA.Y());       // vzeroC has negative pseudorapidity and is taken as subevent A
1064     Double_t Qxa(vB.X());       // vzeroA has positive pseudorapidity and is taken as subevent B
1065     Double_t Qya(vB.Y());
1066     // init some values for the corrections
1067
1068     // default values for vector a (VZEROA)
1069     Double_t Qxamean(0);
1070     Double_t Qxarms(1);
1071     Double_t Qyamean(0);
1072     Double_t Qyarms(1);
1073     // default values for vector b (VZEROC)
1074     Double_t Qxcmean(0);
1075     Double_t Qxcrms(1);
1076     Double_t Qycmean(0);
1077     Double_t Qycrms(1);
1078     // note that defaults are chosen such that for n!=2||n!=3 re-centering is a null-operation
1079
1080     if( n == 2) {       // second order symmetry
1081         Qxamean = fMeanQ[fVZEROcentralityBin][1][0];
1082         Qxarms  = fWidthQ[fVZEROcentralityBin][1][0];
1083         Qyamean = fMeanQ[fVZEROcentralityBin][1][1];
1084         Qyarms  = fWidthQ[fVZEROcentralityBin][1][1];
1085
1086         Qxcmean = fMeanQ[fVZEROcentralityBin][0][0];
1087         Qxcrms  = fWidthQ[fVZEROcentralityBin][0][0];
1088         Qycmean = fMeanQ[fVZEROcentralityBin][0][1];
1089         Qycrms  = fWidthQ[fVZEROcentralityBin][0][1];
1090     } else if (n == 3) {        // third order symmetry
1091         Qxamean = fMeanQv3[fVZEROcentralityBin][1][0];
1092         Qxarms  = fWidthQv3[fVZEROcentralityBin][1][0];
1093         Qyamean = fMeanQv3[fVZEROcentralityBin][1][1];
1094         Qyarms  = fWidthQv3[fVZEROcentralityBin][1][1];
1095
1096         Qxcmean = fMeanQv3[fVZEROcentralityBin][0][0];
1097         Qxcrms  = fWidthQv3[fVZEROcentralityBin][0][0];
1098         Qycmean = fMeanQv3[fVZEROcentralityBin][0][1];
1099         Qycrms  = fWidthQv3[fVZEROcentralityBin][0][1];
1100     }
1101 //    printf(" > n %i Qx %.2f Qxrms %.2f \n", n, Qxamean, Qxarms);
1102     // do the correction
1103     Double_t QxaCor = (Qxa - Qxamean)/Qxarms;
1104     Double_t QyaCor = (Qya - Qyamean)/Qyarms;
1105     Double_t QxcCor = (Qxc - Qxcmean)/Qxcrms;
1106     Double_t QycCor = (Qyc - Qycmean)/Qycrms;
1107     // update the vector
1108     vA.Set(QxcCor, QycCor);
1109     vB.Set(QxaCor, QyaCor);
1110   } else if (fApplyRecentering == 2011) { // 11h style recentering
1111
1112     Double_t QxaCor = 0.;
1113     Double_t QyaCor = 0.;
1114     Double_t QxcCor = 0.;
1115     Double_t QycCor = 0.;
1116
1117     if(fEvent && fEvent->GetEventplane()) {
1118       fEvent->GetEventplane()->CalculateVZEROEventPlane(fEvent, 8, n, QxaCor, QyaCor);
1119       fEvent->GetEventplane()->CalculateVZEROEventPlane(fEvent, 9, n, QxcCor, QycCor);
1120
1121       // set the new q-vectors (which in this case means REPLACING)
1122       vA.Set(QxcCor, QycCor);
1123       vB.Set(QxaCor, QyaCor);
1124     } // if for some reason the info from the event header is not available, the AliFlowTrackCuts object
1125       // has provided the equalized multiplicity info so this should still be relatively safe
1126   }
1127 }
1128 //_____________________________________________________________________________
1129 void AliFlowEvent::SetVZEROCalibrationForTrackCuts(AliFlowTrackCuts* cuts) {
1130     // open calibration info, copied from AliAnalyisTaskVnV0.cxx
1131     fEvent = cuts->GetEvent();
1132     if(!fEvent) return; // coverity. we need to know the event to get the runnumber and centrlaity
1133     // get the vzero centrality percentile (cc dependent calibration)
1134     Float_t v0Centr(fEvent->GetCentrality()->GetCentralityPercentile("V0M"));
1135     if(v0Centr < 5) fVZEROcentralityBin = 0;
1136     else if(v0Centr < 10) fVZEROcentralityBin = 1;
1137     else if(v0Centr < 20) fVZEROcentralityBin = 2;
1138     else if(v0Centr < 30) fVZEROcentralityBin = 3;
1139     else if(v0Centr < 40) fVZEROcentralityBin = 4;
1140     else if(v0Centr < 50) fVZEROcentralityBin = 5;
1141     else if(v0Centr < 60) fVZEROcentralityBin = 6;
1142     else if(v0Centr < 70) fVZEROcentralityBin = 7;
1143     else fVZEROcentralityBin = 8;
1144
1145     // if this event is from the same run as the previous event
1146     // we can use the cached calibration values, no need to re-open the
1147     // aodb file, else cache the new run
1148     Int_t run(fEvent->GetRunNumber());
1149     if(fCachedRun == run) return;
1150     else fCachedRun = run;
1151
1152     TFile *foadb = TFile::Open("$ALICE_ROOT/OADB/PWGCF/VZERO/VZEROcalibEP.root");
1153     if(!foadb){
1154         printf("OADB file $ALICE_ROOT/OADB/PWGCF/VZERO/VZEROcalibEP.root cannot be opened, CALIBRATION FAILED !");
1155         return;
1156     }
1157
1158     AliOADBContainer *cont = (AliOADBContainer*) foadb->Get("hMultV0BefCorr");
1159     if(!cont){
1160         printf("OADB object hMultV0BefCorr is not available in the file\n");
1161         return;
1162     }
1163     if(!(cont->GetObject(run))){
1164         // if the multiplicity correction cannot be found for the specified run,
1165         // loop over the 11h runs to see if it's 11h data
1166         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};
1167         for(Int_t r(0); r < 176; r++) {
1168             if(run == runs11h[r]) {
1169                 printf(" > run has been identified as 11h < \n");
1170  if(cuts->GetVZEROgainEqualizationPerRing()) {
1171                     // enable or disable rings through the weights, weight 1. is enabled, 0. is disabled
1172                     // start with the vzero c rings (segments 0 through 31)
1173                     (cuts->GetUseVZERORing(0)) ? cuts->SetVZEROCpol(0, 1.) : cuts->SetVZEROCpol(0, 0.);
1174                     (cuts->GetUseVZERORing(1)) ? cuts->SetVZEROCpol(1, 1.) : cuts->SetVZEROCpol(1, 0.);
1175                     (cuts->GetUseVZERORing(2)) ? cuts->SetVZEROCpol(2, 1.) : cuts->SetVZEROCpol(2, 0.);
1176                     (cuts->GetUseVZERORing(3)) ? cuts->SetVZEROCpol(3, 1.) : cuts->SetVZEROCpol(3, 0.);
1177                     // same for vzero a
1178                     (cuts->GetUseVZERORing(4)) ? cuts->SetVZEROApol(0, 1.) : cuts->SetVZEROApol(0, 0.);
1179                     (cuts->GetUseVZERORing(5)) ? cuts->SetVZEROApol(1, 1.) : cuts->SetVZEROApol(1, 0.);
1180                     (cuts->GetUseVZERORing(6)) ? cuts->SetVZEROApol(2, 1.) : cuts->SetVZEROApol(2, 0.);
1181                     (cuts->GetUseVZERORing(7)) ? cuts->SetVZEROApol(3, 1.) : cuts->SetVZEROApol(3, 0.);
1182                 } else {
1183                     // else enable all rings, which is also default
1184                     for(Int_t i(0); i < 4; i++) cuts->SetVZEROCpol(i, 1.);
1185                     for(Int_t i(0); i < 4; i++) cuts->SetVZEROApol(i, 1.);
1186                 }
1187                 // pass a NULL pointer to the track cuts object, the NULL pointer will identify 11h runs
1188                 cuts->SetVZEROgainEqualisation(NULL);
1189                 fApplyRecentering = 2011;
1190                 return; // the rest of the steps are not necessary
1191             }
1192         }
1193         // the run has not been identified as lhc11h data, so we assume a template calibration
1194         printf("OADB object hMultVZEROBefCorr is not available for run %i (used default run 137366)\n",run);
1195         run = 137366;
1196     }
1197     printf(" > run has been identified as 10h < \n");
1198     // step 1) get the proper multiplicity weights from the vzero signal
1199     TProfile* fMultVZERO = ((TH2F *) cont->GetObject(run))->ProfileX();
1200
1201     TF1 *fpol0 = new TF1("fpol0","pol0");
1202     if(cuts->GetVZEROgainEqualizationPerRing()) {
1203         // do the calibration per ring
1204         // start with the vzero c rings (segments 0 through 31)
1205         fMultVZERO->Fit(fpol0, "", "", 0, 8);
1206         (cuts->GetUseVZERORing(0)) ? cuts->SetVZEROCpol(0, fpol0->GetParameter(0)) : cuts->SetVZEROCpol(0, 0.);
1207         fMultVZERO->Fit(fpol0, "", "", 8, 16);
1208         (cuts->GetUseVZERORing(1)) ? cuts->SetVZEROCpol(1, fpol0->GetParameter(0)) : cuts->SetVZEROCpol(1, 0.);
1209         fMultVZERO->Fit(fpol0, "", "", 16, 24);
1210         (cuts->GetUseVZERORing(2)) ? cuts->SetVZEROCpol(2, fpol0->GetParameter(0)) : cuts->SetVZEROCpol(2, 0.);
1211         fMultVZERO->Fit(fpol0, "", "", 24, 32);
1212         (cuts->GetUseVZERORing(3)) ? cuts->SetVZEROCpol(3, fpol0->GetParameter(0)) : cuts->SetVZEROCpol(3, 0.);
1213         // same thing for vero A
1214         fMultVZERO->Fit(fpol0, "", "", 32, 40);
1215         (cuts->GetUseVZERORing(4)) ? cuts->SetVZEROApol(0, fpol0->GetParameter(0)) : cuts->SetVZEROApol(0, 0.);
1216         fMultVZERO->Fit(fpol0, "", "", 40, 48);
1217         (cuts->GetUseVZERORing(5)) ? cuts->SetVZEROApol(1, fpol0->GetParameter(0)) : cuts->SetVZEROApol(1, 0.);
1218         fMultVZERO->Fit(fpol0, "", "", 48, 56);
1219         (cuts->GetUseVZERORing(6)) ? cuts->SetVZEROApol(2, fpol0->GetParameter(0)) : cuts->SetVZEROApol(2, 0.);
1220         fMultVZERO->Fit(fpol0, "", "", 56, 64);
1221         (cuts->GetUseVZERORing(7)) ? cuts->SetVZEROApol(3, fpol0->GetParameter(0)) : cuts->SetVZEROApol(3, 0.);
1222     } else {
1223         // do the calibration in one go. the calibration will still be
1224         // stored per ring, but each ring has the same weight now
1225        fMultVZERO->Fit(fpol0,"","",0,31);
1226        for(Int_t i(0); i < 4; i++) cuts->SetVZEROCpol(i, fpol0->GetParameter(0));
1227        fMultVZERO->Fit(fpol0,"","",32,64);
1228        for(Int_t i(0); i < 4; i++) cuts->SetVZEROApol(i, fpol0->GetParameter(0));
1229     }
1230     // the parameters to weigh the vzero track cuts have been extracted now,
1231     // so we can pass them to the current track cuts obect
1232     cuts->SetVZEROgainEqualisation(fMultVZERO);       // passed as a TH1
1233
1234     // step 2) reweight the q-vectors that will be  called by flow methods which use
1235     // subevents
1236     // underlying assumption is that subevent a uses VZEROA
1237     // and subevent b uses VZEROC
1238     for(Int_t iside=0;iside<2;iside++){
1239         for(Int_t icoord=0;icoord<2;icoord++){
1240             for(Int_t i=0;i  < 9;i++){
1241                 char namecont[100];
1242                 if(iside==0 && icoord==0)
1243                   snprintf(namecont,100,"hQxc2_%i",i);
1244                 else if(iside==1 && icoord==0)
1245                   snprintf(namecont,100,"hQxa2_%i",i);
1246                 else if(iside==0 && icoord==1)
1247                   snprintf(namecont,100,"hQyc2_%i",i);
1248                 else if(iside==1 && icoord==1)
1249                   snprintf(namecont,100,"hQya2_%i",i);
1250
1251                 cont = (AliOADBContainer*) foadb->Get(namecont);
1252                 if(!cont){
1253                     printf("OADB object %s is not available in the file\n",namecont);
1254                     return;
1255                 }
1256
1257                 if(!(cont->GetObject(run))){
1258                     printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
1259                     run = 137366;
1260                 }
1261
1262                 // after grabbing all the info, set the CORRECTION TERMS to
1263                 // the 2nd and 3rd order qsub-vectors
1264                 // we do this here for all centralities, so that subsequent events
1265                 // can grab the correction from these cached values
1266                 fMeanQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
1267                 fWidthQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
1268
1269                 //for v3
1270                 if(iside==0 && icoord==0)
1271                   snprintf(namecont,100,"hQxc3_%i",i);
1272                 else if(iside==1 && icoord==0)
1273                   snprintf(namecont,100,"hQxa3_%i",i);
1274                 else if(iside==0 && icoord==1)
1275                   snprintf(namecont,100,"hQyc3_%i",i);
1276                 else if(iside==1 && icoord==1)
1277                   snprintf(namecont,100,"hQya3_%i",i);
1278
1279                 cont = (AliOADBContainer*) foadb->Get(namecont);
1280                 if(!cont){
1281                     printf("OADB object %s is not available in the file\n",namecont);
1282                     return;
1283                 }
1284
1285                 if(!(cont->GetObject(run))){
1286                     printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
1287                     run = 137366;
1288                 }
1289                 fMeanQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
1290                 fWidthQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
1291
1292             }
1293         }
1294     }
1295     // set the recentering style (might be switched back to -1 if recentering is disabeled)
1296     fApplyRecentering = 2010;
1297 }
1298 //_____________________________________________________________________________
1299
1300 void AliFlowEvent::ClearFast()
1301 {
1302   //clear the event without releasing any memory
1303   //note that cached run number of recentering settigns are not clear
1304   //(see AliFlowEvent::ClearCachedRun() )
1305   AliFlowEventSimple::ClearFast();
1306 }
1307
1308 //_____________________________________________________________________________
1309
1310 void AliFlowEvent::ClearCachedRun()
1311 {
1312     //clear the cached run (not in clear fast as cache needs to be persistent in most cases )
1313   fCachedRun=0;
1314   fApplyRecentering=0;
1315 }