1 /*************************************************************************
2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
20 //#define AliFlowAnalysisWithLYZEventPlane_cxx
22 #include "Riostream.h" //needed as include
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
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"
45 // AliFlowAnalysisWithLYZEventPlane:
47 // Class to do flow analysis with the event plane from the LYZ method
49 // author: N. van der Kolk (kolk@nikhef.nl)
52 ClassImp(AliFlowAnalysisWithLYZEventPlane)
54 //-----------------------------------------------------------------------
56 AliFlowAnalysisWithLYZEventPlane::AliFlowAnalysisWithLYZEventPlane():
59 fSecondReDtheta(NULL),
60 fSecondImDtheta(NULL),
66 fHistQsumforChi(NULL),
69 fHistDeltaPhihere(NULL),
75 fCommonHistsRes(NULL),
82 fQsum = new TVector2(); // flow vector sum
84 fHistList = new TList();
85 fSecondRunList = new TList();
90 //-----------------------------------------------------------------------
93 AliFlowAnalysisWithLYZEventPlane::~AliFlowAnalysisWithLYZEventPlane()
98 delete fSecondRunList;
102 //-----------------------------------------------------------------------
103 void AliFlowAnalysisWithLYZEventPlane::Init() {
105 //Initialise all histograms
106 cout<<"---Analysis with Lee Yang Zeros Event Plane Method---"<<endl;
109 if (fSecondRunList) {
110 fSecondReDtheta = (TProfile*)fSecondRunList->FindObject("Second_FlowPro_ReDtheta_LYZ");
111 fHistList->Add(fSecondReDtheta);
113 fSecondImDtheta = (TProfile*)fSecondRunList->FindObject("Second_FlowPro_ImDtheta_LYZ");
114 fHistList->Add(fSecondImDtheta);
116 fFirstr0theta = (TProfile*)fSecondRunList->FindObject("First_FlowPro_r0theta_LYZ");
117 fHistList->Add(fFirstr0theta);
120 if (!fSecondReDtheta) {cout<<"fSecondReDtheta is NULL!"<<endl; }
121 if (!fSecondImDtheta) {cout<<"fSecondImDtheta is NULL!"<<endl; }
122 if (!fFirstr0theta) {cout<<"fFirstr0theta is NULL!"<<endl; }
126 fCommonHists = new AliFlowCommonHist("AliFlowCommonHistLYZEP");
127 fHistList->Add(fCommonHists);
129 fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsLYZEP");
130 fHistList->Add(fCommonHistsRes);
132 Int_t iNbinsPt = AliFlowCommonConstants::GetNbinsPt();
133 Double_t dPtMin = AliFlowCommonConstants::GetPtMin();
134 Double_t dPtMax = AliFlowCommonConstants::GetPtMax();
136 fHistProFlow = new TProfile("FlowPro_VPt_LYZEP","FlowPro_VPt_LYZEP",iNbinsPt,dPtMin,dPtMax);
137 fHistProFlow->SetXTitle("Pt");
138 fHistProFlow->SetYTitle("v2 (%)");
139 fHistList->Add(fHistProFlow);
141 fHistProWr = new TProfile("FlowPro_Wr_LYZEP","FlowPro_Wr_LYZEP",100,0.,0.25);
142 fHistProWr->SetXTitle("Q");
143 fHistProWr->SetYTitle("Wr");
144 fHistList->Add(fHistProWr);
146 fHistQsumforChi = new TH1F("Flow_QsumforChi_LYZEP","Flow_QsumforChi_LYZEP",3,-1.,2.);
147 fHistQsumforChi->SetXTitle("Qsum.X , Qsum.Y, Q2sum");
148 fHistQsumforChi->SetYTitle("value");
149 fHistList->Add(fHistQsumforChi);
151 fHistDeltaPhi = new TH1F("Flow_DeltaPhi_LYZEP","Flow_DeltaPhi_LYZEP",100,0.,3.14);
152 fHistDeltaPhi->SetXTitle("DeltaPhi");
153 fHistDeltaPhi->SetYTitle("Counts");
154 fHistList->Add(fHistDeltaPhi);
156 fHistPhiLYZ = new TH1F("Flow_PhiLYZ_LYZEP","Flow_PhiLYZ_LYZEP",100,0.,3.14);
157 fHistPhiLYZ->SetXTitle("Phi from LYZ");
158 fHistPhiLYZ->SetYTitle("Counts");
159 fHistList->Add(fHistPhiLYZ);
161 fHistPhiEP = new TH1F("Flow_PhiEP_LYZEP","Flow_PhiEP_LYZEP",100,0.,3.14);
162 fHistPhiEP->SetXTitle("Phi from EP");
163 fHistPhiEP->SetYTitle("Counts");
164 fHistList->Add(fHistPhiEP);
166 fEventNumber = 0; //set number of events to zero
170 //-----------------------------------------------------------------------
172 void AliFlowAnalysisWithLYZEventPlane::Make(AliFlowEventSimple* anEvent, AliFlowLYZEventPlane* aLYZEP) {
174 //Get the event plane and weight for each event
177 //fill control histograms
178 fCommonHists->FillControlHistograms(anEvent);
180 //get the Q vector from the FlowEvent
181 AliFlowVector vQ = anEvent->GetQ();
182 if (vQ.X()== 0. && vQ.Y()== 0. ) { cout<<"Q vector is NULL!"<<endl; }
183 //Weight with the multiplicity
186 if (vQ.GetMult()!=0.) {
187 dQX = vQ.X()/vQ.GetMult();
188 dQY = vQ.Y()/vQ.GetMult();
189 } else {cerr<<"vQ.GetMult() is zero!"<<endl; }
191 //cout<<"vQ("<<dQX<<","<<dQY<<")"<<endl;
193 //for chi calculation:
195 fHistQsumforChi->SetBinContent(1,fQsum->X());
196 fHistQsumforChi->SetBinContent(2,fQsum->Y());
198 fHistQsumforChi->SetBinContent(3,fQ2sum);
199 //cout<<"fQ2sum = "<<fQ2sum<<endl;
201 //call AliFlowLYZEventPlane::CalculateRPandW() here!
202 aLYZEP->CalculateRPandW(vQ);
204 Double_t dWR = aLYZEP->GetWR();
205 Double_t dRP = aLYZEP->GetPsi();
207 //fHistProWr->Fill(vQ.Mod(),dWR); //this weight is always positive
208 fHistPhiLYZ->Fill(dRP);
210 //plot difference between event plane from EP-method and LYZ-method
211 Double_t dRPEP = vQ.Phi()/2; //gives distribution from (0 to pi)
212 //Double_t dRPEP = 0.5*TMath::ATan2(vQ.Y(),vQ.X()); //gives distribution from (-pi/2 to pi/2)
213 //cout<<"dRPEP = "<< dRPEP <<endl;
214 fHistPhiEP->Fill(dRPEP);
216 Double_t dDeltaPhi = dRPEP - dRP;
217 if (dDeltaPhi < 0.) { dDeltaPhi += TMath::Pi(); } //to shift distribution from (-pi/2 to pi/2) to (0 to pi)
218 //cout<<"dDeltaPhi = "<<dDeltaPhi<<endl;
219 fHistDeltaPhi->Fill(dDeltaPhi);
222 Double_t dLow = TMath::Pi()/4.;
223 Double_t dHigh = 3.*(TMath::Pi()/4.);
224 if ((dDeltaPhi > dLow) && (dDeltaPhi < dHigh)){
225 dRP -= (TMath::Pi()/2);
227 cerr<<"*** dRP modified ***"<<endl;
229 fHistProWr->Fill(vQ.Mod(),dWR); //corrected weight
232 //loop over the tracks of the event
233 Int_t iNumberOfTracks = anEvent->NumberOfTracks();
234 for (Int_t i=0;i<iNumberOfTracks;i++)
236 AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
238 if (pTrack->UseForDifferentialFlow()) {
239 Double_t dPhi = pTrack->Phi();
240 //if (dPhi<0.) fPhi+=2*TMath::Pi();
242 Double_t dv2 = dWR * TMath::Cos(2*(dPhi-dRP));
243 Double_t dPt = pTrack->Pt();
245 fHistProFlow->Fill(dPt,100*dv2);
251 cout<<"@@@@@ "<<fEventNumber<<" events processed"<<endl;
255 //--------------------------------------------------------------------
256 void AliFlowAnalysisWithLYZEventPlane::Finish() {
258 //plot histograms etc.
259 cout<<"AliFlowAnalysisWithLYZEventPlane::Finish()"<<endl;
262 Double_t dJ01 = 2.405;
263 Int_t iNtheta = AliFlowLYZConstants::kTheta;
264 Int_t iNbinsPt = AliFlowCommonConstants::GetNbinsPt();
265 //set the event number
266 SetEventNumber((int)fCommonHists->GetHistMultOrig()->GetEntries());
267 //cout<<"number of events processed is "<<fEventNumber<<endl;
269 //set the sum of Q vectors
270 fQsum->Set(fHistQsumforChi->GetBinContent(1),fHistQsumforChi->GetBinContent(2));
271 SetQ2sum(fHistQsumforChi->GetBinContent(3));
273 //calculate dV the mean of dVtheta
274 Double_t dVtheta = 0;
276 for (Int_t theta=0;theta<iNtheta;theta++) {
277 Double_t dR0 = fFirstr0theta->GetBinContent(theta+1);
278 if (dR0!=0.) { dVtheta = dJ01/dR0 ;}
284 Double_t dSigma2 = 0;
286 if (fEventNumber!=0) {
287 *fQsum /= fEventNumber;
288 //cerr<<"fQsum->X() = "<<fQsum->X()<<endl;
289 //cerr<<"fQsum->Y() = "<<fQsum->Y()<<endl;
290 fQ2sum /= fEventNumber;
291 //cerr<<"fEventNumber = "<<fEventNumber<<endl;
292 //cerr<<"fQ2sum = "<<fQ2sum<<endl;
293 dSigma2 = fQ2sum - TMath::Power(fQsum->X(),2.) - TMath::Power(fQsum->Y(),2.) - TMath::Power(dV,2.); //BP eq. 62
294 //cerr<<"dSigma2"<<dSigma2<<endl;
295 if (dSigma2>0) dChi = dV/TMath::Sqrt(dSigma2);
297 fCommonHistsRes->FillChi(dChi);
298 cerr<<"dV = "<<dV<<" and chi = "<<dChi<<endl;
301 for(Int_t b=0;b<iNbinsPt;b++){
302 Double_t dv2pro = 0.;
303 Double_t dErr2difcomb = 0.;
304 Double_t dErrdifcomb = 0.;
306 dv2pro = fHistProFlow->GetBinContent(b);
308 for (Int_t theta=0;theta<iNtheta;theta++) {
309 Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
310 Int_t iNprime = TMath::Nint(fHistProFlow->GetBinEntries(b));
311 //cerr<<"iNprime = "<<iNprime<<endl;
313 Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
315 Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
317 dErr2difcomb += (TMath::Cos(dTheta)/(4*iNprime*TMath::BesselJ1(dJ01)*
318 TMath::BesselJ1(dJ01)))*
319 ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
320 (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
322 //else { cout<<"iNprime = 0"<<endl; }
325 if (dErr2difcomb!=0.) {
326 dErr2difcomb /= iNtheta;
327 dErrdifcomb = TMath::Sqrt(dErr2difcomb)*100;
328 //cerr<<"dErrdifcomb = "<<dErrdifcomb<<endl;
330 else {dErrdifcomb = 0.; }
333 fCommonHistsRes->FillDifferentialFlow(b, dv2pro, dErrdifcomb);
337 cout << "Profile Hist missing" << endl;
343 cout<<".....finished"<<endl;