PWG/FLOW/Base/AliFlowEventSimple.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   AliFlowEventSimple: A simple event
  18   for flow analysis
  19
  20   origin: Naomi van der Kolk (kolk@nikhef.nl)
  21           Ante Bilandzic     (anteb@nikhef.nl)
  22           Raimond Snellings  (Raimond.Snellings@nikhef.nl)
  23   mods:   Mikolaj Krzewicki  (mikolaj.krzewicki@cern.ch)
  24 *****************************************************************/
  25
  26 #include "Riostream.h"
  27 #include "TObjArray.h"
  28 #include "TFile.h"
  29 #include "TList.h"
  30 #include "TTree.h"
  31 #include "TParticle.h"
  32 #include "TMath.h"
  33 #include "TH1F.h"
  34 #include "TH1D.h"
  35 #include "TF1.h"
  36 #include "TProfile.h"
  37 #include "TParameter.h"
  38 #include "TBrowser.h"
  39 #include "AliFlowVector.h"
  40 #include "AliFlowTrackSimple.h"
  41 #include "AliFlowTrackSimpleCuts.h"
  42 #include "AliFlowEventSimple.h"
  43 #include "TRandom.h"
  44
  45 using std::cout;
  46 using std::endl;
  47 using std::cerr;
  48 ClassImp(AliFlowEventSimple)
  49
  50 //-----------------------------------------------------------------------
  51 AliFlowEventSimple::AliFlowEventSimple():
  52   fTrackCollection(NULL),
  53   fReferenceMultiplicity(0),
  54   fNumberOfTracks(0),
  55   fNumberOfRPs(0),
  56   fMCReactionPlaneAngle(0.),
  57   fMCReactionPlaneAngleIsSet(kFALSE),
  58   fAfterBurnerPrecision(0.001),
  59   fUserModified(kFALSE),
  60   fNumberOfTracksWrap(NULL),
  61   fNumberOfRPsWrap(NULL),
  62   fMCReactionPlaneAngleWrap(NULL)
  63 {
  64   cout << "AliFlowEventSimple: Default constructor to be used only by root for io" << endl;
  65 }
  66
  67 //-----------------------------------------------------------------------
  68 AliFlowEventSimple::AliFlowEventSimple( Int_t n,
  69                                         ConstructionMethod method,
  70                                         TF1* ptDist,
  71                                         Double_t phiMin,
  72                                         Double_t phiMax,
  73                                         Double_t etaMin,
  74                                         Double_t etaMax):
  75   fTrackCollection(new TObjArray(n)),
  76   fReferenceMultiplicity(0),
  77   fNumberOfTracks(0),
  78   fNumberOfRPs(0),
  79   fMCReactionPlaneAngle(0.),
  80   fMCReactionPlaneAngleIsSet(kFALSE),
  81   fAfterBurnerPrecision(0.001),
  82   fUserModified(kFALSE),
  83   fNumberOfTracksWrap(NULL),
  84   fNumberOfRPsWrap(NULL),
  85   fMCReactionPlaneAngleWrap(NULL)
  86 {
  87   //ctor
  88   // if second argument is set to AliFlowEventSimple::kGenerate
  89   // it generates n random tracks with given Pt distribution
  90   // (a sane default is provided), phi and eta are uniform
  91
  92   if (method==kGenerate)
  93     Generate(n,ptDist,phiMin,phiMax,etaMin,etaMax);
  94 }
  95
  96 //-----------------------------------------------------------------------
  97 AliFlowEventSimple::AliFlowEventSimple(const AliFlowEventSimple& anEvent):
  98   TObject(anEvent),
  99   fTrackCollection((TObjArray*)(anEvent.fTrackCollection)->Clone()),
 100   fReferenceMultiplicity(anEvent.fReferenceMultiplicity),
 101   fNumberOfTracks(anEvent.fNumberOfTracks),
 102   fNumberOfRPs(anEvent.fNumberOfRPs),
 103   fMCReactionPlaneAngle(anEvent.fMCReactionPlaneAngle),
 104   fMCReactionPlaneAngleIsSet(anEvent.fMCReactionPlaneAngleIsSet),
 105   fAfterBurnerPrecision(anEvent.fAfterBurnerPrecision),
 106   fUserModified(anEvent.fUserModified),
 107   fNumberOfTracksWrap(anEvent.fNumberOfTracksWrap),
 108   fNumberOfRPsWrap(anEvent.fNumberOfRPsWrap),
 109   fMCReactionPlaneAngleWrap(anEvent.fMCReactionPlaneAngleWrap)
 110 {
 111   //copy constructor
 112 }
 113
 114 //-----------------------------------------------------------------------
 115 AliFlowEventSimple& AliFlowEventSimple::operator=(const AliFlowEventSimple& anEvent)
 116 {
 117   //assignment operator
 118   if (&anEvent==this) return *this; //check self-assignment
 119   if (fTrackCollection) fTrackCollection->Delete();
 120   delete fTrackCollection;
 121   fTrackCollection = (TObjArray*)(anEvent.fTrackCollection)->Clone(); //deep copy
 122   fReferenceMultiplicity = anEvent.fReferenceMultiplicity;
 123   fNumberOfTracks = anEvent.fNumberOfTracks;
 124   fNumberOfRPs = anEvent.fNumberOfRPs;
 125   fMCReactionPlaneAngle = anEvent.fMCReactionPlaneAngle;
 126   fMCReactionPlaneAngleIsSet = anEvent.fMCReactionPlaneAngleIsSet;
 127   fAfterBurnerPrecision = anEvent.fAfterBurnerPrecision;
 128   fUserModified=anEvent.fUserModified;
 129   fNumberOfTracksWrap = anEvent.fNumberOfTracksWrap;
 130   fNumberOfRPsWrap = anEvent.fNumberOfRPsWrap;
 131   fMCReactionPlaneAngleWrap=anEvent.fMCReactionPlaneAngleWrap;
 132   return *this;
 133 }
 134
 135 //-----------------------------------------------------------------------
 136 AliFlowEventSimple::~AliFlowEventSimple()
 137 {
 138   //destructor
 139   if (fTrackCollection) fTrackCollection->Delete();
 140   delete fTrackCollection;
 141   delete fNumberOfTracksWrap;
 142   delete fNumberOfRPsWrap;
 143   delete fMCReactionPlaneAngleWrap;
 144 }
 145
 146 //-----------------------------------------------------------------------
 147 void AliFlowEventSimple::Generate(Int_t nParticles,
 148                                   TF1* ptDist,
 149                                   Double_t phiMin,
 150                                   Double_t phiMax,
 151                                   Double_t etaMin,
 152                                   Double_t etaMax)
 153 {
 154   //generate nParticles random tracks uniform in phi and eta
 155   //according to the specified pt distribution
 156   if (!ptDist)
 157   {
 158     static TF1 ptdistribution("ptSpectra","x*TMath::Exp(-pow(0.13957*0.13957+x*x,0.5)/0.4)",0.1,10.);
 159     ptDist=&ptdistribution;
 160   }
 161
 162   for (Int_t i=0; i<nParticles; i++)
 163   {
 164     AliFlowTrackSimple* track = new AliFlowTrackSimple();
 165     track->SetPhi( gRandom->Uniform(phiMin,phiMax) );
 166     track->SetEta( gRandom->Uniform(etaMin,etaMax) );
 167     track->SetPt( ptDist->GetRandom() );
 168     track->SetCharge( (gRandom->Uniform()-0.5<0)?-1:1 );
 169     AddTrack(track);
 170   }
 171   fMCReactionPlaneAngle=gRandom->Uniform(0.0,TMath::TwoPi());
 172   fMCReactionPlaneAngleIsSet=kTRUE;
 173   SetUserModified();
 174 }
 175
 176 //-----------------------------------------------------------------------
 177 AliFlowTrackSimple* AliFlowEventSimple::GetTrack(Int_t i)
 178 {
 179   //get track i from collection
 180   if (i>=fNumberOfTracks) return NULL;
 181   AliFlowTrackSimple* pTrack = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i)) ;
 182   return pTrack;
 183 }
 184
 185 //-----------------------------------------------------------------------
 186 void AliFlowEventSimple::AddTrack( AliFlowTrackSimple* track )
 187 {
 188   //add a track, delete the old one if necessary
 189   if (fNumberOfTracks < fTrackCollection->GetEntriesFast())
 190   {
 191     TObject* o = fTrackCollection->At(fNumberOfTracks);
 192     if (o) delete o;
 193   }
 194   fTrackCollection->AddAtAndExpand(track,fNumberOfTracks);
 195   fNumberOfTracks++;
 196 }
 197
 198 //-----------------------------------------------------------------------
 199 AliFlowVector AliFlowEventSimple::GetQ(Int_t n, TList *weightsList, Bool_t usePhiWeights, Bool_t usePtWeights, Bool_t useEtaWeights)
 200 {
 201   // calculate Q-vector in harmonic n without weights (default harmonic n=2)
 202   Double_t dQX = 0.;
 203   Double_t dQY = 0.;
 204   AliFlowVector vQ;
 205   vQ.Set(0.,0.);
 206
 207   Int_t iOrder = n;
 208   Double_t sumOfWeights = 0.;
 209   Double_t dPhi = 0.;
 210   Double_t dPt = 0.;
 211   Double_t dEta = 0.;
 212   Double_t dWeight = 1.;
 213
 214   AliFlowTrackSimple* pTrack = NULL;
 215
 216   Int_t nBinsPhi = 0;
 217   Double_t dBinWidthPt = 0.;
 218   Double_t dPtMin = 0.;
 219   Double_t dBinWidthEta = 0.;
 220   Double_t dEtaMin = 0.;
 221
 222   Double_t wPhi = 1.; // weight Phi
 223   Double_t wPt = 1.;  // weight Pt
 224   Double_t wEta = 1.; // weight Eta
 225
 226   TH1F *phiWeights = NULL;
 227   TH1D *ptWeights  = NULL;
 228   TH1D *etaWeights = NULL;
 229
 230   if(weightsList)
 231   {
 232     if(usePhiWeights)
 233     {
 234       phiWeights = dynamic_cast<TH1F *>(weightsList->FindObject("phi_weights"));
 235       if(phiWeights) nBinsPhi = phiWeights->GetNbinsX();
 236     }
 237     if(usePtWeights)
 238     {
 239       ptWeights = dynamic_cast<TH1D *>(weightsList->FindObject("pt_weights"));
 240       if(ptWeights)
 241       {
 242         dBinWidthPt = ptWeights->GetBinWidth(1); // assuming that all bins have the same width
 243         dPtMin = (ptWeights->GetXaxis())->GetXmin();
 244       }
 245     }
 246     if(useEtaWeights)
 247     {
 248       etaWeights = dynamic_cast<TH1D *>(weightsList->FindObject("eta_weights"));
 249       if(etaWeights)
 250       {
 251         dBinWidthEta = etaWeights->GetBinWidth(1); // assuming that all bins have the same width
 252         dEtaMin = (etaWeights->GetXaxis())->GetXmin();
 253       }
 254     }
 255   } // end of if(weightsList)
 256
 257   // loop over tracks
 258   for(Int_t i=0; i<fNumberOfTracks; i++)
 259   {
 260     pTrack = (AliFlowTrackSimple*)fTrackCollection->At(i);
 261     if(pTrack)
 262     {
 263       if(pTrack->InRPSelection())
 264       {
 265         dPhi = pTrack->Phi();
 266         dPt  = pTrack->Pt();
 267         dEta = pTrack->Eta();
 268               dWeight = pTrack->Weight();
 269
 270         // determine Phi weight: (to be improved, I should here only access it + the treatment of gaps in the if statement)
 271         if(phiWeights && nBinsPhi)
 272         {
 273           wPhi = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*nBinsPhi/TMath::TwoPi())));
 274         }
 275         // determine v'(pt) weight:
 276         if(ptWeights && dBinWidthPt)
 277         {
 278           wPt=ptWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-dPtMin)/dBinWidthPt)));
 279         }
 280         // determine v'(eta) weight:
 281         if(etaWeights && dBinWidthEta)
 282         {
 283           wEta=etaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-dEtaMin)/dBinWidthEta)));
 284         }
 285
 286         // building up the weighted Q-vector:
 287         dQX += dWeight*wPhi*wPt*wEta*TMath::Cos(iOrder*dPhi);
 288         dQY += dWeight*wPhi*wPt*wEta*TMath::Sin(iOrder*dPhi);
 289
 290         // weighted multiplicity:
 291         sumOfWeights += dWeight*wPhi*wPt*wEta;
 292
 293       } // end of if (pTrack->InRPSelection())
 294     } // end of if (pTrack)
 295     else
 296     {
 297       cerr << "no particle!!!"<<endl;
 298     }
 299   } // loop over particles
 300
 301   vQ.Set(dQX,dQY);
 302   vQ.SetMult(sumOfWeights);
 303
 304   return vQ;
 305
 306 }
 307
 308 //-----------------------------------------------------------------------
 309 void AliFlowEventSimple::Get2Qsub(AliFlowVector* Qarray, Int_t n, TList *weightsList, Bool_t usePhiWeights, Bool_t usePtWeights, Bool_t useEtaWeights)
 310 {
 311
 312   // calculate Q-vector in harmonic n without weights (default harmonic n=2)
 313   Double_t dQX = 0.;
 314   Double_t dQY = 0.;
 315
 316   Int_t iOrder = n;
 317   Double_t sumOfWeights = 0.;
 318   Double_t dPhi = 0.;
 319   Double_t dPt  = 0.;
 320   Double_t dEta = 0.;
 321   Double_t dWeight = 1.;
 322
 323   AliFlowTrackSimple* pTrack = NULL;
 324
 325   Int_t    iNbinsPhiSub0 = 0;
 326   Int_t    iNbinsPhiSub1 = 0;
 327   Double_t dBinWidthPt = 0.;
 328   Double_t dPtMin      = 0.;
 329   Double_t dBinWidthEta= 0.;
 330   Double_t dEtaMin     = 0.;
 331
 332   Double_t dWphi = 1.;  // weight Phi
 333   Double_t dWpt  = 1.;  // weight Pt
 334   Double_t dWeta = 1.;  // weight Eta
 335
 336   TH1F* phiWeightsSub0 = NULL;
 337   TH1F* phiWeightsSub1 = NULL;
 338   TH1D* ptWeights  = NULL;
 339   TH1D* etaWeights = NULL;
 340
 341   if(weightsList)
 342   {
 343     if(usePhiWeights)
 344     {
 345       phiWeightsSub0 = dynamic_cast<TH1F *>(weightsList->FindObject("phi_weights_sub0"));
 346       if(phiWeightsSub0) {
 347         iNbinsPhiSub0 = phiWeightsSub0->GetNbinsX();
 348       }
 349       phiWeightsSub1 = dynamic_cast<TH1F *>(weightsList->FindObject("phi_weights_sub1"));
 350       if(phiWeightsSub1) {
 351         iNbinsPhiSub1 = phiWeightsSub1->GetNbinsX();
 352       }
 353     }
 354     if(usePtWeights)
 355     {
 356       ptWeights = dynamic_cast<TH1D *>(weightsList->FindObject("pt_weights"));
 357       if(ptWeights)
 358       {
 359         dBinWidthPt = ptWeights->GetBinWidth(1); // assuming that all bins have the same width
 360         dPtMin = (ptWeights->GetXaxis())->GetXmin();
 361       }
 362     }
 363     if(useEtaWeights)
 364     {
 365       etaWeights = dynamic_cast<TH1D *>(weightsList->FindObject("eta_weights"));
 366       if(etaWeights)
 367       {
 368         dBinWidthEta = etaWeights->GetBinWidth(1); // assuming that all bins have the same width
 369         dEtaMin = (etaWeights->GetXaxis())->GetXmin();
 370       }
 371     }
 372   } // end of if(weightsList)
 373
 374   //loop over the two subevents
 375   for (Int_t s=0; s<2; s++)
 376   {
 377     // loop over tracks
 378     for(Int_t i=0; i<fNumberOfTracks; i++)
 379     {
 380       pTrack = (AliFlowTrackSimple*)fTrackCollection->At(i);
 381       if(pTrack)
 382       {
 383         if(pTrack->InRPSelection())
 384         {
 385           if (pTrack->InSubevent(s))
 386           {
 387             dPhi    = pTrack->Phi();
 388             dPt     = pTrack->Pt();
 389             dEta    = pTrack->Eta();
 390             dWeight = pTrack->Weight();
 391
 392             // determine Phi weight: (to be improved, I should here only access it + the treatment of gaps in the if statement)
 393             //subevent 0
 394             if(s == 0)  {
 395               if(phiWeightsSub0 && iNbinsPhiSub0)  {
 396                 Int_t phiBin = 1+(Int_t)(TMath::Floor(dPhi*iNbinsPhiSub0/TMath::TwoPi()));
 397                 //use the phi value at the center of the bin
 398                 dPhi  = phiWeightsSub0->GetBinCenter(phiBin);
 399                 dWphi = phiWeightsSub0->GetBinContent(phiBin);
 400               }
 401             }
 402             //subevent 1
 403             else if (s == 1) {
 404               if(phiWeightsSub1 && iNbinsPhiSub1) {
 405                 Int_t phiBin = 1+(Int_t)(TMath::Floor(dPhi*iNbinsPhiSub1/TMath::TwoPi()));
 406                 //use the phi value at the center of the bin
 407                 dPhi  = phiWeightsSub1->GetBinCenter(phiBin);
 408                 dWphi = phiWeightsSub1->GetBinContent(phiBin);
 409               }
 410             }
 411
 412             // determine v'(pt) weight:
 413             if(ptWeights && dBinWidthPt)
 414             {
 415               dWpt=ptWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-dPtMin)/dBinWidthPt)));
 416             }
 417
 418             // determine v'(eta) weight:
 419             if(etaWeights && dBinWidthEta)
 420             {
 421               dWeta=etaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-dEtaMin)/dBinWidthEta)));
 422             }
 423
 424             // building up the weighted Q-vector:
 425             dQX += dWeight*dWphi*dWpt*dWeta*TMath::Cos(iOrder*dPhi);
 426             dQY += dWeight*dWphi*dWpt*dWeta*TMath::Sin(iOrder*dPhi);
 427
 428             // weighted multiplicity:
 429             sumOfWeights+=dWeight*dWphi*dWpt*dWeta;
 430
 431           } // end of subevent
 432         } // end of if (pTrack->InRPSelection())
 433       } // end of if (pTrack)
 434       else
 435       {
 436         cerr << "no particle!!!"<<endl;
 437       }
 438     } // loop over particles
 439     Qarray[s].Set(dQX,dQY);
 440     Qarray[s].SetMult(sumOfWeights);
 441     //reset
 442     sumOfWeights = 0.;
 443     dQX = 0.;
 444     dQY = 0.;
 445   }
 446
 447 }
 448
 449
 450 //-----------------------------------------------------------------------
 451 void AliFlowEventSimple::Print(Option_t *option) const
 452 {
 453   //   -*-*-*-*-*Print some global quantities for this histogram collection class *-*-*-*-*-*-*-*
 454   //             ===============================================
 455   //   printf( "TH1.Print Name  = %s, Entries= %d, Total sum= %g\n",GetName(),Int_t(fEntries),GetSumOfWeights());
 456   printf( "Class.Print Name = %s, #tracks= %d, Number of RPs= %d, MC EventPlaneAngle= %f\n",
 457           GetName(),fNumberOfTracks, fNumberOfRPs, fMCReactionPlaneAngle );
 458
 459   TString optionstr(option);
 460   if (!optionstr.Contains("all")) return;
 461   if (fTrackCollection)
 462   {
 463     fTrackCollection->Print(option);
 464   }
 465   else
 466   {
 467     printf( "Empty track collection \n");
 468   }
 469 }
 470
 471 //-----------------------------------------------------------------------
 472 void AliFlowEventSimple::Browse(TBrowser *b)
 473 {
 474   //browse in TBrowser
 475   if (!b) return;
 476   if (!fNumberOfTracksWrap)
 477   {
 478     fNumberOfTracksWrap = new TParameter<int>("fNumberOfTracks", fNumberOfTracks);
 479     b->Add(fNumberOfTracksWrap);
 480   }
 481   if (!fNumberOfRPsWrap)
 482   {
 483     fNumberOfRPsWrap = new TParameter<int>("fNumberOfRPs", fNumberOfRPs);
 484     b->Add(fNumberOfRPsWrap);
 485   }
 486   if (!fMCReactionPlaneAngleWrap)
 487   {
 488     fMCReactionPlaneAngleWrap = new TParameter<double>(" fMCReactionPlaneAngle",  fMCReactionPlaneAngle);
 489     b->Add( fMCReactionPlaneAngleWrap);
 490   }
 491   if (fTrackCollection) b->Add(fTrackCollection,"AliFlowTracksSimple");
 492 }
 493
 494 //-----------------------------------------------------------------------
 495 AliFlowEventSimple::AliFlowEventSimple( TTree* inputTree,
 496                                         const AliFlowTrackSimpleCuts* rpCuts,
 497                                         const AliFlowTrackSimpleCuts* poiCuts):
 498   fTrackCollection(NULL),
 499   fReferenceMultiplicity(0),
 500   fNumberOfTracks(0),
 501   fNumberOfRPs(0),
 502   fMCReactionPlaneAngle(0.),
 503   fMCReactionPlaneAngleIsSet(kFALSE),
 504   fAfterBurnerPrecision(0.001),
 505   fUserModified(kFALSE),
 506   fNumberOfTracksWrap(NULL),
 507   fNumberOfRPsWrap(NULL),
 508   fMCReactionPlaneAngleWrap(NULL)
 509 {
 510   //constructor, fills the event from a TTree of kinematic.root files
 511   //applies RP and POI cuts, tags the tracks
 512
 513   Int_t numberOfInputTracks = inputTree->GetEntries() ;
 514   fTrackCollection = new TObjArray(numberOfInputTracks/2);
 515
 516   TParticle* pParticle = new TParticle();
 517   inputTree->SetBranchAddress("Particles",&pParticle);
 518
 519   Int_t iSelParticlesPOI = 0;
 520
 521   for (Int_t i=0; i<numberOfInputTracks; i++)
 522   {
 523     inputTree->GetEntry(i);   //get input particle
 524
 525     if (!pParticle) continue; //no particle
 526     if (!pParticle->IsPrimary()) continue;
 527
 528     Bool_t rpOK = rpCuts->PassesCuts(pParticle);
 529     Bool_t poiOK = poiCuts->PassesCuts(pParticle);
 530
 531     if (rpOK || poiOK)
 532     {
 533       AliFlowTrackSimple* pTrack = new AliFlowTrackSimple(pParticle);
 534
 535       //marking the particles used for int. flow:
 536       if(rpOK)
 537       {
 538         pTrack->SetForRPSelection(kTRUE);
 539         fNumberOfRPs++;
 540       }
 541       //marking the particles used for diff. flow:
 542       if(poiOK)
 543       {
 544         pTrack->SetForPOISelection(kTRUE);
 545         iSelParticlesPOI++;
 546       }
 547       //adding a particles which were used either for int. or diff. flow to the list
 548       AddTrack(pTrack);
 549     }
 550   }//for i
 551   delete pParticle;
 552 }
 553
 554 //_____________________________________________________________________________
 555 void AliFlowEventSimple::CloneTracks(Int_t n)
 556 {
 557   //clone every track n times to add non-flow
 558   if (n<=0) return; //no use to clone stuff zero or less times
 559   Int_t ntracks = fNumberOfTracks;
 560   fTrackCollection->Expand((n+1)*fNumberOfTracks);
 561   for (Int_t i=0; i<n; i++)
 562   {
 563     for (Int_t itrack=0; itrack<ntracks; itrack++)
 564     {
 565       AliFlowTrackSimple* track = dynamic_cast<AliFlowTrackSimple*>(fTrackCollection->At(itrack));
 566       if (!track) continue;
 567       AddTrack(static_cast<AliFlowTrackSimple*>(track->Clone()));
 568     }
 569   }
 570   SetUserModified();
 571 }
 572
 573 //_____________________________________________________________________________
 574 void AliFlowEventSimple::ResolutionPt(Double_t res)
 575 {
 576   //smear pt of all tracks by gaussian with sigma=res
 577   for (Int_t i=0; i<fNumberOfTracks; i++)
 578   {
 579     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 580     if (track) track->ResolutionPt(res);
 581   }
 582   SetUserModified();
 583 }
 584
 585 //_____________________________________________________________________________
 586 void AliFlowEventSimple::TagSubeventsInEta( Double_t etaMinA,
 587                                             Double_t etaMaxA,
 588                                             Double_t etaMinB,
 589                                             Double_t etaMaxB )
 590 {
 591   //Flag two subevents in given eta ranges
 592   for (Int_t i=0; i<fNumberOfTracks; i++)
 593   {
 594     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 595     if (!track) continue;
 596     track->ResetSubEventTags();
 597     Double_t eta=track->Eta();
 598     if (eta >= etaMinA && eta <= etaMaxA) track->SetForSubevent(0);
 599     if (eta >= etaMinB && eta <= etaMaxB) track->SetForSubevent(1);
 600   }
 601 }
 602
 603 //_____________________________________________________________________________
 604 void AliFlowEventSimple::TagSubeventsByCharge()
 605 {
 606   //Flag two subevents in given eta ranges
 607   for (Int_t i=0; i<fNumberOfTracks; i++)
 608   {
 609     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 610     if (!track) continue;
 611     track->ResetSubEventTags();
 612     Int_t charge=track->Charge();
 613     if (charge<0) track->SetForSubevent(0);
 614     if (charge>0) track->SetForSubevent(1);
 615   }
 616 }
 617
 618 //_____________________________________________________________________________
 619 void AliFlowEventSimple::AddV1( Double_t v1 )
 620 {
 621   //add v2 to all tracks wrt the reaction plane angle
 622   for (Int_t i=0; i<fNumberOfTracks; i++)
 623   {
 624     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 625     if (track) track->AddV1(v1, fMCReactionPlaneAngle, fAfterBurnerPrecision);
 626   }
 627   SetUserModified();
 628 }
 629
 630 //_____________________________________________________________________________
 631 void AliFlowEventSimple::AddV2( Double_t v2 )
 632 {
 633   //add v2 to all tracks wrt the reaction plane angle
 634   for (Int_t i=0; i<fNumberOfTracks; i++)
 635   {
 636     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 637     if (track) track->AddV2(v2, fMCReactionPlaneAngle, fAfterBurnerPrecision);
 638   }
 639   SetUserModified();
 640 }
 641
 642 //_____________________________________________________________________________
 643 void AliFlowEventSimple::AddV3( Double_t v3 )
 644 {
 645   //add v3 to all tracks wrt the reaction plane angle
 646   for (Int_t i=0; i<fNumberOfTracks; i++)
 647   {
 648     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 649     if (track) track->AddV3(v3, fMCReactionPlaneAngle, fAfterBurnerPrecision);
 650   }
 651   SetUserModified();
 652 }
 653
 654 //_____________________________________________________________________________
 655 void AliFlowEventSimple::AddV4( Double_t v4 )
 656 {
 657   //add v4 to all tracks wrt the reaction plane angle
 658   for (Int_t i=0; i<fNumberOfTracks; i++)
 659   {
 660     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 661     if (track) track->AddV4(v4, fMCReactionPlaneAngle, fAfterBurnerPrecision);
 662   }
 663   SetUserModified();
 664 }
 665
 666 //_____________________________________________________________________________
 667 void AliFlowEventSimple::AddV5( Double_t v5 )
 668 {
 669   //add v4 to all tracks wrt the reaction plane angle
 670   for (Int_t i=0; i<fNumberOfTracks; i++)
 671   {
 672     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 673     if (track) track->AddV5(v5, fMCReactionPlaneAngle, fAfterBurnerPrecision);
 674   }
 675   SetUserModified();
 676 }
 677
 678 //_____________________________________________________________________________
 679 void AliFlowEventSimple::AddFlow( Double_t v1, Double_t v2, Double_t v3, Double_t v4, Double_t v5,
 680                                   Double_t rp1, Double_t rp2, Double_t rp3, Double_t rp4, Double_t rp5 )
 681 {
 682   //add flow to all tracks wrt the reaction plane angle, for all harmonic separate angle
 683   for (Int_t i=0; i<fNumberOfTracks; i++)
 684   {
 685     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 686     if (track) track->AddFlow(v1,v2,v3,v4,v5,rp1,rp2,rp3,rp4,rp5,fAfterBurnerPrecision);
 687   }
 688   SetUserModified();
 689 }
 690
 691 //_____________________________________________________________________________
 692 void AliFlowEventSimple::AddFlow( Double_t v1, Double_t v2, Double_t v3, Double_t v4, Double_t v5 )
 693 {
 694   //add flow to all tracks wrt the reaction plane angle
 695   for (Int_t i=0; i<fNumberOfTracks; i++)
 696   {
 697     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 698     if (track) track->AddFlow(v1,v2,v3,v4,v5,fMCReactionPlaneAngle, fAfterBurnerPrecision);
 699   }
 700   SetUserModified();
 701 }
 702
 703 //_____________________________________________________________________________
 704 void AliFlowEventSimple::AddV2( TF1* ptDepV2 )
 705 {
 706   //add v2 to all tracks wrt the reaction plane angle
 707   for (Int_t i=0; i<fNumberOfTracks; i++)
 708   {
 709     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 710     if (!track) continue;
 711     Double_t v2 = ptDepV2->Eval(track->Pt());
 712     track->AddV2(v2, fMCReactionPlaneAngle, fAfterBurnerPrecision);
 713   }
 714   SetUserModified();
 715 }
 716
 717 //_____________________________________________________________________________
 718 void AliFlowEventSimple::TagRP( const AliFlowTrackSimpleCuts* cuts )
 719 {
 720   //tag tracks as reference particles (RPs)
 721   for (Int_t i=0; i<fNumberOfTracks; i++)
 722   {
 723     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 724     if (!track) continue;
 725     Bool_t pass=cuts->PassesCuts(track);
 726     Bool_t rpTrack=track->InRPSelection();
 727     if (pass)
 728     {
 729       if (!rpTrack) fNumberOfRPs++; //only increase if not already tagged
 730     }
 731     else
 732     {
 733       if (rpTrack) fNumberOfRPs--; //only decrease if detagging
 734     }
 735     track->SetForRPSelection(pass);
 736   }
 737 }
 738
 739 //_____________________________________________________________________________
 740 void AliFlowEventSimple::TagPOI( const AliFlowTrackSimpleCuts* cuts )
 741 {
 742   //tag tracks as particles of interest (POIs)
 743   for (Int_t i=0; i<fNumberOfTracks; i++)
 744   {
 745     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 746     if (!track) continue;
 747     Bool_t pass=cuts->PassesCuts(track);
 748     track->SetForPOISelection(pass);
 749   }
 750 }
 751
 752 //_____________________________________________________________________________
 753 void AliFlowEventSimple::DefineDeadZone( Double_t etaMin,
 754                                          Double_t etaMax,
 755                                          Double_t phiMin,
 756                                          Double_t phiMax )
 757 {
 758   //mark tracks in given eta-phi region as dead
 759   //by resetting the flow bits
 760   for (Int_t i=0; i<fNumberOfTracks; i++)
 761   {
 762     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 763     Double_t eta = track->Eta();
 764     Double_t phi = track->Phi();
 765     if (eta>etaMin && eta<etaMax && phi>phiMin && phi<phiMax)
 766     {
 767       if (track->InRPSelection()) fNumberOfRPs--;
 768       track->ResetFlowTags();
 769     }
 770   }
 771 }
 772
 773 //_____________________________________________________________________________
 774 Int_t AliFlowEventSimple::CleanUpDeadTracks()
 775 {
 776   //remove tracks that have no flow tags set and cleanup the container
 777   //returns number of cleaned tracks
 778   Int_t ncleaned=0;
 779   for (Int_t i=0; i<fNumberOfTracks; i++)
 780   {
 781     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 782     if (!track) continue;
 783     if (track->IsDead()) {delete track;track=NULL;ncleaned++;}
 784   }
 785   fTrackCollection->Compress(); //clean up empty slots
 786   fNumberOfTracks-=ncleaned; //update number of tracks
 787   return ncleaned;
 788 }
 789
 790 //_____________________________________________________________________________
 791 TF1* AliFlowEventSimple::SimplePtDepV2()
 792 {
 793   //return a standard pt dependent v2 formula, user has to clean up!
 794   return new TF1("StandardPtDepV2","((x<1.0)*(0.05/1.0)*x+(x>=1.0)*0.05)");
 795 }
 796
 797 //_____________________________________________________________________________
 798 TF1* AliFlowEventSimple::SimplePtSpectrum()
 799 {
 800   //return a standard pt spectrum, user has to clean up!
 801   return new TF1("StandardPtSpectrum","x*TMath::Exp(-pow(0.13957*0.13957+x*x,0.5)/0.4)",0.1,10.);
 802 }
 803
 804 //_____________________________________________________________________________
 805 void AliFlowEventSimple::ClearFast()
 806 {
 807   //clear the counter without deleting allocated objects so they can be reused
 808   fReferenceMultiplicity = 0;
 809   fNumberOfTracks = 0;
 810   fNumberOfRPs = 0;
 811   fMCReactionPlaneAngle = 0.0;
 812   fMCReactionPlaneAngleIsSet = kFALSE;
 813   fAfterBurnerPrecision = 0.001;
 814   fUserModified = kFALSE;
 815 }