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