1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
17 #include "Riostream.h" //needed as include
18 #include "AliFlowCommonConstants.h" //needed as include
19 #include "AliFlowLYZConstants.h" //needed as include
20 #include "AliFlowAnalysisWithLeeYangZeros.h"
21 #include "AliFlowLYZHist1.h" //needed as include
22 #include "AliFlowLYZHist2.h" //needed as include
23 #include "AliFlowCommonHist.h" //needed as include
24 #include "AliFlowCommonHistResults.h" //needed as include
25 #include "AliFlowEventSimple.h" //needed as include
26 #include "AliFlowTrackSimple.h" //needed as include
30 #include "TMath.h" //needed as include
31 #include "TFile.h" //needed as include
40 //Description: Maker to analyze Flow using the LeeYangZeros method
41 // Equation numbers are from Big Paper (BP): Nucl. Phys. A 727, 373 (2003)
42 // Practical Guide (PG): J. Phys. G: Nucl. Part. Phys. 30, S1213 (2004)
43 // Adapted from StFlowLeeYangZerosMaker.cxx
44 // by Markus Oldenberg and Art Poskanzer, LBNL
45 // with advice from Jean-Yves Ollitrault and Nicolas Borghini
47 //Author: Naomi van der Kolk (kolk@nikhef.nl)
51 ClassImp(AliFlowAnalysisWithLeeYangZeros)
53 //-----------------------------------------------------------------------
55 AliFlowAnalysisWithLeeYangZeros::AliFlowAnalysisWithLeeYangZeros():
68 fHistProR0theta(NULL),
69 fHistProReDenom(NULL),
70 fHistProImDenom(NULL),
71 fHistProReDtheta(NULL),
72 fHistProImDtheta(NULL),
73 fHistQsumforChi(NULL),
79 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::AliFlowAnalysisWithLeeYangZeros default constructor****"<<endl;
81 fHistList = new TList();
82 fFirstRunList = new TList();
84 for(Int_t i = 0;i<5;i++)
90 fQsum = new TVector2();
94 //-----------------------------------------------------------------------
97 AliFlowAnalysisWithLeeYangZeros::~AliFlowAnalysisWithLeeYangZeros()
100 if (fDebug) cout<<"****~AliFlowAnalysisWithLeeYangZeros****"<<endl;
103 delete fFirstRunList;
107 //-----------------------------------------------------------------------
109 void AliFlowAnalysisWithLeeYangZeros::WriteHistograms(TString* outputFileName)
111 //store the final results in output .root file
112 TFile *output = new TFile(outputFileName->Data(),"RECREATE");
114 output->mkdir("cobjLYZ1","cobjLYZ1");
115 output->cd("cobjLYZ1");
118 output->mkdir("cobjLYZ2","cobjLYZ2");
119 output->cd("cobjLYZ2");
125 //-----------------------------------------------------------------------
127 Bool_t AliFlowAnalysisWithLeeYangZeros::Init()
130 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::Init()****"<<endl;
133 Int_t iNtheta = AliFlowLYZConstants::kTheta;
134 Int_t iNbinsPt = AliFlowCommonConstants::GetNbinsPt();
135 Int_t iNbinsEta = AliFlowCommonConstants::GetNbinsEta();
137 Double_t dPtMin = AliFlowCommonConstants::GetPtMin();
138 Double_t dPtMax = AliFlowCommonConstants::GetPtMax();
139 Double_t dEtaMin = AliFlowCommonConstants::GetEtaMin();
140 Double_t dEtaMax = AliFlowCommonConstants::GetEtaMax();
142 //for control histograms
143 if (fFirstRun){ fCommonHists = new AliFlowCommonHist("AliFlowCommonHistLYZ1");
144 fHistList->Add(fCommonHists);
145 fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsLYZ1");
146 fHistList->Add(fCommonHistsRes);
148 else { fCommonHists = new AliFlowCommonHist("AliFlowCommonHistLYZ2");
149 fHistList->Add(fCommonHists);
150 fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsLYZ2");
151 fHistList->Add(fCommonHistsRes);
154 fHistQsumforChi = new TH1F("Flow_QsumforChi_LYZ","Flow_QsumforChi_LYZ",3,-1.,2.);
155 fHistQsumforChi->SetXTitle("Qsum.X , Qsum.Y, Q2sum");
156 fHistQsumforChi->SetYTitle("value");
157 fHistList->Add(fHistQsumforChi);
159 //for first loop over events
161 fHistProR0theta = new TProfile("First_FlowPro_r0theta_LYZ","First_FlowPro_r0theta_LYZ",iNtheta,-0.5,iNtheta-0.5);
162 fHistProR0theta->SetXTitle("#theta");
163 fHistProR0theta->SetYTitle("r_{0}^{#theta}");
164 fHistList->Add(fHistProR0theta);
166 fHistProVtheta = new TProfile("First_FlowPro_Vtheta_LYZ","First_FlowPro_Vtheta_LYZ",iNtheta,-0.5,iNtheta-0.5);
167 fHistProVtheta->SetXTitle("#theta");
168 fHistProVtheta->SetYTitle("V_{n}^{#theta}");
169 fHistList->Add(fHistProVtheta);
171 //class AliFlowLYZHist1 defines the histograms: fHistProGtheta, fHistProReGtheta, fHistProImGtheta, fHistProR0theta
172 for (Int_t theta=0;theta<iNtheta;theta++) {
173 TString name = "AliFlowLYZHist1_";
175 fHist1[theta]=new AliFlowLYZHist1(theta, name);
176 fHistList->Add(fHist1[theta]);
180 //for second loop over events
182 fHistProReDenom = new TProfile("Second_FlowPro_ReDenom_LYZ","Second_FlowPro_ReDenom_LYZ" , iNtheta, -0.5, iNtheta-0.5);
183 fHistProReDenom->SetXTitle("#theta");
184 fHistProReDenom->SetYTitle("Re(Q^{#theta}e^{ir_{0}^{#theta}Q^{#theta}})");
185 fHistList->Add(fHistProReDenom);
187 fHistProImDenom = new TProfile("Second_FlowPro_ImDenom_LYZ","Second_FlowPro_ImDenom_LYZ" , iNtheta, -0.5, iNtheta-0.5);
188 fHistProImDenom->SetXTitle("#theta");
189 fHistProImDenom->SetYTitle("Im(Q^{#theta}e^{ir_{0}^{#theta}Q^{#theta}})");
190 fHistList->Add(fHistProImDenom);
192 fHistProVeta = new TProfile("Second_FlowPro_Veta_LYZ","Second_FlowPro_Veta_LYZ",iNbinsEta,dEtaMin,dEtaMax);
193 fHistProVeta->SetXTitle("rapidity");
194 fHistProVeta->SetYTitle("v (%)");
195 fHistList->Add(fHistProVeta);
197 fHistProVPt = new TProfile("Second_FlowPro_VPt_LYZ","Second_FlowPro_VPt_LYZ",iNbinsPt,dPtMin,dPtMax);
198 fHistProVPt->SetXTitle("Pt");
199 fHistProVPt->SetYTitle("v (%)");
200 fHistList->Add(fHistProVPt);
202 fHistProReDtheta = new TProfile("Second_FlowPro_ReDtheta_LYZ","Second_FlowPro_ReDtheta_LYZ",iNtheta, -0.5, iNtheta-0.5);
203 fHistProReDtheta->SetXTitle("#theta");
204 fHistProReDtheta->SetYTitle("Re(D^{#theta})");
205 fHistList->Add(fHistProReDtheta);
207 fHistProImDtheta = new TProfile("Second_FlowPro_ImDtheta_LYZ","Second_FlowPro_ImDtheta_LYZ",iNtheta, -0.5, iNtheta-0.5);
208 fHistProImDtheta->SetXTitle("#theta");
209 fHistProImDtheta->SetYTitle("Im(D^{#theta})");
210 fHistList->Add(fHistProImDtheta);
212 //class AliFlowLYZHist2 defines the histograms:
213 for (Int_t theta=0;theta<iNtheta;theta++) {
214 TString name = "AliFlowLYZHist2_";
216 fHist2[theta]=new AliFlowLYZHist2(theta,name);
217 fHistList->Add(fHist2[theta]);
220 //read histogram fHistProR0theta from the first run list
222 fHistProR0theta = (TProfile*)fFirstRunList->FindObject("First_FlowPro_r0theta_LYZ");
223 if (!fHistProR0theta) {cout<<"fHistProR0theta has a NULL pointer!"<<endl;}
224 fHistList->Add(fHistProR0theta);
225 } else { cout<<"list is NULL pointer!"<<endl; }
230 if (fDebug) cout<<"****Histograms initialised****"<<endl;
232 fEventNumber = 0; //set event counter to zero
237 //-----------------------------------------------------------------------
239 Bool_t AliFlowAnalysisWithLeeYangZeros::Make(AliFlowEventSimple* anEvent)
242 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::Make()****"<<endl;
244 //get tracks from event
247 fCommonHists->FillControlHistograms(anEvent);
248 FillFromFlowEvent(anEvent);
251 fCommonHists->FillControlHistograms(anEvent);
252 SecondFillFromFlowEvent(anEvent);
257 cout<<"##### FlowLeeYangZero: Stack pointer null"<<endl;
260 cout<<"^^^^read event "<<fEventNumber<<endl;
268 //-----------------------------------------------------------------------
269 Bool_t AliFlowAnalysisWithLeeYangZeros::Finish()
272 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::Finish()****"<<endl;
274 //define variables for both runs
275 Double_t dJ01 = 2.405;
276 Int_t iNtheta = AliFlowLYZConstants::kTheta;
277 //set the event number
278 SetEventNumber((int)fCommonHists->GetHistMultOrig()->GetEntries());
279 //cout<<"number of events processed is "<<fEventNumber<<endl;
281 //set the sum of Q vectors
282 fQsum->Set(fHistQsumforChi->GetBinContent(1),fHistQsumforChi->GetBinContent(2));
283 SetQ2sum(fHistQsumforChi->GetBinContent(3));
287 Double_t dVtheta = 0;
290 for (Int_t theta=0;theta<iNtheta;theta++)
292 //get the first minimum r0
293 fHist1[theta]->FillGtheta();
294 dR0 = fHist1[theta]->GetR0();
296 //calculate integrated flow
297 if (dR0!=0) dVtheta = dJ01/dR0;
298 //for estimating systematic error resulting from d0
299 Double_t dBinsize = (AliFlowLYZConstants::fgMax)/(AliFlowLYZConstants::kNbins);
300 Double_t dVplus = dJ01/(dR0+dBinsize);
301 Double_t dVmin = dJ01/(dR0-dBinsize);
303 Double_t dvplus = dVplus;
304 Double_t dvmin = dVmin;
305 cout<<"dv = "<<dv<<" and dvplus = "<<dvplus<< " and dvmin = "<<dvmin<<endl;
307 //fill the histograms
308 fHistProR0theta->Fill(theta,dR0);
309 fHistProVtheta->Fill(theta,dVtheta);
310 //get average value of fVtheta = fV
312 } //end of loop over theta
314 //get average value of fVtheta = fV
318 Double_t dSigma2 = 0;
320 if (fEventNumber!=0) {
321 *fQsum /= fEventNumber;
322 fQ2sum /= fEventNumber;
323 dSigma2 = fQ2sum - TMath::Power(fQsum->X(),2.) - TMath::Power(fQsum->Y(),2.) - TMath::Power(dV,2.); //BP eq. 62
324 if (dSigma2>0) dChi = dV/TMath::Sqrt(dSigma2);
326 fCommonHistsRes->FillChi(dChi);
327 cout<<"dV = "<<dV<<" and chi = "<<dChi<<endl;
330 // recalculate statistical errors on integrated flow
331 //combining 5 theta angles to 1 relative error BP eq. 89
332 Double_t dRelErr2comb = 0.;
333 Int_t iEvts = fEventNumber;
335 for (Int_t theta=0;theta<iNtheta;theta++){
336 Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
337 Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
339 Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
341 dRelErr2comb += (1/(2*iEvts*(dJ01*dJ01)*TMath::BesselJ1(dJ01)*
342 TMath::BesselJ1(dJ01)))*
343 (dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2)) +
344 dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2)));
346 dRelErr2comb /= iNtheta;
348 Double_t dRelErrcomb = TMath::Sqrt(dRelErr2comb);
350 //copy content of profile into TH1D and add error
351 Double_t dv2pro = dV; //in the case that fv is equal to fV
352 Double_t dv2Err = dv2pro*dRelErrcomb ;
353 cout<<"dv2pro +- dv2Err = "<<dv2pro<<" +- "<<dv2Err<<endl;
354 fCommonHistsRes->FillIntegratedFlow(dv2pro, dv2Err);
357 if (fDebug) cout<<"****histograms filled****"<<endl;
360 fEventNumber =0; //set to zero for second round over events
366 //calculate differential flow
368 TComplex cDenom, cNumer, cDtheta;
370 TComplex i = TComplex::I();
371 Double_t dBesselRatio[3] = {1., 1.202, 2.69};
372 Int_t iNbinsPt = AliFlowCommonConstants::GetNbinsPt();
373 Int_t iNbinsEta = AliFlowCommonConstants::GetNbinsEta();
375 Double_t dEta, dPt, dReRatio, dVeta, dVPt;
378 Double_t dVtheta = 0.;
380 Double_t dReDenom = 0.;
381 Double_t dImDenom = 0.;
382 for (Int_t theta=0;theta<iNtheta;theta++) { //loop over theta
383 if (fHistProR0theta) {
384 dR0 = fHistProR0theta->GetBinContent(theta+1);
385 if (fDebug) cerr<<"dR0 = "<<dR0<<endl;
386 if (dR0!=0) dVtheta = dJ01/dR0;
388 // BP Eq. 9 -> Vn^theta = j01/r0^theta
389 if (fHistProReDenom && fHistProImDenom) {
390 dReDenom = fHistProReDenom->GetBinContent(theta+1);
391 dImDenom = fHistProImDenom->GetBinContent(theta+1);
394 cout << "Hist pointer fDenom in file does not exist" <<endl;
399 cout << "Hist pointer R0theta in file does not exist" <<endl;
401 //} //loop over theta
403 cDenom(dReDenom,dImDenom);
405 //for new method and use by others (only with the sum generating function):
407 dR0 = fHistProR0theta->GetBinContent(theta+1);
408 cDtheta = dR0*cDenom/dJ01;
409 fHistProReDtheta->Fill(theta,cDtheta.Re());
410 fHistProImDtheta->Fill(theta,cDtheta.Im());
413 cDenom *= TComplex::Power(i, m-1);
414 //cerr<<"TComplex::Power(i, m-1) = "<<TComplex::Power(i, m-1).Rho()<<endl; //checked ok
416 //v as a function of eta
417 for (Int_t be=1;be<=iNbinsEta;be++) {
418 dEta = fHist2[theta]->GetBinCenter(be);
419 cNumer = fHist2[theta]->GetNumerEta(be);
420 if (cNumer.Rho()==0) {
421 if (fDebug) cerr<<"WARNING: modulus of cNumer is zero in Finish()"<<endl;
424 else if (cDenom.Rho()==0) {
425 if (fDebug) cerr<<"WARNING: modulus of cDenom is zero"<<endl;
429 //if ( j==1 && theta==0) cerr<<"modulus of cNumer = "<<cNumer.Rho() <<endl; //always a number smaller than 1, or 0.
430 dReRatio = (cNumer/cDenom).Re();
433 dVeta = dBesselRatio[m-1]*dReRatio*dVtheta*100.; //BP eq. 12
434 //if ( j==1 && theta==0) cerr<<"eta = "<<dEta<<" cerr::dReRatio for eta = "<<dReRatio<<" cerr::dVeta for eta = "<<dVeta<<endl;
436 fHistProVeta->Fill(dEta,dVeta);
437 } //loop over bins be
439 //v as a function of Pt
440 for (Int_t bp=1;bp<=iNbinsPt;bp++) {
441 dPt = fHist2[theta]->GetBinCenterPt(bp);
442 cNumer = fHist2[theta]->GetNumerPt(bp);
443 if (cNumer.Rho()==0) {
444 if (fDebug) cerr<<"modulus of cNumer is zero"<<endl;
447 else if (cDenom.Rho()==0) {
448 if (fDebug) cerr<<"modulus of cDenom is zero"<<endl;
452 //if ( j==1 && theta==0) cerr<<"modulus of fNumer = "<<fNumer.Rho() <<endl; //always a number smaller than 1, or 0.
453 dReRatio = (cNumer/cDenom).Re();
456 dVPt = dBesselRatio[m-1]*dReRatio*dVtheta*100.; //BP eq. 12
458 fHistProVPt->Fill(dPt,dVPt);
459 } //loop over bins bp
461 }//end of loop over theta
463 //calculate the average of fVtheta = fV
466 //sigma2 and chi (for statistical error calculations)
467 Double_t dSigma2 = 0;
469 if (fEventNumber!=0) {
470 *fQsum /= fEventNumber;
471 //cerr<<"fQsum.X() = "<<fQsum.X()<<endl;
472 //cerr<<"fQsum.Y() = "<<fQsum.Y()<<endl;
473 fQ2sum /= fEventNumber;
474 //cout<<"fQ2sum = "<<fQ2sum<<endl;
475 dSigma2 = fQ2sum - TMath::Power(fQsum->X(),2.) - TMath::Power(fQsum->Y(),2.) - TMath::Power(dV,2.); //BP eq. 62
476 if (dSigma2>0) dChi = dV/TMath::Sqrt(dSigma2);
478 fCommonHistsRes->FillChi(dChi);
479 cout<<"dV = "<<dV<<" and chi = "<<dChi<<endl;
482 //copy content of profile into TH1D and add error
483 for(Int_t b=0;b<iNbinsPt;b++) {
484 Double_t dv2pro = fHistProVPt->GetBinContent(b);
485 Double_t dNprime = fCommonHists->GetEntriesInPtBin(b);
486 Double_t dErrdifcomb = 0.; //set error to zero
487 Double_t dErr2difcomb = 0.; //set error to zero
490 for (Int_t theta=0;theta<iNtheta;theta++) {
491 Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
492 Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
494 Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
496 dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)*
497 TMath::BesselJ1(dJ01)))*
498 ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
499 (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
503 if (dErr2difcomb!=0.) {
504 dErr2difcomb /= iNtheta;
505 dErrdifcomb = TMath::Sqrt(dErr2difcomb)*100;
506 //cerr<<"dErrdifcomb = "<<dErrdifcomb<<endl;
508 else {dErrdifcomb = 0.;}
511 fCommonHistsRes->FillDifferentialFlow(b, dv2pro, dErrdifcomb);
516 cout<<"----LYZ analysis finished....----"<<endl<<endl;
522 //-----------------------------------------------------------------------
524 Bool_t AliFlowAnalysisWithLeeYangZeros::FillFromFlowEvent(AliFlowEventSimple* anEvent)
526 // Get event quantities from AliFlowEvent for all particles
528 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::FillFromFlowEvent()****"<<endl;
531 cout<<"##### FlowLeeYangZero: FlowEvent pointer null"<<endl;
536 TComplex cExpo, cGtheta, cGthetaNew, cZ;
537 Int_t iNtheta = AliFlowLYZConstants::kTheta;
538 Int_t iNbins = AliFlowLYZConstants::kNbins;
542 Double_t dOrder = 2.;
545 AliFlowVector vQ = anEvent->GetQ();
546 //weight by the multiplicity
549 if (vQ.GetMult() != 0) {
550 dQX = vQ.X()/vQ.GetMult();
551 dQY = vQ.Y()/vQ.GetMult();
555 //for chi calculation:
557 fHistQsumforChi->SetBinContent(1,fQsum->X());
558 fHistQsumforChi->SetBinContent(2,fQsum->Y());
560 fHistQsumforChi->SetBinContent(3,fQ2sum);
561 //cerr<<"fQ2sum = "<<fQ2sum<<endl;
563 for (Int_t theta=0;theta<iNtheta;theta++)
565 Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi()/dOrder;
567 //calculate dQtheta = cos(dOrder*(fPhi-dTheta);the projection of the Q vector on the reference direction dTheta
568 Double_t dQtheta = GetQtheta(vQ, dTheta);
570 for (Int_t bin=1;bin<=iNbins;bin++)
572 Double_t dR = fHist1[theta]->GetBinCenter(bin); //bincentre of bins in histogram //FIXED???
573 //if (theta == 0) cerr<<"cerr::dR = "<<dR<<endl;
576 //calculate the sum generating function
577 cExpo(0.,dR*dQtheta); //Re=0 ; Im=dR*dQtheta
578 cGtheta = TComplex::Exp(cExpo);
582 //calculate the product generating function
583 cGtheta = GetGrtheta(anEvent, dR, dTheta); //make this function
584 if (cGtheta.Rho2() > 100.) break;
587 fHist1[theta]->Fill(dR,cGtheta); //fill real and imaginary part of cGtheta
598 //-----------------------------------------------------------------------
599 Bool_t AliFlowAnalysisWithLeeYangZeros::SecondFillFromFlowEvent(AliFlowEventSimple* anEvent)
601 //for differential flow
603 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::SecondFillFromFlowEvent()****"<<endl;
606 cout<<"##### FlowLeeYangZero: FlowEvent pointer null"<<endl;
611 TComplex cExpo, cDenom, cNumer,cCosTermComplex;
612 Double_t dPhi, dEta, dPt;
616 Double_t dOrder = 2.;
617 Int_t iNtheta = AliFlowLYZConstants::kTheta;
621 AliFlowVector vQ = anEvent->GetQ();
622 //weight by the multiplicity
625 if (vQ.GetMult() != 0) {
626 dQX = vQ.X()/vQ.GetMult();
627 dQY = vQ.Y()/vQ.GetMult();
631 //for chi calculation:
633 fHistQsumforChi->SetBinContent(1,fQsum->X());
634 fHistQsumforChi->SetBinContent(2,fQsum->Y());
636 fHistQsumforChi->SetBinContent(3,fQ2sum);
638 for (Int_t theta=0;theta<iNtheta;theta++)
640 Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi()/dOrder;
642 //calculate dQtheta = cos(dOrder*(dPhi-dTheta);the projection of the Q vector on the reference direction dTheta
643 Double_t dQtheta = GetQtheta(vQ, dTheta);
644 //cerr<<"dQtheta for fdenom = "<<dQtheta<<endl;
646 //denominator for differential v
647 if (fHistProR0theta) {
648 dR0 = fHistProR0theta->GetBinContent(theta+1);
651 cout <<"pointer fHistProR0theta does not exist" << endl;
653 //cerr<<"dR0 = "<<dR0 <<endl;
655 if (fUseSum) //sum generating function
657 cExpo(0.,dR0*dQtheta);
658 cDenom = dQtheta*(TComplex::Exp(cExpo)); //BP eq 12
659 //loop over tracks in event
660 Int_t iNumberOfTracks = anEvent->NumberOfTracks();
661 for (Int_t i=0;i<iNumberOfTracks;i++) {
662 AliFlowTrackSimple* pTrack = anEvent->GetTrack(i);
664 if (pTrack->UseForDifferentialFlow()) {
665 dEta = pTrack->Eta();
667 dPhi = pTrack->Phi();
668 dCosTerm = cos(m*dOrder*(dPhi-dTheta));
669 //cerr<<"dCosTerm = "<<dCosTerm <<endl;
670 cNumer = dCosTerm*(TComplex::Exp(cExpo));
671 if (cNumer.Rho()==0) {cerr<<"WARNING: modulus of cNumer is zero in SecondFillFromFlowEvent"<<endl;}
672 if (fDebug) cerr<<"modulus of cNumer is "<<cNumer.Rho()<<endl;
674 fHist2[theta]->Fill(dEta,dPt,cNumer);
677 cout << "fHist2 pointer mising" <<endl;
681 else {cerr << "no particle!!!"<<endl;}
685 else //product generating function
687 cDenom = GetDiffFlow(anEvent, dR0, theta);
690 if (fHistProReDenom && fHistProImDenom) {
691 fHistProReDenom->Fill(theta,cDenom.Re()); //fill the real part of fDenom
692 fHistProImDenom->Fill(theta,cDenom.Im()); //fill the imaginary part of fDenom
695 cout << "Pointers to cDenom mising" << endl;
698 }//end of loop over theta
704 //-----------------------------------------------------------------------
705 Double_t AliFlowAnalysisWithLeeYangZeros::GetQtheta(AliFlowVector aQ, Double_t aTheta)
707 //calculate Qtheta. Qtheta is the sum over all particles of cos(dOrder*(dPhi-dTheta)) BP eq. 3
708 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::GetQtheta()****"<<endl;
710 Double_t dQtheta = 0.;
711 Double_t dOrder = 2.;
713 dQtheta = aQ.X()*cos(dOrder*aTheta)+aQ.Y()*sin(dOrder*aTheta);
720 //-----------------------------------------------------------------------
721 TComplex AliFlowAnalysisWithLeeYangZeros::GetGrtheta(AliFlowEventSimple* anEvent, Double_t aR, Double_t aTheta)
723 // Product Generating Function for LeeYangZeros method
724 // PG Eq. 3 (J. Phys. G Nucl. Part. Phys 30 S1213 (2004))
726 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::GetGrtheta()****"<<endl;
729 TComplex cG = TComplex::One();
730 Double_t dOrder = 2.;
733 Int_t iNumberOfTracks = anEvent->NumberOfTracks();
735 for (Int_t i=0;i<iNumberOfTracks;i++) //loop over tracks in event
737 AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
739 if (pTrack->UseForIntegratedFlow()) {
740 Double_t dPhi = pTrack->Phi();
741 Double_t dGIm = aR * dWgt*cos(dOrder*(dPhi - aTheta));
742 TComplex cGi(1., dGIm);
743 cG *= cGi; //product over all tracks
746 else {cerr << "no particle pointer !!!"<<endl;}
754 //-----------------------------------------------------------------------
755 TComplex AliFlowAnalysisWithLeeYangZeros::GetDiffFlow(AliFlowEventSimple* anEvent, Double_t aR0, Int_t theta)
757 // Sum for the denominator for diff. flow for the Product Generating Function for LeeYangZeros method
758 // PG Eq. 9 (J. Phys. G Nucl. Part. Phys 30 S1213 (2004))
759 // Also for v1 mixed harmonics: DF Eq. 5
760 // It is the deriverative of Grtheta at r0 divided by Grtheta at r0
762 if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::GetGrtheta()****"<<endl;
764 TComplex cG = TComplex::One();
765 TComplex cdGr0(0.,0.);
766 Double_t dOrder = 2.;
769 Int_t iNumberOfTracks = anEvent->NumberOfTracks();
771 Int_t iNtheta = AliFlowLYZConstants::kTheta;
772 Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi()/dOrder;
774 for (Int_t i=0;i<iNumberOfTracks;i++) //loop over tracks in event
776 AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
778 if (pTrack->UseForDifferentialFlow()) {
779 Double_t dPhi = pTrack->Phi();
780 Double_t dCosTerm = dWgt*cos(dOrder*(dPhi - dTheta));
782 Double_t dGIm = aR0 * dCosTerm;
783 TComplex cGi(1., dGIm);
784 cG *= cGi; //product over all tracks
786 TComplex cCosTermComplex(1., aR0*dCosTerm);
787 cdGr0 +=(dCosTerm / cCosTermComplex); //sum over all tracks
790 else {cerr << "no particle!!!"<<endl;}
794 for (Int_t i=0;i<iNumberOfTracks;i++)
796 AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
798 if (pTrack->UseForDifferentialFlow()) {
799 Double_t dEta = pTrack->Eta();
800 Double_t dPt = pTrack->Pt();
801 Double_t dPhi = pTrack->Phi();
802 Double_t dCosTerm = cos(dOrder*(dPhi-dTheta));
803 TComplex cCosTermComplex(1.,aR0*dCosTerm);
804 TComplex cNumer = cG*dCosTerm/cCosTermComplex; //PG Eq. 9
805 fHist2[theta]->Fill(dEta,dPt,cNumer);
808 else {cerr << "no particle pointer!!!"<<endl;}
811 TComplex cDenom = cG*cdGr0;
816 //-----------------------------------------------------------------------