PWG/FLOW/Tasks/AliFlowEvent.cxx

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