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