overloaded write method
[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 }
275 cout<<"^^^^read event "<<fEventNumber<<endl;
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;
288
289 //define variables for both runs
882ffd6a 290 Double_t dJ01 = 2.405;
291 Int_t iNtheta = AliFlowLYZConstants::kTheta;
9d062fe3 292 //set the event number
293 SetEventNumber((int)fCommonHists->GetHistMultOrig()->GetEntries());
1eda9b5e 294 //cout<<"number of events processed is "<<fEventNumber<<endl;
9d062fe3 295
296 //set the sum of Q vectors
297 fQsum->Set(fHistQsumforChi->GetBinContent(1),fHistQsumforChi->GetBinContent(2));
298 SetQ2sum(fHistQsumforChi->GetBinContent(3));
299
f1d945a1 300 if (fFirstRun){
882ffd6a 301 Double_t dR0 = 0;
302 Double_t dVtheta = 0;
303 Double_t dv = 0;
304 Double_t dV = 0;
305 for (Int_t theta=0;theta<iNtheta;theta++)
f1d945a1 306 {
307 //get the first minimum r0
308 fHist1[theta]->FillGtheta();
882ffd6a 309 dR0 = fHist1[theta]->GetR0();
f1d945a1 310
311 //calculate integrated flow
882ffd6a 312 if (dR0!=0) dVtheta = dJ01/dR0;
313 //for estimating systematic error resulting from d0
314 Double_t dBinsize = (AliFlowLYZConstants::fgMax)/(AliFlowLYZConstants::kNbins);
315 Double_t dVplus = dJ01/(dR0+dBinsize);
316 Double_t dVmin = dJ01/(dR0-dBinsize);
317 dv = dVtheta;
318 Double_t dvplus = dVplus;
319 Double_t dvmin = dVmin;
320 cout<<"dv = "<<dv<<" and dvplus = "<<dvplus<< " and dvmin = "<<dvmin<<endl;
f1d945a1 321
322 //fill the histograms
882ffd6a 323 fHistProR0theta->Fill(theta,dR0);
324 fHistProVtheta->Fill(theta,dVtheta);
f1d945a1 325 //get average value of fVtheta = fV
882ffd6a 326 dV += dVtheta;
f1d945a1 327 } //end of loop over theta
328
329 //get average value of fVtheta = fV
882ffd6a 330 dV /=iNtheta;
f1d945a1 331
332 //sigma2 and chi
882ffd6a 333 Double_t dSigma2 = 0;
334 Double_t dChi= 0;
f1d945a1 335 if (fEventNumber!=0) {
0092f3c2 336 *fQsum /= fEventNumber;
f1d945a1 337 fQ2sum /= fEventNumber;
882ffd6a 338 dSigma2 = fQ2sum - TMath::Power(fQsum->X(),2.) - TMath::Power(fQsum->Y(),2.) - TMath::Power(dV,2.); //BP eq. 62
339 if (dSigma2>0) dChi = dV/TMath::Sqrt(dSigma2);
340 else dChi = -1.;
341 fCommonHistsRes->FillChi(dChi);
342 cout<<"dV = "<<dV<<" and chi = "<<dChi<<endl;
f1d945a1 343 }
344
345 // recalculate statistical errors on integrated flow
346 //combining 5 theta angles to 1 relative error BP eq. 89
882ffd6a 347 Double_t dRelErr2comb = 0.;
9d062fe3 348 Int_t iEvts = fEventNumber;
1eda9b5e 349 if (iEvts!=0) {
350 for (Int_t theta=0;theta<iNtheta;theta++){
351 Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
352 Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
882ffd6a 353 TMath::Cos(dTheta));
1eda9b5e 354 Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
882ffd6a 355 TMath::Cos(dTheta));
1eda9b5e 356 dRelErr2comb += (1/(2*iEvts*(dJ01*dJ01)*TMath::BesselJ1(dJ01)*
882ffd6a 357 TMath::BesselJ1(dJ01)))*
1eda9b5e 358 (dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2)) +
359 dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2)));
360 }
361 dRelErr2comb /= iNtheta;
f1d945a1 362 }
882ffd6a 363 Double_t dRelErrcomb = TMath::Sqrt(dRelErr2comb);
f1d945a1 364
365 //copy content of profile into TH1D and add error
882ffd6a 366 Double_t dv2pro = dV; //in the case that fv is equal to fV
367 Double_t dv2Err = dv2pro*dRelErrcomb ;
368 cout<<"dv2pro +- dv2Err = "<<dv2pro<<" +- "<<dv2Err<<endl;
369 fCommonHistsRes->FillIntegratedFlow(dv2pro, dv2Err);
f1d945a1 370
371
372 if (fDebug) cout<<"****histograms filled****"<<endl;
373
f1d945a1 374 return kTRUE;
375 fEventNumber =0; //set to zero for second round over events
376 } //firstrun
377
378
379 else { //second run
380
381 //calculate differential flow
382 //declare variables
882ffd6a 383 TComplex cDenom, cNumer, cDtheta;
f1d945a1 384 Int_t m = 1;
385 TComplex i = TComplex::I();
882ffd6a 386 Double_t dBesselRatio[3] = {1., 1.202, 2.69};
387 Int_t iNbinsPt = AliFlowCommonConstants::GetNbinsPt();
388 Int_t iNbinsEta = AliFlowCommonConstants::GetNbinsEta();
f1d945a1 389
882ffd6a 390 Double_t dEta, dPt, dReRatio, dVeta, dVPt;
9d062fe3 391
882ffd6a 392 Double_t dR0 = 0.;
393 Double_t dVtheta = 0.;
394 Double_t dV = 0.;
395 Double_t dReDenom = 0.;
396 Double_t dImDenom = 0.;
397 for (Int_t theta=0;theta<iNtheta;theta++) { //loop over theta
f1d945a1 398 if (fHistProR0theta) {
882ffd6a 399 dR0 = fHistProR0theta->GetBinContent(theta+1);
400 if (fDebug) cerr<<"dR0 = "<<dR0<<endl;
401 if (dR0!=0) dVtheta = dJ01/dR0;
402 dV += dVtheta;
f1d945a1 403 // BP Eq. 9 -> Vn^theta = j01/r0^theta
404 if (fHistProReDenom && fHistProImDenom) {
882ffd6a 405 dReDenom = fHistProReDenom->GetBinContent(theta+1);
406 dImDenom = fHistProImDenom->GetBinContent(theta+1);
f1d945a1 407 }
408 else {
409 cout << "Hist pointer fDenom in file does not exist" <<endl;
410 }
411
412 }
413 else {
414 cout << "Hist pointer R0theta in file does not exist" <<endl;
415 }
416 //} //loop over theta
417
882ffd6a 418 cDenom(dReDenom,dImDenom);
9d062fe3 419
f1d945a1 420 //for new method and use by others (only with the sum generating function):
421 if (fUseSum) {
882ffd6a 422 dR0 = fHistProR0theta->GetBinContent(theta+1);
423 cDtheta = dR0*cDenom/dJ01;
424 fHistProReDtheta->Fill(theta,cDtheta.Re());
425 fHistProImDtheta->Fill(theta,cDtheta.Im());
f1d945a1 426 }
427
882ffd6a 428 cDenom *= TComplex::Power(i, m-1);
f1d945a1 429 //cerr<<"TComplex::Power(i, m-1) = "<<TComplex::Power(i, m-1).Rho()<<endl; //checked ok
430
431 //v as a function of eta
882ffd6a 432 for (Int_t be=1;be<=iNbinsEta;be++) {
433 dEta = fHist2[theta]->GetBinCenter(be);
434 cNumer = fHist2[theta]->GetNumerEta(be);
435 if (cNumer.Rho()==0) {
436 if (fDebug) cerr<<"WARNING: modulus of cNumer is zero in Finish()"<<endl;
437 dReRatio = 0;
f1d945a1 438 }
882ffd6a 439 else if (cDenom.Rho()==0) {
440 if (fDebug) cerr<<"WARNING: modulus of cDenom is zero"<<endl;
441 dReRatio = 0;
f1d945a1 442 }
443 else {
882ffd6a 444 //if ( j==1 && theta==0) cerr<<"modulus of cNumer = "<<cNumer.Rho() <<endl; //always a number smaller than 1, or 0.
445 dReRatio = (cNumer/cDenom).Re();
f1d945a1 446 }
447
0828acf6 448 dVeta = dBesselRatio[m-1]*dReRatio*dVtheta; //BP eq. 12
882ffd6a 449 //if ( j==1 && theta==0) cerr<<"eta = "<<dEta<<" cerr::dReRatio for eta = "<<dReRatio<<" cerr::dVeta for eta = "<<dVeta<<endl;
f1d945a1 450
882ffd6a 451 fHistProVeta->Fill(dEta,dVeta);
f1d945a1 452 } //loop over bins be
453
454 //v as a function of Pt
882ffd6a 455 for (Int_t bp=1;bp<=iNbinsPt;bp++) {
456 dPt = fHist2[theta]->GetBinCenterPt(bp);
457 cNumer = fHist2[theta]->GetNumerPt(bp);
458 if (cNumer.Rho()==0) {
459 if (fDebug) cerr<<"modulus of cNumer is zero"<<endl;
460 dReRatio = 0;
f1d945a1 461 }
882ffd6a 462 else if (cDenom.Rho()==0) {
463 if (fDebug) cerr<<"modulus of cDenom is zero"<<endl;
464 dReRatio = 0;
f1d945a1 465 }
466 else {
467 //if ( j==1 && theta==0) cerr<<"modulus of fNumer = "<<fNumer.Rho() <<endl; //always a number smaller than 1, or 0.
882ffd6a 468 dReRatio = (cNumer/cDenom).Re();
f1d945a1 469 }
470
0828acf6 471 dVPt = dBesselRatio[m-1]*dReRatio*dVtheta; //BP eq. 12
f1d945a1 472
882ffd6a 473 fHistProVPt->Fill(dPt,dVPt);
f1d945a1 474 } //loop over bins bp
475
476 }//end of loop over theta
477
478 //calculate the average of fVtheta = fV
882ffd6a 479 dV /= iNtheta;
f1d945a1 480
481 //sigma2 and chi (for statistical error calculations)
882ffd6a 482 Double_t dSigma2 = 0;
483 Double_t dChi= 0;
f1d945a1 484 if (fEventNumber!=0) {
0092f3c2 485 *fQsum /= fEventNumber;
f1d945a1 486 //cerr<<"fQsum.X() = "<<fQsum.X()<<endl;
487 //cerr<<"fQsum.Y() = "<<fQsum.Y()<<endl;
488 fQ2sum /= fEventNumber;
489 //cout<<"fQ2sum = "<<fQ2sum<<endl;
882ffd6a 490 dSigma2 = fQ2sum - TMath::Power(fQsum->X(),2.) - TMath::Power(fQsum->Y(),2.) - TMath::Power(dV,2.); //BP eq. 62
491 if (dSigma2>0) dChi = dV/TMath::Sqrt(dSigma2);
492 else dChi = -1.;
493 fCommonHistsRes->FillChi(dChi);
494 cout<<"dV = "<<dV<<" and chi = "<<dChi<<endl;
f1d945a1 495 }
496
497 //copy content of profile into TH1D and add error
882ffd6a 498 for(Int_t b=0;b<iNbinsPt;b++) {
499 Double_t dv2pro = fHistProVPt->GetBinContent(b);
500 Double_t dNprime = fCommonHists->GetEntriesInPtBin(b);
501 Double_t dErrdifcomb = 0.; //set error to zero
502 Double_t dErr2difcomb = 0.; //set error to zero
f1d945a1 503 //calculate error
882ffd6a 504 if (dNprime!=0.) {
505 for (Int_t theta=0;theta<iNtheta;theta++) {
506 Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
507 Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
508 TMath::Cos(dTheta));
509 Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
510 TMath::Cos(dTheta));
511 dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)*
512 TMath::BesselJ1(dJ01)))*
513 ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
514 (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
f1d945a1 515 } //loop over theta
516 }
517
882ffd6a 518 if (dErr2difcomb!=0.) {
519 dErr2difcomb /= iNtheta;
0828acf6 520 dErrdifcomb = TMath::Sqrt(dErr2difcomb);
882ffd6a 521 //cerr<<"dErrdifcomb = "<<dErrdifcomb<<endl;
f1d945a1 522 }
882ffd6a 523 else {dErrdifcomb = 0.;}
f1d945a1 524
525 //fill TH1D
882ffd6a 526 fCommonHistsRes->FillDifferentialFlow(b, dv2pro, dErrdifcomb);
f1d945a1 527 } //loop over bins b
9d062fe3 528
f1d945a1 529 } //secondrun
530
531 cout<<"----LYZ analysis finished....----"<<endl<<endl;
532
533 return kTRUE;
534}
535
536
537//-----------------------------------------------------------------------
538
0092f3c2 539 Bool_t AliFlowAnalysisWithLeeYangZeros::FillFromFlowEvent(AliFlowEventSimple* anEvent)
f1d945a1 540{
541 // Get event quantities from AliFlowEvent for all particles
542
543 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::FillFromFlowEvent()****"<<endl;
544
0092f3c2 545 if (!anEvent){
f1d945a1 546 cout<<"##### FlowLeeYangZero: FlowEvent pointer null"<<endl;
547 return kFALSE;
548 }
549
550 //define variables
882ffd6a 551 TComplex cExpo, cGtheta, cGthetaNew, cZ;
552 Int_t iNtheta = AliFlowLYZConstants::kTheta;
553 Int_t iNbins = AliFlowLYZConstants::kNbins;
f1d945a1 554
555
556 //calculate flow
882ffd6a 557 Double_t dOrder = 2.;
f1d945a1 558
559 //get the Q vector
882ffd6a 560 AliFlowVector vQ = anEvent->GetQ();
f1d945a1 561 //weight by the multiplicity
882ffd6a 562 Double_t dQX = 0;
563 Double_t dQY = 0;
564 if (vQ.GetMult() != 0) {
565 dQX = vQ.X()/vQ.GetMult();
566 dQY = vQ.Y()/vQ.GetMult();
f2ec07bf 567 }
882ffd6a 568 vQ.Set(dQX,dQY);
9d062fe3 569
f1d945a1 570 //for chi calculation:
882ffd6a 571 *fQsum += vQ;
9d062fe3 572 fHistQsumforChi->SetBinContent(1,fQsum->X());
573 fHistQsumforChi->SetBinContent(2,fQsum->Y());
882ffd6a 574 fQ2sum += vQ.Mod2();
9d062fe3 575 fHistQsumforChi->SetBinContent(3,fQ2sum);
f1d945a1 576 //cerr<<"fQ2sum = "<<fQ2sum<<endl;
577
882ffd6a 578 for (Int_t theta=0;theta<iNtheta;theta++)
f1d945a1 579 {
882ffd6a 580 Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi()/dOrder;
f1d945a1 581
882ffd6a 582 //calculate dQtheta = cos(dOrder*(fPhi-dTheta);the projection of the Q vector on the reference direction dTheta
583 Double_t dQtheta = GetQtheta(vQ, dTheta);
f1d945a1 584
882ffd6a 585 for (Int_t bin=1;bin<=iNbins;bin++)
f1d945a1 586 {
882ffd6a 587 Double_t dR = fHist1[theta]->GetBinCenter(bin); //bincentre of bins in histogram //FIXED???
588 //if (theta == 0) cerr<<"cerr::dR = "<<dR<<endl;
f1d945a1 589 if (fUseSum)
590 {
591 //calculate the sum generating function
882ffd6a 592 cExpo(0.,dR*dQtheta); //Re=0 ; Im=dR*dQtheta
593 cGtheta = TComplex::Exp(cExpo);
f1d945a1 594 }
595 else
596 {
597 //calculate the product generating function
882ffd6a 598 cGtheta = GetGrtheta(anEvent, dR, dTheta); //make this function
599 if (cGtheta.Rho2() > 100.) break;
f1d945a1 600 }
601
882ffd6a 602 fHist1[theta]->Fill(dR,cGtheta); //fill real and imaginary part of cGtheta
f1d945a1 603
604 } //loop over bins
605 } //loop over theta
606
607
608 return kTRUE;
609
610
611}
612
613 //-----------------------------------------------------------------------
0092f3c2 614 Bool_t AliFlowAnalysisWithLeeYangZeros::SecondFillFromFlowEvent(AliFlowEventSimple* anEvent)
f1d945a1 615{
616 //for differential flow
617
618 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::SecondFillFromFlowEvent()****"<<endl;
619
0092f3c2 620 if (!anEvent){
f1d945a1 621 cout<<"##### FlowLeeYangZero: FlowEvent pointer null"<<endl;
622 return kFALSE;
623 }
624
625 //define variables
882ffd6a 626 TComplex cExpo, cDenom, cNumer,cCosTermComplex;
627 Double_t dPhi, dEta, dPt;
628 Double_t dR0 = 0.;
629 Double_t dCosTerm;
f1d945a1 630 Double_t m = 1.;
882ffd6a 631 Double_t dOrder = 2.;
632 Int_t iNtheta = AliFlowLYZConstants::kTheta;
f1d945a1 633
882ffd6a 634
f1d945a1 635 //get the Q vector
882ffd6a 636 AliFlowVector vQ = anEvent->GetQ();
f1d945a1 637 //weight by the multiplicity
882ffd6a 638 Double_t dQX = 0.;
639 Double_t dQY = 0.;
640 if (vQ.GetMult() != 0) {
641 dQX = vQ.X()/vQ.GetMult();
642 dQY = vQ.Y()/vQ.GetMult();
f2ec07bf 643 }
9d062fe3 644 vQ.Set(dQX,dQY);
645
f1d945a1 646 //for chi calculation:
882ffd6a 647 *fQsum += vQ;
9d062fe3 648 fHistQsumforChi->SetBinContent(1,fQsum->X());
649 fHistQsumforChi->SetBinContent(2,fQsum->Y());
882ffd6a 650 fQ2sum += vQ.Mod2();
9d062fe3 651 fHistQsumforChi->SetBinContent(3,fQ2sum);
f1d945a1 652
882ffd6a 653 for (Int_t theta=0;theta<iNtheta;theta++)
f1d945a1 654 {
882ffd6a 655 Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi()/dOrder;
f1d945a1 656
882ffd6a 657 //calculate dQtheta = cos(dOrder*(dPhi-dTheta);the projection of the Q vector on the reference direction dTheta
658 Double_t dQtheta = GetQtheta(vQ, dTheta);
659 //cerr<<"dQtheta for fdenom = "<<dQtheta<<endl;
f1d945a1 660
661 //denominator for differential v
662 if (fHistProR0theta) {
882ffd6a 663 dR0 = fHistProR0theta->GetBinContent(theta+1);
f1d945a1 664 }
665 else {
88e00a8a 666 cout <<"pointer fHistProR0theta does not exist" << endl;
f1d945a1 667 }
882ffd6a 668 //cerr<<"dR0 = "<<dR0 <<endl;
f1d945a1 669
670 if (fUseSum) //sum generating function
671 {
882ffd6a 672 cExpo(0.,dR0*dQtheta);
673 cDenom = dQtheta*(TComplex::Exp(cExpo)); //BP eq 12
f1d945a1 674 //loop over tracks in event
882ffd6a 675 Int_t iNumberOfTracks = anEvent->NumberOfTracks();
676 for (Int_t i=0;i<iNumberOfTracks;i++) {
677 AliFlowTrackSimple* pTrack = anEvent->GetTrack(i);
678 if (pTrack) {
679 if (pTrack->UseForDifferentialFlow()) {
680 dEta = pTrack->Eta();
681 dPt = pTrack->Pt();
682 dPhi = pTrack->Phi();
683 dCosTerm = cos(m*dOrder*(dPhi-dTheta));
684 //cerr<<"dCosTerm = "<<dCosTerm <<endl;
685 cNumer = dCosTerm*(TComplex::Exp(cExpo));
686 if (cNumer.Rho()==0) {cerr<<"WARNING: modulus of cNumer is zero in SecondFillFromFlowEvent"<<endl;}
687 if (fDebug) cerr<<"modulus of cNumer is "<<cNumer.Rho()<<endl;
f1d945a1 688 if (fHist2[theta]) {
882ffd6a 689 fHist2[theta]->Fill(dEta,dPt,cNumer);
f1d945a1 690 }
691 else {
692 cout << "fHist2 pointer mising" <<endl;
693 }
694 }
695 } //if particle
696 else {cerr << "no particle!!!"<<endl;}
697 } //loop over tracks
698
699 } //sum
700 else //product generating function
701 {
882ffd6a 702 cDenom = GetDiffFlow(anEvent, dR0, theta);
f1d945a1 703
704 }//product
705 if (fHistProReDenom && fHistProImDenom) {
882ffd6a 706 fHistProReDenom->Fill(theta,cDenom.Re()); //fill the real part of fDenom
707 fHistProImDenom->Fill(theta,cDenom.Im()); //fill the imaginary part of fDenom
f1d945a1 708 }
709 else {
882ffd6a 710 cout << "Pointers to cDenom mising" << endl;
f1d945a1 711 }
712
713 }//end of loop over theta
714
715 return kTRUE;
716
717
718}
719 //-----------------------------------------------------------------------
882ffd6a 720 Double_t AliFlowAnalysisWithLeeYangZeros::GetQtheta(AliFlowVector aQ, Double_t aTheta)
f1d945a1 721{
882ffd6a 722 //calculate Qtheta. Qtheta is the sum over all particles of cos(dOrder*(dPhi-dTheta)) BP eq. 3
f1d945a1 723 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::GetQtheta()****"<<endl;
724
882ffd6a 725 Double_t dQtheta = 0.;
726 Double_t dOrder = 2.;
f1d945a1 727
882ffd6a 728 dQtheta = aQ.X()*cos(dOrder*aTheta)+aQ.Y()*sin(dOrder*aTheta);
f1d945a1 729
882ffd6a 730 return dQtheta;
f1d945a1 731
732}
733
734
735//-----------------------------------------------------------------------
882ffd6a 736TComplex AliFlowAnalysisWithLeeYangZeros::GetGrtheta(AliFlowEventSimple* anEvent, Double_t aR, Double_t aTheta)
f1d945a1 737{
738 // Product Generating Function for LeeYangZeros method
739 // PG Eq. 3 (J. Phys. G Nucl. Part. Phys 30 S1213 (2004))
740
741 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::GetGrtheta()****"<<endl;
742
743
882ffd6a 744 TComplex cG = TComplex::One();
745 Double_t dOrder = 2.;
746 Double_t dWgt = 1.;
f1d945a1 747
882ffd6a 748 Int_t iNumberOfTracks = anEvent->NumberOfTracks();
f1d945a1 749
882ffd6a 750 for (Int_t i=0;i<iNumberOfTracks;i++) //loop over tracks in event
f1d945a1 751 {
882ffd6a 752 AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
753 if (pTrack){
754 if (pTrack->UseForIntegratedFlow()) {
755 Double_t dPhi = pTrack->Phi();
756 Double_t dGIm = aR * dWgt*cos(dOrder*(dPhi - aTheta));
757 TComplex cGi(1., dGIm);
758 cG *= cGi; //product over all tracks
f1d945a1 759 }
760 }
761 else {cerr << "no particle pointer !!!"<<endl;}
762 }//loop over tracks
763
882ffd6a 764 return cG;
f1d945a1 765
766}
767
768
769//-----------------------------------------------------------------------
882ffd6a 770TComplex AliFlowAnalysisWithLeeYangZeros::GetDiffFlow(AliFlowEventSimple* anEvent, Double_t aR0, Int_t theta)
f1d945a1 771{
772 // Sum for the denominator for diff. flow for the Product Generating Function for LeeYangZeros method
773 // PG Eq. 9 (J. Phys. G Nucl. Part. Phys 30 S1213 (2004))
774 // Also for v1 mixed harmonics: DF Eq. 5
775 // It is the deriverative of Grtheta at r0 divided by Grtheta at r0
776
777 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::GetGrtheta()****"<<endl;
778
882ffd6a 779 TComplex cG = TComplex::One();
780 TComplex cdGr0(0.,0.);
781 Double_t dOrder = 2.;
782 Double_t dWgt = 1.;
f1d945a1 783
882ffd6a 784 Int_t iNumberOfTracks = anEvent->NumberOfTracks();
f1d945a1 785
882ffd6a 786 Int_t iNtheta = AliFlowLYZConstants::kTheta;
787 Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi()/dOrder;
f1d945a1 788
882ffd6a 789 for (Int_t i=0;i<iNumberOfTracks;i++) //loop over tracks in event
f1d945a1 790 {
882ffd6a 791 AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
792 if (pTrack){
793 if (pTrack->UseForDifferentialFlow()) {
794 Double_t dPhi = pTrack->Phi();
795 Double_t dCosTerm = dWgt*cos(dOrder*(dPhi - dTheta));
f1d945a1 796 //GetGr0theta
882ffd6a 797 Double_t dGIm = aR0 * dCosTerm;
798 TComplex cGi(1., dGIm);
799 cG *= cGi; //product over all tracks
f1d945a1 800 //GetdGr0theta
882ffd6a 801 TComplex cCosTermComplex(1., aR0*dCosTerm);
802 cdGr0 +=(dCosTerm / cCosTermComplex); //sum over all tracks
f1d945a1 803 }
804 } //if particle
805 else {cerr << "no particle!!!"<<endl;}
806 }//loop over tracks
807
808
882ffd6a 809 for (Int_t i=0;i<iNumberOfTracks;i++)
f1d945a1 810 {
882ffd6a 811 AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
812 if (pTrack){
813 if (pTrack->UseForDifferentialFlow()) {
814 Double_t dEta = pTrack->Eta();
815 Double_t dPt = pTrack->Pt();
816 Double_t dPhi = pTrack->Phi();
817 Double_t dCosTerm = cos(dOrder*(dPhi-dTheta));
818 TComplex cCosTermComplex(1.,aR0*dCosTerm);
819 TComplex cNumer = cG*dCosTerm/cCosTermComplex; //PG Eq. 9
820 fHist2[theta]->Fill(dEta,dPt,cNumer);
f1d945a1 821 }
822 } //if particle
823 else {cerr << "no particle pointer!!!"<<endl;}
824 }//loop over tracks
825
882ffd6a 826 TComplex cDenom = cG*cdGr0;
827 return cDenom;
f1d945a1 828
829}
830
831//-----------------------------------------------------------------------
832