name change Int/Diff RP/POI
[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 #include "Riostream.h"
17 #include "TObjArray.h"
18 #include "TFile.h"
19 #include "TList.h"
20 #include "TMath.h"
21 #include "TH1F.h"
22 #include "TH1D.h"
23 #include "TProfile.h"
24 #include "TBrowser.h"
25 #include "AliFlowVector.h"
26 #include "AliFlowTrackSimple.h"
27 #include "AliFlowEventSimple.h"
28
29 /**************************************
30  * AliFlowEventSimple: A simple event *
31  *  for flow analysis                 * 
32  *                                    * 
33  * authors: Naomi van der Kolk        *
34  *           (kolk@nikhef.nl)         *  
35  *          Ante Bilandzic            *
36  *           (anteb@nikhef.nl)        * 
37  * ***********************************/
38  
39 ClassImp(AliFlowEventSimple)
40
41 //-----------------------------------------------------------------------
42
43 AliFlowEventSimple::AliFlowEventSimple():
44   fTrackCollection(NULL),
45   fNumberOfTracks(0),
46   fEventNSelTracksRP(0),
47   fMCReactionPlaneAngle(0.),
48   fNumberOfTracksWrap(NULL),
49   fEventNSelTracksRPWrap(NULL),
50   fMCReactionPlaneAngleWrap(NULL)
51 {
52   cout << "AliFlowEventSimple: Default constructor to be used only by root for io" << endl;
53 }
54
55 //-----------------------------------------------------------------------
56
57 AliFlowEventSimple::AliFlowEventSimple(Int_t aLenght):
58   fTrackCollection(NULL),
59   fNumberOfTracks(0),
60   fEventNSelTracksRP(0),
61   fMCReactionPlaneAngle(0.),
62   fNumberOfTracksWrap(NULL),
63   fEventNSelTracksRPWrap(NULL),
64   fMCReactionPlaneAngleWrap(NULL)
65 {
66   //constructor 
67   fTrackCollection =  new TObjArray(aLenght) ;
68 }
69
70 //-----------------------------------------------------------------------
71
72 AliFlowEventSimple::AliFlowEventSimple(const AliFlowEventSimple& anEvent):
73   TObject(),
74   fTrackCollection(anEvent.fTrackCollection),
75   fNumberOfTracks(anEvent.fNumberOfTracks),
76   fEventNSelTracksRP(anEvent.fEventNSelTracksRP),
77   fMCReactionPlaneAngle(anEvent.fMCReactionPlaneAngle),
78   fNumberOfTracksWrap(anEvent.fNumberOfTracksWrap),
79   fEventNSelTracksRPWrap(anEvent.fEventNSelTracksRPWrap),
80   fMCReactionPlaneAngleWrap(anEvent.fMCReactionPlaneAngleWrap)
81 {
82   //copy constructor 
83 }
84
85 //-----------------------------------------------------------------------
86
87 AliFlowEventSimple& AliFlowEventSimple::operator=(const AliFlowEventSimple& anEvent)
88 {
89   *fTrackCollection = *anEvent.fTrackCollection ;
90   fNumberOfTracks = anEvent.fNumberOfTracks;
91   fEventNSelTracksRP = anEvent.fEventNSelTracksRP;
92   fMCReactionPlaneAngle = anEvent.fMCReactionPlaneAngle;
93   fNumberOfTracksWrap = anEvent.fNumberOfTracksWrap; 
94   fEventNSelTracksRPWrap = anEvent.fEventNSelTracksRPWrap;
95   fMCReactionPlaneAngleWrap=anEvent.fMCReactionPlaneAngleWrap;
96
97   return *this;
98 }
99
100 //----------------------------------------------------------------------- 
101
102 AliFlowEventSimple::~AliFlowEventSimple()
103 {
104   //destructor
105   if (fTrackCollection) fTrackCollection->Delete(); delete fTrackCollection;
106   if (fNumberOfTracksWrap) delete fNumberOfTracksWrap;
107   if (fEventNSelTracksRPWrap) delete fEventNSelTracksRPWrap;
108   if (fMCReactionPlaneAngleWrap) delete fMCReactionPlaneAngleWrap;
109 }
110
111 //----------------------------------------------------------------------- 
112
113 AliFlowTrackSimple* AliFlowEventSimple::GetTrack(Int_t i)
114 {
115   //get track i from collection
116   AliFlowTrackSimple* pTrack = (AliFlowTrackSimple*)TrackCollection()->At(i) ;
117   return pTrack;
118 }
119
120 //-----------------------------------------------------------------------   
121 AliFlowVector AliFlowEventSimple::GetQ(Int_t n, TList *weightsList, Bool_t usePhiWeights, Bool_t usePtWeights, Bool_t useEtaWeights) 
122 {
123   // calculate Q-vector in harmonic n without weights (default harmonic n=2)  
124   Double_t dQX = 0.;
125   Double_t dQY = 0.;
126   AliFlowVector vQ;
127   vQ.Set(0.,0.);
128   
129   Int_t iOrder = n;
130   Double_t iUsedTracks = 0;
131   Double_t dPhi=0.;
132   Double_t dPt=0.;
133   Double_t dEta=0.;
134   
135   AliFlowTrackSimple* pTrack = NULL;
136  
137   Int_t nBinsPhi=0; 
138   Double_t dBinWidthPt=0.;
139   Double_t dPtMin=0.;
140   Double_t dBinWidthEta=0.;
141   Double_t dEtaMin=0.;
142  
143   Double_t wPhi=1.; // weight Phi  
144   Double_t wPt=1.;  // weight Pt 
145   Double_t wEta=1.; // weight Eta 
146   
147   TH1F *phiWeights = NULL;
148   TH1D *ptWeights  = NULL;
149   TH1D *etaWeights = NULL;
150
151   Double_t dSumOfWeightsToPower2 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 2)
152   Double_t dSumOfWeightsToPower3 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 3)
153   Double_t dSumOfWeightsToPower4 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 4)
154   Double_t dSumOfWeightsToPower5 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 5)
155   Double_t dSumOfWeightsToPower6 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 6)
156   Double_t dSumOfWeightsToPower7 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 7)
157   Double_t dSumOfWeightsToPower8 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 8) 
158
159   if(weightsList)
160   {
161    if(usePhiWeights)
162    {
163     phiWeights = dynamic_cast<TH1F *>(weightsList->FindObject("phi_weights"));
164     if(phiWeights) nBinsPhi = phiWeights->GetNbinsX();
165    }          
166    if(usePtWeights)
167    {
168     ptWeights = dynamic_cast<TH1D *>(weightsList->FindObject("pt_weights"));
169     if(ptWeights)
170     {
171      dBinWidthPt = ptWeights->GetBinWidth(1); // assuming that all bins have the same width
172      dPtMin = (ptWeights->GetXaxis())->GetXmin();
173     } 
174    }       
175    if(useEtaWeights)
176    {
177     etaWeights = dynamic_cast<TH1D *>(weightsList->FindObject("eta_weights"));
178     if(etaWeights)
179     {
180      dBinWidthEta = etaWeights->GetBinWidth(1); // assuming that all bins have the same width
181      dEtaMin = (etaWeights->GetXaxis())->GetXmin();
182     } 
183    }          
184   } // end of if(weightsList)
185   
186   // loop over tracks    
187   for(Int_t i=0;i<fNumberOfTracks;i++)                               
188   {
189    pTrack = (AliFlowTrackSimple*)TrackCollection()->At(i); 
190    if(pTrack)
191    {
192     if(pTrack->InRPSelection()) 
193     {
194      dPhi = pTrack->Phi();
195      dPt  = pTrack->Pt();
196      dEta = pTrack->Eta();
197      
198      // determine Phi weight: (to be improved, I should here only access it + the treatment of gaps in the if statement)
199      if(phiWeights && nBinsPhi)
200      {
201       wPhi = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*nBinsPhi/TMath::TwoPi())));
202      }
203      // determine v'(pt) weight:    
204      if(ptWeights && dBinWidthPt)
205      {
206       wPt=ptWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-dPtMin)/dBinWidthPt))); 
207      }            
208      // determine v'(eta) weight:    
209      if(etaWeights && dBinWidthEta)
210      {
211       wEta=etaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-dEtaMin)/dBinWidthEta))); 
212      } 
213
214      // building up the weighted Q-vector:       
215      dQX += wPhi*wPt*wEta*TMath::Cos(iOrder*dPhi);
216      dQY += wPhi*wPt*wEta*TMath::Sin(iOrder*dPhi); 
217     
218      // weighted multiplicity:
219      iUsedTracks+=wPhi*wPt*wEta;
220     
221      // weights raised to various powers are summed up:
222      dSumOfWeightsToPower2+=pow(wPhi*wPt*wEta, 2); 
223      dSumOfWeightsToPower3+=pow(wPhi*wPt*wEta, 3); 
224      dSumOfWeightsToPower4+=pow(wPhi*wPt*wEta, 4); 
225      dSumOfWeightsToPower5+=pow(wPhi*wPt*wEta, 5); 
226      dSumOfWeightsToPower6+=pow(wPhi*wPt*wEta, 6); 
227      dSumOfWeightsToPower7+=pow(wPhi*wPt*wEta, 7); 
228      dSumOfWeightsToPower8+=pow(wPhi*wPt*wEta, 8); 
229      
230     } // end of if (pTrack->InRPSelection())
231    } // end of if (pTrack)
232    else {cerr << "no particle!!!"<<endl;}
233   } // loop over particles
234     
235   vQ.Set(dQX,dQY);
236   vQ.SetMult(iUsedTracks);
237   vQ.SetSumOfWeightsToPower2(dSumOfWeightsToPower2);
238   vQ.SetSumOfWeightsToPower3(dSumOfWeightsToPower3);
239   vQ.SetSumOfWeightsToPower4(dSumOfWeightsToPower4);
240   vQ.SetSumOfWeightsToPower5(dSumOfWeightsToPower5);
241   vQ.SetSumOfWeightsToPower6(dSumOfWeightsToPower6);
242   vQ.SetSumOfWeightsToPower7(dSumOfWeightsToPower7);
243   vQ.SetSumOfWeightsToPower8(dSumOfWeightsToPower8);
244
245   return vQ;
246   
247 }
248
249 //-----------------------------------------------------------------------   
250 AliFlowVector AliFlowEventSimple::GetQsub(Double_t etaMin, Double_t etaMax, Int_t n, TList *weightsList, Bool_t usePhiWeights, Bool_t usePtWeights, Bool_t useEtaWeights) 
251 {
252   //for eta subevents
253   
254   // calculate Q-vector in harmonic n without weights (default harmonic n=2)  
255   Double_t dQX = 0.;
256   Double_t dQY = 0.;
257   AliFlowVector vQ;
258   vQ.Set(0.,0.);
259   
260   Int_t iOrder = n;
261   Double_t iUsedTracks = 0;
262   Double_t dPhi = 0.;
263   Double_t dPt  = 0.;
264   Double_t dEta = 0.;
265   
266   AliFlowTrackSimple* pTrack = NULL;
267  
268   Int_t    nBinsPhi    = 0; 
269   Double_t dBinWidthPt = 0.;
270   Double_t dPtMin      = 0.;
271   Double_t dBinWidthEta= 0.;
272   Double_t dEtaMin     = 0.;
273  
274   Double_t wPhi = 1.;  // weight Phi  
275   Double_t wPt  = 1.;  // weight Pt 
276   Double_t wEta = 1.;  // weight Eta 
277   
278   TH1F *phiWeights = NULL;
279   TH1D *ptWeights  = NULL;
280   TH1D *etaWeights = NULL;
281
282   Double_t dSumOfWeightsToPower2 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 2)
283   Double_t dSumOfWeightsToPower3 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 3)
284   Double_t dSumOfWeightsToPower4 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 4)
285   Double_t dSumOfWeightsToPower5 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 5)
286   Double_t dSumOfWeightsToPower6 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 6)
287   Double_t dSumOfWeightsToPower7 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 7)
288   Double_t dSumOfWeightsToPower8 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 8) 
289
290   if(weightsList)
291     {
292       if(usePhiWeights)
293         {
294     phiWeights = dynamic_cast<TH1F *>(weightsList->FindObject("phi_weights"));
295     if(phiWeights) nBinsPhi = phiWeights->GetNbinsX();
296    }          
297    if(usePtWeights)
298    {
299     ptWeights = dynamic_cast<TH1D *>(weightsList->FindObject("pt_weights"));
300     if(ptWeights)
301     {
302      dBinWidthPt = ptWeights->GetBinWidth(1); // assuming that all bins have the same width
303      dPtMin = (ptWeights->GetXaxis())->GetXmin();
304     } 
305    }       
306    if(useEtaWeights)
307    {
308     etaWeights = dynamic_cast<TH1D *>(weightsList->FindObject("eta_weights"));
309     if(etaWeights)
310     {
311      dBinWidthEta = etaWeights->GetBinWidth(1); // assuming that all bins have the same width
312      dEtaMin = (etaWeights->GetXaxis())->GetXmin();
313     } 
314    }          
315   } // end of if(weightsList)
316   
317   // loop over tracks    
318   for(Int_t i=0;i<fNumberOfTracks;i++)                               
319   {
320    pTrack = (AliFlowTrackSimple*)TrackCollection()->At(i); 
321    if(pTrack)
322    {
323     if(pTrack->InRPSelection())
324     {
325      dPhi = pTrack->Phi();
326      dPt  = pTrack->Pt();
327      dEta = pTrack->Eta();
328      if (dEta>etaMin && dEta<etaMax) {
329        // determine Phi weight: (to be improved, I should here only access it + the treatment of gaps in the if statement)
330        if(phiWeights && nBinsPhi)
331          {
332            wPhi = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*nBinsPhi/TMath::TwoPi())));
333          }
334        // determine v'(pt) weight:    
335        if(ptWeights && dBinWidthPt)
336          {
337            wPt=ptWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-dPtMin)/dBinWidthPt))); 
338          }            
339        // determine v'(eta) weight:    
340        if(etaWeights && dBinWidthEta)
341          {
342            wEta=etaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-dEtaMin)/dBinWidthEta))); 
343          } 
344
345        // building up the weighted Q-vector:       
346        dQX += wPhi*wPt*wEta*TMath::Cos(iOrder*dPhi);
347        dQY += wPhi*wPt*wEta*TMath::Sin(iOrder*dPhi); 
348     
349        // weighted multiplicity:
350        iUsedTracks+=wPhi*wPt*wEta;
351     
352        // weights raised to various powers are summed up:
353        dSumOfWeightsToPower2+=pow(wPhi*wPt*wEta, 2); 
354        dSumOfWeightsToPower3+=pow(wPhi*wPt*wEta, 3); 
355        dSumOfWeightsToPower4+=pow(wPhi*wPt*wEta, 4); 
356        dSumOfWeightsToPower5+=pow(wPhi*wPt*wEta, 5); 
357        dSumOfWeightsToPower6+=pow(wPhi*wPt*wEta, 6); 
358        dSumOfWeightsToPower7+=pow(wPhi*wPt*wEta, 7); 
359        dSumOfWeightsToPower8+=pow(wPhi*wPt*wEta, 8); 
360      } // end of if dEta in eta range
361     } // end of if (pTrack->InRPSelection())
362    } // end of if (pTrack)
363    else {cerr << "no particle!!!"<<endl;}
364   } // loop over particles
365     
366   vQ.Set(dQX,dQY);
367   vQ.SetMult(iUsedTracks);
368   vQ.SetSumOfWeightsToPower2(dSumOfWeightsToPower2);
369   vQ.SetSumOfWeightsToPower3(dSumOfWeightsToPower3);
370   vQ.SetSumOfWeightsToPower4(dSumOfWeightsToPower4);
371   vQ.SetSumOfWeightsToPower5(dSumOfWeightsToPower5);
372   vQ.SetSumOfWeightsToPower6(dSumOfWeightsToPower6);
373   vQ.SetSumOfWeightsToPower7(dSumOfWeightsToPower7);
374   vQ.SetSumOfWeightsToPower8(dSumOfWeightsToPower8);
375
376   return vQ;
377   
378 }
379
380
381 //----------------------------------------------------------------------- 
382 void AliFlowEventSimple::Print(Option_t *option) const
383 {
384   //   -*-*-*-*-*Print some global quantities for this histogram collection class *-*-*-*-*-*-*-*
385   //             ===============================================
386   //   printf( "TH1.Print Name  = %s, Entries= %d, Total sum= %g\n",GetName(),Int_t(fEntries),GetSumOfWeights());
387   printf( "Class.Print Name = %s, Total number of tracks= %d, Number of selected tracks= %d, MC EventPlaneAngle= %f",
388           GetName(),fNumberOfTracks, fEventNSelTracksRP, fMCReactionPlaneAngle );
389
390   if (fTrackCollection) {  
391     fTrackCollection->Print(option);
392   }
393   else {
394     printf( "Empty track collection \n");
395   }
396 }
397
398 //----------------------------------------------------------------------- 
399  void AliFlowEventSimple::Browse(TBrowser *b)
400 {
401   if (!b) return;
402   if (!fNumberOfTracksWrap) {
403     fNumberOfTracksWrap = new TParameter<int>("fNumberOfTracks", fNumberOfTracks);
404     b->Add(fNumberOfTracksWrap);
405   }
406   if (!fEventNSelTracksRPWrap) {
407     fEventNSelTracksRPWrap = new TParameter<int>("fEventNSelTracksRP", fEventNSelTracksRP);
408     b->Add(fEventNSelTracksRPWrap);
409   }
410   if (!fMCReactionPlaneAngleWrap) {
411     fMCReactionPlaneAngleWrap = new TParameter<double>(" fMCReactionPlaneAngle",  fMCReactionPlaneAngle);
412     b->Add( fMCReactionPlaneAngleWrap);
413   }
414   if (fTrackCollection) b->Add(fTrackCollection,"AliFlowTracksSimple");
415 }
416
417