cumulant tasks now use exchange container FlowEvent Task fills reaction plane
[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   fEventNSelTracksIntFlow(0),
47   fMCReactionPlaneAngle(0.),
48   fNumberOfTracksWrap(NULL),
49   fEventNSelTracksIntFlowWrap(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   fEventNSelTracksIntFlow(0),
61   fMCReactionPlaneAngle(0.),
62   fNumberOfTracksWrap(NULL),
63   fEventNSelTracksIntFlowWrap(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   fEventNSelTracksIntFlow(anEvent.fEventNSelTracksIntFlow),
77   fMCReactionPlaneAngle(anEvent.fMCReactionPlaneAngle),
78   fNumberOfTracksWrap(anEvent.fNumberOfTracksWrap),
79   fEventNSelTracksIntFlowWrap(anEvent.fEventNSelTracksIntFlowWrap),
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   fEventNSelTracksIntFlow = anEvent.fEventNSelTracksIntFlow;
92   fMCReactionPlaneAngle = anEvent.fMCReactionPlaneAngle;
93   fNumberOfTracksWrap = anEvent.fNumberOfTracksWrap; 
94   fEventNSelTracksIntFlowWrap = anEvent.fEventNSelTracksIntFlowWrap;
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 (fEventNSelTracksIntFlowWrap) delete fEventNSelTracksIntFlowWrap;
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->UseForIntegratedFlow()) 
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->UseForIntegratedFlow())
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 void AliFlowEventSimple::Print(Option_t *option) const
251 {
252   //   -*-*-*-*-*Print some global quantities for this histogram collection class *-*-*-*-*-*-*-*
253   //             ===============================================
254   //   printf( "TH1.Print Name  = %s, Entries= %d, Total sum= %g\n",GetName(),Int_t(fEntries),GetSumOfWeights());
255   printf( "Class.Print Name = %s, Total number of tracks= %d, Number of selected tracks= %d, MC EventPlaneAngle= %f",
256           GetName(),fNumberOfTracks, fEventNSelTracksIntFlow, fMCReactionPlaneAngle );
257
258   if (fTrackCollection) {  
259     fTrackCollection->Print(option);
260   }
261   else {
262     printf( "Empty track collection \n");
263   }
264 }
265
266 //----------------------------------------------------------------------- 
267  void AliFlowEventSimple::Browse(TBrowser *b)
268 {
269   if (!b) return;
270   if (!fNumberOfTracksWrap) {
271     fNumberOfTracksWrap = new TParameter<int>("fNumberOfTracks", fNumberOfTracks);
272     b->Add(fNumberOfTracksWrap);
273   }
274   if (!fEventNSelTracksIntFlowWrap) {
275     fEventNSelTracksIntFlowWrap = new TParameter<int>("fEventNSelTracksIntFlow", fEventNSelTracksIntFlow);
276     b->Add(fEventNSelTracksIntFlowWrap);
277   }
278   if (!fMCReactionPlaneAngleWrap) {
279     fMCReactionPlaneAngleWrap = new TParameter<double>(" fMCReactionPlaneAngle",  fMCReactionPlaneAngle);
280     b->Add( fMCReactionPlaneAngleWrap);
281   }
282   if (fTrackCollection) b->Add(fTrackCollection,"AliFlowTracksSimple");
283 }
284
285
286
287
288
289
290
291
292
293
294