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