small mods for separate task
[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 "AliFlowVector.h"
25 #include "AliFlowTrackSimple.h"
26 #include "AliFlowEventSimple.h"
27
28 /**************************************
29  * AliFlowEventSimple: A simple event *
30  *  for flow analysis                 * 
31  *                                    * 
32  * authors: Naomi van der Kolk        *
33  *           (kolk@nikhef.nl)         *  
34  *          Ante Bilandzic            *
35  *           (anteb@nikhef.nl)        * 
36  * ***********************************/
37  
38 ClassImp(AliFlowEventSimple)
39
40 //-----------------------------------------------------------------------
41
42 AliFlowEventSimple::AliFlowEventSimple():
43   fTrackCollection(NULL),
44   fNumberOfTracks(0),
45   fEventNSelTracksIntFlow(0)
46 {
47   cout << "AliFlowEventSimple: Default constructor to be used only by root for io" << endl;
48 }
49
50 //-----------------------------------------------------------------------
51
52 AliFlowEventSimple::AliFlowEventSimple(Int_t aLenght):
53   fTrackCollection(NULL),
54   fNumberOfTracks(0),
55   fEventNSelTracksIntFlow(0)
56 {
57   //constructor 
58   fTrackCollection =  new TObjArray(aLenght) ;
59 }
60
61 //-----------------------------------------------------------------------
62
63 AliFlowEventSimple::AliFlowEventSimple(const AliFlowEventSimple& anEvent):
64   TObject(),
65   fTrackCollection(anEvent.fTrackCollection),
66   fNumberOfTracks(anEvent.fNumberOfTracks),
67   fEventNSelTracksIntFlow(anEvent.fEventNSelTracksIntFlow)
68 {
69   //copy constructor 
70 }
71
72 //-----------------------------------------------------------------------
73
74 AliFlowEventSimple& AliFlowEventSimple::operator=(const AliFlowEventSimple& anEvent)
75 {
76   *fTrackCollection = *anEvent.fTrackCollection ;
77   fNumberOfTracks = anEvent.fNumberOfTracks;
78   fEventNSelTracksIntFlow = anEvent.fEventNSelTracksIntFlow;
79   
80   return *this;
81 }
82
83 //----------------------------------------------------------------------- 
84
85 AliFlowEventSimple::~AliFlowEventSimple()
86 {
87   //destructor
88   if (fTrackCollection) {
89     fTrackCollection->Delete() ; delete fTrackCollection ;
90   }
91   else { 
92     cout << "AliFlowEventSimple: Warning trying to delete track collections NULL pointer" << endl; 
93   }
94 }
95
96 //----------------------------------------------------------------------- 
97
98 AliFlowTrackSimple* AliFlowEventSimple::GetTrack(Int_t i)
99 {
100   //get track i from collection
101   AliFlowTrackSimple* pTrack = (AliFlowTrackSimple*)TrackCollection()->At(i) ;
102   return pTrack;
103 }
104
105 //-----------------------------------------------------------------------   
106 AliFlowVector AliFlowEventSimple::GetQ(Int_t n, TList *weightsList, Bool_t usePhiWeights, Bool_t usePtWeights, Bool_t useEtaWeights) 
107 {
108   // calculate Q-vector in harmonic n without weights (default harmonic n=2)  
109   Double_t dQX = 0.;
110   Double_t dQY = 0.;
111   AliFlowVector vQ;
112   vQ.Set(0.,0.);
113   
114   Int_t iOrder = n;
115   Double_t iUsedTracks = 0;
116   Double_t dPhi=0.;
117   Double_t dPt=0.;
118   Double_t dEta=0.;
119   
120   AliFlowTrackSimple* pTrack = NULL;
121  
122   Int_t nBinsPhi=0; 
123   Double_t dBinWidthPt=0.;
124   Double_t dPtMin=0.;
125   Double_t dBinWidthEta=0.;
126   Double_t dEtaMin=0.;
127  
128   Double_t wPhi=1.; // weight Phi  
129   Double_t wPt=1.;  // weight Pt 
130   Double_t wEta=1.; // weight Eta 
131   
132   TH1F *phiWeights = NULL;
133   TH1D *ptWeights  = NULL;
134   TH1D *etaWeights = NULL;
135
136   Double_t dSumOfWeightsToPower2 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 2)
137   Double_t dSumOfWeightsToPower3 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 3)
138   Double_t dSumOfWeightsToPower4 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 4)
139   Double_t dSumOfWeightsToPower5 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 5)
140   Double_t dSumOfWeightsToPower6 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 6)
141   Double_t dSumOfWeightsToPower7 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 7)
142   Double_t dSumOfWeightsToPower8 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 8) 
143
144   if(weightsList)
145   {
146    if(usePhiWeights)
147    {
148     phiWeights = dynamic_cast<TH1F *>(weightsList->FindObject("phi_weights"));
149     if(phiWeights) nBinsPhi = phiWeights->GetNbinsX();
150    }          
151    if(usePtWeights)
152    {
153     ptWeights = dynamic_cast<TH1D *>(weightsList->FindObject("pt_weights"));
154     if(ptWeights)
155     {
156      dBinWidthPt = ptWeights->GetBinWidth(1); // assuming that all bins have the same width
157      dPtMin = (ptWeights->GetXaxis())->GetXmin();
158     } 
159    }       
160    if(useEtaWeights)
161    {
162     etaWeights = dynamic_cast<TH1D *>(weightsList->FindObject("eta_weights"));
163     if(etaWeights)
164     {
165      dBinWidthEta = etaWeights->GetBinWidth(1); // assuming that all bins have the same width
166      dEtaMin = (etaWeights->GetXaxis())->GetXmin();
167     } 
168    }          
169   } // end of if(weightsList)
170   
171   // loop over tracks    
172   for(Int_t i=0;i<fNumberOfTracks;i++)                               
173   {
174    pTrack = (AliFlowTrackSimple*)TrackCollection()->At(i); 
175    if(pTrack)
176    {
177     if(pTrack->UseForIntegratedFlow()) 
178     {
179      dPhi = pTrack->Phi();
180      dPt  = pTrack->Pt();
181      dEta = pTrack->Eta();
182      
183      // determine Phi weight: (to be improved, I should here only access it + the treatment of gaps in the if statement)
184      if(phiWeights && nBinsPhi)
185      {
186       wPhi = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*nBinsPhi/TMath::TwoPi())));
187      }
188      // determine v'(pt) weight:    
189      if(ptWeights && dBinWidthPt)
190      {
191       wPt=ptWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-dPtMin)/dBinWidthPt))); 
192      }            
193      // determine v'(eta) weight:    
194      if(etaWeights && dBinWidthEta)
195      {
196       wEta=etaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-dEtaMin)/dBinWidthEta))); 
197      } 
198
199      // building up the weighted Q-vector:       
200      dQX += wPhi*wPt*wEta*TMath::Cos(iOrder*dPhi);
201      dQY += wPhi*wPt*wEta*TMath::Sin(iOrder*dPhi); 
202     
203      // weighted multiplicity:
204      iUsedTracks+=wPhi*wPt*wEta;
205     
206      // weights raised to various powers are summed up:
207      dSumOfWeightsToPower2+=pow(wPhi*wPt*wEta, 2); 
208      dSumOfWeightsToPower3+=pow(wPhi*wPt*wEta, 3); 
209      dSumOfWeightsToPower4+=pow(wPhi*wPt*wEta, 4); 
210      dSumOfWeightsToPower5+=pow(wPhi*wPt*wEta, 5); 
211      dSumOfWeightsToPower6+=pow(wPhi*wPt*wEta, 6); 
212      dSumOfWeightsToPower7+=pow(wPhi*wPt*wEta, 7); 
213      dSumOfWeightsToPower8+=pow(wPhi*wPt*wEta, 8); 
214      
215     } // end of if (pTrack->UseForIntegratedFlow())
216    } // end of if (pTrack)
217    else {cerr << "no particle!!!"<<endl;}
218   } // loop over particles
219     
220   vQ.Set(dQX,dQY);
221   vQ.SetMult(iUsedTracks);
222   vQ.SetSumOfWeightsToPower2(dSumOfWeightsToPower2);
223   vQ.SetSumOfWeightsToPower3(dSumOfWeightsToPower3);
224   vQ.SetSumOfWeightsToPower4(dSumOfWeightsToPower4);
225   vQ.SetSumOfWeightsToPower5(dSumOfWeightsToPower5);
226   vQ.SetSumOfWeightsToPower6(dSumOfWeightsToPower6);
227   vQ.SetSumOfWeightsToPower7(dSumOfWeightsToPower7);
228   vQ.SetSumOfWeightsToPower8(dSumOfWeightsToPower8);
229
230   return vQ;
231   
232 }
233
234
235
236
237
238
239
240
241
242