overloaded write method
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowAnalysisWithLYZEventPlane.cxx
1 /*************************************************************************
2 * Copyright(c) 1998-2008, 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 $Log$
18 */ 
19
20 //#define AliFlowAnalysisWithLYZEventPlane_cxx
21  
22 #include "Riostream.h"  //needed as include
23
24 #include "TFile.h"
25 #include "TList.h"
26 #include "TComplex.h"   //needed as include
27 #include "TCanvas.h"   //needed as include
28 #include "TLegend.h"   //needed as include
29 #include "TProfile.h"  //needed as include
30 #include "TVector2.h"
31
32 class TH1F;
33
34 #include "AliFlowLYZConstants.h"    //needed as include
35 #include "AliFlowCommonConstants.h" //needed as include
36 #include "AliFlowEventSimple.h"
37 #include "AliFlowTrackSimple.h"
38 #include "AliFlowCommonHist.h"
39 #include "AliFlowCommonHistResults.h"
40 #include "AliFlowLYZEventPlane.h"
41 #include "AliFlowAnalysisWithLYZEventPlane.h"
42
43 class AliFlowVector;
44
45 // AliFlowAnalysisWithLYZEventPlane:
46 //
47 // Class to do flow analysis with the event plane from the LYZ method
48 //
49 // author: N. van der Kolk (kolk@nikhef.nl)
50
51
52 ClassImp(AliFlowAnalysisWithLYZEventPlane)
53
54   //-----------------------------------------------------------------------
55  
56  AliFlowAnalysisWithLYZEventPlane::AliFlowAnalysisWithLYZEventPlane():
57    fHistList(NULL),
58    fSecondRunList(NULL),
59    fSecondReDtheta(NULL),
60    fSecondImDtheta(NULL),
61    fFirstr0theta(NULL),
62    fHistProFlow(NULL),
63    fHistProFlow2(NULL),
64    fHistProWr(NULL),
65    fHistProWrCorr(NULL),
66    fHistQsumforChi(NULL),
67    fHistDeltaPhi(NULL),
68    fHistDeltaPhi2(NULL),
69    fHistDeltaPhihere(NULL),
70    fHistPhiEP(NULL),
71    fHistPhiEPhere(NULL),
72    fHistPhiLYZ(NULL),
73    fHistPhiLYZ2(NULL),
74    fCommonHists(NULL),
75    fCommonHistsRes(NULL),
76    fEventNumber(0),
77    fQsum(NULL),
78    fQ2sum(0)
79 {
80
81   // Constructor.
82   fQsum = new TVector2();        // flow vector sum
83
84   fHistList = new TList();
85   fSecondRunList = new TList();
86 }
87
88  
89
90  //-----------------------------------------------------------------------
91
92
93  AliFlowAnalysisWithLYZEventPlane::~AliFlowAnalysisWithLYZEventPlane() 
94  {
95    //destructor
96    delete fQsum;
97    delete fHistList;
98    delete fSecondRunList;
99  }
100  
101
102 //-----------------------------------------------------------------------
103
104 void AliFlowAnalysisWithLYZEventPlane::WriteHistograms(TString* outputFileName)
105 {
106  //store the final results in output .root file
107
108   TFile *output = new TFile(outputFileName->Data(),"RECREATE");
109   output->WriteObject(fHistList, "cobjLYZEP","SingleKey");
110   delete output;
111 }
112
113 //-----------------------------------------------------------------------
114
115 void AliFlowAnalysisWithLYZEventPlane::WriteHistograms(TString outputFileName)
116 {
117  //store the final results in output .root file
118
119   TFile *output = new TFile(outputFileName.Data(),"RECREATE");
120   output->WriteObject(fHistList, "cobjLYZEP","SingleKey");
121   delete output;
122 }
123
124 //-----------------------------------------------------------------------
125 void AliFlowAnalysisWithLYZEventPlane::Init() {
126
127   //Initialise all histograms
128   cout<<"---Analysis with Lee Yang Zeros Event Plane Method---"<<endl;
129
130   //input histograms
131   if (fSecondRunList) {
132     fSecondReDtheta = (TProfile*)fSecondRunList->FindObject("Second_FlowPro_ReDtheta_LYZ");
133     fHistList->Add(fSecondReDtheta);
134
135     fSecondImDtheta = (TProfile*)fSecondRunList->FindObject("Second_FlowPro_ImDtheta_LYZ");
136     fHistList->Add(fSecondImDtheta);
137     
138     fFirstr0theta = (TProfile*)fSecondRunList->FindObject("First_FlowPro_r0theta_LYZ");
139     fHistList->Add(fFirstr0theta);
140
141     //warnings
142     if (!fSecondReDtheta) {cout<<"fSecondReDtheta is NULL!"<<endl; }
143     if (!fSecondImDtheta) {cout<<"fSecondImDtheta is NULL!"<<endl; }
144     if (!fFirstr0theta)   {cout<<"fFirstr0theta is NULL!"<<endl; }
145
146   }
147
148   fCommonHists = new AliFlowCommonHist("AliFlowCommonHistLYZEP");
149   fHistList->Add(fCommonHists);
150   
151   fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsLYZEP"); 
152   fHistList->Add(fCommonHistsRes); 
153     
154   Int_t iNbinsPt = AliFlowCommonConstants::GetNbinsPt();
155   Double_t  dPtMin = AliFlowCommonConstants::GetPtMin();             
156   Double_t  dPtMax = AliFlowCommonConstants::GetPtMax();
157
158   fHistProFlow = new TProfile("FlowPro_VPt_LYZEP","FlowPro_VPt_LYZEP",iNbinsPt,dPtMin,dPtMax);
159   fHistProFlow->SetXTitle("Pt");
160   fHistProFlow->SetYTitle("v2 (%)");
161   fHistList->Add(fHistProFlow);
162   
163   fHistProWr = new TProfile("FlowPro_Wr_LYZEP","FlowPro_Wr_LYZEP",100,0.,0.25);
164   fHistProWr->SetXTitle("Q");
165   fHistProWr->SetYTitle("Wr");
166   fHistList->Add(fHistProWr);
167
168   fHistQsumforChi = new TH1F("Flow_QsumforChi_LYZEP","Flow_QsumforChi_LYZEP",3,-1.,2.);
169   fHistQsumforChi->SetXTitle("Qsum.X , Qsum.Y, Q2sum");
170   fHistQsumforChi->SetYTitle("value");
171   fHistList->Add(fHistQsumforChi);
172
173   fHistDeltaPhi = new TH1F("Flow_DeltaPhi_LYZEP","Flow_DeltaPhi_LYZEP",100,0.,3.14);
174   fHistDeltaPhi->SetXTitle("DeltaPhi");
175   fHistDeltaPhi->SetYTitle("Counts");
176   fHistList->Add(fHistDeltaPhi);
177
178   fHistPhiLYZ = new TH1F("Flow_PhiLYZ_LYZEP","Flow_PhiLYZ_LYZEP",100,0.,3.14);
179   fHistPhiLYZ->SetXTitle("Phi from LYZ");
180   fHistPhiLYZ->SetYTitle("Counts");
181   fHistList->Add(fHistPhiLYZ);
182
183   fHistPhiEP = new TH1F("Flow_PhiEP_LYZEP","Flow_PhiEP_LYZEP",100,0.,3.14);
184   fHistPhiEP->SetXTitle("Phi from EP");
185   fHistPhiEP->SetYTitle("Counts");
186   fHistList->Add(fHistPhiEP);
187
188   fEventNumber = 0;  //set number of events to zero
189       
190
191  
192 //-----------------------------------------------------------------------
193  
194 void AliFlowAnalysisWithLYZEventPlane::Make(AliFlowEventSimple* anEvent, AliFlowLYZEventPlane* aLYZEP) {
195   
196   //Get the event plane and weight for each event
197   if (anEvent) {
198          
199     //fill control histograms     
200     fCommonHists->FillControlHistograms(anEvent);
201
202     //get the Q vector from the FlowEvent
203     AliFlowVector vQ = anEvent->GetQ(); 
204     if (vQ.X()== 0. && vQ.Y()== 0. ) { cout<<"Q vector is NULL!"<<endl; }
205     //Weight with the multiplicity
206     Double_t dQX = 0.;
207     Double_t dQY = 0.;
208     if (vQ.GetMult()!=0.) {
209       dQX = vQ.X()/vQ.GetMult();
210       dQY = vQ.Y()/vQ.GetMult();
211     } else {cerr<<"vQ.GetMult() is zero!"<<endl; }
212     vQ.Set(dQX,dQY);
213     //cout<<"vQ("<<dQX<<","<<dQY<<")"<<endl;
214
215     //for chi calculation:
216     *fQsum += vQ;
217     fHistQsumforChi->SetBinContent(1,fQsum->X());
218     fHistQsumforChi->SetBinContent(2,fQsum->Y());
219     fQ2sum += vQ.Mod2();
220     fHistQsumforChi->SetBinContent(3,fQ2sum);
221     //cout<<"fQ2sum = "<<fQ2sum<<endl;
222
223     //call AliFlowLYZEventPlane::CalculateRPandW() here!
224     aLYZEP->CalculateRPandW(vQ);
225
226     Double_t dWR = aLYZEP->GetWR();     
227     Double_t dRP = aLYZEP->GetPsi();
228
229     //fHistProWr->Fill(vQ.Mod(),dWR); //this weight is always positive
230     fHistPhiLYZ->Fill(dRP);   
231     
232     //plot difference between event plane from EP-method and LYZ-method
233     Double_t dRPEP = vQ.Phi()/2;                              //gives distribution from (0 to pi)
234     //Double_t dRPEP = 0.5*TMath::ATan2(vQ.Y(),vQ.X());       //gives distribution from (-pi/2 to pi/2)
235     //cout<<"dRPEP = "<< dRPEP <<endl;
236     fHistPhiEP->Fill(dRPEP);
237
238     Double_t dDeltaPhi = dRPEP - dRP;
239     if (dDeltaPhi < 0.) { dDeltaPhi += TMath::Pi(); }        //to shift distribution from (-pi/2 to pi/2) to (0 to pi)
240     //cout<<"dDeltaPhi = "<<dDeltaPhi<<endl;
241     fHistDeltaPhi->Fill(dDeltaPhi); 
242
243     //Flip sign of WR
244     Double_t dLow = TMath::Pi()/4.;
245     Double_t dHigh = 3.*(TMath::Pi()/4.);
246     if ((dDeltaPhi > dLow) && (dDeltaPhi < dHigh)){
247       dRP -= (TMath::Pi()/2);
248       dWR = -dWR;
249       cerr<<"*** dRP modified ***"<<endl;
250     }
251     fHistProWr->Fill(vQ.Mod(),dWR); //corrected weight
252        
253     //calculate flow
254     //loop over the tracks of the event
255     Int_t iNumberOfTracks = anEvent->NumberOfTracks(); 
256     for (Int_t i=0;i<iNumberOfTracks;i++) 
257       {
258         AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ; 
259         if (pTrack){
260           if (pTrack->UseForDifferentialFlow()) {
261             Double_t dPhi = pTrack->Phi();
262             //if (dPhi<0.) fPhi+=2*TMath::Pi();
263             //calculate flow v2:
264             Double_t dv2 = dWR * TMath::Cos(2*(dPhi-dRP));
265             Double_t dPt = pTrack->Pt();
266             //fill histogram
267             fHistProFlow->Fill(dPt,100*dv2);  
268           }  
269         }//track selected
270       }//loop over tracks
271           
272     fEventNumber++;
273     cout<<"@@@@@ "<<fEventNumber<<" events processed"<<endl;
274   }
275 }
276
277   //--------------------------------------------------------------------    
278 void AliFlowAnalysisWithLYZEventPlane::Finish() {
279    
280   //plot histograms etc. 
281   cout<<"AliFlowAnalysisWithLYZEventPlane::Finish()"<<endl;
282   
283   //constants:
284   Double_t  dJ01 = 2.405; 
285   Int_t iNtheta = AliFlowLYZConstants::kTheta;
286   Int_t iNbinsPt = AliFlowCommonConstants::GetNbinsPt();
287   //set the event number
288   SetEventNumber((int)fCommonHists->GetHistMultOrig()->GetEntries());
289   //cout<<"number of events processed is "<<fEventNumber<<endl;
290   
291   //set the sum of Q vectors
292   fQsum->Set(fHistQsumforChi->GetBinContent(1),fHistQsumforChi->GetBinContent(2));
293   SetQ2sum(fHistQsumforChi->GetBinContent(3));  
294
295   //calculate dV the mean of dVtheta
296   Double_t  dVtheta = 0; 
297   Double_t  dV = 0; 
298   for (Int_t theta=0;theta<iNtheta;theta++)     {
299     Double_t dR0 = fFirstr0theta->GetBinContent(theta+1); 
300     if (dR0!=0.) { dVtheta = dJ01/dR0 ;}
301     dV += dVtheta;
302   }
303   dV /= iNtheta;
304
305   //calculate dChi 
306   Double_t  dSigma2 = 0;
307   Double_t  dChi= 0;
308   if (fEventNumber!=0) {
309     *fQsum /= fEventNumber;
310     //cerr<<"fQsum->X() = "<<fQsum->X()<<endl;
311     //cerr<<"fQsum->Y() = "<<fQsum->Y()<<endl;
312     fQ2sum /= fEventNumber;
313     //cerr<<"fEventNumber = "<<fEventNumber<<endl;
314     //cerr<<"fQ2sum = "<<fQ2sum<<endl;
315     dSigma2 = fQ2sum - TMath::Power(fQsum->X(),2.) - TMath::Power(fQsum->Y(),2.) - TMath::Power(dV,2.);  //BP eq. 62
316     //cerr<<"dSigma2"<<dSigma2<<endl;
317     if (dSigma2>0) dChi = dV/TMath::Sqrt(dSigma2);
318     else dChi = -1.;
319     fCommonHistsRes->FillChi(dChi);
320     cerr<<"dV = "<<dV<<" and chi = "<<dChi<<endl;
321   }
322   
323   for(Int_t b=0;b<iNbinsPt;b++){
324     Double_t dv2pro = 0.;
325     Double_t dErr2difcomb = 0.;   
326     Double_t dErrdifcomb = 0.;
327     if(fHistProFlow) {
328       dv2pro = fHistProFlow->GetBinContent(b);
329       //calculate error
330       for (Int_t theta=0;theta<iNtheta;theta++) {
331         Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi(); 
332         Int_t iNprime = TMath::Nint(fHistProFlow->GetBinEntries(b));
333         //cerr<<"iNprime = "<<iNprime<<endl;
334         if (iNprime!=0) { 
335           Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
336                                            TMath::Cos(dTheta));
337           Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
338                                           TMath::Cos(dTheta));
339           dErr2difcomb += (TMath::Cos(dTheta)/(4*iNprime*TMath::BesselJ1(dJ01)*
340                                                  TMath::BesselJ1(dJ01)))*
341             ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) - 
342              (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
343         } //if !=0
344         //else { cout<<"iNprime = 0"<<endl; }
345       } //loop over theta
346       
347       if (dErr2difcomb!=0.) {
348         dErr2difcomb /= iNtheta;
349         dErrdifcomb = TMath::Sqrt(dErr2difcomb)*100;
350         //cerr<<"dErrdifcomb = "<<dErrdifcomb<<endl;
351       }
352       else {dErrdifcomb = 0.; }
353
354       //fill TH1D
355       fCommonHistsRes->FillDifferentialFlow(b, dv2pro, dErrdifcomb); 
356     
357     } //if fHistProFLow
358     else  {
359       cout << "Profile Hist missing" << endl;
360       break;
361     }
362     
363   } //loop over b
364
365   cout<<".....finished"<<endl;
366  }
367