Coding convention fixes
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowCommon / AliFlowAnalysisWithScalarProduct.cxx
CommitLineData
8d312f00 1/*************************************************************************
2* Copyright(c) 1998-2008, 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
16#define AliFlowAnalysisWithScalarProduct_cxx
17
929098e4 18#include "Riostream.h"
19#include "TFile.h"
d7eb18ec 20#include "TList.h"
8d312f00 21#include "TMath.h"
22#include "TProfile.h"
23#include "TVector2.h"
929098e4 24#include "TH1F.h"
8d312f00 25
929098e4 26#include "AliFlowCommonConstants.h"
8d312f00 27#include "AliFlowEventSimple.h"
28#include "AliFlowTrackSimple.h"
29#include "AliFlowCommonHist.h"
395fadba 30#include "AliFlowCommonHistResults.h"
8d312f00 31#include "AliFlowAnalysisWithScalarProduct.h"
929098e4 32#include "AliFlowVector.h"
8d312f00 33
12dc476e 34//////////////////////////////////////////////////////////////////////////////
8d312f00 35// AliFlowAnalysisWithScalarProduct:
36// Description:
37// Maker to analyze Flow with the Scalar product method.
38//
12dc476e 39// authors: N. van der Kolk (kolk@nikhef.nl), A. Bilandzic (anteb@nikhef.nl)
40//////////////////////////////////////////////////////////////////////////////
8d312f00 41
42ClassImp(AliFlowAnalysisWithScalarProduct)
43
44 //-----------------------------------------------------------------------
12dc476e 45 AliFlowAnalysisWithScalarProduct::AliFlowAnalysisWithScalarProduct():
e35ddff0 46 fEventNumber(0),
47 fDebug(kFALSE),
29195b69 48 fWeightsList(NULL),
49 fUsePhiWeights(kFALSE),
50 fPhiWeights(NULL),
e35ddff0 51 fHistList(NULL),
395fadba 52 fHistProUQetaRP(NULL),
53 fHistProUQetaPOI(NULL),
54 fHistProUQPtRP(NULL),
55 fHistProUQPtPOI(NULL),
e4cecfa0 56 fHistProQaQb(NULL),
12dc476e 57 fHistSumOfLinearWeights(NULL),
58 fHistSumOfQuadraticWeights(NULL),
59 fHistProUQQaQbPtRP(NULL),
60 fHistProUQQaQbEtaRP(NULL),
61 fHistProUQQaQbPtPOI(NULL),
62 fHistProUQQaQbEtaPOI(NULL),
395fadba 63 fCommonHists(NULL),
64 fCommonHistsRes(NULL)
8d312f00 65{
8d312f00 66 // Constructor.
29195b69 67 fWeightsList = new TList();
e35ddff0 68 fHistList = new TList();
12dc476e 69
70 // Initialize arrays:
71 for(Int_t i=0;i<3;i++)
72 {
73 fHistSumOfWeightsPtRP[i] = NULL;
74 fHistSumOfWeightsEtaRP[i] = NULL;
75 fHistSumOfWeightsPtPOI[i] = NULL;
76 fHistSumOfWeightsEtaPOI[i] = NULL;
77 }
8d312f00 78}
79 //-----------------------------------------------------------------------
8d312f00 80 AliFlowAnalysisWithScalarProduct::~AliFlowAnalysisWithScalarProduct()
81 {
82 //destructor
29195b69 83 delete fWeightsList;
e35ddff0 84 delete fHistList;
8d312f00 85 }
86
1dfa3c16 87//-----------------------------------------------------------------------
1dfa3c16 88void AliFlowAnalysisWithScalarProduct::WriteHistograms(TString* outputFileName)
89{
90 //store the final results in output .root file
91
92 TFile *output = new TFile(outputFileName->Data(),"RECREATE");
0fe80f88 93 //output->WriteObject(fHistList, "cobjSP","SingleKey");
94 fHistList->SetName("cobjSP");
9455e15e 95 fHistList->SetOwner(kTRUE);
0fe80f88 96 fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
1dfa3c16 97 delete output;
98}
99
b0fda271 100//-----------------------------------------------------------------------
b0fda271 101void AliFlowAnalysisWithScalarProduct::WriteHistograms(TString outputFileName)
102{
103 //store the final results in output .root file
104
105 TFile *output = new TFile(outputFileName.Data(),"RECREATE");
0fe80f88 106 //output->WriteObject(fHistList, "cobjSP","SingleKey");
107 fHistList->SetName("cobjSP");
9455e15e 108 fHistList->SetOwner(kTRUE);
0fe80f88 109 fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
b0fda271 110 delete output;
111}
112
ad87ae62 113//-----------------------------------------------------------------------
ad87ae62 114void AliFlowAnalysisWithScalarProduct::WriteHistograms(TDirectoryFile *outputFileName)
115{
116 //store the final results in output .root file
117 fHistList->SetName("cobjSP");
118 fHistList->SetOwner(kTRUE);
119 outputFileName->Add(fHistList);
120 outputFileName->Write(outputFileName->GetName(), TObject::kSingleKey);
121}
122
8d312f00 123//-----------------------------------------------------------------------
124void AliFlowAnalysisWithScalarProduct::Init() {
125
126 //Define all histograms
d7eb18ec 127 cout<<"---Analysis with the Scalar Product Method--- Init"<<endl;
8d312f00 128
489d5531 129 //save old value and prevent histograms from being added to directory
130 //to avoid name clashes in case multiple analaysis objects are used
131 //in an analysis
132 Bool_t oldHistAddStatus = TH1::AddDirectoryStatus();
133 TH1::AddDirectory(kFALSE);
134
bee2efdc 135 Int_t iNbinsPt = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
136 Double_t dPtMin = AliFlowCommonConstants::GetMaster()->GetPtMin();
137 Double_t dPtMax = AliFlowCommonConstants::GetMaster()->GetPtMax();
138 Int_t iNbinsEta = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
139 Double_t dEtaMin = AliFlowCommonConstants::GetMaster()->GetEtaMin();
140 Double_t dEtaMax = AliFlowCommonConstants::GetMaster()->GetEtaMax();
8d312f00 141
12dc476e 142 fHistProUQetaRP = new TProfile("Flow_UQetaRP_SP","Flow_UQetaRP_SP",iNbinsEta,dEtaMin,dEtaMax,"s");
395fadba 143 fHistProUQetaRP->SetXTitle("{eta}");
144 fHistProUQetaRP->SetYTitle("<uQ>");
145 fHistList->Add(fHistProUQetaRP);
146
12dc476e 147 fHistProUQetaPOI = new TProfile("Flow_UQetaPOI_SP","Flow_UQetaPOI_SP",iNbinsEta,dEtaMin,dEtaMax,"s");
395fadba 148 fHistProUQetaPOI->SetXTitle("{eta}");
149 fHistProUQetaPOI->SetYTitle("<uQ>");
150 fHistList->Add(fHistProUQetaPOI);
151
12dc476e 152 fHistProUQPtRP = new TProfile("Flow_UQPtRP_SP","Flow_UQPtRP_SP",iNbinsPt,dPtMin,dPtMax,"s");
395fadba 153 fHistProUQPtRP->SetXTitle("p_t (GeV)");
154 fHistProUQPtRP->SetYTitle("<uQ>");
155 fHistList->Add(fHistProUQPtRP);
156
12dc476e 157 fHistProUQPtPOI = new TProfile("Flow_UQPtPOI_SP","Flow_UQPtPOI_SP",iNbinsPt,dPtMin,dPtMax,"s");
395fadba 158 fHistProUQPtPOI->SetXTitle("p_t (GeV)");
159 fHistProUQPtPOI->SetYTitle("<uQ>");
160 fHistList->Add(fHistProUQPtPOI);
161
12dc476e 162 fHistProQaQb = new TProfile("Flow_QaQb_SP","Flow_QaQb_SP", 1, 0.5, 1.5,"s");
e4cecfa0 163 fHistProQaQb->SetYTitle("<QaQb>");
164 fHistList->Add(fHistProQaQb);
8d312f00 165
12dc476e 166 fHistSumOfLinearWeights = new TH1D("Flow_SumOfLinearWeights_SP","Flow_SumOfLinearWeights_SP",1,-0.5, 0.5);
167 fHistSumOfLinearWeights -> SetYTitle("sum (*)");
168 fHistSumOfLinearWeights -> SetXTitle("sum (Ma*Mb)");
169 fHistList->Add(fHistSumOfLinearWeights);
170
171 fHistSumOfQuadraticWeights = new TH1D("Flow_SumOfQuadraticWeights_SP","Flow_SumOfQuadraticWeights_SP",1,-0.5, 0.5);
172 fHistSumOfQuadraticWeights -> SetYTitle("sum (*)");
173 fHistSumOfQuadraticWeights -> SetXTitle("sum (Ma*Mb)^2");
174 fHistList->Add(fHistSumOfQuadraticWeights);
175
176 fHistProUQQaQbPtRP = new TProfile("Flow_UQQaQbPtRP_SP","Flow_UQQaQbPtRP_SP",iNbinsPt,dPtMin,dPtMax);
177 fHistProUQQaQbPtRP -> SetYTitle("<*>");
178 fHistProUQQaQbPtRP -> SetXTitle("<Qu QaQb>");
179 fHistList->Add(fHistProUQQaQbPtRP);
180
181 fHistProUQQaQbEtaRP = new TProfile("Flow_UQQaQbEtaRP_SP","Flow_UQQaQbEtaRP_SP",iNbinsEta,dEtaMin,dEtaMax);
182 fHistProUQQaQbEtaRP -> SetYTitle("<*>");
183 fHistProUQQaQbEtaRP -> SetXTitle("<Qu QaQb>");
184 fHistList->Add(fHistProUQQaQbEtaRP);
185
186 fHistProUQQaQbPtPOI = new TProfile("Flow_UQQaQbPtPOI_SP","Flow_UQQaQbPtPOI_SP",iNbinsPt,dPtMin,dPtMax);
187 fHistProUQQaQbPtPOI -> SetYTitle("<*>");
188 fHistProUQQaQbPtPOI -> SetXTitle("<Qu QaQb>");
189 fHistList->Add(fHistProUQQaQbPtPOI);
190
191 fHistProUQQaQbEtaPOI = new TProfile("Flow_UQQaQbEtaPOI_SP","Flow_UQQaQbEtaPOI_SP",iNbinsEta,dEtaMin,dEtaMax);
192 fHistProUQQaQbEtaPOI -> SetYTitle("<*>");
193 fHistProUQQaQbEtaPOI -> SetXTitle("<Qu QaQb>");
194 fHistList->Add(fHistProUQQaQbEtaPOI);
195
196 TString weightFlag[3] = {"w_Qu_","w_Qu^2_","w_QuQaQb_"};
197 for(Int_t i=0;i<3;i++)
198 {
cff1bec4 199 fHistSumOfWeightsPtRP[i] = new TH1D(Form("Flow_SumOfWeights%sPtRP_SP",weightFlag[i].Data()),
200 Form("Flow_SumOfWeights%sPtRP_SP",weightFlag[i].Data()),iNbinsPt,dPtMin,dPtMax);
12dc476e 201 fHistSumOfWeightsPtRP[i] -> SetYTitle("sum (*)");
202 fHistSumOfWeightsPtRP[i] -> SetXTitle("p_{T}");
203 fHistList->Add(fHistSumOfWeightsPtRP[i]);
204
cff1bec4 205 fHistSumOfWeightsEtaRP[i] = new TH1D(Form("Flow_SumOfWeights%sEtaRP_SP",weightFlag[i].Data()),
206 Form("Flow_SumOfWeights%sEtaRP_SP",weightFlag[i].Data()),iNbinsEta,dEtaMin,dEtaMax);
12dc476e 207 fHistSumOfWeightsEtaRP[i] -> SetYTitle("sum (*)");
208 fHistSumOfWeightsEtaRP[i] -> SetXTitle("#eta");
209 fHistList->Add(fHistSumOfWeightsEtaRP[i]);
210
cff1bec4 211 fHistSumOfWeightsPtPOI[i] = new TH1D(Form("Flow_SumOfWeights%sPtPOI_SP",weightFlag[i].Data()),
212 Form("Flow_SumOfWeights%sPtPOI_SP",weightFlag[i].Data()),iNbinsPt,dPtMin,dPtMax);
12dc476e 213 fHistSumOfWeightsPtPOI[i] -> SetYTitle("sum (*)");
214 fHistSumOfWeightsPtPOI[i] -> SetXTitle("p_{T}");
215 fHistList->Add(fHistSumOfWeightsPtPOI[i]);
216
cff1bec4 217 fHistSumOfWeightsEtaPOI[i] = new TH1D(Form("Flow_SumOfWeights%sEtaPOI_SP",weightFlag[i].Data()),
218 Form("Flow_SumOfWeights%sEtaPOI_SP",weightFlag[i].Data()),iNbinsEta,dEtaMin,dEtaMax);
12dc476e 219 fHistSumOfWeightsEtaPOI[i] -> SetYTitle("sum (*)");
220 fHistSumOfWeightsEtaPOI[i] -> SetXTitle("#eta");
221 fHistList->Add(fHistSumOfWeightsEtaPOI[i]);
222 }
223
04f6283b 224 fCommonHists = new AliFlowCommonHist("AliFlowCommonHistSP");
099e1753 225 fHistList->Add(fCommonHists);
395fadba 226 fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsSP");
227 fHistList->Add(fCommonHistsRes);
228
29195b69 229 //weights
230 if(fUsePhiWeights) {
231 if(!fWeightsList) {
232 cout<<"WARNING: fWeightsList is NULL in the Scalar Product method."<<endl;
233 exit(0);
234 }
235 if(fWeightsList->FindObject("phi_weights")) {
236 fPhiWeights = dynamic_cast<TH1F*>(fWeightsList->FindObject("phi_weights"));
237 fHistList->Add(fPhiWeights);
238 } else {
239 cout<<"WARNING: histogram with phi weights is not accessible in Scalar Product"<<endl;
240 exit(0);
241 }
242 } // end of if(fUsePhiWeights)
243
e2d51347 244 fEventNumber = 0; //set number of events to zero
489d5531 245
246 TH1::AddDirectory(oldHistAddStatus);
e2d51347 247}
248
8d312f00 249//-----------------------------------------------------------------------
8232a5ec 250void AliFlowAnalysisWithScalarProduct::Make(AliFlowEventSimple* anEvent) {
8d312f00 251
489d5531 252 //Calculate flow based on <QaQb/MaMb> = <v^2>
12dc476e 253
254 //Fill histograms
8232a5ec 255 if (anEvent) {
8d312f00 256
cbbaf54a 257 //get Q vectors for the eta-subevents
b125a454 258 AliFlowVector* vQarray = new AliFlowVector[2];
29195b69 259 if (fUsePhiWeights) {
260 anEvent->Get2Qsub(vQarray,2,fWeightsList,kTRUE);
261 } else {
262 anEvent->Get2Qsub(vQarray);
263 }
12dc476e 264 //subevent a
b125a454 265 AliFlowVector vQa = vQarray[0];
12dc476e 266 //subevent b
b125a454 267 AliFlowVector vQb = vQarray[1];
12dc476e 268
7a01f4a7 269 //check that the subevents are not empty:
12dc476e 270 Double_t dMa = vQa.GetMult();
12dc476e 271 Double_t dMb = vQb.GetMult();
7a01f4a7 272 if (dMa != 0 && dMb != 0) {
273
274 //fill control histograms
275 fCommonHists->FillControlHistograms(anEvent);
276
277 //get total Q vector = the sum of subevent a and subevent b
278 AliFlowVector vQ = vQa + vQb;
279 //weight the Q vectors for the subevents by the multiplicity
280 //Note: Weight Q only in the particle loop when it is clear if it should be (m-1) or M
281 Double_t dQXa = vQa.X()/dMa;
282 Double_t dQYa = vQa.Y()/dMa;
283 vQa.Set(dQXa,dQYa);
284
285 Double_t dQXb = vQb.X()/dMb;
286 Double_t dQYb = vQb.Y()/dMb;
287 vQb.Set(dQXb,dQYb);
12dc476e 288
7a01f4a7 289 //scalar product of the two subevents
290 Double_t dQaQb = (vQa*vQb);
291 fHistProQaQb -> Fill(1.,dQaQb,dMa*dMb); //Fill (QaQb/MaMb) with weight (MaMb).
292 //needed for the error calculation:
293 fHistSumOfLinearWeights -> Fill(0.,dMa*dMb);
294 fHistSumOfQuadraticWeights -> Fill(0.,pow(dMa*dMb,2.));
295
296 //loop over the tracks of the event
297 AliFlowTrackSimple* pTrack = NULL;
298 Int_t iNumberOfTracks = anEvent->NumberOfTracks();
299 Double_t dMq = vQ.GetMult();
300
301 for (Int_t i=0;i<iNumberOfTracks;i++)
302 {
303 pTrack = anEvent->GetTrack(i) ;
304 if (pTrack){
305 Double_t dPhi = pTrack->Phi();
306 //set default phi weight to 1
307 Double_t dW = 1.;
308 //phi weight of pTrack
309 if(fUsePhiWeights && fPhiWeights) {
310 Int_t iNBinsPhi = fPhiWeights->GetNbinsX();
311 dW = fPhiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*iNBinsPhi/TMath::TwoPi())));
312 //bin = 1 + value*nbins/range
313 //TMath::Floor rounds to the lower integer
314 }
12dc476e 315
7a01f4a7 316 //calculate vU
317 TVector2 vU;
318 Double_t dUX = TMath::Cos(2*dPhi);
319 Double_t dUY = TMath::Sin(2*dPhi);
320 vU.Set(dUX,dUY);
321 Double_t dModulus = vU.Mod();
322 if (dModulus!=0.) vU.Set(dUX/dModulus,dUY/dModulus); // make length 1
323 else cerr<<"dModulus is zero!"<<endl;
12dc476e 324
7a01f4a7 325 //redefine the Q vector and devide by its multiplicity
326 TVector2 vQm;
327 Double_t dQmX = 0.;
328 Double_t dQmY = 0.;
329 //subtract particle from the flowvector if used to define it
330 if (pTrack->InSubevent(0) || pTrack->InSubevent(1)) {
331 dQmX = (vQ.X() - dW*dUX)/(dMq-1);
332 dQmY = (vQ.Y() - dW*dUY)/(dMq-1);
333
334 vQm.Set(dQmX,dQmY);
335
336 //dUQ = scalar product of vU and vQm
337 Double_t dUQ = (vU * vQm);
338 Double_t dPt = pTrack->Pt();
339 Double_t dEta = pTrack->Eta();
340 //fill the profile histograms
341 if (pTrack->InRPSelection()) {
342 fHistProUQetaRP -> Fill(dEta,dUQ,(dMq-1)); //Fill (Qu/(Mq-1)) with weight (Mq-1)
343 //needed for the error calculation:
344 fHistProUQQaQbEtaRP -> Fill(dEta,dUQ*dQaQb,(dMq-1)*dMa*dMb); //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb
345 fHistProUQPtRP -> Fill(dPt,dUQ,(dMq-1)); //Fill (Qu/(Mq-1)) with weight (Mq-1)
346 fHistProUQQaQbPtRP -> Fill(dPt,dUQ*dQaQb,(dMq-1)*dMa*dMb); //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb
347
348 fHistSumOfWeightsEtaRP[0]->Fill(dEta,(dMq-1)); // sum of Mq-1
349 fHistSumOfWeightsEtaRP[1]->Fill(dEta,pow((dMq-1),2.));// sum of (Mq-1)^2
350 fHistSumOfWeightsEtaRP[2]->Fill(dEta,(dMq-1)*dMa*dMb);// sum of (Mq-1)*MaMb
351 fHistSumOfWeightsPtRP[0]->Fill(dPt,(dMq-1)); // sum of Mq-1
352 fHistSumOfWeightsPtRP[1]->Fill(dPt,pow((dMq-1),2.)); // sum of (Mq-1)^2
353 fHistSumOfWeightsPtRP[2]->Fill(dPt,(dMq-1)*dMa*dMb); // sum of (Mq-1)*MaMb
354 }
355 if (pTrack->InPOISelection()) {
356 fHistProUQetaPOI -> Fill(dEta,dUQ,(dMq-1));//Fill (Qu/(Mq-1)) with weight (Mq-1)
357 //needed for the error calculation:
358 fHistProUQQaQbEtaPOI -> Fill(dEta,dUQ*dQaQb,(dMq-1)*dMa*dMb); //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb
359 fHistProUQPtPOI -> Fill(dPt,dUQ,(dMq-1)); //Fill (Qu/(Mq-1)) with weight (Mq-1)
360 fHistProUQQaQbPtPOI -> Fill(dPt,dUQ*dQaQb,(dMq-1)*dMa*dMb); //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb
361
362 fHistSumOfWeightsEtaPOI[0]->Fill(dEta,(dMq-1)); // sum of Mq-1
363 fHistSumOfWeightsEtaPOI[1]->Fill(dEta,pow((dMq-1),2.));// sum of (Mq-1)^2
364 fHistSumOfWeightsEtaPOI[2]->Fill(dEta,(dMq-1)*dMa*dMb);// sum of (Mq-1)*MaMb
365 fHistSumOfWeightsPtPOI[0]->Fill(dPt,(dMq-1)); // sum of Mq-1
366 fHistSumOfWeightsPtPOI[1]->Fill(dPt,pow((dMq-1),2.)); // sum of (Mq-1)^2
367 fHistSumOfWeightsPtPOI[2]->Fill(dPt,(dMq-1)*dMa*dMb); // sum of (Mq-1)*MaMb
368 }
369
370 } else { //do not subtract the particle from the flowvector
371 dQmX = vQ.X()/dMq;
372 dQmY = vQ.Y()/dMq;
373
374 vQm.Set(dQmX,dQmY);
12dc476e 375
7a01f4a7 376 //dUQ = scalar product of vU and vQm
377 Double_t dUQ = (vU * vQm);
378 Double_t dPt = pTrack->Pt();
379 Double_t dEta = pTrack->Eta();
380 //fill the profile histograms
381 if (pTrack->InRPSelection()) {
382 fHistProUQetaRP -> Fill(dEta,dUQ,dMq); //Fill (Qu/Mq) with weight Mq
383 //needed for the error calculation:
384 fHistProUQQaQbEtaRP -> Fill(dEta,dUQ*dQaQb,dMq*dMa*dMb); //Fill [Qu/Mq]*[QaQb/MaMb] with weight Mq*MaMb
385 fHistProUQPtRP -> Fill(dPt,dUQ,dMq); //Fill (Qu/Mq) with weight Mq
386 fHistProUQQaQbPtRP -> Fill(dPt,dUQ*dQaQb,dMq*dMa*dMb); //Fill [Qu/Mq]*[QaQb/MaMb] with weight Mq*MaMb
387
388 fHistSumOfWeightsEtaRP[0]->Fill(dEta,dMq); // sum of Mq
389 fHistSumOfWeightsEtaRP[1]->Fill(dEta,pow(dMq,2.));// sum of Mq^2
390 fHistSumOfWeightsEtaRP[2]->Fill(dEta,dMq*dMa*dMb);// sum of Mq*MaMb
391 fHistSumOfWeightsPtRP[0]->Fill(dPt,dMq); // sum of Mq
392 fHistSumOfWeightsPtRP[1]->Fill(dPt,pow(dMq,2.)); // sum of Mq^2
393 fHistSumOfWeightsPtRP[2]->Fill(dPt,dMq*dMa*dMb); // sum of Mq*MaMb
394 }
395 if (pTrack->InPOISelection()) {
396 fHistProUQetaPOI -> Fill(dEta,dUQ,dMq); //Fill (Qu/Mq) with weight Mq
397 //needed for the error calculation:
398 fHistProUQQaQbEtaPOI -> Fill(dEta,dUQ*dQaQb,dMq*dMa*dMb); //Fill [Qu/Mq]*[QaQb/MaMb] with weight Mq*MaMb
399 fHistProUQPtPOI -> Fill(dPt,dUQ,dMq); //Fill (Qu/Mq) with weight Mq
400 fHistProUQQaQbPtPOI -> Fill(dPt,dUQ*dQaQb,dMq*dMa*dMb); //Fill [Qu/Mq]*[QaQb/MaMb] with weight Mq*MaMb
401
402 fHistSumOfWeightsEtaPOI[0]->Fill(dEta,dMq); // sum of Mq
403 fHistSumOfWeightsEtaPOI[1]->Fill(dEta,pow(dMq,2.));// sum of Mq^2
404 fHistSumOfWeightsEtaPOI[2]->Fill(dEta,dMq*dMa*dMb);// sum of Mq*MaMb
405 fHistSumOfWeightsPtPOI[0]->Fill(dPt,dMq); // sum of Mq
406 fHistSumOfWeightsPtPOI[1]->Fill(dPt,pow(dMq,2.)); // sum of Mq^2
407 fHistSumOfWeightsPtPOI[2]->Fill(dPt,dMq*dMa*dMb); // sum of Mq*MaMb
408 }
409 }//track not in subevents
410
411 }//track
12dc476e 412
7a01f4a7 413 }//loop over tracks
414
415 fEventNumber++;
416 // cout<<"@@@@@ "<<fEventNumber<<" events processed"<<endl;
12dc476e 417
7a01f4a7 418 }// subevents not empty
b125a454 419 delete [] vQarray;
7a01f4a7 420 } //event
421}//end of Make()
8d312f00 422
12dc476e 423//--------------------------------------------------------------------
fd46c3dd 424void AliFlowAnalysisWithScalarProduct::GetOutputHistograms(TList *outputListHistos){
425
426 //get pointers to all output histograms (called before Finish())
12dc476e 427
fd46c3dd 428 if (outputListHistos) {
429 //Get the common histograms from the output list
430 AliFlowCommonHist *pCommonHist = dynamic_cast<AliFlowCommonHist*>
431 (outputListHistos->FindObject("AliFlowCommonHistSP"));
432 AliFlowCommonHistResults *pCommonHistResults = dynamic_cast<AliFlowCommonHistResults*>
433 (outputListHistos->FindObject("AliFlowCommonHistResultsSP"));
434 TProfile* pHistProQaQb = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_QaQb_SP"));
fd46c3dd 435 TProfile* pHistProUQetaRP = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQetaRP_SP"));
436 TProfile* pHistProUQetaPOI = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQetaPOI_SP"));
437 TProfile* pHistProUQPtRP = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQPtRP_SP"));
438 TProfile* pHistProUQPtPOI = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQPtPOI_SP"));
12dc476e 439 TH1D* pHistSumOfLinearWeights = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_SumOfLinearWeights_SP"));
440 TH1D* pHistSumOfQuadraticWeights = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_SumOfQuadraticWeights_SP"));
441 TProfile* pHistProUQQaQbPtRP = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQQaQbPtRP_SP"));
442 TProfile* pHistProUQQaQbEtaRP = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQQaQbEtaRP_SP"));
443 TProfile* pHistProUQQaQbPtPOI = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQQaQbPtPOI_SP"));
444 TProfile* pHistProUQQaQbEtaPOI = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQQaQbEtaPOI_SP"));
cff1bec4 445
12dc476e 446 TString weightFlag[3] = {"w_Qu_","w_Qu^2_","w_QuQaQb_"};
447 TH1D* pHistSumOfWeightsPtRP[3] = {NULL};
448 TH1D* pHistSumOfWeightsEtaRP[3] = {NULL};
449 TH1D* pHistSumOfWeightsPtPOI[3] = {NULL};
450 TH1D* pHistSumOfWeightsEtaPOI[3] = {NULL};
451 for(Int_t i=0;i<3;i++)
452 {
cff1bec4 453 pHistSumOfWeightsPtRP[i] = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("Flow_SumOfWeights%sPtRP_SP",weightFlag[i].Data())));
454 pHistSumOfWeightsEtaRP[i] = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("Flow_SumOfWeights%sEtaRP_SP",weightFlag[i].Data())));
455 pHistSumOfWeightsPtPOI[i] = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("Flow_SumOfWeights%sPtPOI_SP",weightFlag[i].Data())));
456 pHistSumOfWeightsEtaPOI[i] = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("Flow_SumOfWeights%sEtaPOI_SP",weightFlag[i].Data())));
12dc476e 457 }
458
459 //pass the pointers to the task
460 if (pCommonHist && pCommonHistResults && pHistProQaQb &&
461 pHistProUQetaRP && pHistProUQetaPOI && pHistProUQPtRP && pHistProUQPtPOI
462 && pHistSumOfLinearWeights && pHistSumOfQuadraticWeights && pHistProUQQaQbPtRP
463 && pHistProUQQaQbEtaRP && pHistProUQQaQbPtPOI && pHistProUQQaQbEtaPOI) {
fd46c3dd 464 this -> SetCommonHists(pCommonHist);
465 this -> SetCommonHistsRes(pCommonHistResults);
466 this -> SetHistProQaQb(pHistProQaQb);
fd46c3dd 467 this -> SetHistProUQetaRP(pHistProUQetaRP);
468 this -> SetHistProUQetaPOI(pHistProUQetaPOI);
469 this -> SetHistProUQPtRP(pHistProUQPtRP);
470 this -> SetHistProUQPtPOI(pHistProUQPtPOI);
12dc476e 471 this -> SetHistProUQQaQbPtRP(pHistProUQQaQbPtRP);
472 this -> SetHistProUQQaQbEtaRP(pHistProUQQaQbEtaRP);
473 this -> SetHistProUQQaQbPtPOI(pHistProUQQaQbPtPOI);
474 this -> SetHistProUQQaQbEtaPOI(pHistProUQQaQbEtaPOI);
475 this -> SetHistSumOfLinearWeights(pHistSumOfLinearWeights);
476 this -> SetHistSumOfQuadraticWeights(pHistSumOfQuadraticWeights);
cff1bec4 477 }
12dc476e 478
479 for(Int_t i=0;i<3;i++)
480 {
481 if(pHistSumOfWeightsPtRP[i]) this -> SetHistSumOfWeightsPtRP(pHistSumOfWeightsPtRP[i],i);
cff1bec4 482 if(pHistSumOfWeightsEtaRP[i]) this -> SetHistSumOfWeightsEtaRP(pHistSumOfWeightsEtaRP[i],i);
483 if(pHistSumOfWeightsPtPOI[i]) this -> SetHistSumOfWeightsPtPOI(pHistSumOfWeightsPtPOI[i],i);
484 if(pHistSumOfWeightsEtaPOI[i]) this -> SetHistSumOfWeightsEtaPOI(pHistSumOfWeightsEtaPOI[i],i);
12dc476e 485 }
12dc476e 486 } // end of if(outputListHistos)
fd46c3dd 487}
488
489//--------------------------------------------------------------------
8d312f00 490void AliFlowAnalysisWithScalarProduct::Finish() {
491
489d5531 492 //calculate flow and fill the AliFlowCommonHistResults
395fadba 493 if (fDebug) cout<<"AliFlowAnalysisWithScalarProduct::Finish()"<<endl;
494
495 cout<<"*************************************"<<endl;
496 cout<<"*************************************"<<endl;
497 cout<<" Integrated flow from "<<endl;
498 cout<<" Scalar product "<<endl;
499 cout<<endl;
e4cecfa0 500
bee2efdc 501 Int_t iNbinsPt = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
502 Int_t iNbinsEta = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
12dc476e 503
12dc476e 504 //Calculate reference flow (noname)
505 //----------------------------------
506 //weighted average over (QaQb/MaMb) with weight (MaMb)
507 Double_t dQaQb = fHistProQaQb->GetBinContent(1);
524798bf 508 Double_t dV = 0.;
509 if(dQaQb>=0.)
510 {
511 dV = TMath::Sqrt(dQaQb);
512 }
12dc476e 513 //statistical error of dQaQb:
514 // statistical error = term1 * spread * term2:
515 // term1 = sqrt{sum_{i=1}^{N} w^2}/(sum_{i=1}^{N} w)
516 // term2 = 1/sqrt(1-term1^2)
517 Double_t dSpreadQaQb = fHistProQaQb->GetBinError(1);
518 Double_t dSumOfLinearWeights = fHistSumOfLinearWeights->GetBinContent(1);
519 Double_t dSumOfQuadraticWeights = fHistSumOfQuadraticWeights->GetBinContent(1);
520 Double_t dTerm1 = 0.;
521 Double_t dTerm2 = 0.;
522 if(dSumOfLinearWeights) {
523 dTerm1 = pow(dSumOfQuadraticWeights,0.5)/dSumOfLinearWeights;
524 }
525 if(1.-pow(dTerm1,2.) > 0.) {
526 dTerm2 = 1./pow(1-pow(dTerm1,2.),0.5);
527 }
528 Double_t dStatErrorQaQb = dTerm1 * dSpreadQaQb * dTerm2;
529 //calculate the statistical error
530 Double_t dVerr = 0.;
531 if(dQaQb > 0.) {
532 dVerr = (1./(2.*pow(dQaQb,0.5)))*dStatErrorQaQb;
533 }
534 fCommonHistsRes->FillIntegratedFlow(dV,dVerr);
535 cout<<"dV = "<<dV<<" +- "<<dVerr<<endl;
536
537 //Calculate differential flow and integrated flow (RP, POI)
538 //---------------------------------------------------------
12dc476e 539 //v as a function of eta for RP selection
a9de3704 540 for(Int_t b=1;b<iNbinsEta+1;b++) {
12dc476e 541 Double_t duQpro = fHistProUQetaRP->GetBinContent(b);
542 Double_t dv2pro = 0.;
543 if (dV!=0.) { dv2pro = duQpro/dV; }
544 //calculate the statistical error
545 Double_t dv2ProErr = CalculateStatisticalError(b, dStatErrorQaQb, fHistProUQetaRP, fHistProUQQaQbEtaRP, fHistSumOfWeightsEtaRP);
12dc476e 546 //fill TH1D
cff1bec4 547 fCommonHistsRes->FillDifferentialFlowEtaRP(b, dv2pro, dv2ProErr);
12dc476e 548 } //loop over bins b
549
80f72ed6 550
12dc476e 551 //v as a function of eta for POI selection
a9de3704 552 for(Int_t b=1;b<iNbinsEta+1;b++) {
12dc476e 553 Double_t duQpro = fHistProUQetaPOI->GetBinContent(b);
554 Double_t dv2pro = 0.;
555 if (dV!=0.) { dv2pro = duQpro/dV; }
556 //calculate the statistical error
557 Double_t dv2ProErr = CalculateStatisticalError(b, dStatErrorQaQb, fHistProUQetaPOI, fHistProUQQaQbEtaPOI, fHistSumOfWeightsEtaPOI);
558
559 //fill TH1D
560 fCommonHistsRes->FillDifferentialFlowEtaPOI(b, dv2pro, dv2ProErr);
561 } //loop over bins b
562
80f72ed6 563
12dc476e 564 //v as a function of Pt for RP selection
565 TH1F* fHistPtRP = fCommonHists->GetHistPtRP(); //for calculating integrated flow
566 Double_t dVRP = 0.;
567 Double_t dSumRP = 0.;
568 Double_t dErrVRP =0.;
569
a9de3704 570 for(Int_t b=1;b<iNbinsPt+1;b++) {
12dc476e 571 Double_t duQpro = fHistProUQPtRP->GetBinContent(b);
572 Double_t dv2pro = 0.;
573 if (dV!=0.) { dv2pro = duQpro/dV; }
574 //calculate the statistical error
575 Double_t dv2ProErr = CalculateStatisticalError(b, dStatErrorQaQb, fHistProUQPtRP, fHistProUQQaQbPtRP, fHistSumOfWeightsPtRP);
576
577 //fill TH1D
578 fCommonHistsRes->FillDifferentialFlowPtRP(b, dv2pro, dv2ProErr);
579
580 //calculate integrated flow for RP selection
581 if (fHistPtRP){
582 Double_t dYieldPt = fHistPtRP->GetBinContent(b);
583 dVRP += dv2pro*dYieldPt;
584 dSumRP +=dYieldPt;
585 dErrVRP += dYieldPt*dYieldPt*dv2ProErr*dv2ProErr;
586 } else { cout<<"fHistPtRP is NULL"<<endl; }
587 } //loop over bins b
588
589 if (dSumRP != 0.) {
590 dVRP /= dSumRP; //the pt distribution should be normalised
591 dErrVRP /= (dSumRP*dSumRP);
592 dErrVRP = TMath::Sqrt(dErrVRP);
80f72ed6 593 }
12dc476e 594 fCommonHistsRes->FillIntegratedFlowRP(dVRP,dErrVRP);
595 cout<<"dV(RP) = "<<dVRP<<" +- "<<dErrVRP<<endl;
596
80f72ed6 597
12dc476e 598 //v as a function of Pt for POI selection
599 TH1F* fHistPtPOI = fCommonHists->GetHistPtPOI(); //for calculating integrated flow
600 Double_t dVPOI = 0.;
601 Double_t dSumPOI = 0.;
602 Double_t dErrVPOI =0.;
603
a9de3704 604 for(Int_t b=1;b<iNbinsPt+1;b++) {
12dc476e 605 Double_t duQpro = fHistProUQPtPOI->GetBinContent(b);
606 Double_t dv2pro = 0.;
607 if (dV!=0.) { dv2pro = duQpro/dV; }
608 //calculate the statistical error
609 Double_t dv2ProErr = CalculateStatisticalError(b, dStatErrorQaQb, fHistProUQPtPOI, fHistProUQQaQbPtPOI, fHistSumOfWeightsPtPOI);
610
611 //fill TH1D
612 fCommonHistsRes->FillDifferentialFlowPtPOI(b, dv2pro, dv2ProErr);
395fadba 613
12dc476e 614 //calculate integrated flow for POI selection
615 if (fHistPtPOI){
616 Double_t dYieldPt = fHistPtPOI->GetBinContent(b);
617 dVPOI += dv2pro*dYieldPt;
618 dSumPOI +=dYieldPt;
619 dErrVPOI += dYieldPt*dYieldPt*dv2ProErr*dv2ProErr;
620 } else { cout<<"fHistPtPOI is NULL"<<endl; }
621 } //loop over bins b
395fadba 622
12dc476e 623 if (dSumPOI != 0.) {
489d5531 624 dVPOI /= dSumPOI; //the pt distribution should be normalised
12dc476e 625 dErrVPOI /= (dSumPOI*dSumPOI);
626 dErrVPOI = TMath::Sqrt(dErrVPOI);
627 }
628 fCommonHistsRes->FillIntegratedFlowPOI(dVPOI,dErrVPOI);
629 cout<<"dV(POI) = "<<dVPOI<<" +- "<<dErrVPOI<<endl;
e4cecfa0 630
395fadba 631 cout<<endl;
632 cout<<"*************************************"<<endl;
633 cout<<"*************************************"<<endl;
12dc476e 634
395fadba 635 //cout<<".....finished"<<endl;
12dc476e 636}
8232a5ec 637
638
12dc476e 639//--------------------------------------------------------------------
640Double_t AliFlowAnalysisWithScalarProduct::CalculateStatisticalError(Int_t b, Double_t aStatErrorQaQb, TProfile* pHistProUQ, TProfile* pHistProUQQaQb, TH1D** pHistSumOfWeights) {
12dc476e 641 //calculate the statistical error for differential flow for bin b
642 Double_t duQproSpread = pHistProUQ->GetBinError(b);
643 Double_t sumOfMq = pHistSumOfWeights[0]->GetBinContent(b);
644 Double_t sumOfMqSquared = pHistSumOfWeights[1]->GetBinContent(b);
645 Double_t dTerm1 = 0.;
646 Double_t dTerm2 = 0.;
647 if(sumOfMq) {
648 dTerm1 = (pow(sumOfMqSquared,0.5)/sumOfMq);
649 }
650 if(1.-pow(dTerm1,2.)>0.) {
651 dTerm2 = 1./pow(1.-pow(dTerm1,2.),0.5);
652 }
653 Double_t duQproErr = dTerm1*duQproSpread*dTerm2;
12dc476e 654 // covariances:
655 Double_t dTerm1Cov = pHistSumOfWeights[2]->GetBinContent(b);
656 Double_t dTerm2Cov = fHistSumOfLinearWeights->GetBinContent(1);
657 Double_t dTerm3Cov = sumOfMq;
658 Double_t dWeightedCovariance = 0.;
659 if(dTerm2Cov*dTerm3Cov>0.) {
660 Double_t dDenominator = 1.-dTerm1Cov/(dTerm2Cov*dTerm3Cov);
661 Double_t dPrefactor = dTerm1Cov/(dTerm2Cov*dTerm3Cov);
662 if(dDenominator!=0) {
663 Double_t dCovariance = (pHistProUQQaQb->GetBinContent(b)-fHistProQaQb->GetBinContent(1)*pHistProUQ->GetBinContent(b))/dDenominator;
664 dWeightedCovariance = dCovariance*dPrefactor;
665 }
666 }
cff1bec4 667
12dc476e 668 Double_t dv2ProErr = 0.; // final statitical error:
669 if(fHistProQaQb->GetBinContent(1)>0.) {
670 Double_t dv2ProErrorSquared = (1./4.)*pow(fHistProQaQb->GetBinContent(1),-3.)*
671 (pow(pHistProUQ->GetBinContent(b),2.)*pow(aStatErrorQaQb,2.)
672 + 4.*pow(fHistProQaQb->GetBinContent(1),2.)*pow(duQproErr,2.)
673 - 4.*fHistProQaQb->GetBinContent(1)*pHistProUQ->GetBinContent(b)*dWeightedCovariance);
674 if(dv2ProErrorSquared>0.) dv2ProErr = pow(dv2ProErrorSquared,0.5);
675 }
676
677 return dv2ProErr;
678}
679//--------------------------------------------------------------------