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