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