AliPhysics  6cf2591 (6cf2591)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliFlowAnalysisWithSimpleSP.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 ALIFLOWANALYSISWITHSIMPLESP_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 
36 using std::cout;
37 using std::endl;
38 
39 
41 
42  //-----------------------------------------------------------------------
44  fDebug(0),
45  fMinimalBook(kFALSE),
46  fUsePhiWeights(0),
47  fApplyCorrectionForNUA(0),
48  fHarmonic(2),
49  fWeights(kFALSE),
50  fScaling(kTRUE),
51  fNormalizationType(1),
52  fV0SanityCheck(0),
53  fExternalResolution(-999.),
54  fExternalResErr(0.),
55  fTotalQvector(3),
56  fWeightsList(NULL),
57  fHistList(NULL),
58  fHistProConfig(NULL),
59  fHistProQaQbNorm(NULL),
60  fHistSumOfWeights(NULL),
61  fHistProNUAq(NULL),
62  fCentralityWeight(1.),
63  fHistProQNorm(NULL),
64  fHistProQaQb(NULL),
65  fHistProQaQbM(NULL),
66  fHistMaMb(NULL),
67  fHistQNormQaQbNorm(NULL),
68  fHistQaNormMa(NULL),
69  fHistQbNormMb(NULL),
70  fResolution(NULL),
71  fHistQaQb(NULL),
72  fHistQaQc(NULL),
73  fHistQbQc(NULL),
74  fHistQaQbCos(NULL),
75  fHistNumberOfSubtractedDaughters(NULL),
76  fCommonHists(NULL),
77  fCommonHistsuQ(NULL),
78  fCommonHistsRes(NULL)
79 {
80  //ctor
81  for(int i=0; i!=2; ++i) {
82  fPhiWeightsSub[i] = NULL;
83  for(int j=0; j!=2; ++j) {
84  fHistProUQ[i][j] = NULL;
85  fHistProUQQaQb[i][j] = NULL;
86  fHistProNUAu[i][j][0] = NULL;
87  fHistProNUAu[i][j][1] = NULL;
88  for(int k=0; k!=3; ++k)
89  fHistSumOfWeightsu[i][j][k] = NULL;
90  }
91  }
92 }
93 //-----------------------------------------------------------------------
95 {
96  //destructor
97  delete fWeightsList;
98  delete fHistList;
99 }
100 //-----------------------------------------------------------------------
102  //Define all histograms
103  //printf("---Analysis with the Scalar Product Method--- Init\n");
104  //printf("--- fNormalizationType %d ---\n", fNormalizationType);
105  //save old value and prevent histograms from being added to directory
106  //to avoid name clashes in case multiple analaysis objects are used
107  //in an analysis
108  Bool_t oldHistAddStatus = TH1::AddDirectoryStatus();
109  TH1::AddDirectory(kFALSE);
110 
111  fHistList = new TList();
112  fHistList->SetName("cobjSP");
113  fHistList->SetOwner();
114 
115  TList *uQRelated = new TList();
116  uQRelated->SetName("uQ");
117  uQRelated->SetOwner();
118 
119  TList *nuaRelated = new TList();
120  nuaRelated->SetName("NUA");
121  nuaRelated->SetOwner();
122 
123  TList *errorRelated = new TList();
124  errorRelated->SetName("error");
125  errorRelated->SetOwner();
126 
127  TList *tQARelated = new TList();
128  tQARelated->SetName("QA");
129  tQARelated->SetOwner();
130 
131  fCommonHists = new AliFlowCommonHist("AliFlowCommonHist_SP","AliFlowCommonHist",fMinimalBook);
132  (fCommonHists->GetHarmonic())->Fill(0.5,fHarmonic); // store harmonic
133  fHistList->Add(fCommonHists);
134 
135  // fCommonHistsuQ = new AliFlowCommonHist("AliFlowCommonHist_uQ");
136  // (fCommonHistsuQ->GetHarmonic())->Fill(0.5,fHarmonic); // store harmonic
137  // fHistList->Add(fCommonHistsuQ);
138 
139  fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResults_SP","",fHarmonic);
141 
142  fHistProConfig = new TProfile("FlowPro_Flags_SP","Flow_Flags_SP",4,0.5,4.5,"s");
143  fHistProConfig->GetXaxis()->SetBinLabel(1,"fApplyCorrectionForNUA");
144  fHistProConfig->GetXaxis()->SetBinLabel(2,"fNormalizationType");
145  fHistProConfig->GetXaxis()->SetBinLabel(3,"fUsePhiWeights");
146  fHistProConfig->GetXaxis()->SetBinLabel(4,"fHarmonic");
150  fHistProConfig->Fill(4,fHarmonic);
152 
153  fHistProQaQbNorm = new TProfile("FlowPro_QaQbNorm_SP","FlowPro_QaQbNorm_SP", 10000, -1000, 1000);
154  fHistProQaQbNorm->SetYTitle("<QaQb/Na/Nb>");
155  errorRelated->Add(fHistProQaQbNorm);
156 
157  fHistProNUAq = new TProfile("FlowPro_NUAq_SP","FlowPro_NUAq_SP", 6, 0.5, 6.5,"s");
158  fHistProNUAq->GetXaxis()->SetBinLabel( 1,"<<sin(#Phi_{a})>>");
159  fHistProNUAq->GetXaxis()->SetBinLabel( 2,"<<cos(#Phi_{a})>>");
160  fHistProNUAq->GetXaxis()->SetBinLabel( 3,"<<sin(#Phi_{b})>>");
161  fHistProNUAq->GetXaxis()->SetBinLabel( 4,"<<cos(#Phi_{b})>>");
162  fHistProNUAq->GetXaxis()->SetBinLabel( 5,"<<sin(#Phi_{t})>>");
163  fHistProNUAq->GetXaxis()->SetBinLabel( 6,"<<cos(#Phi_{t})>>");
164  nuaRelated->Add(fHistProNUAq);
165 
166  fHistSumOfWeights = new TH1D("Flow_SumOfWeights_SP","Flow_SumOfWeights_SP",2,0.5, 2.5);
167  fHistSumOfWeights->GetXaxis()->SetBinLabel( 1,"#Sigma Na*Nb");
168  fHistSumOfWeights->GetXaxis()->SetBinLabel( 2,"#Sigma (Na*Nb)^2");
169  errorRelated->Add(fHistSumOfWeights);
170 
171  TString sPOI[2] = {"RP","POI"}; // backward compatibility
172  TString sEta[2] = {"Pt","eta"}; // backward compatibility
173  TString sTitle[2] = {"p_{T} [GeV]","#eta"};
174  TString sWeights[3] = {"uQ","uQuQ","uQQaQb"};
175  Int_t iNbins[2];
176  Double_t dMin[2], dMax[2];
183  for(Int_t iPOI=0; iPOI!=2; ++iPOI)
184  for(Int_t iSpace=0; iSpace!=2; ++iSpace) {
185  // uQ
186  fHistProUQ[iPOI][iSpace] = new TProfile( Form( "FlowPro_UQ_%s%s_SP", sEta[iSpace].Data(), sPOI[iPOI].Data() ),
187  Form( "FlowPro_UQ%s%s_SP", sEta[iSpace].Data(), sPOI[iPOI].Data() ),
188  iNbins[iSpace], dMin[iSpace], dMax[iSpace], "s");
189  fHistProUQ[iPOI][iSpace]->SetXTitle(sTitle[iSpace].Data());
190  fHistProUQ[iPOI][iSpace]->SetYTitle("<uQ>");
191  uQRelated->Add(fHistProUQ[iPOI][iSpace]);
192 
193  // NUAu
194  fHistProNUAu[iPOI][iSpace][0] = new TProfile( Form("FlowPro_NUAu_%s%s_IM_SP", sEta[iSpace].Data(), sPOI[iPOI].Data() ),
195  Form("FlowPro_NUAu_%s%s_IM_SP", sEta[iSpace].Data(), sPOI[iPOI].Data() ),
196  iNbins[iSpace], dMin[iSpace], dMax[iSpace]);
197  fHistProNUAu[iPOI][iSpace][0]->SetXTitle(sTitle[iSpace].Data());
198  nuaRelated->Add(fHistProNUAu[iPOI][iSpace][0]);
199  fHistProNUAu[iPOI][iSpace][1] = new TProfile( Form("FlowPro_NUAu_%s%s_RE_SP", sEta[iSpace].Data(), sPOI[iPOI].Data() ),
200  Form("FlowPro_NUAu_%s%s_RE_SP", sEta[iSpace].Data(), sPOI[iPOI].Data() ),
201  iNbins[iSpace], dMin[iSpace], dMax[iSpace]);
202  fHistProNUAu[iPOI][iSpace][1]->SetXTitle(sTitle[iSpace].Data());
203  nuaRelated->Add(fHistProNUAu[iPOI][iSpace][1]);
204 
205  // uQ QaQb
206  fHistProUQQaQb[iPOI][iSpace] = new TProfile( Form("FlowPro_UQQaQb_%s%s_SP", sEta[iSpace].Data(), sPOI[iPOI].Data() ),
207  Form("FlowPro_UQQaQb_%s%s_SP", sEta[iSpace].Data(), sPOI[iPOI].Data() ),
208  iNbins[iSpace], dMin[iSpace], dMax[iSpace]);
209  fHistProUQQaQb[iPOI][iSpace]->SetXTitle(sTitle[iSpace].Data());
210  fHistProUQQaQb[iPOI][iSpace]-> SetYTitle("<Qu QaQb>");
211  errorRelated->Add(fHistProUQQaQb[iPOI][iSpace]);
212 
213  // uWeights
214  for(Int_t i=0; i!=3; ++i) {
215  fHistSumOfWeightsu[iPOI][iSpace][i] = new TH1D( Form("Flow_SumOfWeights_%s%s_%s_SP",sWeights[i].Data(),sPOI[iPOI].Data(),sEta[iSpace].Data()),
216  Form("Flow_SumOfWeights_%s%s_%s_SP",sWeights[i].Data(),sPOI[iPOI].Data(),sEta[iSpace].Data()),
217  iNbins[iSpace], dMin[iSpace], dMax[iSpace]);
218  fHistSumOfWeightsu[iPOI][iSpace][i]->SetYTitle(sWeights[i].Data());
219  fHistSumOfWeightsu[iPOI][iSpace][i]->SetXTitle(sTitle[iSpace].Data());
220  errorRelated->Add(fHistSumOfWeightsu[iPOI][iSpace][i]);
221  }
222  }
223  //weights
224  if(fUsePhiWeights) {
225  if(!fWeightsList) {
226  printf( "WARNING: fWeightsList is NULL in the Scalar Product method.\n" );
227  exit(0);
228  }
229  fPhiWeightsSub[0] = dynamic_cast<TH1F*>
230  (fWeightsList->FindObject("phi_weights_sub0"));
231  if(!fPhiWeightsSub[0]) {
232  printf( "WARNING: phi_weights_sub0 not found in the Scalar Product method.\n" );
233  exit(0);
234  }
235  nuaRelated->Add( fPhiWeightsSub[0] );
236  fPhiWeightsSub[1] = dynamic_cast<TH1F*>
237  (fWeightsList->FindObject("phi_weights_sub1"));
238  if(!fPhiWeightsSub[1]) {
239  printf( "WARNING: phi_weights_sub1 not found in the Scalar Product method.\n" );
240  exit(0);
241  }
242  nuaRelated->Add( fPhiWeightsSub[1] );
243  }
244 
245  if(!fMinimalBook) {
246  fHistProQNorm = new TProfile("FlowPro_QNorm_SP","FlowPro_QNorm_SP", 1,0.5,1.5,"s");
247  fHistProQNorm->SetYTitle("<|Qa+Qb|>");
248  tQARelated->Add(fHistProQNorm);
249 
250  fHistProQaQb = new TH1D("FlowPro_QaQb_SP","FlowPro_QaQb_SP", 10000,-100,100);
251  fHistProQaQb->SetYTitle("<QaQb>");
252  fHistProQaQb->StatOverflows(kTRUE);
253  tQARelated->Add(fHistProQaQb);
254 
255  fHistProQaQbM = new TH1D("FlowPro_QaQbvsM_SP","FlowPro_QaQbvsM_SP",1000,0.0,10000);
256  fHistProQaQbM->SetYTitle("<QaQb>");
257  fHistProQaQbM->SetXTitle("M");
258  fHistProQaQbM->Sumw2();
259  tQARelated->Add(fHistProQaQbM);
260 
261  fHistMaMb = new TH2D("Flow_MavsMb_SP","Flow_MavsMb_SP",100,0.,100.,100,0.,100.);
262  fHistMaMb->SetYTitle("Ma");
263  fHistMaMb->SetXTitle("Mb");
264  tQARelated->Add(fHistMaMb);
265 
266  fHistQNormQaQbNorm = new TH2D("Flow_QNormvsQaQbNorm_SP","Flow_QNormvsQaQbNorm_SP",88,-1.1,1.1,22,0.,1.1);
267  fHistQNormQaQbNorm->SetYTitle("|Q/Mq|");
268  fHistQNormQaQbNorm->SetXTitle("QaQb/MaMb");
269  tQARelated->Add(fHistQNormQaQbNorm);
270 
271  fHistQaNormMa = new TH2D("Flow_QaNormvsMa_SP","Flow_QaNormvsMa_SP",100,0.,100.,22,0.,1.1);
272  fHistQaNormMa->SetYTitle("|Qa/Ma|");
273  fHistQaNormMa->SetXTitle("Ma");
274  tQARelated->Add(fHistQaNormMa);
275 
276  fHistQbNormMb = new TH2D("Flow_QbNormvsMb_SP","Flow_QbNormvsMb_SP",100,0.,100.,22,0.,1.1);
277  fHistQbNormMb->SetYTitle("|Qb/Mb|");
278  fHistQbNormMb->SetXTitle("Mb");
279  tQARelated->Add(fHistQbNormMb);
280 
281  fResolution = new TH1D("Flow_resolution_SP","Flow_resolution_SP",100,-1.0,1.0);
282  fResolution->SetYTitle("dN/d(Cos2(#phi_a - #phi_b))");
283  fResolution->SetXTitle("Cos2(#phi_a - #phi_b)");
284  tQARelated->Add(fResolution);
285 
286  fHistQaQb = new TH1D("Flow_QaQb_SP","Flow_QaQb_SP",20000,-1000.,1000.);
287  fHistQaQb->SetYTitle("dN/dQaQb");
288  fHistQaQb->SetXTitle("dQaQb");
289  fHistQaQb->StatOverflows(kTRUE);
290  tQARelated->Add(fHistQaQb);
291 
292  fHistQaQc = new TH1D("Flow_QaQc_SP","Flow_QaQc_SP",20000,-1000.,1000.);
293  fHistQaQc->SetYTitle("dN/dQaQc");
294  fHistQaQc->SetXTitle("dQaQc");
295  fHistQaQc->StatOverflows(kTRUE);
296  tQARelated->Add(fHistQaQc);
297 
298  fHistQbQc = new TH1D("Flow_QbQc_SP","Flow_QbQc_SP",20000,-1000.,1000.);
299  fHistQbQc->SetYTitle("dN/dQbQc");
300  fHistQbQc->SetXTitle("dQbQc");
301  fHistQbQc->StatOverflows(kTRUE);
302  tQARelated->Add(fHistQbQc);
303 
304  fHistQaQbCos = new TH1D("Flow_QaQbCos_SP","Flow_QaQbCos_SP",63,0.,TMath::Pi());
305  fHistQaQbCos->SetYTitle("dN/d#phi");
306  fHistQaQbCos->SetXTitle("#phi");
307  tQARelated->Add(fHistQaQbCos);
308 
309  fHistNumberOfSubtractedDaughters = new TH1I("NumberOfSubtractedDaughtersInQ",";daughters;counts;Number of daughters subtracted from Q",5,0.,5.);
310  tQARelated->Add(fHistNumberOfSubtractedDaughters);
311  }
312 
313  fHistList->Add(uQRelated);
314  fHistList->Add(nuaRelated);
315  fHistList->Add(errorRelated);
316  fHistList->Add(tQARelated);
317 
318  TH1::AddDirectory(oldHistAddStatus);
319 }
320 
321 //-----------------------------------------------------------------------
323  // Scalar Product method
324  if (!anEvent) return; // for coverity
325 
326  fCentralityWeight = anEvent->GetPsi5();
327  // Get Q vectors for the subevents
328  AliFlowVector* vQarray = new AliFlowVector[2];
329  if (fUsePhiWeights)
330  anEvent->Get2Qsub(vQarray,fHarmonic,fWeightsList,kTRUE);
331  else
332  anEvent->Get2Qsub(vQarray,fHarmonic);
333  // Subevent a
334  AliFlowVector vQa = vQarray[0];
335  // Subevent b
336  AliFlowVector vQb = vQarray[1];
337  delete [] vQarray;
338 
339  // multiplicity here corresponds to the V0 equalized multiplicity
340  Double_t dMa = vQa.GetMult();
341  if( dMa < 2 ) return;
342  Double_t dMb = vQb.GetMult();
343  if( dMb < 2 ) return;
344  //fill control histograms
346 
347  //Normalizing: weight the Q vectors for the subevents
348  Double_t dNa = dMa;//fNormalizationType ? dMa: vQa.Mod(); // SP corresponds to true
349  Double_t dNb = dMb;//fNormalizationType ? dMb: vQb.Mod(); // SP corresponds to true
350  Double_t dWa = 1.;//fNormalizationType ? dMa: 1; // SP corresponds to true
351  Double_t dWb = 1.;//NormalizationType ? dMb: 1; // SP corresponds to true
352 
353  // in the original SP method, event weights correspond to multiplicity
354  // for the V0 this does not necessarily make sense as the eq mult
355  // does not scale directly with charged track mult / centrality
356 
357  // this task does not support EP style running, as this is
358  // ambiguous at best
359 
360  //Scalar product of the two subevents vectors
361  Double_t dQaQb = (vQa*vQb);
362 
363  // 01 10 11 <=== fTotalQVector
364  // Q ?= Qa or Qb or QaQb
365  AliFlowVector vQm;
366  vQm.Set(0.0,0.0);
367  AliFlowVector vQTPC;
368  vQTPC.Set(anEvent->GetPsi2(), anEvent->GetPsi3());
369 
370  Double_t dNq=0;
371  if( (fTotalQvector%2)>0 ) { // 01 or 11
372  vQm += vQa;
373  dNq += dMa;
374  }
375  if( fTotalQvector>1 ) { // 10 or 11
376  vQm += vQb;
377  dNq += dMb;
378  }
379  Double_t dWq = 1.;//fNormalizationType ? dNq: 1; // SP corresponds to true
380 
381 
382  // if we dont normalize the q-vectors to the multiplicity, just use 1 here
383  if(!fScaling) {
384  dNa = 1.;
385  dNb = 1.;
386  }
387 
388  // this profile stores the scalar product of the two
389  // subevent q-vectors (the denominator or 'resolution' of the
390  // measurement)
391  // fHistProQaQbNorm->Fill(1., dQaQb/dNa/dNb);
392 
393  //loop over the tracks of the event
394  AliFlowTrackSimple* pTrack = NULL;
395  Int_t iNumberOfTracks = anEvent->NumberOfTracks();
396  Double_t dMq = 0;
397  for (Int_t i=0;i<iNumberOfTracks;i++) {
398  // so this is a track loop ...
399  pTrack = anEvent->GetTrack(i) ;
400  if (!pTrack) continue;
401  Double_t dPhi = pTrack->Phi();
402  Double_t dPt = pTrack->Pt();
403  Double_t dEta = pTrack->Eta();
404 
405  // make sure it's not a vzero track
406  if(TMath::Abs(dEta) > 0.9) continue;
407 
408  //calculate vU
409  // this is the track q-vecotr (small u)
410  TVector2 vU;
411  Double_t dUX = TMath::Cos(fHarmonic*dPhi);
412  Double_t dUY = TMath::Sin(fHarmonic*dPhi);
413  vU.Set(dUX,dUY);
414 
415  // 01 10 11 <=== fTotalQVector
416  // Q ?= Qa or Qb or QaQb
417  // this will be the vector for hte scalar product itself
418  vQm.Set(0.0,0.0); //start the loop fresh
419  dMq=0;
420  if( (fTotalQvector%2)>0 ) { // 01 or 11
421  vQm += vQa;
422  dMq += dMa;
423  }
424  if( fTotalQvector>1 ) { // 10 or 11
425  vQm += vQb;
426  dMq += dMb;
427  }
428 
429  dNq = dMq;//fNormalizationType ? dMq : vQm.Mod();
430  dWq = 1;//fNormalizationType ? dMq : 1;
431 
432  // this little guy is the enumerator of the final equation
433  Double_t dUQ = vU*vQm;
434  // if we dont want scaling, disable it here
435  if(!fScaling) dNq = 1;
436  //fill the profile histograms
437  for(Int_t iPOI=0; iPOI!=2; ++iPOI) {
438  if( (iPOI==0)&&(!pTrack->InRPSelection()) )
439  continue;
440  if( (iPOI==1)&&(!pTrack->InPOISelection(fPOItype)) )
441  continue;
442  fHistProUQ[iPOI][0]->Fill(dPt ,fCentralityWeight*dUQ/dNq); //Fill (uQ/Nq') with weight (Nq')
443  fHistProUQ[iPOI][1]->Fill(dEta,fCentralityWeight*dUQ/dNq); //Fill (uQ/Nq') with weight (Nq')
444  //fHistProUQQaQb[iPOI][0]-> Fill(dPt,(dUQ*dUQ/dNq)/(dQaQb/dNa/dNb)); //Fill [Qu/Nq']*[QaQb/NaNb] with weight (Nq')NaNb
445  }
446  }//loop over tracks
447  //Filling QA (for compatibility with previous version)
448  if(!fMinimalBook) {
449  fHistQaQb->Fill(vQa*vQb/dNa/dNb, fCentralityWeight);
450  double qc = anEvent->GetPsi2();
451 
452  fHistQaQc->Fill(vQa*vQTPC, fCentralityWeight);
453  fHistQbQc->Fill(vQb*vQTPC, fCentralityWeight);
454 
455  }
456 
457 
458 }
459 
460 //--------------------------------------------------------------------
462  //get pointers to all output histograms (called before Finish())
463  fHistList = outputListHistos;
464 
465  fCommonHists = (AliFlowCommonHist*) fHistList->FindObject("AliFlowCommonHist_SP");
466  // fCommonHistsuQ = (AliFlowCommonHist*) fHistList->FindObject("AliFlowCommonHist_uQ");
467  fCommonHistsRes = (AliFlowCommonHistResults*) fHistList->FindObject("AliFlowCommonHistResults_SP");
468  fHistProConfig = (TProfile*) fHistList->FindObject("FlowPro_Flags_SP");
469  if(!fHistProConfig) printf("Error loading fHistProConfig\n");
470  TList *uQ = (TList*) fHistList->FindObject("uQ");
471  TList *nua = (TList*) fHistList->FindObject("NUA");
472  TList *error = (TList*) fHistList->FindObject("error");
473  TList* tQARelated = (TList*) fHistList->FindObject("QA");
474 
475  fHistProQaQbNorm = (TProfile*) error->FindObject("FlowPro_QaQbNorm_SP");
476  if(!fHistProQaQbNorm) printf("Error loading fHistProQaQbNorm\n");
477  fHistProNUAq = (TProfile*) nua->FindObject("FlowPro_NUAq_SP");
478  if(!fHistProNUAq) printf("Error loading fHistProNUAq\n");
479  fHistSumOfWeights = (TH1D*) error->FindObject("Flow_SumOfWeights_SP");
480  if(!fHistSumOfWeights) printf("Error loading fHistSumOfWeights\n");
481 
482  TString sPOI[2] = {"RP","POI"};
483  TString sEta[2] = {"Pt","eta"};
484  TString sWeights[3] = {"uQ","uQuQ","uQQaQb"};
485  for(Int_t iPOI=0; iPOI!=2; ++iPOI) for(Int_t iSpace=0; iSpace!=2; ++iSpace) {
486  fHistProUQ[iPOI][iSpace] = (TProfile*) uQ->FindObject( Form( "FlowPro_UQ_%s%s_SP", sEta[iSpace].Data(), sPOI[iPOI].Data() ) );
487  if(!fHistProUQ[iPOI][iSpace]) printf("Error loading fHistProUQ[%d][%d]\n",iPOI,iSpace);
488  fHistProNUAu[iPOI][iSpace][0] = (TProfile*) nua->FindObject( Form("FlowPro_NUAu_%s%s_IM_SP", sEta[iSpace].Data(), sPOI[iPOI].Data() ) );
489  if(!fHistProNUAu[iPOI][iSpace][0]) printf("Error loading fHistProNUAu[%d][%d][0]\n",iPOI,iSpace);
490  fHistProNUAu[iPOI][iSpace][1] = (TProfile*) nua->FindObject( Form("FlowPro_NUAu_%s%s_RE_SP", sEta[iSpace].Data(), sPOI[iPOI].Data() ) );
491  if(!fHistProNUAu[iPOI][iSpace][1]) printf("Error loading fHistProNUAu[%d][%d][1]\n",iPOI,iSpace);
492  fHistProUQQaQb[iPOI][iSpace] = (TProfile*) error->FindObject( Form("FlowPro_UQQaQb_%s%s_SP", sEta[iSpace].Data(), sPOI[iPOI].Data() ) );
493  for(Int_t i=0; i!=3; ++i){
494  fHistSumOfWeightsu[iPOI][iSpace][i] = (TH1D*) error->FindObject( Form("Flow_SumOfWeights_%s%s_%s_SP",sWeights[i].Data(),sPOI[iPOI].Data(),sEta[iSpace].Data()) );
495  if(!fHistSumOfWeightsu[iPOI][iSpace][i]) printf("Error loading fHistSumOfWeightsu[%d][%d][%d]\n",iPOI,iSpace,i);
496  }
497  }
498  if(fHistProConfig) {
499  fApplyCorrectionForNUA = (Int_t) fHistProConfig->GetBinContent(1);
500  fNormalizationType = (Int_t) fHistProConfig->GetBinContent(2);
501  fUsePhiWeights = (Int_t) fHistProConfig->GetBinContent(3);
502  fHarmonic = (Int_t) fHistProConfig->GetBinContent(4);
503  }
504  fHistQaQb = (TH1D*)tQARelated->FindObject("Flow_QaQb_SP");
505  fHistQaQc = (TH1D*)tQARelated->FindObject("Flow_QaQc_SP");
506  fHistQbQc = (TH1D*)tQARelated->FindObject("Flow_QbQc_SP");
507 
508 
509 }
510 
511 
512 
513 //--------------------------------------------------------------------
515  //calculate flow and fill the AliFlowCommonHistResults
516  printf("AliFlowAnalysisWithSimpleSP::Finish()\n");
517 
518  // access harmonic:
519  fApplyCorrectionForNUA = 0;//(Int_t)(fHistProConfig->GetBinContent(1));
520  fNormalizationType = 1;//(Int_t)(fHistProConfig->GetBinContent(2));
521  fHarmonic = (Int_t)(fHistProConfig->GetBinContent(4));
522 
523  printf("*************************************\n");
524  printf("*************************************\n");
525  printf(" Integrated flow from SIMPLE \n");
526  printf(" Scalar Product \n\n");
527 
528  //Calculate reference flow
529  //----------------------------------
530  //weighted average over (QaQb/NaNb) with weight (NaNb)
531 
532 
533  if(!fHistQaQb) {
534  printf(" PANIC: run with full booking !");
535  return;
536  }
537  Double_t dEntriesQaQb = fHistQaQb->GetEntries();
538  if( dEntriesQaQb < 1 )
539  return;
540  //fHistQaQb->GetXaxis()->SetRangeUser(-10.,10.);
541  //fHistQaQC->GetXaxis()->SetRangeUser(-1000., 1000.);
542  //fHistQbQc->GetXaxis()->SetRangeUser(-1000., 1000.);
543 
544 
545 
546 
547  Double_t dQaQb = fHistQaQb->GetMean();
548  Double_t dQaQc = fHistQaQc->GetMean();
549  Double_t dQbQc = fHistQbQc->GetMean();
550 
551  if(A) {
552  // now it's the resolution of A (which is vzero C)
553  dQaQb = (dQaQb*dQaQc)/dQbQc;
554  } else {
555 
556  // else we select the resolution of B (which is V0A)
557  dQaQb = (dQaQb*dQbQc)/dQaQc;
558  }
559 
560  if( dQaQb < 0 )
561  return;
562  Double_t dSpreadQaQb = fHistQaQb->GetRMS();
563 
564  // this is the `resolution'
565  if(dQaQb <= .0 ) {
566  printf(" Panic! the average of QaQb <= 0! Probably you need to run on more events !\n");
567  printf(" \tusing dummy value 1 to avoid segfault \n");
568  dQaQb = 1.;
569  dSpreadQaQb = 1.;
570  }
571  Double_t dV = TMath::Sqrt(dQaQb);
572 
573  printf("ResSub (sqrt of scalar product of sub-event qvectors) = %f\n", dV );
574  printf("fTotalQvector %d \n",fTotalQvector);
575  // we just take the spread as the uncertainty
576 
577  Double_t dStatErrorQaQb = dSpreadQaQb;
578 
579  Double_t dVerr = 0.;
580  if(dQaQb > 0.) dVerr = (1./(2.*pow(dQaQb,0.5)))*dStatErrorQaQb;
582  printf("v%d(subevents) = %f +- %f\n",fHarmonic,dV,dVerr);
583 
584  Int_t iNbins[2];
587 
588  //Calculate differential flow and integrated flow (RP, POI)
589  //---------------------------------------------------------
590  //v as a function of eta for RP selection
591  for(Int_t iRFPorPOI=0; iRFPorPOI != 2; ++iRFPorPOI)
592  for(Int_t iPTorETA=0; iPTorETA != 2; ++iPTorETA)
593  for(Int_t b=1; b != iNbins[iPTorETA]+1; ++b) {
594  Double_t duQpro = fHistProUQ[iRFPorPOI][iPTorETA]->GetBinContent(b);
595  Double_t dv2pro = -999.;
596  if( TMath::Abs(dV!=0.) && fExternalResolution < 0 ) { dv2pro = duQpro/dV; }
597  else if(fExternalResolution > 0) { dv2pro = duQpro / fExternalResolution; }
598  Double_t dv2ProErr = fHistProUQ[iRFPorPOI][iPTorETA]->GetBinError(b);
599  if(fHistProUQ[iRFPorPOI][iPTorETA]->GetBinEntries(b) > 1) dv2ProErr/=TMath::Sqrt(fHistProUQ[iRFPorPOI][iPTorETA]->GetBinEntries(b));
600  if(fExternalResErr > 0) {
601  // do default error propagation for a fraction where
602  Double_t a(0), b(0), t(0);
603  if(dv2ProErr > 0) a = duQpro*duQpro/(dv2ProErr*dv2ProErr); // squared relative error
605  t = duQpro*fExternalResolution/(dv2ProErr*fExternalResErr);
606  t = dv2pro*dv2pro*(a + b -2.*t);
607  if(t>0) dv2ProErr = TMath::Sqrt(t);
608  }
609 
610  if( (iRFPorPOI==0)&&(iPTorETA==0) )
611  fCommonHistsRes->FillDifferentialFlowPtRP( b, dv2pro, dv2ProErr);
612  if( (iRFPorPOI==0)&&(iPTorETA==1) )
613  fCommonHistsRes->FillDifferentialFlowEtaRP( b, dv2pro, dv2ProErr);
614  if( (iRFPorPOI==1)&&(iPTorETA==0) )
615  fCommonHistsRes->FillDifferentialFlowPtPOI( b, dv2pro, dv2ProErr);
616  if( (iRFPorPOI==1)&&(iPTorETA==1) )
617  fCommonHistsRes->FillDifferentialFlowEtaPOI(b, dv2pro, dv2ProErr);
618  }
619 
620  printf("\n");
621  printf("*************************************\n");
622  printf("*************************************\n");
623 }
624 
625 //-----------------------------------------------------------------------
626 void AliFlowAnalysisWithSimpleSP::WriteHistograms(TDirectoryFile *outputFileName) const {
627  //store the final results in output .root file
628  outputFileName->Add(fHistList);
629  outputFileName->Write(outputFileName->GetName(), TObject::kSingleKey);
630 }
631 
632 //--------------------------------------------------------------------
AliFlowCommonHistResults * fCommonHistsRes
double Double_t
Definition: External.C:58
AliFlowTrackSimple * GetTrack(Int_t i)
Double_t GetMult() const
Definition: AliFlowVector.h:46
virtual void Finish()
Bool_t FillDifferentialFlowPtRP(Int_t aBin, Double_t av, Double_t anError)
void Make(AliFlowEventSimple *anEvent)
Double_t GetPsi3() const
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
Definition: External.C:228
Definition: External.C:212
static AliFlowCommonConstants * GetMaster()
void GetOutputHistograms(TList *outputListHistos)
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)
Double_t GetPsi2() const
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
void WriteHistograms(TDirectoryFile *outputFileName) const
Bool_t InPOISelection(Int_t poiType=1) const
Double_t Eta() const
bool Bool_t
Definition: External.C:53
TProfile * GetHarmonic()
ClassImp(AliFlowAnalysisWithSimpleSP) AliFlowAnalysisWithSimpleSP
Double_t GetPsi5() const
Bool_t FillDifferentialFlowEtaRP(Int_t aBin, Double_t av, Double_t anError)
Int_t NumberOfTracks() const