]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG/FLOW/Base/AliFlowAnalysisWithScalarProduct.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
AOD handeling, added mass to the flowtracks, new task to reuse flowevent
[u/mrichter/AliRoot.git] / PWG / FLOW / Base / 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
02914cfa 16#define ALIFLOWANALYSISWITHSCALARPRODUCT_CXX
8d312f00 17
929098e4 18#include "TFile.h"
d7eb18ec 19#include "TList.h"
8d312f00 20#include "TMath.h"
02914cfa 21#include "TString.h"
8d312f00 22#include "TProfile.h"
23#include "TVector2.h"
a554b203 24#include "TH1D.h"
ad1b7114 25#include "TH1F.h"
a554b203 26#include "TH2D.h"
8d312f00 27
929098e4 28#include "AliFlowCommonConstants.h"
8d312f00 29#include "AliFlowEventSimple.h"
a93de0f0 30#include "AliFlowVector.h"
8d312f00 31#include "AliFlowTrackSimple.h"
32#include "AliFlowCommonHist.h"
395fadba 33#include "AliFlowCommonHistResults.h"
8d312f00 34#include "AliFlowAnalysisWithScalarProduct.h"
8d312f00 35
12dc476e 36//////////////////////////////////////////////////////////////////////////////
8d312f00 37// AliFlowAnalysisWithScalarProduct:
02914cfa 38// Description: Maker to analyze Flow from the Scalar Product method.
39// authors: Naomi van del Kolk (kolk@nikhef.nl)
40// Ante Bilandzic (anteb@nikhef.nl)
41// mods: Carlos Perez (cperez@nikhef.nl)
12dc476e 42//////////////////////////////////////////////////////////////////////////////
8d312f00 43
44ClassImp(AliFlowAnalysisWithScalarProduct)
45
02914cfa 46//-----------------------------------------------------------------------
47AliFlowAnalysisWithScalarProduct::AliFlowAnalysisWithScalarProduct():
48fDebug(0),
49fUsePhiWeights(0),
50fApplyCorrectionForNUA(0),
51fHarmonic(2),
52fNormalizationType(1),
53fTotalQvector(3),
54fWeightsList(NULL),
55fHistList(NULL),
56fHistProConfig(NULL),
57fHistProQaQbNorm(NULL),
58fHistSumOfWeights(NULL),
59fHistProNUAq(NULL),
60fHistProQNorm(NULL),
61fHistProQaQb(NULL),
62fHistProQaQbM(NULL),
63fHistMaMb(NULL),
64fHistQNormQaQbNorm(NULL),
65fHistQaNormMa(NULL),
66fHistQbNormMb(NULL),
67fResolution(NULL),
68fHistQaQb(NULL),
69fHistQaQbCos(NULL),
70fCommonHists(NULL),
71fCommonHistsuQ(NULL),
72fCommonHistsRes(NULL)
8d312f00 73{
02914cfa 74 //ctor
75 for(int i=0; i!=2; ++i) {
76 fPhiWeightsSub[i] = NULL;
77 for(int j=0; j!=2; ++j) {
78 fHistProUQ[i][j] = NULL;
79 fHistProUQQaQb[i][j] = NULL;
80 fHistProNUAu[i][j][0] = NULL;
81 fHistProNUAu[i][j][1] = NULL;
82 for(int k=0; k!=3; ++k)
83 fHistSumOfWeightsu[i][j][k] = NULL;
29ecccee 84 }
29ecccee 85 }
8d312f00 86}
1dfa3c16 87//-----------------------------------------------------------------------
02914cfa 88AliFlowAnalysisWithScalarProduct::~AliFlowAnalysisWithScalarProduct()
1dfa3c16 89{
02914cfa 90 //destructor
91 delete fWeightsList;
92 delete fHistList;
1dfa3c16 93}
8d312f00 94//-----------------------------------------------------------------------
95void AliFlowAnalysisWithScalarProduct::Init() {
8d312f00 96 //Define all histograms
02914cfa 97 printf("---Analysis with the Scalar Product Method--- Init\n");
98 printf("--- fNormalizationType %d ---\n", fNormalizationType);
489d5531 99 //save old value and prevent histograms from being added to directory
100 //to avoid name clashes in case multiple analaysis objects are used
101 //in an analysis
102 Bool_t oldHistAddStatus = TH1::AddDirectoryStatus();
103 TH1::AddDirectory(kFALSE);
104
02914cfa 105 fHistList = new TList();
106 fHistList->SetName("cobjSP");
107 fHistList->SetOwner();
108
109 TList *uQRelated = new TList();
110 uQRelated->SetName("uQ");
111 uQRelated->SetOwner();
112
113 TList *nuaRelated = new TList();
114 nuaRelated->SetName("NUA");
115 nuaRelated->SetOwner();
116
117 TList *errorRelated = new TList();
118 errorRelated->SetName("error");
119 errorRelated->SetOwner();
120
9812922b 121 TList *tQARelated = new TList();
122 tQARelated->SetName("QA");
123 tQARelated->SetOwner();
02914cfa 124
ad1b7114 125 fCommonHists = new AliFlowCommonHist("AliFlowCommonHist_SP");
02914cfa 126 (fCommonHists->GetHarmonic())->Fill(0.5,fHarmonic); // store harmonic
127 fHistList->Add(fCommonHists);
128
ad1b7114 129 fCommonHistsuQ = new AliFlowCommonHist("AliFlowCommonHist_uQ");
02914cfa 130 (fCommonHistsuQ->GetHarmonic())->Fill(0.5,fHarmonic); // store harmonic
131 fHistList->Add(fCommonHistsuQ);
132
ad1b7114 133 fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResults_SP","",fHarmonic);
02914cfa 134 fHistList->Add(fCommonHistsRes);
135
136 fHistProConfig = new TProfile("FlowPro_Flags_SP","Flow_Flags_SP",4,0.5,4.5,"s");
137 fHistProConfig->GetXaxis()->SetBinLabel(1,"fApplyCorrectionForNUA");
138 fHistProConfig->GetXaxis()->SetBinLabel(2,"fNormalizationType");
139 fHistProConfig->GetXaxis()->SetBinLabel(3,"fUsePhiWeights");
140 fHistProConfig->GetXaxis()->SetBinLabel(4,"fHarmonic");
141 fHistProConfig->Fill(1,fApplyCorrectionForNUA);
142 fHistProConfig->Fill(2,fNormalizationType);
143 fHistProConfig->Fill(3,fUsePhiWeights);
144 fHistProConfig->Fill(4,fHarmonic);
145 fHistList->Add(fHistProConfig);
381342ec 146
a554b203 147 fHistProQaQbNorm = new TProfile("FlowPro_QaQbNorm_SP","FlowPro_QaQbNorm_SP", 1, 0.5, 1.5,"s");
02914cfa 148 fHistProQaQbNorm->SetYTitle("<QaQb/Na/Nb>");
149 errorRelated->Add(fHistProQaQbNorm);
150
ad1b7114 151 fHistProNUAq = new TProfile("FlowPro_NUAq_SP","FlowPro_NUAq_SP", 6, 0.5, 6.5,"s");
02914cfa 152 fHistProNUAq->GetXaxis()->SetBinLabel( 1,"<<sin(#Phi_{a})>>");
153 fHistProNUAq->GetXaxis()->SetBinLabel( 2,"<<cos(#Phi_{a})>>");
154 fHistProNUAq->GetXaxis()->SetBinLabel( 3,"<<sin(#Phi_{b})>>");
155 fHistProNUAq->GetXaxis()->SetBinLabel( 4,"<<cos(#Phi_{b})>>");
156 fHistProNUAq->GetXaxis()->SetBinLabel( 5,"<<sin(#Phi_{t})>>");
157 fHistProNUAq->GetXaxis()->SetBinLabel( 6,"<<cos(#Phi_{t})>>");
158 nuaRelated->Add(fHistProNUAq);
159
160 fHistSumOfWeights = new TH1D("Flow_SumOfWeights_SP","Flow_SumOfWeights_SP",2,0.5, 2.5);
161 fHistSumOfWeights->GetXaxis()->SetBinLabel( 1,"#Sigma Na*Nb");
162 fHistSumOfWeights->GetXaxis()->SetBinLabel( 2,"#Sigma (Na*Nb)^2");
163 errorRelated->Add(fHistSumOfWeights);
164
165 TString sPOI[2] = {"RP","POI"}; // backward compatibility
166 TString sEta[2] = {"Pt","eta"}; // backward compatibility
167 TString sTitle[2] = {"p_{T} [GeV]","#eta"};
168 TString sWeights[3] = {"uQ","uQuQ","uQQaQb"};
169 Int_t iNbins[2];
170 Double_t dMin[2], dMax[2];
171 iNbins[0] = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
172 iNbins[1] = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
173 dMin[0] = AliFlowCommonConstants::GetMaster()->GetPtMin();
174 dMin[1] = AliFlowCommonConstants::GetMaster()->GetEtaMin();
175 dMax[0] = AliFlowCommonConstants::GetMaster()->GetPtMax();
176 dMax[1] = AliFlowCommonConstants::GetMaster()->GetEtaMax();
177 for(Int_t iPOI=0; iPOI!=2; ++iPOI) for(Int_t iSpace=0; iSpace!=2; ++iSpace) {
178 // uQ
ad1b7114 179 fHistProUQ[iPOI][iSpace] = new TProfile( Form( "FlowPro_UQ_%s%s_SP", sEta[iSpace].Data(), sPOI[iPOI].Data() ),
02914cfa 180 Form( "FlowPro_UQ%s%s_SP", sEta[iSpace].Data(), sPOI[iPOI].Data() ),
181 iNbins[iSpace], dMin[iSpace], dMax[iSpace], "s");
182 fHistProUQ[iPOI][iSpace]->SetXTitle(sTitle[iSpace].Data());
183 fHistProUQ[iPOI][iSpace]->SetYTitle("<uQ>");
184 uQRelated->Add(fHistProUQ[iPOI][iSpace]);
185
186 // NUAu
187 fHistProNUAu[iPOI][iSpace][0] = new TProfile( Form("FlowPro_NUAu_%s%s_IM_SP", sEta[iSpace].Data(), sPOI[iPOI].Data() ),
188 Form("FlowPro_NUAu_%s%s_IM_SP", sEta[iSpace].Data(), sPOI[iPOI].Data() ),
189 iNbins[iSpace], dMin[iSpace], dMax[iSpace]);
190 fHistProNUAu[iPOI][iSpace][0]->SetXTitle(sTitle[iSpace].Data());
191 nuaRelated->Add(fHistProNUAu[iPOI][iSpace][0]);
192 fHistProNUAu[iPOI][iSpace][1] = new TProfile( Form("FlowPro_NUAu_%s%s_RE_SP", sEta[iSpace].Data(), sPOI[iPOI].Data() ),
193 Form("FlowPro_NUAu_%s%s_RE_SP", sEta[iSpace].Data(), sPOI[iPOI].Data() ),
194 iNbins[iSpace], dMin[iSpace], dMax[iSpace]);
195 fHistProNUAu[iPOI][iSpace][1]->SetXTitle(sTitle[iSpace].Data());
196 nuaRelated->Add(fHistProNUAu[iPOI][iSpace][1]);
197
198 // uQ QaQb
199 fHistProUQQaQb[iPOI][iSpace] = new TProfile( Form("FlowPro_UQQaQb_%s%s_SP", sEta[iSpace].Data(), sPOI[iPOI].Data() ),
200 Form("FlowPro_UQQaQb_%s%s_SP", sEta[iSpace].Data(), sPOI[iPOI].Data() ),
201 iNbins[iSpace], dMin[iSpace], dMax[iSpace]);
202 fHistProUQQaQb[iPOI][iSpace]->SetXTitle(sTitle[iSpace].Data());
203 fHistProUQQaQb[iPOI][iSpace]-> SetYTitle("<Qu QaQb>");
204 errorRelated->Add(fHistProUQQaQb[iPOI][iSpace]);
205
206 // uWeights
207 for(Int_t i=0; i!=3; ++i) {
208 fHistSumOfWeightsu[iPOI][iSpace][i] = new TH1D( Form("Flow_SumOfWeights_%s%s_%s_SP",sWeights[i].Data(),sPOI[iPOI].Data(),sEta[iSpace].Data()),
209 Form("Flow_SumOfWeights_%s%s_%s_SP",sWeights[i].Data(),sPOI[iPOI].Data(),sEta[iSpace].Data()),
210 iNbins[iSpace], dMin[iSpace], dMax[iSpace]);
211 fHistSumOfWeightsu[iPOI][iSpace][i]->SetYTitle(sWeights[i].Data());
212 fHistSumOfWeightsu[iPOI][iSpace][i]->SetXTitle(sTitle[iSpace].Data());
213 errorRelated->Add(fHistSumOfWeightsu[iPOI][iSpace][i]);
214 }
12dc476e 215 }
29195b69 216 //weights
217 if(fUsePhiWeights) {
218 if(!fWeightsList) {
02914cfa 219 printf( "WARNING: fWeightsList is NULL in the Scalar Product method.\n" );
29195b69 220 exit(0);
221 }
02914cfa 222 fPhiWeightsSub[0] = dynamic_cast<TH1F*>
223 (fWeightsList->FindObject("phi_weights_sub0"));
224 if(!fPhiWeightsSub[0]) {
225 printf( "WARNING: phi_weights_sub0 not found in the Scalar Product method.\n" );
226 exit(0);
04c6b875 227 }
02914cfa 228 nuaRelated->Add( fPhiWeightsSub[0] );
229 fPhiWeightsSub[1] = dynamic_cast<TH1F*>
230 (fWeightsList->FindObject("phi_weights_sub1"));
231 if(!fPhiWeightsSub[1]) {
232 printf( "WARNING: phi_weights_sub1 not found in the Scalar Product method.\n" );
233 exit(0);
234 }
235 nuaRelated->Add( fPhiWeightsSub[1] );
236 }
04c6b875 237
e2d51347 238
02914cfa 239 fHistProQNorm = new TProfile("FlowPro_QNorm_SP","FlowPro_QNorm_SP", 1,0.5,1.5,"s");
240 fHistProQNorm->SetYTitle("<|Qa+Qb|>");
9812922b 241 tQARelated->Add(fHistProQNorm);
8d312f00 242
02914cfa 243 fHistProQaQb = new TProfile("FlowPro_QaQb_SP","FlowPro_QaQb_SP", 1,0.5,1.5,"s");
244 fHistProQaQb->SetYTitle("<QaQb>");
9812922b 245 tQARelated->Add(fHistProQaQb);
02914cfa 246
247 fHistProQaQbM = new TProfile("FlowPro_QaQbvsM_SP","FlowPro_QaQbvsM_SP",1000,0.0,10000);
248 fHistProQaQbM->SetYTitle("<QaQb>");
249 fHistProQaQbM->SetXTitle("M");
250 fHistProQaQbM->Sumw2();
9812922b 251 tQARelated->Add(fHistProQaQbM);
02914cfa 252
253 fHistMaMb = new TH2D("Flow_MavsMb_SP","Flow_MavsMb_SP",100,0.,100.,100,0.,100.);
254 fHistMaMb->SetYTitle("Ma");
255 fHistMaMb->SetXTitle("Mb");
9812922b 256 tQARelated->Add(fHistMaMb);
02914cfa 257
258 fHistQNormQaQbNorm = new TH2D("Flow_QNormvsQaQbNorm_SP","Flow_QNormvsQaQbNorm_SP",88,-1.1,1.1,22,0.,1.1);
259 fHistQNormQaQbNorm->SetYTitle("|Q/Mq|");
260 fHistQNormQaQbNorm->SetXTitle("QaQb/MaMb");
9812922b 261 tQARelated->Add(fHistQNormQaQbNorm);
02914cfa 262
263 fHistQaNormMa = new TH2D("Flow_QaNormvsMa_SP","Flow_QaNormvsMa_SP",100,0.,100.,22,0.,1.1);
264 fHistQaNormMa->SetYTitle("|Qa/Ma|");
265 fHistQaNormMa->SetXTitle("Ma");
9812922b 266 tQARelated->Add(fHistQaNormMa);
02914cfa 267
268 fHistQbNormMb = new TH2D("Flow_QbNormvsMb_SP","Flow_QbNormvsMb_SP",100,0.,100.,22,0.,1.1);
269 fHistQbNormMb->SetYTitle("|Qb/Mb|");
270 fHistQbNormMb->SetXTitle("Mb");
9812922b 271 tQARelated->Add(fHistQbNormMb);
02914cfa 272
273 fResolution = new TH1D("Flow_resolution_SP","Flow_resolution_SP",100,-1.0,1.0);
274 fResolution->SetYTitle("dN/d(Cos2(#phi_a - #phi_b))");
275 fResolution->SetXTitle("Cos2(#phi_a - #phi_b)");
9812922b 276 tQARelated->Add(fResolution);
1b302085 277
02914cfa 278 fHistQaQb = new TH1D("Flow_QaQb_SP","Flow_QaQb_SP",200,-100.,100.);
279 fHistQaQb->SetYTitle("dN/dQaQb");
280 fHistQaQb->SetXTitle("dQaQb");
9812922b 281 tQARelated->Add(fHistQaQb);
1b302085 282
02914cfa 283 fHistQaQbCos = new TH1D("Flow_QaQbCos_SP","Flow_QaQbCos_SP",63,0.,TMath::Pi());
284 fHistQaQbCos->SetYTitle("dN/d#phi");
285 fHistQaQbCos->SetXTitle("#phi");
9812922b 286 tQARelated->Add(fHistQaQbCos);
1b302085 287
02914cfa 288 fHistList->Add(uQRelated);
289 fHistList->Add(nuaRelated);
290 fHistList->Add(errorRelated);
9812922b 291 fHistList->Add(tQARelated);
1b302085 292
02914cfa 293 TH1::AddDirectory(oldHistAddStatus);
1b302085 294}
295
296//-----------------------------------------------------------------------
02914cfa 297void AliFlowAnalysisWithScalarProduct::Make(AliFlowEventSimple* anEvent) {
298 // Scalar Product method
299 if (!anEvent) return; // for coverity
300
301 // Get Q vectors for the subevents
302 AliFlowVector* vQarray = new AliFlowVector[2];
303 if (fUsePhiWeights)
304 anEvent->Get2Qsub(vQarray,fHarmonic,fWeightsList,kTRUE);
305 else
306 anEvent->Get2Qsub(vQarray,fHarmonic);
307 // Subevent a
308 AliFlowVector vQa = vQarray[0];
309 // Subevent b
310 AliFlowVector vQb = vQarray[1];
311 delete [] vQarray;
312
313 Double_t dMa = vQa.GetMult();
314 if( dMa < 2 ) return;
315 Double_t dMb = vQb.GetMult();
316 if( dMb < 2 ) return;
317 //fill control histograms
318 if (fUsePhiWeights) {
319 fCommonHists->FillControlHistograms(anEvent,fWeightsList,kTRUE);
320 } else {
321 fCommonHists->FillControlHistograms(anEvent);
322 }
1b302085 323
02914cfa 324 //Normalizing: weight the Q vectors for the subevents
325 Double_t dNa = fNormalizationType ? dMa: vQa.Mod(); // SP corresponds to true
326 Double_t dNb = fNormalizationType ? dMb: vQb.Mod(); // SP corresponds to true
327 Double_t dWa = fNormalizationType ? dMa: 1; // SP corresponds to true
328 Double_t dWb = fNormalizationType ? dMb: 1; // SP corresponds to true
12dc476e 329
02914cfa 330 //Scalar product of the two subevents vectors
331 Double_t dQaQb = (vQa*vQb);
8d312f00 332
02914cfa 333 //printf("==============\n");
334 //printf("vQa: { %f, %f }\n",vQa.X(),vQa.Y());
335 //printf("QaQb/|Qa||Qb|: %f\n",dQaQb/vQa.Mod()/vQb.Mod());
336 //printf("QaQb/|Ma||Mb|: %f\n",dQaQb/dMa/dMb);
337
338 // 01 10 11 <=== fTotalQVector
339 // Q ?= Qa or Qb or QaQb
340 AliFlowVector vQm;
341 vQm.Set(0.0,0.0);
342 Double_t dNq=0;
343 if( (fTotalQvector%2)>0 ) { // 01 or 11
344 vQm += vQa;
345 dNq += dMa;
346 }
347 if( fTotalQvector>1 ) { // 10 or 11
348 vQm += vQb;
349 dNq += dMb;
350 }
351 Double_t dWq = fNormalizationType ? dNq: 1; // SP corresponds to true
352 dNq = fNormalizationType ? dNq: vQm.Mod(); // SP corresponds to true
353
354 //fill some EP control histograms
355 fHistProQaQbNorm->Fill(1.,dQaQb/dNa/dNb,dWa*dWb); //Fill (QaQb/NaNb) with weight (WaWb).
356 //needed for the error calculation:
357 fHistSumOfWeights -> Fill(1.,dWa*dWb);
358 fHistSumOfWeights -> Fill(2.,pow(dWa*dWb,2.));
359 //needed for correcting non-uniform acceptance:
360 fHistProNUAq->Fill(1.,vQa.Y()/dNa,dWa); // to get <<sin(phi_a)>>
361 fHistProNUAq->Fill(2.,vQa.X()/dNa,dWa); // to get <<cos(phi_a)>>
362 fHistProNUAq->Fill(3.,vQb.Y()/dNb,dWb); // to get <<sin(phi_b)>>
363 fHistProNUAq->Fill(4.,vQb.X()/dNb,dWb); // to get <<cos(phi_b)>>
364 fHistProNUAq->Fill(5.,vQm.Y()/dNq,dWq);
365 fHistProNUAq->Fill(6.,vQm.X()/dNq,dWq);
366
367 //loop over the tracks of the event
368 AliFlowTrackSimple* pTrack = NULL;
369 Int_t iNumberOfTracks = anEvent->NumberOfTracks();
370 for (Int_t i=0;i<iNumberOfTracks;i++) {
371 pTrack = anEvent->GetTrack(i) ;
372 if (!pTrack) continue;
373 Double_t dPhi = pTrack->Phi();
374 Double_t dPt = pTrack->Pt();
375 Double_t dEta = pTrack->Eta();
376
377 //calculate vU
378 TVector2 vU;
379 Double_t dUX = TMath::Cos(fHarmonic*dPhi);
380 Double_t dUY = TMath::Sin(fHarmonic*dPhi);
381 vU.Set(dUX,dUY);
382
383 // 01 10 11 <=== fTotalQVector
384 // Q ?= Qa or Qb or QaQb
385 vQm.Set(0.0,0.0); //start the loop fresh
386 Double_t dMq=0;
387 if( (fTotalQvector%2)>0 ) { // 01 or 11
388 vQm += vQa;
389 dMq += dMa;
29195b69 390 }
02914cfa 391 if( fTotalQvector>1 ) { // 10 or 11
392 vQm += vQb;
393 dMq += dMb;
1b302085 394 }
1b302085 395
02914cfa 396 //remove track if in subevent
397 for(Int_t inSubEvent=0; inSubEvent!=2; ++inSubEvent) {
398 if( !pTrack->InSubevent( inSubEvent ) )
399 continue;
400 if(inSubEvent==0)
401 if( (fTotalQvector%2)!=1 )
402 continue;
403 if(inSubEvent==1)
404 if( fTotalQvector>0 )
405 continue;
406 Double_t dW=1;
407 Double_t dPhiCenter = dPhi;
408 //if phi weights are used
409 if(fUsePhiWeights && fPhiWeightsSub[inSubEvent]) {
410 Int_t iNBinsPhiSub = fPhiWeightsSub[inSubEvent]->GetNbinsX();
411 Int_t phiBin = 1+(Int_t)(TMath::Floor(dPhi*iNBinsPhiSub/TMath::TwoPi()));
412 dW = fPhiWeightsSub[inSubEvent]->GetBinContent(phiBin);
413 dPhiCenter = fPhiWeightsSub[inSubEvent]->GetBinCenter(phiBin);
1b302085 414 }
02914cfa 415 Double_t dQmX = vQm.X() - dW*(pTrack->Weight())* TMath::Cos(fHarmonic*dPhiCenter);
416 Double_t dQmY = vQm.Y() - dW*(pTrack->Weight())* TMath::Sin(fHarmonic*dPhiCenter);
417 vQm.Set(dQmX,dQmY);
418 dMq = dMq-dW*pTrack->Weight();
419 }
420 dNq = fNormalizationType ? dMq : vQm.Mod();
421 dWq = fNormalizationType ? dMq : 1;
422
423 //Filling QA (for compatibility with previous version)
424 fHistProQaQb->Fill(1,vQa*vQb,1);
425 fHistProQaQbM->Fill(anEvent->GetNumberOfRPs()+0.5,(vQa*vQb)/dMa/dMb,dMa*dMb);
426 fHistQaQbCos->Fill(TMath::ACos((vQa*vQb)/vQa.Mod()/vQb.Mod()));
427 fResolution->Fill( TMath::Cos( vQa.Phi()-vQb.Phi() ) );
428 fHistQaQb->Fill(vQa*vQb);
429 fHistMaMb->Fill(dMb,dMa);
430 fHistProQNorm->Fill(1,vQm.Mod()/dMq,dMq);
431 fHistQNormQaQbNorm->Fill((vQa*vQb)/dMa/dMb,vQm.Mod()/dMq);
432 fHistQaNormMa->Fill(dMa,vQa.Mod()/dMa);
433 fHistQbNormMb->Fill(dMb,vQb.Mod()/dMb);
434
435 Double_t dUQ = vU*vQm;
436
437 //fill the profile histograms
438 for(Int_t iPOI=0; iPOI!=2; ++iPOI) {
439 if( (iPOI==0)&&(!pTrack->InRPSelection()) )
440 continue;
441 if( (iPOI==1)&&(!pTrack->InPOISelection()) )
442 continue;
443 fHistProUQ[iPOI][0]->Fill(dPt ,dUQ/dNq,dWq); //Fill (uQ/Nq') with weight (Nq')
444 fHistProUQ[iPOI][1]->Fill(dEta,dUQ/dNq,dWq); //Fill (uQ/Nq') with weight (Nq')
445 //needed for the error calculation:
446 fHistProUQQaQb[iPOI][0]-> Fill(dPt ,dUQ/dNq*dQaQb/dNa/dNb,dWq*dWa*dWb); //Fill [Qu/Nq']*[QaQb/NaNb] with weight (Nq')NaNb
447 fHistProUQQaQb[iPOI][1]-> Fill(dEta,dUQ/dNq*dQaQb/dNa/dNb,dWq*dWa*dWb); //Fill [Qu/Nq']*[QaQb/NaNb] with weight (Nq')NaNb
448 fHistSumOfWeightsu[iPOI][0][0]->Fill(dPt ,dWq); // sum of Nq'
449 fHistSumOfWeightsu[iPOI][0][1]->Fill(dPt ,pow(dWq,2.));// sum of Nq'^2
450 fHistSumOfWeightsu[iPOI][0][2]->Fill(dPt ,dWq*dWa*dWb);// sum of Nq'*Na*Nb
451 fHistSumOfWeightsu[iPOI][1][0]->Fill(dEta,dWq); // sum of Nq'
452 fHistSumOfWeightsu[iPOI][1][1]->Fill(dEta,pow(dWq,2.));// sum of Nq'^2
453 fHistSumOfWeightsu[iPOI][1][2]->Fill(dEta,dNq*dWa*dWb);// sum of Nq'*Na*Nb
454 //NUA:
455 fHistProNUAu[iPOI][0][0]->Fill(dPt,dUY,1.); //sin u
456 fHistProNUAu[iPOI][0][1]->Fill(dPt,dUX,1.); //cos u
457 fHistProNUAu[iPOI][1][0]->Fill(dEta,dUY,1.); //sin u
458 fHistProNUAu[iPOI][1][1]->Fill(dEta,dUX,1.); //cos u
459 }
460 }//loop over tracks
1b302085 461
02914cfa 462}
8d312f00 463
12dc476e 464//--------------------------------------------------------------------
fd46c3dd 465void AliFlowAnalysisWithScalarProduct::GetOutputHistograms(TList *outputListHistos){
fd46c3dd 466 //get pointers to all output histograms (called before Finish())
02914cfa 467 fHistList = outputListHistos;
468
469 fCommonHists = (AliFlowCommonHist*) fHistList->FindObject("AliFlowCommonHist_SP");
470 fCommonHistsuQ = (AliFlowCommonHist*) fHistList->FindObject("AliFlowCommonHist_uQ");
471 fCommonHistsRes = (AliFlowCommonHistResults*) fHistList->FindObject("AliFlowCommonHistResults_SP");
472 fHistProConfig = (TProfile*) fHistList->FindObject("FlowPro_Flags_SP");
473 if(!fHistProConfig) printf("Error loading fHistProConfig\n");
474 TList *uQ = (TList*) fHistList->FindObject("uQ");
475 TList *nua = (TList*) fHistList->FindObject("NUA");
476 TList *error = (TList*) fHistList->FindObject("error");
477
478 fHistProQaQbNorm = (TProfile*) error->FindObject("FlowPro_QaQbNorm_SP");
479 if(!fHistProQaQbNorm) printf("Error loading fHistProQaQbNorm\n");
ad1b7114 480 fHistProNUAq = (TProfile*) nua->FindObject("FlowPro_NUAq_SP");
02914cfa 481 if(!fHistProNUAq) printf("Error loading fHistProNUAq\n");
482 fHistSumOfWeights = (TH1D*) error->FindObject("Flow_SumOfWeights_SP");
483 if(!fHistSumOfWeights) printf("Error loading fHistSumOfWeights\n");
484
485 TString sPOI[2] = {"RP","POI"};
486 TString sEta[2] = {"Pt","eta"};
487 TString sWeights[3] = {"uQ","uQuQ","uQQaQb"};
02914cfa 488 for(Int_t iPOI=0; iPOI!=2; ++iPOI) for(Int_t iSpace=0; iSpace!=2; ++iSpace) {
489 fHistProUQ[iPOI][iSpace] = (TProfile*) uQ->FindObject( Form( "FlowPro_UQ_%s%s_SP", sEta[iSpace].Data(), sPOI[iPOI].Data() ) );
490 if(!fHistProUQ[iPOI][iSpace]) printf("Error loading fHistProUQ[%d][%d]\n",iPOI,iSpace);
491 fHistProNUAu[iPOI][iSpace][0] = (TProfile*) nua->FindObject( Form("FlowPro_NUAu_%s%s_IM_SP", sEta[iSpace].Data(), sPOI[iPOI].Data() ) );
492 if(!fHistProNUAu[iPOI][iSpace][0]) printf("Error loading fHistProNUAu[%d][%d][0]\n",iPOI,iSpace);
493 fHistProNUAu[iPOI][iSpace][1] = (TProfile*) nua->FindObject( Form("FlowPro_NUAu_%s%s_RE_SP", sEta[iSpace].Data(), sPOI[iPOI].Data() ) );
494 if(!fHistProNUAu[iPOI][iSpace][1]) printf("Error loading fHistProNUAu[%d][%d][1]\n",iPOI,iSpace);
495 fHistProUQQaQb[iPOI][iSpace] = (TProfile*) error->FindObject( Form("FlowPro_UQQaQb_%s%s_SP", sEta[iSpace].Data(), sPOI[iPOI].Data() ) );
496 for(Int_t i=0; i!=3; ++i){
497 fHistSumOfWeightsu[iPOI][iSpace][i] = (TH1D*) error->FindObject( Form("Flow_SumOfWeights_%s%s_%s_SP",sWeights[i].Data(),sPOI[iPOI].Data(),sEta[iSpace].Data()) );
498 if(!fHistSumOfWeightsu[iPOI][iSpace][i]) printf("Error loading fHistSumOfWeightsu[%d][%d][%d]\n",iPOI,iSpace,i);
499 }
500 }
0006c329 501 if(fHistProConfig) {
502 fApplyCorrectionForNUA = (Int_t) fHistProConfig->GetBinContent(1);
503 fNormalizationType = (Int_t) fHistProConfig->GetBinContent(2);
504 fUsePhiWeights = (Int_t) fHistProConfig->GetBinContent(3);
505 fHarmonic = (Int_t) fHistProConfig->GetBinContent(4);
506 }
fd46c3dd 507}
508
509//--------------------------------------------------------------------
8d312f00 510void AliFlowAnalysisWithScalarProduct::Finish() {
489d5531 511 //calculate flow and fill the AliFlowCommonHistResults
02914cfa 512 printf("AliFlowAnalysisWithScalarProduct::Finish()\n");
29ecccee 513
50fe33aa 514 // access harmonic:
9812922b 515 fApplyCorrectionForNUA = fHistProConfig->GetBinContent(1);
516 fNormalizationType = fHistProConfig->GetBinContent(2);
517 fHarmonic = fHistProConfig->GetBinContent(4);
50fe33aa 518
02914cfa 519 printf("*************************************\n");
520 printf("*************************************\n");
521 printf(" Integrated flow from \n");
ad1b7114 522 printf(" Scalar Product \n\n");
523 if(!fNormalizationType)
524 printf(" (BehaveAsEP) \n\n");
e4cecfa0 525
02914cfa 526 //Calculate reference flow
12dc476e 527 //----------------------------------
02914cfa 528 //weighted average over (QaQb/NaNb) with weight (NaNb)
529 Double_t dEntriesQaQb = fHistProQaQbNorm->GetEntries();
530 if( dEntriesQaQb < 1 )
531 return;
a554b203 532 Double_t dQaQb = fHistProQaQbNorm->GetBinContent(1);
02914cfa 533 if( dQaQb < 0 )
534 return;
a93de0f0 535 Double_t dSpreadQaQb = fHistProQaQbNorm->GetBinError(1);
02914cfa 536
537 //NUA qq:
538 Double_t dImQa = fHistProNUAq->GetBinContent(1);
539 Double_t dReQa = fHistProNUAq->GetBinContent(2);
540 Double_t dImQb = fHistProNUAq->GetBinContent(3);
541 Double_t dReQb = fHistProNUAq->GetBinContent(4);
29ecccee 542 if(fApplyCorrectionForNUA)
02914cfa 543 dQaQb = dQaQb - dImQa*dImQb - dReQa*dReQb;
544 printf("QaQb = %f +- %f\n", dQaQb, (dSpreadQaQb/TMath::Sqrt(dEntriesQaQb)) );
545 Double_t dV = TMath::Sqrt(dQaQb);
546
547 printf("ResSub = %f\n", dV );
548 printf("fTotalQvector %d \n",fTotalQvector);
549 if(!fNormalizationType) {
ad1b7114 550 if(fTotalQvector>2) {
9812922b 551 dV = ComputeResolution( TMath::Sqrt2()*FindXi(dV,1e-6) );
ad1b7114 552 printf("An estimate of the event plane resolution is: %f\n", dV );
553 }
524798bf 554 }
02914cfa 555 printf("ResTot = %f\n", dV );
12dc476e 556 //statistical error of dQaQb:
557 // statistical error = term1 * spread * term2:
558 // term1 = sqrt{sum_{i=1}^{N} w^2}/(sum_{i=1}^{N} w)
559 // term2 = 1/sqrt(1-term1^2)
02914cfa 560 Double_t dSumOfLinearWeights = fHistSumOfWeights->GetBinContent(1);
561 Double_t dSumOfQuadraticWeights = fHistSumOfWeights->GetBinContent(2);
12dc476e 562 Double_t dTerm1 = 0.;
563 Double_t dTerm2 = 0.;
02914cfa 564 if(dSumOfLinearWeights)
12dc476e 565 dTerm1 = pow(dSumOfQuadraticWeights,0.5)/dSumOfLinearWeights;
02914cfa 566 if(1.-pow(dTerm1,2.) > 0.)
12dc476e 567 dTerm2 = 1./pow(1-pow(dTerm1,2.),0.5);
12dc476e 568 Double_t dStatErrorQaQb = dTerm1 * dSpreadQaQb * dTerm2;
12dc476e 569 Double_t dVerr = 0.;
02914cfa 570 if(dQaQb > 0.)
12dc476e 571 dVerr = (1./(2.*pow(dQaQb,0.5)))*dStatErrorQaQb;
02914cfa 572 fCommonHistsRes->FillIntegratedFlow(dV,dVerr);
573 printf("v%d(subevents) = %f +- %f\n",fHarmonic,dV,dVerr);
574
575 Int_t iNbins[2];
576 iNbins[0] = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
577 iNbins[1] = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
578
12dc476e 579 //Calculate differential flow and integrated flow (RP, POI)
580 //---------------------------------------------------------
12dc476e 581 //v as a function of eta for RP selection
02914cfa 582 for(Int_t iRFPorPOI=0; iRFPorPOI != 2; ++iRFPorPOI)
583 for(Int_t iPTorETA=0; iPTorETA != 2; ++iPTorETA)
584 for(Int_t b=1; b != iNbins[iPTorETA]+1; ++b) {
585 Double_t duQpro = fHistProUQ[iRFPorPOI][iPTorETA]->GetBinContent(b);
586 if(fApplyCorrectionForNUA)
a93de0f0 587 duQpro = duQpro
02914cfa 588 - fHistProNUAu[iRFPorPOI][iPTorETA][1]->GetBinContent(b)*fHistProNUAq->GetBinContent(6)
589 - fHistProNUAu[iRFPorPOI][iPTorETA][0]->GetBinContent(b)*fHistProNUAq->GetBinContent(5);
a554b203 590 Double_t dv2pro = -999.;
02914cfa 591 if( TMath::Abs(dV!=0.) ) { dv2pro = duQpro/dV; }
12dc476e 592 //calculate the statistical error
02914cfa 593 Double_t dv2ProErr = CalculateStatisticalError(iRFPorPOI, iPTorETA, b, dStatErrorQaQb);
594 if( (iRFPorPOI==0)&&(iPTorETA==0) )
595 fCommonHistsRes->FillDifferentialFlowPtRP( b, dv2pro, dv2ProErr);
596 if( (iRFPorPOI==0)&&(iPTorETA==1) )
597 fCommonHistsRes->FillDifferentialFlowEtaRP( b, dv2pro, dv2ProErr);
598 if( (iRFPorPOI==1)&&(iPTorETA==0) )
599 fCommonHistsRes->FillDifferentialFlowPtPOI( b, dv2pro, dv2ProErr);
600 if( (iRFPorPOI==1)&&(iPTorETA==1) )
601 fCommonHistsRes->FillDifferentialFlowEtaPOI(b, dv2pro, dv2ProErr);
602 //printf("POI %d | PT %d >>> %f +- %f\n",iRFPorPOI,iPTorETA,dv2pro,dv2ProErr);
80f72ed6 603 }
12dc476e 604
02914cfa 605 printf("\n");
606 printf("*************************************\n");
607 printf("*************************************\n");
12dc476e 608}
8232a5ec 609
02914cfa 610//-----------------------------------------------------------------------
9812922b 611void AliFlowAnalysisWithScalarProduct::WriteHistograms(TDirectoryFile *outputFileName) const {
02914cfa 612 //store the final results in output .root file
613 outputFileName->Add(fHistList);
614 outputFileName->Write(outputFileName->GetName(), TObject::kSingleKey);
615}
8232a5ec 616
12dc476e 617//--------------------------------------------------------------------
9812922b 618Double_t AliFlowAnalysisWithScalarProduct::CalculateStatisticalError(Int_t iRFPorPOI, Int_t iPTorETA, Int_t b, Double_t aStatErrorQaQb) const {
12dc476e 619 //calculate the statistical error for differential flow for bin b
02914cfa 620 Double_t duQproSpread = fHistProUQ[iRFPorPOI][iPTorETA]->GetBinError(b);
621 Double_t sumOfMq = fHistSumOfWeightsu[iRFPorPOI][iPTorETA][0]->GetBinContent(b);
622 Double_t sumOfMqSquared = fHistSumOfWeightsu[iRFPorPOI][iPTorETA][1]->GetBinContent(b);
623 Double_t dQaQb = fHistProQaQbNorm->GetBinContent(1);
624
ccad6081 625 //non-isotropic terms:
02914cfa 626 if(fApplyCorrectionForNUA) {
627 Double_t dImQa = fHistProNUAq->GetBinContent(1); // <<sin(phi_a)>>
628 Double_t dReQa = fHistProNUAq->GetBinContent(2); // <<cos(phi_a)>>
629 Double_t dImQb = fHistProNUAq->GetBinContent(3); // <<sin(phi_b)>>
630 Double_t dReQb = fHistProNUAq->GetBinContent(4); // <<cos(phi_b)>>
631 dQaQb = dQaQb - dImQa*dImQb - dReQa*dReQb;
632 }
633
12dc476e 634 Double_t dTerm1 = 0.;
635 Double_t dTerm2 = 0.;
636 if(sumOfMq) {
637 dTerm1 = (pow(sumOfMqSquared,0.5)/sumOfMq);
638 }
639 if(1.-pow(dTerm1,2.)>0.) {
640 dTerm2 = 1./pow(1.-pow(dTerm1,2.),0.5);
641 }
642 Double_t duQproErr = dTerm1*duQproSpread*dTerm2;
12dc476e 643 // covariances:
02914cfa 644 Double_t dTerm1Cov = fHistSumOfWeightsu[iRFPorPOI][iPTorETA][2]->GetBinContent(b);
645 Double_t dTerm2Cov = fHistSumOfWeights->GetBinContent(1);
12dc476e 646 Double_t dTerm3Cov = sumOfMq;
647 Double_t dWeightedCovariance = 0.;
648 if(dTerm2Cov*dTerm3Cov>0.) {
649 Double_t dDenominator = 1.-dTerm1Cov/(dTerm2Cov*dTerm3Cov);
650 Double_t dPrefactor = dTerm1Cov/(dTerm2Cov*dTerm3Cov);
651 if(dDenominator!=0) {
02914cfa 652 Double_t dCovariance = ( fHistProUQQaQb[iRFPorPOI][iPTorETA]->GetBinContent(b)-dQaQb*fHistProUQ[iRFPorPOI][iPTorETA]->GetBinContent(b))/dDenominator;
12dc476e 653 dWeightedCovariance = dCovariance*dPrefactor;
654 }
655 }
656 Double_t dv2ProErr = 0.; // final statitical error:
a554b203 657 if(dQaQb>0.) {
658 Double_t dv2ProErrorSquared = (1./4.)*pow(dQaQb,-3.)*
02914cfa 659 (pow(fHistProUQ[iRFPorPOI][iPTorETA]->GetBinContent(b),2.)*pow(aStatErrorQaQb,2.)
a554b203 660 + 4.*pow(dQaQb,2.)*pow(duQproErr,2.)
02914cfa 661 - 4.*dQaQb*fHistProUQ[iRFPorPOI][iPTorETA]->GetBinContent(b)*dWeightedCovariance);
12dc476e 662 if(dv2ProErrorSquared>0.) dv2ProErr = pow(dv2ProErrorSquared,0.5);
02914cfa 663 }
12dc476e 664 return dv2ProErr;
665}
29ecccee 666
9812922b 667Double_t AliFlowAnalysisWithScalarProduct::ComputeResolution( Double_t x ) const {
668 // Computes resolution for Event Plane method
02914cfa 669 if(x > 51.3) {
670 printf("Warning: Estimation of total resolution might be WRONG. Please check!");
671 return 0.99981;
672 }
673 Double_t a = x*x/4;
674 Double_t b = TMath::Exp(-a)*TMath::BesselI0(a)+TMath::Exp(-a)*TMath::BesselI1(a);
675 return TMath::Sqrt(TMath::PiOver2())/2*x*b;
676}
29ecccee 677
9812922b 678Double_t AliFlowAnalysisWithScalarProduct::FindXi( Double_t res, Double_t prec ) const {
679 // Computes x(res) for Event Plane method
02914cfa 680 if(res > 0.99981) {
681 printf("Warning: Resolution for subEvent is high. You reached the precision limit.");
682 return 51.3;
683 }
684 int nSteps =0;
685 Double_t xtmp=0, xmin=0, xmax=51.3, rtmp=0, delta=2*prec;
686 while( delta > prec ) {
687 xtmp = 0.5*(xmin+xmax);
9812922b 688 rtmp = ComputeResolution(xtmp);
02914cfa 689 delta = TMath::Abs( res-rtmp );
690 if(rtmp>res) xmax = xtmp;
691 if(rtmp<res) xmin = xtmp;
692 nSteps++;
693 }
694 return xtmp;
695}
29ecccee 696