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 AliFlowTrackSimple* AliFlowEventSimple::MakeNewTrack()
 200 {
 201    AliFlowTrackSimple *t=dynamic_cast<AliFlowTrackSimple *>(fTrackCollection->RemoveAt(fNumberOfTracks));
 202    if( !t ) {  // If there was no track at the end of the list then create a new track
 203       t=new AliFlowTrackSimple();
 204    }
 205
 206    return t;
 207 }
 208
 209 //-----------------------------------------------------------------------
 210 AliFlowVector AliFlowEventSimple::GetQ(Int_t n, TList *weightsList, Bool_t usePhiWeights, Bool_t usePtWeights, Bool_t useEtaWeights)
 211 {
 212   // calculate Q-vector in harmonic n without weights (default harmonic n=2)
 213   Double_t dQX = 0.;
 214   Double_t dQY = 0.;
 215   AliFlowVector vQ;
 216   vQ.Set(0.,0.);
 217
 218   Int_t iOrder = n;
 219   Double_t sumOfWeights = 0.;
 220   Double_t dPhi = 0.;
 221   Double_t dPt = 0.;
 222   Double_t dEta = 0.;
 223   Double_t dWeight = 1.;
 224
 225   AliFlowTrackSimple* pTrack = NULL;
 226
 227   Int_t nBinsPhi = 0;
 228   Double_t dBinWidthPt = 0.;
 229   Double_t dPtMin = 0.;
 230   Double_t dBinWidthEta = 0.;
 231   Double_t dEtaMin = 0.;
 232
 233   Double_t wPhi = 1.; // weight Phi
 234   Double_t wPt = 1.;  // weight Pt
 235   Double_t wEta = 1.; // weight Eta
 236
 237   TH1F *phiWeights = NULL;
 238   TH1D *ptWeights  = NULL;
 239   TH1D *etaWeights = NULL;
 240
 241   if(weightsList)
 242   {
 243     if(usePhiWeights)
 244     {
 245       phiWeights = dynamic_cast<TH1F *>(weightsList->FindObject("phi_weights"));
 246       if(phiWeights) nBinsPhi = phiWeights->GetNbinsX();
 247     }
 248     if(usePtWeights)
 249     {
 250       ptWeights = dynamic_cast<TH1D *>(weightsList->FindObject("pt_weights"));
 251       if(ptWeights)
 252       {
 253         dBinWidthPt = ptWeights->GetBinWidth(1); // assuming that all bins have the same width
 254         dPtMin = (ptWeights->GetXaxis())->GetXmin();
 255       }
 256     }
 257     if(useEtaWeights)
 258     {
 259       etaWeights = dynamic_cast<TH1D *>(weightsList->FindObject("eta_weights"));
 260       if(etaWeights)
 261       {
 262         dBinWidthEta = etaWeights->GetBinWidth(1); // assuming that all bins have the same width
 263         dEtaMin = (etaWeights->GetXaxis())->GetXmin();
 264       }
 265     }
 266   } // end of if(weightsList)
 267
 268   // loop over tracks
 269   for(Int_t i=0; i<fNumberOfTracks; i++)
 270   {
 271     pTrack = (AliFlowTrackSimple*)fTrackCollection->At(i);
 272     if(pTrack)
 273     {
 274       if(pTrack->InRPSelection())
 275       {
 276         dPhi = pTrack->Phi();
 277         dPt  = pTrack->Pt();
 278         dEta = pTrack->Eta();
 279               dWeight = pTrack->Weight();
 280
 281         // determine Phi weight: (to be improved, I should here only access it + the treatment of gaps in the if statement)
 282         if(phiWeights && nBinsPhi)
 283         {
 284           wPhi = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*nBinsPhi/TMath::TwoPi())));
 285         }
 286         // determine v'(pt) weight:
 287         if(ptWeights && dBinWidthPt)
 288         {
 289           wPt=ptWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-dPtMin)/dBinWidthPt)));
 290         }
 291         // determine v'(eta) weight:
 292         if(etaWeights && dBinWidthEta)
 293         {
 294           wEta=etaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-dEtaMin)/dBinWidthEta)));
 295         }
 296
 297         // building up the weighted Q-vector:
 298         dQX += dWeight*wPhi*wPt*wEta*TMath::Cos(iOrder*dPhi);
 299         dQY += dWeight*wPhi*wPt*wEta*TMath::Sin(iOrder*dPhi);
 300
 301         // weighted multiplicity:
 302         sumOfWeights += dWeight*wPhi*wPt*wEta;
 303
 304       } // end of if (pTrack->InRPSelection())
 305     } // end of if (pTrack)
 306     else
 307     {
 308       cerr << "no particle!!!"<<endl;
 309     }
 310   } // loop over particles
 311
 312   vQ.Set(dQX,dQY);
 313   vQ.SetMult(sumOfWeights);
 314
 315   return vQ;
 316
 317 }
 318
 319 //-----------------------------------------------------------------------
 320 void AliFlowEventSimple::Get2Qsub(AliFlowVector* Qarray, Int_t n, TList *weightsList, Bool_t usePhiWeights, Bool_t usePtWeights, Bool_t useEtaWeights)
 321 {
 322
 323   // calculate Q-vector in harmonic n without weights (default harmonic n=2)
 324   Double_t dQX = 0.;
 325   Double_t dQY = 0.;
 326
 327   Int_t iOrder = n;
 328   Double_t sumOfWeights = 0.;
 329   Double_t dPhi = 0.;
 330   Double_t dPt  = 0.;
 331   Double_t dEta = 0.;
 332   Double_t dWeight = 1.;
 333
 334   AliFlowTrackSimple* pTrack = NULL;
 335
 336   Int_t    iNbinsPhiSub0 = 0;
 337   Int_t    iNbinsPhiSub1 = 0;
 338   Double_t dBinWidthPt = 0.;
 339   Double_t dPtMin      = 0.;
 340   Double_t dBinWidthEta= 0.;
 341   Double_t dEtaMin     = 0.;
 342
 343   Double_t dWphi = 1.;  // weight Phi
 344   Double_t dWpt  = 1.;  // weight Pt
 345   Double_t dWeta = 1.;  // weight Eta
 346
 347   TH1F* phiWeightsSub0 = NULL;
 348   TH1F* phiWeightsSub1 = NULL;
 349   TH1D* ptWeights  = NULL;
 350   TH1D* etaWeights = NULL;
 351
 352   if(weightsList)
 353   {
 354     if(usePhiWeights)
 355     {
 356       phiWeightsSub0 = dynamic_cast<TH1F *>(weightsList->FindObject("phi_weights_sub0"));
 357       if(phiWeightsSub0) {
 358         iNbinsPhiSub0 = phiWeightsSub0->GetNbinsX();
 359       }
 360       phiWeightsSub1 = dynamic_cast<TH1F *>(weightsList->FindObject("phi_weights_sub1"));
 361       if(phiWeightsSub1) {
 362         iNbinsPhiSub1 = phiWeightsSub1->GetNbinsX();
 363       }
 364     }
 365     if(usePtWeights)
 366     {
 367       ptWeights = dynamic_cast<TH1D *>(weightsList->FindObject("pt_weights"));
 368       if(ptWeights)
 369       {
 370         dBinWidthPt = ptWeights->GetBinWidth(1); // assuming that all bins have the same width
 371         dPtMin = (ptWeights->GetXaxis())->GetXmin();
 372       }
 373     }
 374     if(useEtaWeights)
 375     {
 376       etaWeights = dynamic_cast<TH1D *>(weightsList->FindObject("eta_weights"));
 377       if(etaWeights)
 378       {
 379         dBinWidthEta = etaWeights->GetBinWidth(1); // assuming that all bins have the same width
 380         dEtaMin = (etaWeights->GetXaxis())->GetXmin();
 381       }
 382     }
 383   } // end of if(weightsList)
 384
 385   //loop over the two subevents
 386   for (Int_t s=0; s<2; s++)
 387   {
 388     // loop over tracks
 389     for(Int_t i=0; i<fNumberOfTracks; i++)
 390     {
 391       pTrack = (AliFlowTrackSimple*)fTrackCollection->At(i);
 392       if(pTrack)
 393       {
 394         if(pTrack->InRPSelection())
 395         {
 396           if (pTrack->InSubevent(s))
 397           {
 398             dPhi    = pTrack->Phi();
 399             dPt     = pTrack->Pt();
 400             dEta    = pTrack->Eta();
 401             dWeight = pTrack->Weight();
 402
 403             // determine Phi weight: (to be improved, I should here only access it + the treatment of gaps in the if statement)
 404             //subevent 0
 405             if(s == 0)  {
 406               if(phiWeightsSub0 && iNbinsPhiSub0)  {
 407                 Int_t phiBin = 1+(Int_t)(TMath::Floor(dPhi*iNbinsPhiSub0/TMath::TwoPi()));
 408                 //use the phi value at the center of the bin
 409                 dPhi  = phiWeightsSub0->GetBinCenter(phiBin);
 410                 dWphi = phiWeightsSub0->GetBinContent(phiBin);
 411               }
 412             }
 413             //subevent 1
 414             else if (s == 1) {
 415               if(phiWeightsSub1 && iNbinsPhiSub1) {
 416                 Int_t phiBin = 1+(Int_t)(TMath::Floor(dPhi*iNbinsPhiSub1/TMath::TwoPi()));
 417                 //use the phi value at the center of the bin
 418                 dPhi  = phiWeightsSub1->GetBinCenter(phiBin);
 419                 dWphi = phiWeightsSub1->GetBinContent(phiBin);
 420               }
 421             }
 422
 423             // determine v'(pt) weight:
 424             if(ptWeights && dBinWidthPt)
 425             {
 426               dWpt=ptWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-dPtMin)/dBinWidthPt)));
 427             }
 428
 429             // determine v'(eta) weight:
 430             if(etaWeights && dBinWidthEta)
 431             {
 432               dWeta=etaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-dEtaMin)/dBinWidthEta)));
 433             }
 434
 435             // building up the weighted Q-vector:
 436             dQX += dWeight*dWphi*dWpt*dWeta*TMath::Cos(iOrder*dPhi);
 437             dQY += dWeight*dWphi*dWpt*dWeta*TMath::Sin(iOrder*dPhi);
 438
 439             // weighted multiplicity:
 440             sumOfWeights+=dWeight*dWphi*dWpt*dWeta;
 441
 442           } // end of subevent
 443         } // end of if (pTrack->InRPSelection())
 444       } // end of if (pTrack)
 445       else
 446       {
 447         cerr << "no particle!!!"<<endl;
 448       }
 449     } // loop over particles
 450     Qarray[s].Set(dQX,dQY);
 451     Qarray[s].SetMult(sumOfWeights);
 452     //reset
 453     sumOfWeights = 0.;
 454     dQX = 0.;
 455     dQY = 0.;
 456   }
 457
 458 }
 459
 460
 461 //-----------------------------------------------------------------------
 462 void AliFlowEventSimple::Print(Option_t *option) const
 463 {
 464   //   -*-*-*-*-*Print some global quantities for this histogram collection class *-*-*-*-*-*-*-*
 465   //             ===============================================
 466   //   printf( "TH1.Print Name  = %s, Entries= %d, Total sum= %g\n",GetName(),Int_t(fEntries),GetSumOfWeights());
 467   printf( "Class.Print Name = %s, #tracks= %d, Number of RPs= %d, MC EventPlaneAngle= %f\n",
 468           GetName(),fNumberOfTracks, fNumberOfRPs, fMCReactionPlaneAngle );
 469
 470   TString optionstr(option);
 471   if (!optionstr.Contains("all")) return;
 472   if (fTrackCollection)
 473   {
 474     fTrackCollection->Print(option);
 475   }
 476   else
 477   {
 478     printf( "Empty track collection \n");
 479   }
 480 }
 481
 482 //-----------------------------------------------------------------------
 483 void AliFlowEventSimple::Browse(TBrowser *b)
 484 {
 485   //browse in TBrowser
 486   if (!b) return;
 487   if (!fNumberOfTracksWrap)
 488   {
 489     fNumberOfTracksWrap = new TParameter<int>("fNumberOfTracks", fNumberOfTracks);
 490     b->Add(fNumberOfTracksWrap);
 491   }
 492   if (!fNumberOfRPsWrap)
 493   {
 494     fNumberOfRPsWrap = new TParameter<int>("fNumberOfRPs", fNumberOfRPs);
 495     b->Add(fNumberOfRPsWrap);
 496   }
 497   if (!fMCReactionPlaneAngleWrap)
 498   {
 499     fMCReactionPlaneAngleWrap = new TParameter<double>(" fMCReactionPlaneAngle",  fMCReactionPlaneAngle);
 500     b->Add( fMCReactionPlaneAngleWrap);
 501   }
 502   if (fTrackCollection) b->Add(fTrackCollection,"AliFlowTracksSimple");
 503 }
 504
 505 //-----------------------------------------------------------------------
 506 AliFlowEventSimple::AliFlowEventSimple( TTree* inputTree,
 507                                         const AliFlowTrackSimpleCuts* rpCuts,
 508                                         const AliFlowTrackSimpleCuts* poiCuts):
 509   fTrackCollection(NULL),
 510   fReferenceMultiplicity(0),
 511   fNumberOfTracks(0),
 512   fNumberOfRPs(0),
 513   fMCReactionPlaneAngle(0.),
 514   fMCReactionPlaneAngleIsSet(kFALSE),
 515   fAfterBurnerPrecision(0.001),
 516   fUserModified(kFALSE),
 517   fNumberOfTracksWrap(NULL),
 518   fNumberOfRPsWrap(NULL),
 519   fMCReactionPlaneAngleWrap(NULL)
 520 {
 521   //constructor, fills the event from a TTree of kinematic.root files
 522   //applies RP and POI cuts, tags the tracks
 523
 524   Int_t numberOfInputTracks = inputTree->GetEntries() ;
 525   fTrackCollection = new TObjArray(numberOfInputTracks/2);
 526
 527   TParticle* pParticle = new TParticle();
 528   inputTree->SetBranchAddress("Particles",&pParticle);
 529
 530   Int_t iSelParticlesPOI = 0;
 531
 532   for (Int_t i=0; i<numberOfInputTracks; i++)
 533   {
 534     inputTree->GetEntry(i);   //get input particle
 535
 536     if (!pParticle) continue; //no particle
 537     if (!pParticle->IsPrimary()) continue;
 538
 539     Bool_t rpOK = rpCuts->PassesCuts(pParticle);
 540     Bool_t poiOK = poiCuts->PassesCuts(pParticle);
 541
 542     if (rpOK || poiOK)
 543     {
 544       AliFlowTrackSimple* pTrack = new AliFlowTrackSimple(pParticle);
 545
 546       //marking the particles used for int. flow:
 547       if(rpOK)
 548       {
 549         pTrack->SetForRPSelection(kTRUE);
 550         fNumberOfRPs++;
 551       }
 552       //marking the particles used for diff. flow:
 553       if(poiOK)
 554       {
 555         pTrack->SetForPOISelection(kTRUE);
 556         iSelParticlesPOI++;
 557       }
 558       //adding a particles which were used either for int. or diff. flow to the list
 559       AddTrack(pTrack);
 560     }
 561   }//for i
 562   delete pParticle;
 563 }
 564
 565 //_____________________________________________________________________________
 566 void AliFlowEventSimple::CloneTracks(Int_t n)
 567 {
 568   //clone every track n times to add non-flow
 569   if (n<=0) return; //no use to clone stuff zero or less times
 570   Int_t ntracks = fNumberOfTracks;
 571   fTrackCollection->Expand((n+1)*fNumberOfTracks);
 572   for (Int_t i=0; i<n; i++)
 573   {
 574     for (Int_t itrack=0; itrack<ntracks; itrack++)
 575     {
 576       AliFlowTrackSimple* track = dynamic_cast<AliFlowTrackSimple*>(fTrackCollection->At(itrack));
 577       if (!track) continue;
 578       AddTrack(static_cast<AliFlowTrackSimple*>(track->Clone()));
 579     }
 580   }
 581   SetUserModified();
 582 }
 583
 584 //_____________________________________________________________________________
 585 void AliFlowEventSimple::ResolutionPt(Double_t res)
 586 {
 587   //smear pt of all tracks by gaussian with sigma=res
 588   for (Int_t i=0; i<fNumberOfTracks; i++)
 589   {
 590     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 591     if (track) track->ResolutionPt(res);
 592   }
 593   SetUserModified();
 594 }
 595
 596 //_____________________________________________________________________________
 597 void AliFlowEventSimple::TagSubeventsInEta( Double_t etaMinA,
 598                                             Double_t etaMaxA,
 599                                             Double_t etaMinB,
 600                                             Double_t etaMaxB )
 601 {
 602   //Flag two subevents in given eta ranges
 603   for (Int_t i=0; i<fNumberOfTracks; i++)
 604   {
 605     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 606     if (!track) continue;
 607     track->ResetSubEventTags();
 608     Double_t eta=track->Eta();
 609     if (eta >= etaMinA && eta <= etaMaxA) track->SetForSubevent(0);
 610     if (eta >= etaMinB && eta <= etaMaxB) track->SetForSubevent(1);
 611   }
 612 }
 613
 614 //_____________________________________________________________________________
 615 void AliFlowEventSimple::TagSubeventsByCharge()
 616 {
 617   //Flag two subevents in given eta ranges
 618   for (Int_t i=0; i<fNumberOfTracks; i++)
 619   {
 620     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 621     if (!track) continue;
 622     track->ResetSubEventTags();
 623     Int_t charge=track->Charge();
 624     if (charge<0) track->SetForSubevent(0);
 625     if (charge>0) track->SetForSubevent(1);
 626   }
 627 }
 628
 629 //_____________________________________________________________________________
 630 void AliFlowEventSimple::AddV1( Double_t v1 )
 631 {
 632   //add v2 to all tracks wrt the reaction plane angle
 633   for (Int_t i=0; i<fNumberOfTracks; i++)
 634   {
 635     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 636     if (track) track->AddV1(v1, fMCReactionPlaneAngle, fAfterBurnerPrecision);
 637   }
 638   SetUserModified();
 639 }
 640
 641 //_____________________________________________________________________________
 642 void AliFlowEventSimple::AddV2( Double_t v2 )
 643 {
 644   //add v2 to all tracks wrt the reaction plane angle
 645   for (Int_t i=0; i<fNumberOfTracks; i++)
 646   {
 647     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 648     if (track) track->AddV2(v2, fMCReactionPlaneAngle, fAfterBurnerPrecision);
 649   }
 650   SetUserModified();
 651 }
 652
 653 //_____________________________________________________________________________
 654 void AliFlowEventSimple::AddV3( Double_t v3 )
 655 {
 656   //add v3 to all tracks wrt the reaction plane angle
 657   for (Int_t i=0; i<fNumberOfTracks; i++)
 658   {
 659     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 660     if (track) track->AddV3(v3, fMCReactionPlaneAngle, fAfterBurnerPrecision);
 661   }
 662   SetUserModified();
 663 }
 664
 665 //_____________________________________________________________________________
 666 void AliFlowEventSimple::AddV4( Double_t v4 )
 667 {
 668   //add v4 to all tracks wrt the reaction plane angle
 669   for (Int_t i=0; i<fNumberOfTracks; i++)
 670   {
 671     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 672     if (track) track->AddV4(v4, fMCReactionPlaneAngle, fAfterBurnerPrecision);
 673   }
 674   SetUserModified();
 675 }
 676
 677 //_____________________________________________________________________________
 678 void AliFlowEventSimple::AddV5( Double_t v5 )
 679 {
 680   //add v4 to all tracks wrt the reaction plane angle
 681   for (Int_t i=0; i<fNumberOfTracks; i++)
 682   {
 683     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 684     if (track) track->AddV5(v5, fMCReactionPlaneAngle, fAfterBurnerPrecision);
 685   }
 686   SetUserModified();
 687 }
 688
 689 //_____________________________________________________________________________
 690 void AliFlowEventSimple::AddFlow( Double_t v1, Double_t v2, Double_t v3, Double_t v4, Double_t v5,
 691                                   Double_t rp1, Double_t rp2, Double_t rp3, Double_t rp4, Double_t rp5 )
 692 {
 693   //add flow to all tracks wrt the reaction plane angle, for all harmonic separate angle
 694   for (Int_t i=0; i<fNumberOfTracks; i++)
 695   {
 696     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 697     if (track) track->AddFlow(v1,v2,v3,v4,v5,rp1,rp2,rp3,rp4,rp5,fAfterBurnerPrecision);
 698   }
 699   SetUserModified();
 700 }
 701
 702 //_____________________________________________________________________________
 703 void AliFlowEventSimple::AddFlow( Double_t v1, Double_t v2, Double_t v3, Double_t v4, Double_t v5 )
 704 {
 705   //add flow to all tracks wrt the reaction plane angle
 706   for (Int_t i=0; i<fNumberOfTracks; i++)
 707   {
 708     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 709     if (track) track->AddFlow(v1,v2,v3,v4,v5,fMCReactionPlaneAngle, fAfterBurnerPrecision);
 710   }
 711   SetUserModified();
 712 }
 713
 714 //_____________________________________________________________________________
 715 void AliFlowEventSimple::AddV2( TF1* ptDepV2 )
 716 {
 717   //add v2 to all tracks wrt the reaction plane angle
 718   for (Int_t i=0; i<fNumberOfTracks; i++)
 719   {
 720     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 721     if (!track) continue;
 722     Double_t v2 = ptDepV2->Eval(track->Pt());
 723     track->AddV2(v2, fMCReactionPlaneAngle, fAfterBurnerPrecision);
 724   }
 725   SetUserModified();
 726 }
 727
 728 //_____________________________________________________________________________
 729 void AliFlowEventSimple::TagRP( const AliFlowTrackSimpleCuts* cuts )
 730 {
 731   //tag tracks as reference particles (RPs)
 732   for (Int_t i=0; i<fNumberOfTracks; i++)
 733   {
 734     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 735     if (!track) continue;
 736     Bool_t pass=cuts->PassesCuts(track);
 737     Bool_t rpTrack=track->InRPSelection();
 738     if (pass)
 739     {
 740       if (!rpTrack) fNumberOfRPs++; //only increase if not already tagged
 741     }
 742     else
 743     {
 744       if (rpTrack) fNumberOfRPs--; //only decrease if detagging
 745     }
 746     track->SetForRPSelection(pass);
 747   }
 748 }
 749
 750 //_____________________________________________________________________________
 751 void AliFlowEventSimple::TagPOI( const AliFlowTrackSimpleCuts* cuts )
 752 {
 753   //tag tracks as particles of interest (POIs)
 754   for (Int_t i=0; i<fNumberOfTracks; i++)
 755   {
 756     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 757     if (!track) continue;
 758     Bool_t pass=cuts->PassesCuts(track);
 759     track->SetForPOISelection(pass);
 760   }
 761 }
 762
 763 //_____________________________________________________________________________
 764 void AliFlowEventSimple::DefineDeadZone( Double_t etaMin,
 765                                          Double_t etaMax,
 766                                          Double_t phiMin,
 767                                          Double_t phiMax )
 768 {
 769   //mark tracks in given eta-phi region as dead
 770   //by resetting the flow bits
 771   for (Int_t i=0; i<fNumberOfTracks; i++)
 772   {
 773     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 774     Double_t eta = track->Eta();
 775     Double_t phi = track->Phi();
 776     if (eta>etaMin && eta<etaMax && phi>phiMin && phi<phiMax)
 777     {
 778       if (track->InRPSelection()) fNumberOfRPs--;
 779       track->ResetFlowTags();
 780     }
 781   }
 782 }
 783
 784 //_____________________________________________________________________________
 785 Int_t AliFlowEventSimple::CleanUpDeadTracks()
 786 {
 787   //remove tracks that have no flow tags set and cleanup the container
 788   //returns number of cleaned tracks
 789   Int_t ncleaned=0;
 790   for (Int_t i=0; i<fNumberOfTracks; i++)
 791   {
 792     AliFlowTrackSimple* track = static_cast<AliFlowTrackSimple*>(fTrackCollection->At(i));
 793     if (!track) continue;
 794     if (track->IsDead()) {delete track;track=NULL;ncleaned++;}
 795   }
 796   fTrackCollection->Compress(); //clean up empty slots
 797   fNumberOfTracks-=ncleaned; //update number of tracks
 798   return ncleaned;
 799 }
 800
 801 //_____________________________________________________________________________
 802 TF1* AliFlowEventSimple::SimplePtDepV2()
 803 {
 804   //return a standard pt dependent v2 formula, user has to clean up!
 805   return new TF1("StandardPtDepV2","((x<1.0)*(0.05/1.0)*x+(x>=1.0)*0.05)");
 806 }
 807
 808 //_____________________________________________________________________________
 809 TF1* AliFlowEventSimple::SimplePtSpectrum()
 810 {
 811   //return a standard pt spectrum, user has to clean up!
 812   return new TF1("StandardPtSpectrum","x*TMath::Exp(-pow(0.13957*0.13957+x*x,0.5)/0.4)",0.1,10.);
 813 }
 814
 815 //_____________________________________________________________________________
 816 void AliFlowEventSimple::ClearFast()
 817 {
 818   //clear the counter without deleting allocated objects so they can be reused
 819   fReferenceMultiplicity = 0;
 820   fNumberOfTracks = 0;
 821   fNumberOfRPs = 0;
 822   fMCReactionPlaneAngle = 0.0;
 823   fMCReactionPlaneAngleIsSet = kFALSE;
 824   fAfterBurnerPrecision = 0.001;
 825   fUserModified = kFALSE;
 826 }