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 //-----------------------------------------------------------------------
104 void AliFlowAnalysisWithLYZEventPlane::WriteHistograms(TString* outputFileName)
106 //store the final results in output .root file
108 TFile *output = new TFile(outputFileName->Data(),"RECREATE");
109 output->WriteObject(fHistList, "cobjLYZEP","SingleKey");
113 //-----------------------------------------------------------------------
114 void AliFlowAnalysisWithLYZEventPlane::Init() {
116 //Initialise all histograms
117 cout<<"---Analysis with Lee Yang Zeros Event Plane Method---"<<endl;
120 if (fSecondRunList) {
121 fSecondReDtheta = (TProfile*)fSecondRunList->FindObject("Second_FlowPro_ReDtheta_LYZ");
122 fHistList->Add(fSecondReDtheta);
124 fSecondImDtheta = (TProfile*)fSecondRunList->FindObject("Second_FlowPro_ImDtheta_LYZ");
125 fHistList->Add(fSecondImDtheta);
127 fFirstr0theta = (TProfile*)fSecondRunList->FindObject("First_FlowPro_r0theta_LYZ");
128 fHistList->Add(fFirstr0theta);
131 if (!fSecondReDtheta) {cout<<"fSecondReDtheta is NULL!"<<endl; }
132 if (!fSecondImDtheta) {cout<<"fSecondImDtheta is NULL!"<<endl; }
133 if (!fFirstr0theta) {cout<<"fFirstr0theta is NULL!"<<endl; }
137 fCommonHists = new AliFlowCommonHist("AliFlowCommonHistLYZEP");
138 fHistList->Add(fCommonHists);
140 fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsLYZEP");
141 fHistList->Add(fCommonHistsRes);
143 Int_t iNbinsPt = AliFlowCommonConstants::GetNbinsPt();
144 Double_t dPtMin = AliFlowCommonConstants::GetPtMin();
145 Double_t dPtMax = AliFlowCommonConstants::GetPtMax();
147 fHistProFlow = new TProfile("FlowPro_VPt_LYZEP","FlowPro_VPt_LYZEP",iNbinsPt,dPtMin,dPtMax);
148 fHistProFlow->SetXTitle("Pt");
149 fHistProFlow->SetYTitle("v2 (%)");
150 fHistList->Add(fHistProFlow);
152 fHistProWr = new TProfile("FlowPro_Wr_LYZEP","FlowPro_Wr_LYZEP",100,0.,0.25);
153 fHistProWr->SetXTitle("Q");
154 fHistProWr->SetYTitle("Wr");
155 fHistList->Add(fHistProWr);
157 fHistQsumforChi = new TH1F("Flow_QsumforChi_LYZEP","Flow_QsumforChi_LYZEP",3,-1.,2.);
158 fHistQsumforChi->SetXTitle("Qsum.X , Qsum.Y, Q2sum");
159 fHistQsumforChi->SetYTitle("value");
160 fHistList->Add(fHistQsumforChi);
162 fHistDeltaPhi = new TH1F("Flow_DeltaPhi_LYZEP","Flow_DeltaPhi_LYZEP",100,0.,3.14);
163 fHistDeltaPhi->SetXTitle("DeltaPhi");
164 fHistDeltaPhi->SetYTitle("Counts");
165 fHistList->Add(fHistDeltaPhi);
167 fHistPhiLYZ = new TH1F("Flow_PhiLYZ_LYZEP","Flow_PhiLYZ_LYZEP",100,0.,3.14);
168 fHistPhiLYZ->SetXTitle("Phi from LYZ");
169 fHistPhiLYZ->SetYTitle("Counts");
170 fHistList->Add(fHistPhiLYZ);
172 fHistPhiEP = new TH1F("Flow_PhiEP_LYZEP","Flow_PhiEP_LYZEP",100,0.,3.14);
173 fHistPhiEP->SetXTitle("Phi from EP");
174 fHistPhiEP->SetYTitle("Counts");
175 fHistList->Add(fHistPhiEP);
177 fEventNumber = 0; //set number of events to zero
181 //-----------------------------------------------------------------------
183 void AliFlowAnalysisWithLYZEventPlane::Make(AliFlowEventSimple* anEvent, AliFlowLYZEventPlane* aLYZEP) {
185 //Get the event plane and weight for each event
188 //fill control histograms
189 fCommonHists->FillControlHistograms(anEvent);
191 //get the Q vector from the FlowEvent
192 AliFlowVector vQ = anEvent->GetQ();
193 if (vQ.X()== 0. && vQ.Y()== 0. ) { cout<<"Q vector is NULL!"<<endl; }
194 //Weight with the multiplicity
197 if (vQ.GetMult()!=0.) {
198 dQX = vQ.X()/vQ.GetMult();
199 dQY = vQ.Y()/vQ.GetMult();
200 } else {cerr<<"vQ.GetMult() is zero!"<<endl; }
202 //cout<<"vQ("<<dQX<<","<<dQY<<")"<<endl;
204 //for chi calculation:
206 fHistQsumforChi->SetBinContent(1,fQsum->X());
207 fHistQsumforChi->SetBinContent(2,fQsum->Y());
209 fHistQsumforChi->SetBinContent(3,fQ2sum);
210 //cout<<"fQ2sum = "<<fQ2sum<<endl;
212 //call AliFlowLYZEventPlane::CalculateRPandW() here!
213 aLYZEP->CalculateRPandW(vQ);
215 Double_t dWR = aLYZEP->GetWR();
216 Double_t dRP = aLYZEP->GetPsi();
218 //fHistProWr->Fill(vQ.Mod(),dWR); //this weight is always positive
219 fHistPhiLYZ->Fill(dRP);
221 //plot difference between event plane from EP-method and LYZ-method
222 Double_t dRPEP = vQ.Phi()/2; //gives distribution from (0 to pi)
223 //Double_t dRPEP = 0.5*TMath::ATan2(vQ.Y(),vQ.X()); //gives distribution from (-pi/2 to pi/2)
224 //cout<<"dRPEP = "<< dRPEP <<endl;
225 fHistPhiEP->Fill(dRPEP);
227 Double_t dDeltaPhi = dRPEP - dRP;
228 if (dDeltaPhi < 0.) { dDeltaPhi += TMath::Pi(); } //to shift distribution from (-pi/2 to pi/2) to (0 to pi)
229 //cout<<"dDeltaPhi = "<<dDeltaPhi<<endl;
230 fHistDeltaPhi->Fill(dDeltaPhi);
233 Double_t dLow = TMath::Pi()/4.;
234 Double_t dHigh = 3.*(TMath::Pi()/4.);
235 if ((dDeltaPhi > dLow) && (dDeltaPhi < dHigh)){
236 dRP -= (TMath::Pi()/2);
238 cerr<<"*** dRP modified ***"<<endl;
240 fHistProWr->Fill(vQ.Mod(),dWR); //corrected weight
243 //loop over the tracks of the event
244 Int_t iNumberOfTracks = anEvent->NumberOfTracks();
245 for (Int_t i=0;i<iNumberOfTracks;i++)
247 AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
249 if (pTrack->UseForDifferentialFlow()) {
250 Double_t dPhi = pTrack->Phi();
251 //if (dPhi<0.) fPhi+=2*TMath::Pi();
253 Double_t dv2 = dWR * TMath::Cos(2*(dPhi-dRP));
254 Double_t dPt = pTrack->Pt();
256 fHistProFlow->Fill(dPt,100*dv2);
262 cout<<"@@@@@ "<<fEventNumber<<" events processed"<<endl;
266 //--------------------------------------------------------------------
267 void AliFlowAnalysisWithLYZEventPlane::Finish() {
269 //plot histograms etc.
270 cout<<"AliFlowAnalysisWithLYZEventPlane::Finish()"<<endl;
273 Double_t dJ01 = 2.405;
274 Int_t iNtheta = AliFlowLYZConstants::kTheta;
275 Int_t iNbinsPt = AliFlowCommonConstants::GetNbinsPt();
276 //set the event number
277 SetEventNumber((int)fCommonHists->GetHistMultOrig()->GetEntries());
278 //cout<<"number of events processed is "<<fEventNumber<<endl;
280 //set the sum of Q vectors
281 fQsum->Set(fHistQsumforChi->GetBinContent(1),fHistQsumforChi->GetBinContent(2));
282 SetQ2sum(fHistQsumforChi->GetBinContent(3));
284 //calculate dV the mean of dVtheta
285 Double_t dVtheta = 0;
287 for (Int_t theta=0;theta<iNtheta;theta++) {
288 Double_t dR0 = fFirstr0theta->GetBinContent(theta+1);
289 if (dR0!=0.) { dVtheta = dJ01/dR0 ;}
295 Double_t dSigma2 = 0;
297 if (fEventNumber!=0) {
298 *fQsum /= fEventNumber;
299 //cerr<<"fQsum->X() = "<<fQsum->X()<<endl;
300 //cerr<<"fQsum->Y() = "<<fQsum->Y()<<endl;
301 fQ2sum /= fEventNumber;
302 //cerr<<"fEventNumber = "<<fEventNumber<<endl;
303 //cerr<<"fQ2sum = "<<fQ2sum<<endl;
304 dSigma2 = fQ2sum - TMath::Power(fQsum->X(),2.) - TMath::Power(fQsum->Y(),2.) - TMath::Power(dV,2.); //BP eq. 62
305 //cerr<<"dSigma2"<<dSigma2<<endl;
306 if (dSigma2>0) dChi = dV/TMath::Sqrt(dSigma2);
308 fCommonHistsRes->FillChi(dChi);
309 cerr<<"dV = "<<dV<<" and chi = "<<dChi<<endl;
312 for(Int_t b=0;b<iNbinsPt;b++){
313 Double_t dv2pro = 0.;
314 Double_t dErr2difcomb = 0.;
315 Double_t dErrdifcomb = 0.;
317 dv2pro = fHistProFlow->GetBinContent(b);
319 for (Int_t theta=0;theta<iNtheta;theta++) {
320 Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
321 Int_t iNprime = TMath::Nint(fHistProFlow->GetBinEntries(b));
322 //cerr<<"iNprime = "<<iNprime<<endl;
324 Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
326 Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
328 dErr2difcomb += (TMath::Cos(dTheta)/(4*iNprime*TMath::BesselJ1(dJ01)*
329 TMath::BesselJ1(dJ01)))*
330 ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
331 (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
333 //else { cout<<"iNprime = 0"<<endl; }
336 if (dErr2difcomb!=0.) {
337 dErr2difcomb /= iNtheta;
338 dErrdifcomb = TMath::Sqrt(dErr2difcomb)*100;
339 //cerr<<"dErrdifcomb = "<<dErrdifcomb<<endl;
341 else {dErrdifcomb = 0.; }
344 fCommonHistsRes->FillDifferentialFlow(b, dv2pro, dErrdifcomb);
348 cout << "Profile Hist missing" << endl;
354 cout<<".....finished"<<endl;