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