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