]>
Commit | Line | Data |
---|---|---|
f1d945a1 | 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 | ||
f1d945a1 | 16 | |
17 | #include "Riostream.h" //needed as include | |
ebaded83 | 18 | #include "TObject.h" //needed as include |
f1d945a1 | 19 | #include "AliFlowCommonConstants.h" //needed as include |
20 | #include "AliFlowLYZConstants.h" //needed as include | |
21 | #include "AliFlowAnalysisWithLeeYangZeros.h" | |
22 | #include "AliFlowLYZHist1.h" //needed as include | |
23 | #include "AliFlowLYZHist2.h" //needed as include | |
24 | #include "AliFlowCommonHist.h" //needed as include | |
25 | #include "AliFlowCommonHistResults.h" //needed as include | |
26 | #include "AliFlowEventSimple.h" //needed as include | |
27 | #include "AliFlowTrackSimple.h" //needed as include | |
28 | ||
29 | class AliFlowVector; | |
30 | ||
31 | #include "TMath.h" //needed as include | |
32 | #include "TFile.h" //needed as include | |
88e00a8a | 33 | #include "TList.h" |
f1d945a1 | 34 | |
35 | class TComplex; | |
36 | class TProfile; | |
37 | class TH1F; | |
38 | class TH1D; | |
39 | ||
40 | ||
41 | //Description: Maker to analyze Flow using the LeeYangZeros method | |
42 | // Equation numbers are from Big Paper (BP): Nucl. Phys. A 727, 373 (2003) | |
43 | // Practical Guide (PG): J. Phys. G: Nucl. Part. Phys. 30, S1213 (2004) | |
44 | // Adapted from StFlowLeeYangZerosMaker.cxx | |
45 | // by Markus Oldenberg and Art Poskanzer, LBNL | |
46 | // with advice from Jean-Yves Ollitrault and Nicolas Borghini | |
47 | // | |
48 | //Author: Naomi van der Kolk (kolk@nikhef.nl) | |
49 | ||
50 | ||
51 | ||
52 | ClassImp(AliFlowAnalysisWithLeeYangZeros) | |
53 | ||
54 | //----------------------------------------------------------------------- | |
55 | ||
56 | AliFlowAnalysisWithLeeYangZeros::AliFlowAnalysisWithLeeYangZeros(): | |
882ffd6a | 57 | fQsum(NULL), |
58 | fQ2sum(0), | |
59 | fEventNumber(0), | |
60 | fFirstRun(kTRUE), | |
61 | fUseSum(kTRUE), | |
62 | fDoubleLoop(kFALSE), | |
63 | fDebug(kFALSE), | |
88e00a8a | 64 | fHistList(NULL), |
9d062fe3 | 65 | fFirstRunList(NULL), |
882ffd6a | 66 | fHistProVtheta(NULL), |
67 | fHistProVeta(NULL), | |
68 | fHistProVPt(NULL), | |
69 | fHistProR0theta(NULL), | |
70 | fHistProReDenom(NULL), | |
71 | fHistProImDenom(NULL), | |
72 | fHistProReDtheta(NULL), | |
73 | fHistProImDtheta(NULL), | |
9d062fe3 | 74 | fHistQsumforChi(NULL), |
882ffd6a | 75 | fCommonHists(NULL), |
76 | fCommonHistsRes(NULL) | |
f1d945a1 | 77 | |
78 | { | |
79 | //default constructor | |
80 | if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::AliFlowAnalysisWithLeeYangZeros default constructor****"<<endl; | |
81 | ||
88e00a8a | 82 | fHistList = new TList(); |
9d062fe3 | 83 | fFirstRunList = new TList(); |
f1d945a1 | 84 | |
85 | for(Int_t i = 0;i<5;i++) | |
86 | { | |
87 | fHist1[i]=0; | |
88 | fHist2[i]=0; | |
89 | } | |
90 | ||
0092f3c2 | 91 | fQsum = new TVector2(); |
f1d945a1 | 92 | |
93 | } | |
94 | ||
f1d945a1 | 95 | //----------------------------------------------------------------------- |
96 | ||
97 | ||
98 | AliFlowAnalysisWithLeeYangZeros::~AliFlowAnalysisWithLeeYangZeros() | |
99 | { | |
100 | //default destructor | |
101 | if (fDebug) cout<<"****~AliFlowAnalysisWithLeeYangZeros****"<<endl; | |
0092f3c2 | 102 | delete fQsum; |
88e00a8a | 103 | delete fHistList; |
9d062fe3 | 104 | delete fFirstRunList; |
105 | ||
f1d945a1 | 106 | } |
107 | ||
7ebf0358 | 108 | //----------------------------------------------------------------------- |
109 | ||
110 | void AliFlowAnalysisWithLeeYangZeros::WriteHistograms(TString* outputFileName) | |
111 | { | |
112 | //store the final results in output .root file | |
ebaded83 | 113 | |
114 | TFile *output = new TFile(outputFileName->Data(),"RECREATE"); | |
115 | if (GetFirstRun()) { | |
116 | output->WriteObject(fHistList, "cobjLYZ1","SingleKey"); | |
117 | } | |
118 | else { | |
119 | output->WriteObject(fHistList, "cobjLYZ2","SingleKey"); | |
120 | } | |
121 | delete output; | |
7ebf0358 | 122 | } |
123 | ||
124 | //----------------------------------------------------------------------- | |
f1d945a1 | 125 | |
b0fda271 | 126 | void AliFlowAnalysisWithLeeYangZeros::WriteHistograms(TString outputFileName) |
127 | { | |
128 | //store the final results in output .root file | |
129 | ||
130 | TFile *output = new TFile(outputFileName.Data(),"RECREATE"); | |
131 | if (GetFirstRun()) { | |
132 | output->WriteObject(fHistList, "cobjLYZ1","SingleKey"); | |
133 | } | |
134 | else { | |
135 | output->WriteObject(fHistList, "cobjLYZ2","SingleKey"); | |
136 | } | |
137 | delete output; | |
138 | } | |
139 | ||
140 | //----------------------------------------------------------------------- | |
141 | ||
f1d945a1 | 142 | Bool_t AliFlowAnalysisWithLeeYangZeros::Init() |
143 | { | |
144 | //init method | |
145 | if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::Init()****"<<endl; | |
146 | ||
f1d945a1 | 147 | // Book histograms |
882ffd6a | 148 | Int_t iNtheta = AliFlowLYZConstants::kTheta; |
149 | Int_t iNbinsPt = AliFlowCommonConstants::GetNbinsPt(); | |
150 | Int_t iNbinsEta = AliFlowCommonConstants::GetNbinsEta(); | |
151 | ||
152 | Double_t dPtMin = AliFlowCommonConstants::GetPtMin(); | |
153 | Double_t dPtMax = AliFlowCommonConstants::GetPtMax(); | |
154 | Double_t dEtaMin = AliFlowCommonConstants::GetEtaMin(); | |
155 | Double_t dEtaMax = AliFlowCommonConstants::GetEtaMax(); | |
9d062fe3 | 156 | |
157 | //for control histograms | |
7ebf0358 | 158 | if (fFirstRun){ fCommonHists = new AliFlowCommonHist("AliFlowCommonHistLYZ1"); |
9d062fe3 | 159 | fHistList->Add(fCommonHists); |
7ebf0358 | 160 | fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsLYZ1"); |
161 | fHistList->Add(fCommonHistsRes); | |
162 | } | |
163 | else { fCommonHists = new AliFlowCommonHist("AliFlowCommonHistLYZ2"); | |
164 | fHistList->Add(fCommonHists); | |
165 | fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsLYZ2"); | |
9d062fe3 | 166 | fHistList->Add(fCommonHistsRes); |
7ebf0358 | 167 | } |
168 | ||
169 | fHistQsumforChi = new TH1F("Flow_QsumforChi_LYZ","Flow_QsumforChi_LYZ",3,-1.,2.); | |
9d062fe3 | 170 | fHistQsumforChi->SetXTitle("Qsum.X , Qsum.Y, Q2sum"); |
171 | fHistQsumforChi->SetYTitle("value"); | |
172 | fHistList->Add(fHistQsumforChi); | |
f1d945a1 | 173 | |
174 | //for first loop over events | |
175 | if (fFirstRun){ | |
882ffd6a | 176 | fHistProR0theta = new TProfile("First_FlowPro_r0theta_LYZ","First_FlowPro_r0theta_LYZ",iNtheta,-0.5,iNtheta-0.5); |
f1d945a1 | 177 | fHistProR0theta->SetXTitle("#theta"); |
178 | fHistProR0theta->SetYTitle("r_{0}^{#theta}"); | |
88e00a8a | 179 | fHistList->Add(fHistProR0theta); |
f1d945a1 | 180 | |
882ffd6a | 181 | fHistProVtheta = new TProfile("First_FlowPro_Vtheta_LYZ","First_FlowPro_Vtheta_LYZ",iNtheta,-0.5,iNtheta-0.5); |
f1d945a1 | 182 | fHistProVtheta->SetXTitle("#theta"); |
183 | fHistProVtheta->SetYTitle("V_{n}^{#theta}"); | |
88e00a8a | 184 | fHistList->Add(fHistProVtheta); |
f1d945a1 | 185 | |
186 | //class AliFlowLYZHist1 defines the histograms: fHistProGtheta, fHistProReGtheta, fHistProImGtheta, fHistProR0theta | |
882ffd6a | 187 | for (Int_t theta=0;theta<iNtheta;theta++) { |
9d062fe3 | 188 | TString name = "AliFlowLYZHist1_"; |
189 | name += theta; | |
190 | fHist1[theta]=new AliFlowLYZHist1(theta, name); | |
191 | fHistList->Add(fHist1[theta]); | |
f1d945a1 | 192 | } |
9d062fe3 | 193 | |
f1d945a1 | 194 | } |
195 | //for second loop over events | |
196 | else { | |
882ffd6a | 197 | fHistProReDenom = new TProfile("Second_FlowPro_ReDenom_LYZ","Second_FlowPro_ReDenom_LYZ" , iNtheta, -0.5, iNtheta-0.5); |
f1d945a1 | 198 | fHistProReDenom->SetXTitle("#theta"); |
199 | fHistProReDenom->SetYTitle("Re(Q^{#theta}e^{ir_{0}^{#theta}Q^{#theta}})"); | |
88e00a8a | 200 | fHistList->Add(fHistProReDenom); |
f1d945a1 | 201 | |
882ffd6a | 202 | fHistProImDenom = new TProfile("Second_FlowPro_ImDenom_LYZ","Second_FlowPro_ImDenom_LYZ" , iNtheta, -0.5, iNtheta-0.5); |
f1d945a1 | 203 | fHistProImDenom->SetXTitle("#theta"); |
204 | fHistProImDenom->SetYTitle("Im(Q^{#theta}e^{ir_{0}^{#theta}Q^{#theta}})"); | |
88e00a8a | 205 | fHistList->Add(fHistProImDenom); |
f1d945a1 | 206 | |
882ffd6a | 207 | fHistProVeta = new TProfile("Second_FlowPro_Veta_LYZ","Second_FlowPro_Veta_LYZ",iNbinsEta,dEtaMin,dEtaMax); |
f1d945a1 | 208 | fHistProVeta->SetXTitle("rapidity"); |
209 | fHistProVeta->SetYTitle("v (%)"); | |
88e00a8a | 210 | fHistList->Add(fHistProVeta); |
f1d945a1 | 211 | |
882ffd6a | 212 | fHistProVPt = new TProfile("Second_FlowPro_VPt_LYZ","Second_FlowPro_VPt_LYZ",iNbinsPt,dPtMin,dPtMax); |
f1d945a1 | 213 | fHistProVPt->SetXTitle("Pt"); |
214 | fHistProVPt->SetYTitle("v (%)"); | |
88e00a8a | 215 | fHistList->Add(fHistProVPt); |
f1d945a1 | 216 | |
882ffd6a | 217 | fHistProReDtheta = new TProfile("Second_FlowPro_ReDtheta_LYZ","Second_FlowPro_ReDtheta_LYZ",iNtheta, -0.5, iNtheta-0.5); |
f1d945a1 | 218 | fHistProReDtheta->SetXTitle("#theta"); |
219 | fHistProReDtheta->SetYTitle("Re(D^{#theta})"); | |
88e00a8a | 220 | fHistList->Add(fHistProReDtheta); |
f1d945a1 | 221 | |
882ffd6a | 222 | fHistProImDtheta = new TProfile("Second_FlowPro_ImDtheta_LYZ","Second_FlowPro_ImDtheta_LYZ",iNtheta, -0.5, iNtheta-0.5); |
f1d945a1 | 223 | fHistProImDtheta->SetXTitle("#theta"); |
224 | fHistProImDtheta->SetYTitle("Im(D^{#theta})"); | |
88e00a8a | 225 | fHistList->Add(fHistProImDtheta); |
f1d945a1 | 226 | |
227 | //class AliFlowLYZHist2 defines the histograms: | |
882ffd6a | 228 | for (Int_t theta=0;theta<iNtheta;theta++) { |
9d062fe3 | 229 | TString name = "AliFlowLYZHist2_"; |
230 | name += theta; | |
231 | fHist2[theta]=new AliFlowLYZHist2(theta,name); | |
232 | fHistList->Add(fHist2[theta]); | |
f1d945a1 | 233 | } |
9d062fe3 | 234 | |
235 | //read histogram fHistProR0theta from the first run list | |
236 | if (fFirstRunList) { | |
237 | fHistProR0theta = (TProfile*)fFirstRunList->FindObject("First_FlowPro_r0theta_LYZ"); | |
88e00a8a | 238 | if (!fHistProR0theta) {cout<<"fHistProR0theta has a NULL pointer!"<<endl;} |
239 | fHistList->Add(fHistProR0theta); | |
9d062fe3 | 240 | } else { cout<<"list is NULL pointer!"<<endl; } |
241 | ||
f1d945a1 | 242 | } |
243 | ||
244 | ||
245 | if (fDebug) cout<<"****Histograms initialised****"<<endl; | |
246 | ||
247 | fEventNumber = 0; //set event counter to zero | |
248 | ||
249 | return kTRUE; | |
250 | } | |
251 | ||
252 | //----------------------------------------------------------------------- | |
253 | ||
0092f3c2 | 254 | Bool_t AliFlowAnalysisWithLeeYangZeros::Make(AliFlowEventSimple* anEvent) |
f1d945a1 | 255 | { |
256 | //make method | |
257 | if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::Make()****"<<endl; | |
258 | ||
259 | //get tracks from event | |
0092f3c2 | 260 | if (anEvent) { |
f1d945a1 | 261 | if (fFirstRun){ |
0092f3c2 | 262 | fCommonHists->FillControlHistograms(anEvent); |
263 | FillFromFlowEvent(anEvent); | |
f1d945a1 | 264 | } |
265 | else { | |
0092f3c2 | 266 | fCommonHists->FillControlHistograms(anEvent); |
267 | SecondFillFromFlowEvent(anEvent); | |
f1d945a1 | 268 | } |
269 | } | |
270 | ||
271 | else { | |
272 | cout<<"##### FlowLeeYangZero: Stack pointer null"<<endl; | |
273 | return kFALSE; | |
274 | } | |
deebd72f | 275 | // cout<<"^^^^read event "<<fEventNumber<<endl; |
f1d945a1 | 276 | fEventNumber++; |
277 | ||
278 | ||
279 | return kTRUE; | |
280 | } | |
281 | ||
282 | ||
283 | //----------------------------------------------------------------------- | |
284 | Bool_t AliFlowAnalysisWithLeeYangZeros::Finish() | |
285 | { | |
286 | //finish method | |
287 | if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::Finish()****"<<endl; | |
7b88dba1 | 288 | |
289 | cout << endl; | |
290 | cout << "******************************************************" << endl; | |
f1d945a1 | 291 | //define variables for both runs |
882ffd6a | 292 | Double_t dJ01 = 2.405; |
293 | Int_t iNtheta = AliFlowLYZConstants::kTheta; | |
9d062fe3 | 294 | //set the event number |
295 | SetEventNumber((int)fCommonHists->GetHistMultOrig()->GetEntries()); | |
1eda9b5e | 296 | //cout<<"number of events processed is "<<fEventNumber<<endl; |
9d062fe3 | 297 | |
298 | //set the sum of Q vectors | |
299 | fQsum->Set(fHistQsumforChi->GetBinContent(1),fHistQsumforChi->GetBinContent(2)); | |
300 | SetQ2sum(fHistQsumforChi->GetBinContent(3)); | |
301 | ||
f1d945a1 | 302 | if (fFirstRun){ |
882ffd6a | 303 | Double_t dR0 = 0; |
304 | Double_t dVtheta = 0; | |
305 | Double_t dv = 0; | |
306 | Double_t dV = 0; | |
307 | for (Int_t theta=0;theta<iNtheta;theta++) | |
f1d945a1 | 308 | { |
309 | //get the first minimum r0 | |
310 | fHist1[theta]->FillGtheta(); | |
882ffd6a | 311 | dR0 = fHist1[theta]->GetR0(); |
f1d945a1 | 312 | |
313 | //calculate integrated flow | |
882ffd6a | 314 | if (dR0!=0) dVtheta = dJ01/dR0; |
315 | //for estimating systematic error resulting from d0 | |
316 | Double_t dBinsize = (AliFlowLYZConstants::fgMax)/(AliFlowLYZConstants::kNbins); | |
317 | Double_t dVplus = dJ01/(dR0+dBinsize); | |
318 | Double_t dVmin = dJ01/(dR0-dBinsize); | |
319 | dv = dVtheta; | |
320 | Double_t dvplus = dVplus; | |
321 | Double_t dvmin = dVmin; | |
322 | cout<<"dv = "<<dv<<" and dvplus = "<<dvplus<< " and dvmin = "<<dvmin<<endl; | |
f1d945a1 | 323 | |
324 | //fill the histograms | |
882ffd6a | 325 | fHistProR0theta->Fill(theta,dR0); |
326 | fHistProVtheta->Fill(theta,dVtheta); | |
f1d945a1 | 327 | //get average value of fVtheta = fV |
882ffd6a | 328 | dV += dVtheta; |
f1d945a1 | 329 | } //end of loop over theta |
330 | ||
331 | //get average value of fVtheta = fV | |
882ffd6a | 332 | dV /=iNtheta; |
f1d945a1 | 333 | |
334 | //sigma2 and chi | |
882ffd6a | 335 | Double_t dSigma2 = 0; |
336 | Double_t dChi= 0; | |
f1d945a1 | 337 | if (fEventNumber!=0) { |
0092f3c2 | 338 | *fQsum /= fEventNumber; |
f1d945a1 | 339 | fQ2sum /= fEventNumber; |
882ffd6a | 340 | dSigma2 = fQ2sum - TMath::Power(fQsum->X(),2.) - TMath::Power(fQsum->Y(),2.) - TMath::Power(dV,2.); //BP eq. 62 |
341 | if (dSigma2>0) dChi = dV/TMath::Sqrt(dSigma2); | |
342 | else dChi = -1.; | |
343 | fCommonHistsRes->FillChi(dChi); | |
344 | cout<<"dV = "<<dV<<" and chi = "<<dChi<<endl; | |
f1d945a1 | 345 | } |
346 | ||
347 | // recalculate statistical errors on integrated flow | |
348 | //combining 5 theta angles to 1 relative error BP eq. 89 | |
882ffd6a | 349 | Double_t dRelErr2comb = 0.; |
9d062fe3 | 350 | Int_t iEvts = fEventNumber; |
1eda9b5e | 351 | if (iEvts!=0) { |
352 | for (Int_t theta=0;theta<iNtheta;theta++){ | |
353 | Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi(); | |
354 | Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)* | |
882ffd6a | 355 | TMath::Cos(dTheta)); |
1eda9b5e | 356 | Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)* |
882ffd6a | 357 | TMath::Cos(dTheta)); |
1eda9b5e | 358 | dRelErr2comb += (1/(2*iEvts*(dJ01*dJ01)*TMath::BesselJ1(dJ01)* |
882ffd6a | 359 | TMath::BesselJ1(dJ01)))* |
1eda9b5e | 360 | (dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2)) + |
361 | dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))); | |
362 | } | |
363 | dRelErr2comb /= iNtheta; | |
f1d945a1 | 364 | } |
882ffd6a | 365 | Double_t dRelErrcomb = TMath::Sqrt(dRelErr2comb); |
f1d945a1 | 366 | |
367 | //copy content of profile into TH1D and add error | |
882ffd6a | 368 | Double_t dv2pro = dV; //in the case that fv is equal to fV |
369 | Double_t dv2Err = dv2pro*dRelErrcomb ; | |
370 | cout<<"dv2pro +- dv2Err = "<<dv2pro<<" +- "<<dv2Err<<endl; | |
371 | fCommonHistsRes->FillIntegratedFlow(dv2pro, dv2Err); | |
f1d945a1 | 372 | |
373 | ||
374 | if (fDebug) cout<<"****histograms filled****"<<endl; | |
375 | ||
f1d945a1 | 376 | return kTRUE; |
377 | fEventNumber =0; //set to zero for second round over events | |
378 | } //firstrun | |
379 | ||
380 | ||
381 | else { //second run | |
382 | ||
383 | //calculate differential flow | |
384 | //declare variables | |
882ffd6a | 385 | TComplex cDenom, cNumer, cDtheta; |
f1d945a1 | 386 | Int_t m = 1; |
387 | TComplex i = TComplex::I(); | |
882ffd6a | 388 | Double_t dBesselRatio[3] = {1., 1.202, 2.69}; |
389 | Int_t iNbinsPt = AliFlowCommonConstants::GetNbinsPt(); | |
390 | Int_t iNbinsEta = AliFlowCommonConstants::GetNbinsEta(); | |
f1d945a1 | 391 | |
882ffd6a | 392 | Double_t dEta, dPt, dReRatio, dVeta, dVPt; |
9d062fe3 | 393 | |
882ffd6a | 394 | Double_t dR0 = 0.; |
395 | Double_t dVtheta = 0.; | |
396 | Double_t dV = 0.; | |
397 | Double_t dReDenom = 0.; | |
398 | Double_t dImDenom = 0.; | |
399 | for (Int_t theta=0;theta<iNtheta;theta++) { //loop over theta | |
f1d945a1 | 400 | if (fHistProR0theta) { |
882ffd6a | 401 | dR0 = fHistProR0theta->GetBinContent(theta+1); |
402 | if (fDebug) cerr<<"dR0 = "<<dR0<<endl; | |
403 | if (dR0!=0) dVtheta = dJ01/dR0; | |
404 | dV += dVtheta; | |
f1d945a1 | 405 | // BP Eq. 9 -> Vn^theta = j01/r0^theta |
406 | if (fHistProReDenom && fHistProImDenom) { | |
882ffd6a | 407 | dReDenom = fHistProReDenom->GetBinContent(theta+1); |
408 | dImDenom = fHistProImDenom->GetBinContent(theta+1); | |
f1d945a1 | 409 | } |
410 | else { | |
411 | cout << "Hist pointer fDenom in file does not exist" <<endl; | |
412 | } | |
413 | ||
414 | } | |
415 | else { | |
416 | cout << "Hist pointer R0theta in file does not exist" <<endl; | |
417 | } | |
418 | //} //loop over theta | |
419 | ||
882ffd6a | 420 | cDenom(dReDenom,dImDenom); |
9d062fe3 | 421 | |
f1d945a1 | 422 | //for new method and use by others (only with the sum generating function): |
423 | if (fUseSum) { | |
882ffd6a | 424 | dR0 = fHistProR0theta->GetBinContent(theta+1); |
425 | cDtheta = dR0*cDenom/dJ01; | |
426 | fHistProReDtheta->Fill(theta,cDtheta.Re()); | |
427 | fHistProImDtheta->Fill(theta,cDtheta.Im()); | |
f1d945a1 | 428 | } |
429 | ||
882ffd6a | 430 | cDenom *= TComplex::Power(i, m-1); |
f1d945a1 | 431 | //cerr<<"TComplex::Power(i, m-1) = "<<TComplex::Power(i, m-1).Rho()<<endl; //checked ok |
432 | ||
433 | //v as a function of eta | |
882ffd6a | 434 | for (Int_t be=1;be<=iNbinsEta;be++) { |
435 | dEta = fHist2[theta]->GetBinCenter(be); | |
436 | cNumer = fHist2[theta]->GetNumerEta(be); | |
437 | if (cNumer.Rho()==0) { | |
438 | if (fDebug) cerr<<"WARNING: modulus of cNumer is zero in Finish()"<<endl; | |
439 | dReRatio = 0; | |
f1d945a1 | 440 | } |
882ffd6a | 441 | else if (cDenom.Rho()==0) { |
442 | if (fDebug) cerr<<"WARNING: modulus of cDenom is zero"<<endl; | |
443 | dReRatio = 0; | |
f1d945a1 | 444 | } |
445 | else { | |
882ffd6a | 446 | //if ( j==1 && theta==0) cerr<<"modulus of cNumer = "<<cNumer.Rho() <<endl; //always a number smaller than 1, or 0. |
447 | dReRatio = (cNumer/cDenom).Re(); | |
f1d945a1 | 448 | } |
449 | ||
0828acf6 | 450 | dVeta = dBesselRatio[m-1]*dReRatio*dVtheta; //BP eq. 12 |
882ffd6a | 451 | //if ( j==1 && theta==0) cerr<<"eta = "<<dEta<<" cerr::dReRatio for eta = "<<dReRatio<<" cerr::dVeta for eta = "<<dVeta<<endl; |
f1d945a1 | 452 | |
882ffd6a | 453 | fHistProVeta->Fill(dEta,dVeta); |
f1d945a1 | 454 | } //loop over bins be |
455 | ||
456 | //v as a function of Pt | |
882ffd6a | 457 | for (Int_t bp=1;bp<=iNbinsPt;bp++) { |
458 | dPt = fHist2[theta]->GetBinCenterPt(bp); | |
459 | cNumer = fHist2[theta]->GetNumerPt(bp); | |
460 | if (cNumer.Rho()==0) { | |
461 | if (fDebug) cerr<<"modulus of cNumer is zero"<<endl; | |
462 | dReRatio = 0; | |
f1d945a1 | 463 | } |
882ffd6a | 464 | else if (cDenom.Rho()==0) { |
465 | if (fDebug) cerr<<"modulus of cDenom is zero"<<endl; | |
466 | dReRatio = 0; | |
f1d945a1 | 467 | } |
468 | else { | |
469 | //if ( j==1 && theta==0) cerr<<"modulus of fNumer = "<<fNumer.Rho() <<endl; //always a number smaller than 1, or 0. | |
882ffd6a | 470 | dReRatio = (cNumer/cDenom).Re(); |
f1d945a1 | 471 | } |
472 | ||
0828acf6 | 473 | dVPt = dBesselRatio[m-1]*dReRatio*dVtheta; //BP eq. 12 |
f1d945a1 | 474 | |
882ffd6a | 475 | fHistProVPt->Fill(dPt,dVPt); |
f1d945a1 | 476 | } //loop over bins bp |
477 | ||
478 | }//end of loop over theta | |
479 | ||
480 | //calculate the average of fVtheta = fV | |
882ffd6a | 481 | dV /= iNtheta; |
f1d945a1 | 482 | |
483 | //sigma2 and chi (for statistical error calculations) | |
882ffd6a | 484 | Double_t dSigma2 = 0; |
485 | Double_t dChi= 0; | |
f1d945a1 | 486 | if (fEventNumber!=0) { |
0092f3c2 | 487 | *fQsum /= fEventNumber; |
f1d945a1 | 488 | //cerr<<"fQsum.X() = "<<fQsum.X()<<endl; |
489 | //cerr<<"fQsum.Y() = "<<fQsum.Y()<<endl; | |
490 | fQ2sum /= fEventNumber; | |
491 | //cout<<"fQ2sum = "<<fQ2sum<<endl; | |
882ffd6a | 492 | dSigma2 = fQ2sum - TMath::Power(fQsum->X(),2.) - TMath::Power(fQsum->Y(),2.) - TMath::Power(dV,2.); //BP eq. 62 |
493 | if (dSigma2>0) dChi = dV/TMath::Sqrt(dSigma2); | |
494 | else dChi = -1.; | |
495 | fCommonHistsRes->FillChi(dChi); | |
496 | cout<<"dV = "<<dV<<" and chi = "<<dChi<<endl; | |
f1d945a1 | 497 | } |
498 | ||
499 | //copy content of profile into TH1D and add error | |
882ffd6a | 500 | for(Int_t b=0;b<iNbinsPt;b++) { |
501 | Double_t dv2pro = fHistProVPt->GetBinContent(b); | |
502 | Double_t dNprime = fCommonHists->GetEntriesInPtBin(b); | |
503 | Double_t dErrdifcomb = 0.; //set error to zero | |
504 | Double_t dErr2difcomb = 0.; //set error to zero | |
f1d945a1 | 505 | //calculate error |
882ffd6a | 506 | if (dNprime!=0.) { |
507 | for (Int_t theta=0;theta<iNtheta;theta++) { | |
508 | Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi(); | |
509 | Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)* | |
510 | TMath::Cos(dTheta)); | |
511 | Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)* | |
512 | TMath::Cos(dTheta)); | |
513 | dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)* | |
514 | TMath::BesselJ1(dJ01)))* | |
515 | ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) - | |
516 | (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2)))); | |
f1d945a1 | 517 | } //loop over theta |
518 | } | |
519 | ||
882ffd6a | 520 | if (dErr2difcomb!=0.) { |
521 | dErr2difcomb /= iNtheta; | |
0828acf6 | 522 | dErrdifcomb = TMath::Sqrt(dErr2difcomb); |
882ffd6a | 523 | //cerr<<"dErrdifcomb = "<<dErrdifcomb<<endl; |
f1d945a1 | 524 | } |
882ffd6a | 525 | else {dErrdifcomb = 0.;} |
f1d945a1 | 526 | |
527 | //fill TH1D | |
882ffd6a | 528 | fCommonHistsRes->FillDifferentialFlow(b, dv2pro, dErrdifcomb); |
f1d945a1 | 529 | } //loop over bins b |
9d062fe3 | 530 | |
f1d945a1 | 531 | } //secondrun |
532 | ||
533 | cout<<"----LYZ analysis finished....----"<<endl<<endl; | |
534 | ||
535 | return kTRUE; | |
536 | } | |
537 | ||
538 | ||
539 | //----------------------------------------------------------------------- | |
540 | ||
0092f3c2 | 541 | Bool_t AliFlowAnalysisWithLeeYangZeros::FillFromFlowEvent(AliFlowEventSimple* anEvent) |
f1d945a1 | 542 | { |
543 | // Get event quantities from AliFlowEvent for all particles | |
544 | ||
545 | if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::FillFromFlowEvent()****"<<endl; | |
546 | ||
0092f3c2 | 547 | if (!anEvent){ |
f1d945a1 | 548 | cout<<"##### FlowLeeYangZero: FlowEvent pointer null"<<endl; |
549 | return kFALSE; | |
550 | } | |
551 | ||
552 | //define variables | |
882ffd6a | 553 | TComplex cExpo, cGtheta, cGthetaNew, cZ; |
554 | Int_t iNtheta = AliFlowLYZConstants::kTheta; | |
555 | Int_t iNbins = AliFlowLYZConstants::kNbins; | |
f1d945a1 | 556 | |
557 | ||
558 | //calculate flow | |
882ffd6a | 559 | Double_t dOrder = 2.; |
f1d945a1 | 560 | |
561 | //get the Q vector | |
882ffd6a | 562 | AliFlowVector vQ = anEvent->GetQ(); |
f1d945a1 | 563 | //weight by the multiplicity |
882ffd6a | 564 | Double_t dQX = 0; |
565 | Double_t dQY = 0; | |
566 | if (vQ.GetMult() != 0) { | |
567 | dQX = vQ.X()/vQ.GetMult(); | |
568 | dQY = vQ.Y()/vQ.GetMult(); | |
f2ec07bf | 569 | } |
882ffd6a | 570 | vQ.Set(dQX,dQY); |
9d062fe3 | 571 | |
f1d945a1 | 572 | //for chi calculation: |
882ffd6a | 573 | *fQsum += vQ; |
9d062fe3 | 574 | fHistQsumforChi->SetBinContent(1,fQsum->X()); |
575 | fHistQsumforChi->SetBinContent(2,fQsum->Y()); | |
882ffd6a | 576 | fQ2sum += vQ.Mod2(); |
9d062fe3 | 577 | fHistQsumforChi->SetBinContent(3,fQ2sum); |
f1d945a1 | 578 | //cerr<<"fQ2sum = "<<fQ2sum<<endl; |
579 | ||
882ffd6a | 580 | for (Int_t theta=0;theta<iNtheta;theta++) |
f1d945a1 | 581 | { |
882ffd6a | 582 | Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi()/dOrder; |
f1d945a1 | 583 | |
882ffd6a | 584 | //calculate dQtheta = cos(dOrder*(fPhi-dTheta);the projection of the Q vector on the reference direction dTheta |
585 | Double_t dQtheta = GetQtheta(vQ, dTheta); | |
f1d945a1 | 586 | |
882ffd6a | 587 | for (Int_t bin=1;bin<=iNbins;bin++) |
f1d945a1 | 588 | { |
882ffd6a | 589 | Double_t dR = fHist1[theta]->GetBinCenter(bin); //bincentre of bins in histogram //FIXED??? |
590 | //if (theta == 0) cerr<<"cerr::dR = "<<dR<<endl; | |
f1d945a1 | 591 | if (fUseSum) |
592 | { | |
593 | //calculate the sum generating function | |
882ffd6a | 594 | cExpo(0.,dR*dQtheta); //Re=0 ; Im=dR*dQtheta |
595 | cGtheta = TComplex::Exp(cExpo); | |
f1d945a1 | 596 | } |
597 | else | |
598 | { | |
599 | //calculate the product generating function | |
882ffd6a | 600 | cGtheta = GetGrtheta(anEvent, dR, dTheta); //make this function |
601 | if (cGtheta.Rho2() > 100.) break; | |
f1d945a1 | 602 | } |
603 | ||
882ffd6a | 604 | fHist1[theta]->Fill(dR,cGtheta); //fill real and imaginary part of cGtheta |
f1d945a1 | 605 | |
606 | } //loop over bins | |
607 | } //loop over theta | |
608 | ||
609 | ||
610 | return kTRUE; | |
611 | ||
612 | ||
613 | } | |
614 | ||
615 | //----------------------------------------------------------------------- | |
0092f3c2 | 616 | Bool_t AliFlowAnalysisWithLeeYangZeros::SecondFillFromFlowEvent(AliFlowEventSimple* anEvent) |
f1d945a1 | 617 | { |
618 | //for differential flow | |
619 | ||
620 | if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::SecondFillFromFlowEvent()****"<<endl; | |
621 | ||
0092f3c2 | 622 | if (!anEvent){ |
f1d945a1 | 623 | cout<<"##### FlowLeeYangZero: FlowEvent pointer null"<<endl; |
624 | return kFALSE; | |
625 | } | |
626 | ||
627 | //define variables | |
882ffd6a | 628 | TComplex cExpo, cDenom, cNumer,cCosTermComplex; |
629 | Double_t dPhi, dEta, dPt; | |
630 | Double_t dR0 = 0.; | |
631 | Double_t dCosTerm; | |
f1d945a1 | 632 | Double_t m = 1.; |
882ffd6a | 633 | Double_t dOrder = 2.; |
634 | Int_t iNtheta = AliFlowLYZConstants::kTheta; | |
f1d945a1 | 635 | |
882ffd6a | 636 | |
f1d945a1 | 637 | //get the Q vector |
882ffd6a | 638 | AliFlowVector vQ = anEvent->GetQ(); |
f1d945a1 | 639 | //weight by the multiplicity |
882ffd6a | 640 | Double_t dQX = 0.; |
641 | Double_t dQY = 0.; | |
642 | if (vQ.GetMult() != 0) { | |
643 | dQX = vQ.X()/vQ.GetMult(); | |
644 | dQY = vQ.Y()/vQ.GetMult(); | |
f2ec07bf | 645 | } |
9d062fe3 | 646 | vQ.Set(dQX,dQY); |
647 | ||
f1d945a1 | 648 | //for chi calculation: |
882ffd6a | 649 | *fQsum += vQ; |
9d062fe3 | 650 | fHistQsumforChi->SetBinContent(1,fQsum->X()); |
651 | fHistQsumforChi->SetBinContent(2,fQsum->Y()); | |
882ffd6a | 652 | fQ2sum += vQ.Mod2(); |
9d062fe3 | 653 | fHistQsumforChi->SetBinContent(3,fQ2sum); |
f1d945a1 | 654 | |
882ffd6a | 655 | for (Int_t theta=0;theta<iNtheta;theta++) |
f1d945a1 | 656 | { |
882ffd6a | 657 | Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi()/dOrder; |
f1d945a1 | 658 | |
882ffd6a | 659 | //calculate dQtheta = cos(dOrder*(dPhi-dTheta);the projection of the Q vector on the reference direction dTheta |
660 | Double_t dQtheta = GetQtheta(vQ, dTheta); | |
661 | //cerr<<"dQtheta for fdenom = "<<dQtheta<<endl; | |
f1d945a1 | 662 | |
663 | //denominator for differential v | |
664 | if (fHistProR0theta) { | |
882ffd6a | 665 | dR0 = fHistProR0theta->GetBinContent(theta+1); |
f1d945a1 | 666 | } |
667 | else { | |
88e00a8a | 668 | cout <<"pointer fHistProR0theta does not exist" << endl; |
f1d945a1 | 669 | } |
882ffd6a | 670 | //cerr<<"dR0 = "<<dR0 <<endl; |
f1d945a1 | 671 | |
672 | if (fUseSum) //sum generating function | |
673 | { | |
882ffd6a | 674 | cExpo(0.,dR0*dQtheta); |
675 | cDenom = dQtheta*(TComplex::Exp(cExpo)); //BP eq 12 | |
f1d945a1 | 676 | //loop over tracks in event |
882ffd6a | 677 | Int_t iNumberOfTracks = anEvent->NumberOfTracks(); |
678 | for (Int_t i=0;i<iNumberOfTracks;i++) { | |
679 | AliFlowTrackSimple* pTrack = anEvent->GetTrack(i); | |
680 | if (pTrack) { | |
681 | if (pTrack->UseForDifferentialFlow()) { | |
682 | dEta = pTrack->Eta(); | |
683 | dPt = pTrack->Pt(); | |
684 | dPhi = pTrack->Phi(); | |
685 | dCosTerm = cos(m*dOrder*(dPhi-dTheta)); | |
686 | //cerr<<"dCosTerm = "<<dCosTerm <<endl; | |
687 | cNumer = dCosTerm*(TComplex::Exp(cExpo)); | |
688 | if (cNumer.Rho()==0) {cerr<<"WARNING: modulus of cNumer is zero in SecondFillFromFlowEvent"<<endl;} | |
689 | if (fDebug) cerr<<"modulus of cNumer is "<<cNumer.Rho()<<endl; | |
f1d945a1 | 690 | if (fHist2[theta]) { |
882ffd6a | 691 | fHist2[theta]->Fill(dEta,dPt,cNumer); |
f1d945a1 | 692 | } |
693 | else { | |
694 | cout << "fHist2 pointer mising" <<endl; | |
695 | } | |
696 | } | |
697 | } //if particle | |
698 | else {cerr << "no particle!!!"<<endl;} | |
699 | } //loop over tracks | |
700 | ||
701 | } //sum | |
702 | else //product generating function | |
703 | { | |
882ffd6a | 704 | cDenom = GetDiffFlow(anEvent, dR0, theta); |
f1d945a1 | 705 | |
706 | }//product | |
707 | if (fHistProReDenom && fHistProImDenom) { | |
882ffd6a | 708 | fHistProReDenom->Fill(theta,cDenom.Re()); //fill the real part of fDenom |
709 | fHistProImDenom->Fill(theta,cDenom.Im()); //fill the imaginary part of fDenom | |
f1d945a1 | 710 | } |
711 | else { | |
882ffd6a | 712 | cout << "Pointers to cDenom mising" << endl; |
f1d945a1 | 713 | } |
714 | ||
715 | }//end of loop over theta | |
716 | ||
717 | return kTRUE; | |
718 | ||
719 | ||
720 | } | |
721 | //----------------------------------------------------------------------- | |
882ffd6a | 722 | Double_t AliFlowAnalysisWithLeeYangZeros::GetQtheta(AliFlowVector aQ, Double_t aTheta) |
f1d945a1 | 723 | { |
882ffd6a | 724 | //calculate Qtheta. Qtheta is the sum over all particles of cos(dOrder*(dPhi-dTheta)) BP eq. 3 |
f1d945a1 | 725 | if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::GetQtheta()****"<<endl; |
726 | ||
882ffd6a | 727 | Double_t dQtheta = 0.; |
728 | Double_t dOrder = 2.; | |
f1d945a1 | 729 | |
882ffd6a | 730 | dQtheta = aQ.X()*cos(dOrder*aTheta)+aQ.Y()*sin(dOrder*aTheta); |
f1d945a1 | 731 | |
882ffd6a | 732 | return dQtheta; |
f1d945a1 | 733 | |
734 | } | |
735 | ||
736 | ||
737 | //----------------------------------------------------------------------- | |
882ffd6a | 738 | TComplex AliFlowAnalysisWithLeeYangZeros::GetGrtheta(AliFlowEventSimple* anEvent, Double_t aR, Double_t aTheta) |
f1d945a1 | 739 | { |
740 | // Product Generating Function for LeeYangZeros method | |
741 | // PG Eq. 3 (J. Phys. G Nucl. Part. Phys 30 S1213 (2004)) | |
742 | ||
743 | if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::GetGrtheta()****"<<endl; | |
744 | ||
745 | ||
882ffd6a | 746 | TComplex cG = TComplex::One(); |
747 | Double_t dOrder = 2.; | |
748 | Double_t dWgt = 1.; | |
f1d945a1 | 749 | |
882ffd6a | 750 | Int_t iNumberOfTracks = anEvent->NumberOfTracks(); |
f1d945a1 | 751 | |
882ffd6a | 752 | for (Int_t i=0;i<iNumberOfTracks;i++) //loop over tracks in event |
f1d945a1 | 753 | { |
882ffd6a | 754 | AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ; |
755 | if (pTrack){ | |
756 | if (pTrack->UseForIntegratedFlow()) { | |
757 | Double_t dPhi = pTrack->Phi(); | |
758 | Double_t dGIm = aR * dWgt*cos(dOrder*(dPhi - aTheta)); | |
759 | TComplex cGi(1., dGIm); | |
760 | cG *= cGi; //product over all tracks | |
f1d945a1 | 761 | } |
762 | } | |
763 | else {cerr << "no particle pointer !!!"<<endl;} | |
764 | }//loop over tracks | |
765 | ||
882ffd6a | 766 | return cG; |
f1d945a1 | 767 | |
768 | } | |
769 | ||
770 | ||
771 | //----------------------------------------------------------------------- | |
882ffd6a | 772 | TComplex AliFlowAnalysisWithLeeYangZeros::GetDiffFlow(AliFlowEventSimple* anEvent, Double_t aR0, Int_t theta) |
f1d945a1 | 773 | { |
774 | // Sum for the denominator for diff. flow for the Product Generating Function for LeeYangZeros method | |
775 | // PG Eq. 9 (J. Phys. G Nucl. Part. Phys 30 S1213 (2004)) | |
776 | // Also for v1 mixed harmonics: DF Eq. 5 | |
777 | // It is the deriverative of Grtheta at r0 divided by Grtheta at r0 | |
778 | ||
779 | if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::GetGrtheta()****"<<endl; | |
780 | ||
882ffd6a | 781 | TComplex cG = TComplex::One(); |
782 | TComplex cdGr0(0.,0.); | |
783 | Double_t dOrder = 2.; | |
784 | Double_t dWgt = 1.; | |
f1d945a1 | 785 | |
882ffd6a | 786 | Int_t iNumberOfTracks = anEvent->NumberOfTracks(); |
f1d945a1 | 787 | |
882ffd6a | 788 | Int_t iNtheta = AliFlowLYZConstants::kTheta; |
789 | Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi()/dOrder; | |
f1d945a1 | 790 | |
882ffd6a | 791 | for (Int_t i=0;i<iNumberOfTracks;i++) //loop over tracks in event |
f1d945a1 | 792 | { |
882ffd6a | 793 | AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ; |
794 | if (pTrack){ | |
795 | if (pTrack->UseForDifferentialFlow()) { | |
796 | Double_t dPhi = pTrack->Phi(); | |
797 | Double_t dCosTerm = dWgt*cos(dOrder*(dPhi - dTheta)); | |
f1d945a1 | 798 | //GetGr0theta |
882ffd6a | 799 | Double_t dGIm = aR0 * dCosTerm; |
800 | TComplex cGi(1., dGIm); | |
801 | cG *= cGi; //product over all tracks | |
f1d945a1 | 802 | //GetdGr0theta |
882ffd6a | 803 | TComplex cCosTermComplex(1., aR0*dCosTerm); |
804 | cdGr0 +=(dCosTerm / cCosTermComplex); //sum over all tracks | |
f1d945a1 | 805 | } |
806 | } //if particle | |
807 | else {cerr << "no particle!!!"<<endl;} | |
808 | }//loop over tracks | |
809 | ||
810 | ||
882ffd6a | 811 | for (Int_t i=0;i<iNumberOfTracks;i++) |
f1d945a1 | 812 | { |
882ffd6a | 813 | AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ; |
814 | if (pTrack){ | |
815 | if (pTrack->UseForDifferentialFlow()) { | |
816 | Double_t dEta = pTrack->Eta(); | |
817 | Double_t dPt = pTrack->Pt(); | |
818 | Double_t dPhi = pTrack->Phi(); | |
819 | Double_t dCosTerm = cos(dOrder*(dPhi-dTheta)); | |
820 | TComplex cCosTermComplex(1.,aR0*dCosTerm); | |
821 | TComplex cNumer = cG*dCosTerm/cCosTermComplex; //PG Eq. 9 | |
822 | fHist2[theta]->Fill(dEta,dPt,cNumer); | |
f1d945a1 | 823 | } |
824 | } //if particle | |
825 | else {cerr << "no particle pointer!!!"<<endl;} | |
826 | }//loop over tracks | |
827 | ||
882ffd6a | 828 | TComplex cDenom = cG*cdGr0; |
829 | return cDenom; | |
f1d945a1 | 830 | |
831 | } | |
832 | ||
833 | //----------------------------------------------------------------------- | |
834 |