cleanup of some cout's
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowAnalysisWithLeeYangZeros.cxx
CommitLineData
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
29class AliFlowVector;
30
31#include "TMath.h" //needed as include
32#include "TFile.h" //needed as include
88e00a8a 33#include "TList.h"
f1d945a1 34
35class TComplex;
36class TProfile;
37class TH1F;
38class 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
52ClassImp(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
110void 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 126void 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 142Bool_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 254Bool_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 738TComplex 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 772TComplex 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