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