added MC Reaction Plane information to FlowEventSimple
[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   fNumberOfTracksWrap(NULL),
48   fEventNSelTracksIntFlowWrap(NULL),
49   fMCReactionPlaneAngleWrap(NULL)
50 {
51   cout << "AliFlowEventSimple: Default constructor to be used only by root for io" << endl;
52 }
53
54 //-----------------------------------------------------------------------
55
56 AliFlowEventSimple::AliFlowEventSimple(Int_t aLenght):
57   fTrackCollection(NULL),
58   fNumberOfTracks(0),
59   fEventNSelTracksIntFlow(0),
60   fNumberOfTracksWrap(NULL),
61   fEventNSelTracksIntFlowWrap(NULL),
62   fMCReactionPlaneAngleWrap(NULL)
63 {
64   //constructor 
65   fTrackCollection =  new TObjArray(aLenght) ;
66 }
67
68 //-----------------------------------------------------------------------
69
70 AliFlowEventSimple::AliFlowEventSimple(const AliFlowEventSimple& anEvent):
71   TObject(),
72   fTrackCollection(anEvent.fTrackCollection),
73   fNumberOfTracks(anEvent.fNumberOfTracks),
74   fEventNSelTracksIntFlow(anEvent.fEventNSelTracksIntFlow),
75   fMCReactionPlaneAngle(anEvent.fMCReactionPlaneAngle),
76   fNumberOfTracksWrap(anEvent.fNumberOfTracksWrap),
77   fEventNSelTracksIntFlowWrap(anEvent.fEventNSelTracksIntFlowWrap),
78   fMCReactionPlaneAngleWrap(anEvent.fMCReactionPlaneAngleWrap)
79 {
80   //copy constructor 
81 }
82
83 //-----------------------------------------------------------------------
84
85 AliFlowEventSimple& AliFlowEventSimple::operator=(const AliFlowEventSimple& anEvent)
86 {
87   *fTrackCollection = *anEvent.fTrackCollection ;
88   fNumberOfTracks = anEvent.fNumberOfTracks;
89   fEventNSelTracksIntFlow = anEvent.fEventNSelTracksIntFlow;
90   fMCReactionPlaneAngle = anEvent.fMCReactionPlaneAngle;
91   fNumberOfTracksWrap = anEvent.fNumberOfTracksWrap; 
92   fEventNSelTracksIntFlowWrap = anEvent.fEventNSelTracksIntFlowWrap;
93   fMCReactionPlaneAngleWrap=anEvent.fMCReactionPlaneAngleWrap;
94
95   return *this;
96 }
97
98 //----------------------------------------------------------------------- 
99
100 AliFlowEventSimple::~AliFlowEventSimple()
101 {
102   //destructor
103   if (fTrackCollection) fTrackCollection->Delete(); delete fTrackCollection;
104   if (fNumberOfTracksWrap) delete fNumberOfTracksWrap;
105   if (fEventNSelTracksIntFlowWrap) delete fEventNSelTracksIntFlowWrap;
106   if (fMCReactionPlaneAngleWrap) delete fMCReactionPlaneAngleWrap;
107 }
108
109 //----------------------------------------------------------------------- 
110
111 AliFlowTrackSimple* AliFlowEventSimple::GetTrack(Int_t i)
112 {
113   //get track i from collection
114   AliFlowTrackSimple* pTrack = (AliFlowTrackSimple*)TrackCollection()->At(i) ;
115   return pTrack;
116 }
117
118 //-----------------------------------------------------------------------   
119 AliFlowVector AliFlowEventSimple::GetQ(Int_t n, TList *weightsList, Bool_t usePhiWeights, Bool_t usePtWeights, Bool_t useEtaWeights) 
120 {
121   // calculate Q-vector in harmonic n without weights (default harmonic n=2)  
122   Double_t dQX = 0.;
123   Double_t dQY = 0.;
124   AliFlowVector vQ;
125   vQ.Set(0.,0.);
126   
127   Int_t iOrder = n;
128   Double_t iUsedTracks = 0;
129   Double_t dPhi=0.;
130   Double_t dPt=0.;
131   Double_t dEta=0.;
132   
133   AliFlowTrackSimple* pTrack = NULL;
134  
135   Int_t nBinsPhi=0; 
136   Double_t dBinWidthPt=0.;
137   Double_t dPtMin=0.;
138   Double_t dBinWidthEta=0.;
139   Double_t dEtaMin=0.;
140  
141   Double_t wPhi=1.; // weight Phi  
142   Double_t wPt=1.;  // weight Pt 
143   Double_t wEta=1.; // weight Eta 
144   
145   TH1F *phiWeights = NULL;
146   TH1D *ptWeights  = NULL;
147   TH1D *etaWeights = NULL;
148
149   Double_t dSumOfWeightsToPower2 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 2)
150   Double_t dSumOfWeightsToPower3 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 3)
151   Double_t dSumOfWeightsToPower4 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 4)
152   Double_t dSumOfWeightsToPower5 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 5)
153   Double_t dSumOfWeightsToPower6 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 6)
154   Double_t dSumOfWeightsToPower7 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 7)
155   Double_t dSumOfWeightsToPower8 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 8) 
156
157   if(weightsList)
158   {
159    if(usePhiWeights)
160    {
161     phiWeights = dynamic_cast<TH1F *>(weightsList->FindObject("phi_weights"));
162     if(phiWeights) nBinsPhi = phiWeights->GetNbinsX();
163    }          
164    if(usePtWeights)
165    {
166     ptWeights = dynamic_cast<TH1D *>(weightsList->FindObject("pt_weights"));
167     if(ptWeights)
168     {
169      dBinWidthPt = ptWeights->GetBinWidth(1); // assuming that all bins have the same width
170      dPtMin = (ptWeights->GetXaxis())->GetXmin();
171     } 
172    }       
173    if(useEtaWeights)
174    {
175     etaWeights = dynamic_cast<TH1D *>(weightsList->FindObject("eta_weights"));
176     if(etaWeights)
177     {
178      dBinWidthEta = etaWeights->GetBinWidth(1); // assuming that all bins have the same width
179      dEtaMin = (etaWeights->GetXaxis())->GetXmin();
180     } 
181    }          
182   } // end of if(weightsList)
183   
184   // loop over tracks    
185   for(Int_t i=0;i<fNumberOfTracks;i++)                               
186   {
187    pTrack = (AliFlowTrackSimple*)TrackCollection()->At(i); 
188    if(pTrack)
189    {
190     if(pTrack->UseForIntegratedFlow()) 
191     {
192      dPhi = pTrack->Phi();
193      dPt  = pTrack->Pt();
194      dEta = pTrack->Eta();
195      
196      // determine Phi weight: (to be improved, I should here only access it + the treatment of gaps in the if statement)
197      if(phiWeights && nBinsPhi)
198      {
199       wPhi = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*nBinsPhi/TMath::TwoPi())));
200      }
201      // determine v'(pt) weight:    
202      if(ptWeights && dBinWidthPt)
203      {
204       wPt=ptWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-dPtMin)/dBinWidthPt))); 
205      }            
206      // determine v'(eta) weight:    
207      if(etaWeights && dBinWidthEta)
208      {
209       wEta=etaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-dEtaMin)/dBinWidthEta))); 
210      } 
211
212      // building up the weighted Q-vector:       
213      dQX += wPhi*wPt*wEta*TMath::Cos(iOrder*dPhi);
214      dQY += wPhi*wPt*wEta*TMath::Sin(iOrder*dPhi); 
215     
216      // weighted multiplicity:
217      iUsedTracks+=wPhi*wPt*wEta;
218     
219      // weights raised to various powers are summed up:
220      dSumOfWeightsToPower2+=pow(wPhi*wPt*wEta, 2); 
221      dSumOfWeightsToPower3+=pow(wPhi*wPt*wEta, 3); 
222      dSumOfWeightsToPower4+=pow(wPhi*wPt*wEta, 4); 
223      dSumOfWeightsToPower5+=pow(wPhi*wPt*wEta, 5); 
224      dSumOfWeightsToPower6+=pow(wPhi*wPt*wEta, 6); 
225      dSumOfWeightsToPower7+=pow(wPhi*wPt*wEta, 7); 
226      dSumOfWeightsToPower8+=pow(wPhi*wPt*wEta, 8); 
227      
228     } // end of if (pTrack->UseForIntegratedFlow())
229    } // end of if (pTrack)
230    else {cerr << "no particle!!!"<<endl;}
231   } // loop over particles
232     
233   vQ.Set(dQX,dQY);
234   vQ.SetMult(iUsedTracks);
235   vQ.SetSumOfWeightsToPower2(dSumOfWeightsToPower2);
236   vQ.SetSumOfWeightsToPower3(dSumOfWeightsToPower3);
237   vQ.SetSumOfWeightsToPower4(dSumOfWeightsToPower4);
238   vQ.SetSumOfWeightsToPower5(dSumOfWeightsToPower5);
239   vQ.SetSumOfWeightsToPower6(dSumOfWeightsToPower6);
240   vQ.SetSumOfWeightsToPower7(dSumOfWeightsToPower7);
241   vQ.SetSumOfWeightsToPower8(dSumOfWeightsToPower8);
242
243   return vQ;
244   
245 }
246
247 //----------------------------------------------------------------------- 
248 void AliFlowEventSimple::Print(Option_t *option) const
249 {
250   //   -*-*-*-*-*Print some global quantities for this histogram collection class *-*-*-*-*-*-*-*
251   //             ===============================================
252   //   printf( "TH1.Print Name  = %s, Entries= %d, Total sum= %g\n",GetName(),Int_t(fEntries),GetSumOfWeights());
253   printf( "Class.Print Name = %s, Total number of tracks= %d, Number of selected tracks= %d, MC EventPlaneAngle= %d",
254           GetName(),fNumberOfTracks, fEventNSelTracksIntFlow, fMCReactionPlaneAngle );
255
256   if (fTrackCollection) {  
257     fTrackCollection->Print(option);
258   }
259   else {
260     printf( "Empty track collection \n");
261   }
262 }
263
264 //----------------------------------------------------------------------- 
265  void AliFlowEventSimple::Browse(TBrowser *b)
266 {
267   if (!b) return;
268   if (!fNumberOfTracksWrap) {
269     fNumberOfTracksWrap = new TParameter<int>("fNumberOfTracks", fNumberOfTracks);
270     b->Add(fNumberOfTracksWrap);
271   }
272   if (!fEventNSelTracksIntFlowWrap) {
273     fEventNSelTracksIntFlowWrap = new TParameter<int>("fEventNSelTracksIntFlow", fEventNSelTracksIntFlow);
274     b->Add(fEventNSelTracksIntFlowWrap);
275   }
276   if (!fMCReactionPlaneAngleWrap) {
277      fMCReactionPlaneAngleWrap = new TParameter<double>(" fMCReactionPlaneAngle",  fMCReactionPlaneAngle);
278     b->Add( fMCReactionPlaneAngleWrap);
279   }
280   if (fTrackCollection) b->Add(fTrackCollection,"AliFlowTracksSimple");
281 }
282
283
284
285
286
287
288
289
290
291
292