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