Add task to handle the muon type in simulations. The class allows to identify, among...
[u/mrichter/AliRoot.git] / PWG / muon / AliAnalysisTaskFlowSingleMu.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2007, 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 /* $Id: AliAnalysisTaskFlowSingleMu.cxx 55545 2012-04-04 07:16:39Z pcrochet $ */
17
18 //-----------------------------------------------------------------------------
19 /// \class AliAnalysisTaskFlowSingleMu
20 /// Analysis task for flow of single muons in the spectrometer.
21 /// The output is a list of histograms and CF containers.
22 /// The macro class can run on AODs or ESDs.
23 /// If Monte Carlo information is present, some basics checks are performed.
24 ///
25 /// \author Diego Stocco
26 //-----------------------------------------------------------------------------
27
28 #define AliAnalysisTaskFlowSingleMu_cxx
29
30 #include "AliAnalysisTaskFlowSingleMu.h"
31
32 // ROOT includes
33 #include "TROOT.h"
34 #include "TH1.h"
35 #include "TH2.h"
36 #include "TGraphErrors.h"
37 #include "TAxis.h"
38 #include "TCanvas.h"
39 #include "TLegend.h"
40 #include "TMath.h"
41 #include "TObjString.h"
42 #include "TObjArray.h"
43 #include "TF1.h"
44 #include "TStyle.h"
45 //#include "TMCProcess.h"
46 #include "TFitResultPtr.h"
47 #include "TFile.h"
48 #include "TArrayI.h"
49 #include "TArrayD.h"
50 #include "TSystem.h"
51 #include "TGraphErrors.h"
52 #include "TLatex.h"
53
54 // STEER includes
55 #include "AliAODEvent.h"
56 #include "AliAODTrack.h"
57 #include "AliAODMCParticle.h"
58 #include "AliMCEvent.h"
59 #include "AliMCParticle.h"
60 #include "AliESDEvent.h"
61 #include "AliESDMuonTrack.h"
62 #include "AliVHeader.h"
63 #include "AliAODMCHeader.h"
64 #include "AliStack.h"
65 #include "AliEventplane.h"
66
67 // ANALYSIS includes
68 #include "AliAnalysisManager.h"
69
70 // CORRFW includes
71 #include "AliCFContainer.h"
72 #include "AliCFGridSparse.h"
73 #include "AliCFEffGrid.h"
74
75 // PWG includes
76 #include "AliVAnalysisMuon.h"
77 #include "AliMergeableCollection.h"
78 #include "AliCounterCollection.h"
79 #include "AliMuonEventCuts.h"
80 #include "AliMuonTrackCuts.h"
81 #include "AliAnalysisMuonUtility.h"
82
83
84 /// \cond CLASSIMP
85 ClassImp(AliAnalysisTaskFlowSingleMu) // Class implementation in ROOT context
86 /// \endcond
87
88 using std::ifstream;
89
90 //________________________________________________________________________
91 AliAnalysisTaskFlowSingleMu::AliAnalysisTaskFlowSingleMu() :
92   AliVAnalysisMuon(),
93   fEPKeys(0x0),
94   fRandom(0x0)
95 {
96   /// Default ctor.
97 }
98
99 //________________________________________________________________________
100 AliAnalysisTaskFlowSingleMu::AliAnalysisTaskFlowSingleMu(const char *name, const AliMuonTrackCuts& cuts) :
101   AliVAnalysisMuon(name, cuts),
102   fEPKeys(0x0),
103   fRandom(0x0)
104 {
105   //
106   /// Constructor.
107   //
108   TString epMethods = "V0A V0C TPC V0Cr1 V0Cr4";
109   fEPKeys = epMethods.Tokenize(" ");
110 }
111
112
113 //________________________________________________________________________
114 AliAnalysisTaskFlowSingleMu::~AliAnalysisTaskFlowSingleMu()
115 {
116   //
117   /// Destructor
118   //
119   delete fEPKeys;
120   delete fRandom;
121 }
122
123 //___________________________________________________________________________
124 void AliAnalysisTaskFlowSingleMu::MyUserCreateOutputObjects()
125 {
126   //
127   /// Create outputs
128   //
129
130   if ( ! fRandom ) fRandom = new TRandom3();
131   // Set seed here, to avoid having the same seed for all jobs in grid
132   fRandom->SetSeed(0);
133
134   Int_t nPtBins = 160;
135   Double_t ptMin = 0., ptMax = 80.;
136   TString ptName("Pt"), ptTitle("p_{t}"), ptUnits("GeV/c");
137   
138   Int_t nEtaBins = 25;
139   Double_t etaMin = -4.5, etaMax = -2.;
140   TString etaName("Eta"), etaTitle("#eta"), etaUnits("");
141   
142   Int_t nPhiBins = 36;
143   Double_t phiMin = 0.; Double_t phiMax = 2.*TMath::Pi();
144   TString phiName("Phi"), phiTitle("#phi"), phiUnits("rad");
145
146   Int_t nDeltaPhiBins = 18;
147   Double_t deltaPhiMin = 0.; Double_t deltaPhiMax = TMath::Pi();
148   TString deltaPhiName("DeltaPhi"), deltaPhiTitle("#phi_{#mu} - #psi_{plane}"), deltaPhiUnits("rad");
149     
150   Int_t nChargeBins = 2;
151   Double_t chargeMin = -2, chargeMax = 2.;
152   TString chargeName("Charge"), chargeTitle("charge"), chargeUnits("e");
153   
154   Int_t nMotherTypeBins = kNtrackSources;
155   Double_t motherTypeMin = -0.5, motherTypeMax = (Double_t)kNtrackSources - 0.5;
156   TString motherType("MotherType"), motherTypeTitle("motherType"), motherTypeUnits("");
157
158   Int_t nbins[kNvars] = {nPtBins, nEtaBins, nPhiBins, nDeltaPhiBins, nChargeBins, nMotherTypeBins};
159   Double_t xmin[kNvars] = {ptMin, etaMin, phiMin, deltaPhiMin, chargeMin, motherTypeMin};
160   Double_t xmax[kNvars] = {ptMax, etaMax, phiMax, deltaPhiMax, chargeMax, motherTypeMax};
161   TString axisTitle[kNvars] = {ptTitle, etaTitle, phiTitle, deltaPhiTitle, chargeTitle, motherTypeTitle};
162   TString axisUnits[kNvars] = {ptUnits, etaUnits, phiUnits, deltaPhiUnits, chargeUnits, motherTypeUnits};
163
164   Int_t iobj=0;
165   TH1* histo = 0x0;
166   TString currTitle = "";
167   for ( Int_t iep=0; iep<fEPKeys->GetEntries(); iep++ ) {
168     TString currEP = fEPKeys->At(iep)->GetName();
169     AliCFContainer* cfContainer = new AliCFContainer(Form("FlowSingleMuContainer%s", currEP.Data()),"Container for tracks",kNsteps,kNvars,nbins);
170   
171     TString currName = "";
172     for ( Int_t idim = 0; idim<kNvars; idim++){
173       currName = Form("%s (%s)", axisTitle[idim].Data(), axisUnits[idim].Data());
174       currName.ReplaceAll("()","");
175
176       cfContainer->SetVarTitle(idim, currName.Data());
177       cfContainer->SetBinLimits(idim, xmin[idim], xmax[idim]);
178     }
179     
180     TString stepTitle[kNsteps] = {"reconstructed", "generated"};
181
182     TAxis* currAxis = 0x0;
183     for (Int_t istep=0; istep<kNsteps; istep++){
184       cfContainer->SetStepTitle(istep, stepTitle[istep].Data());
185       AliCFGridSparse* gridSparse = cfContainer->GetGrid(istep);
186
187       currAxis = gridSparse->GetAxis(kHvarMotherType);
188       for ( Int_t ibin=0; ibin<fSrcKeys->GetEntries(); ibin++ ) {
189         currAxis->SetBinLabel(ibin+1, fSrcKeys->At(ibin)->GetName());
190       }
191     }
192
193     AddObjectToCollection(cfContainer, iobj++);
194
195     Bool_t isTPC = currEP.Contains("TPC");
196     Double_t minPsiPlane = ( isTPC ) ? 0. : -TMath::Pi()/2.;
197     Double_t maxPsiPlane = ( isTPC ) ? TMath::Pi() : TMath::Pi()/2.;
198     
199     histo = new TH1D(Form("hEventPlane%s", currEP.Data()), Form("Event plane distribution: %s", currEP.Data()), 18, minPsiPlane, maxPsiPlane);
200     histo->GetXaxis()->SetTitle(Form("#psi_{%s} (rad)", currEP.Data()));
201     AddObjectToCollection(histo, iobj++);
202
203     for ( Int_t jep=iep+1; jep<fEPKeys->GetEntries(); jep++ ) {
204       TString auxEP = fEPKeys->At(jep)->GetName();
205       currTitle = Form("cos(#psi_{%s} - #psi_{%s})",currEP.Data(),auxEP.Data());
206       histo = new TH1D(Form("hResoSub3%s%s",currEP.Data(),auxEP.Data()), currTitle.Data(), 100, -1.,1.);
207       histo->GetXaxis()->SetTitle(currTitle.Data());
208       AddObjectToCollection(histo, iobj++);
209     }
210
211     currTitle = Form("cos(#psi_{1} - #psi_{2})");
212     histo = new TH1D(Form("hResoSub2%s",currEP.Data()), Form("%s: %s", currEP.Data(), currTitle.Data()), 100, -1.,1.);
213     histo->GetXaxis()->SetTitle(currTitle.Data());
214
215     AddObjectToCollection(histo, iobj++);
216   } // loop on EPs
217
218   fMuonTrackCuts->Print("mask");
219 }
220
221 //________________________________________________________________________
222 void AliAnalysisTaskFlowSingleMu::ProcessEvent(TString physSel, const TObjArray& selectTrigClasses, TString centrality)
223 {
224   //
225   /// Fill output objects
226   //
227
228   //
229   // Base event cuts
230   //
231   //Bool_t isPileupFromSPD = ( fAODEvent && ! fAODEvent->GetTracklets() ) ? InputEvent()->IsPileupFromSPD(3, 0.8, 3., 2., 5.) : InputEvent()->IsPileupFromSPDInMultBins(); // Avoid break when reading Muon AODs (tracklet info is not present and IsPileupFromSPDInMultBins crashes
232   //if ( isPileupFromSPD ) return;
233
234   TArrayD psiPlane(fEPKeys->GetEntriesFast());
235   psiPlane.Reset(-999.);
236   const Double_t kLowestVal = -10.;
237   Int_t harmonic = 2;
238
239   //
240   // Fill EP info
241   //
242   AliEventplane* eventPlane = InputEvent()->GetEventplane();
243   psiPlane[0] = eventPlane->GetEventplane("V0A",InputEvent(),harmonic);
244   psiPlane[1] = eventPlane->GetEventplane("V0C",InputEvent(),harmonic);
245   psiPlane[2] = eventPlane->GetEventplane("Q");
246   if ( psiPlane[2] < 0.) psiPlane[2] = -999.;
247   Double_t qx = 0., qy = 0.;
248   psiPlane[3] = eventPlane->CalculateVZEROEventPlane(InputEvent(),0,harmonic,qx,qy); // V0C ring 1
249   psiPlane[4] = eventPlane->CalculateVZEROEventPlane(InputEvent(),3,harmonic,qx,qy); // V0C ring 3
250   //psiPlane[3] = fRandom->Uniform(-TMath::Pi()/2., TMath::Pi()/2.);
251
252   Bool_t noValidEP = kTRUE;
253   TArrayD psiPlaneCalc(fEPKeys->GetEntriesFast());
254   for ( Int_t iep=0; iep<fEPKeys->GetEntriesFast(); iep++ ) {
255     psiPlaneCalc[iep] = psiPlane[iep];
256     if ( psiPlane[iep] < kLowestVal ) continue;
257     if ( psiPlane[iep] < 0. ) psiPlaneCalc[iep] += TMath::Pi();
258     noValidEP = kFALSE;
259   }
260   if ( noValidEP ) return;
261
262   for ( Int_t iep=0; iep<fEPKeys->GetEntriesFast(); iep++ ) {
263     if ( psiPlane[iep] < kLowestVal ) continue; 
264     for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
265       TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
266       TString currEP = fEPKeys->At(iep)->GetName();
267       ((TH1*)GetMergeableObject(physSel, trigClassName, centrality, Form("hEventPlane%s",currEP.Data())))->Fill(psiPlane[iep]);
268       for ( Int_t jep=iep+1; jep<fEPKeys->GetEntriesFast(); jep++ ) {
269         if ( psiPlane[jep] < kLowestVal ) continue; 
270         ((TH1*)GetMergeableObject(physSel, trigClassName, centrality, Form("hResoSub3%s%s",currEP.Data(), fEPKeys->At(jep)->GetName())))->Fill(TMath::Cos(2.*(psiPlaneCalc[iep]-psiPlaneCalc[jep])));
271       } // loop on auxialiary EP
272       if ( iep == 2 ) {
273         TVector2* qsub1 = eventPlane->GetQsub1();
274         TVector2* qsub2 = eventPlane->GetQsub2();
275         if ( qsub1 && qsub2 )
276           ((TH1*)GetMergeableObject(physSel, trigClassName, centrality, Form("hResoSub2%s",currEP.Data())))->Fill(TMath::Cos(qsub1->Phi()/2.-qsub2->Phi()/2.));
277       }
278       //else ((TH1*)GetMergeableObject(physSel, trigClassName, centrality, Form("hResoSub2%s",currEP.Data())))->Fill(TMath::Cos(2.*psiPlane[iep])); // REMEMBER TO CUT
279     } // loop on trigger
280   } // loop on EP
281
282   //
283   // Fill track info
284   //
285   Double_t containerInput[kNvars];
286   AliVParticle* track = 0x0;
287   Int_t nSteps = MCEvent() ? 2 : 1;
288   for ( Int_t istep = 0; istep<nSteps; ++istep ) {
289     Int_t nTracks = ( istep == kStepReconstructed ) ? AliAnalysisMuonUtility::GetNTracks(InputEvent()) : MCEvent()->GetNumberOfTracks();
290     for (Int_t itrack = 0; itrack < nTracks; itrack++) {
291       track = ( istep == kStepReconstructed ) ? AliAnalysisMuonUtility::GetTrack(itrack,InputEvent()) : MCEvent()->GetTrack(itrack);
292       
293       Bool_t isSelected = ( istep == kStepReconstructed ) ? fMuonTrackCuts->IsSelected(track) : ( TMath::Abs(track->PdgCode()) == 13 );
294       if ( ! isSelected ) continue;
295       
296       Int_t trackSrc = GetParticleType(track);
297  
298
299       containerInput[kHvarPt]         = track->Pt();
300       containerInput[kHvarEta]        = track->Eta();
301       containerInput[kHvarPhi]        = track->Phi();
302       containerInput[kHvarCharge]     = track->Charge()/3.;
303       containerInput[kHvarMotherType] = (Double_t)trackSrc;
304
305       for ( Int_t iep=0; iep<fEPKeys->GetEntriesFast(); iep++ ) {
306         if ( psiPlane[iep] < kLowestVal ) continue; 
307         TString currEP = fEPKeys->At(iep)->GetName();
308         Double_t deltaPhi = containerInput[kHvarPhi] - psiPlane[iep];
309
310         // The difference between phi - psiPlane must be between (0, pi)
311         // Since the reaction plain is symmetric in pi,
312         // we can do in such a way that:
313         // -pi < delta_phi + N*pi < pi
314         // We choose N accrodingly
315         Double_t nPi = TMath::Floor( 1. - deltaPhi / TMath::Pi() );
316         deltaPhi += nPi * TMath::Pi();
317         containerInput[kHvarDeltaPhi] = deltaPhi;
318
319         for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
320           TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
321           if ( istep == kStepReconstructed && ! fMuonTrackCuts->TrackPtCutMatchTrigClass(track, fMuonEventCuts->GetTrigClassPtCutLevel(trigClassName)) ) continue;
322           ((AliCFContainer*)GetMergeableObject(physSel, trigClassName, centrality, Form("FlowSingleMuContainer%s",currEP.Data())))->Fill(containerInput,istep);
323         } // loop on selected trigger classes
324       } // loop on event plane
325     } // loop on tracks
326   } // loop on container steps
327 }
328
329
330 //________________________________________________________________________
331 TArrayD AliAnalysisTaskFlowSingleMu::GetCentralityRange(TString sRange)
332 {
333   //
334   /// Get the range from string
335   //
336   TArrayD centralityRange(2);
337   centralityRange.Reset(-999.);
338   TObjArray* array = sRange.Tokenize("_");
339   if ( array->GetEntries() == 2 ) {
340     for ( Int_t ientry=0; ientry<2; ientry++ ) {
341       centralityRange[ientry] = ((TObjString*)array->At(ientry))->GetString().Atof();
342     }
343   }
344   delete array;
345   return centralityRange;
346 }
347
348 //________________________________________________________________________
349 void AliAnalysisTaskFlowSingleMu::Terminate(Option_t *) {
350   //
351   /// Draw some histograms at the end.
352   //
353
354   AliVAnalysisMuon::Terminate("");
355
356   if ( ! fMergeableCollection ) return;
357   
358   TString physSel = fTerminateOptions->At(0)->GetName();
359   TString trigClassName = fTerminateOptions->At(1)->GetName();
360   TString centralityRange = fTerminateOptions->At(2)->GetName();
361   TString furtherOpt = fTerminateOptions->At(3)->GetName();
362
363   TObjArray* centralityClasses = centralityRange.Tokenize(",");
364   TArrayD centralityArray(centralityClasses->GetEntries()+1);
365   Int_t ib = 0;
366   for ( Int_t icent=0; icent<centralityClasses->GetEntries(); icent++ ) {
367     TArrayD range = GetCentralityRange(centralityClasses->At(icent)->GetName());
368     if ( icent == 0 ) centralityArray.AddAt(range[0],ib++);
369     centralityArray.AddAt(range[1],ib++);
370   }
371   
372   TString selectEP = "";
373   TObjArray resoEParray;
374   resoEParray.SetOwner();
375   
376   TString outFilename = "";
377   TObjArray* optArr = furtherOpt.Tokenize(" ");
378   TString currName = "";
379   TString resoFilename = "";
380   for ( Int_t iopt=0; iopt<optArr->GetEntries(); iopt++ ) {
381     currName = optArr->At(iopt)->GetName();
382     if ( currName.Contains(".root") ) outFilename = currName;
383     else if ( currName.Contains(".txt") ) resoFilename = currName;
384     else {
385       for ( Int_t iep=0; iep<fEPKeys->GetEntries(); iep++ ) {
386         if ( ! currName.CompareTo(fEPKeys->At(iep)->GetName()) ) resoEParray.Add(new TObjString(currName.Data()));
387       }
388     }
389   }
390   delete optArr;
391
392   if ( resoEParray.At(0) ) selectEP = resoEParray.At(0)->GetName();
393
394   furtherOpt.ToUpper();
395
396   Bool_t v2only = furtherOpt.Contains("V2ONLY");
397
398   TAxis* customBins = 0x0;
399   
400   //Double_t myBins[] = {0., 0.5, 1., 1.5, 2., 3., 4., 5., 6., 8., 10.};
401   //Double_t myBins[] = {4., 5., 6., 8., 10.}; // REMEMBER TO CHOOSE
402   //Double_t myBins[] = {6., 8., 10.};
403   Double_t myBins[] = {2., 3., 4., 5., 6., 8., 10.};
404   //Double_t myBins[] = {0., 0.5, 1., 1.5, 2., 3., 4., 10.};
405   //Double_t myBins[] = {2., 3., 4., 5., 6., 8., 10., 80.};
406   Int_t nMyBins = sizeof(myBins)/sizeof(myBins[0])-1;
407   customBins = new TAxis(nMyBins, myBins); // REMEMBER TO CHOOSE
408
409   TCanvas *can = 0x0, *showCan = 0x0;
410   Int_t xshift = 100;
411   Int_t yshift = 20;
412   Int_t igroup1 = 0;
413   Int_t igroup2 = 0;
414   
415   TList outList;
416
417   //
418   /// Read plane resolution file (if present)
419   //
420   TObjArray resoCentrality;
421   resoCentrality.SetOwner();
422   TArrayD resolution[3];
423   for ( Int_t ires=0; ires<3; ires++ ) {
424     resolution[ires].Set(50);
425   }
426
427   Bool_t externalReso = ( ! resoFilename.IsNull() && ! gSystem->AccessPathName(resoFilename.Data()) );
428   Bool_t internalReso = kTRUE;
429   if ( externalReso ) {
430     // If plane resolution file exists, fill resolution
431     TString currLine = "";
432     std::ifstream inFile(resoFilename.Data());
433     if (inFile.is_open()) {
434       ib = 0;
435       while (! inFile.eof() ) {
436         currLine.ReadLine(inFile,kFALSE);
437         if ( currLine.BeginsWith("#") || currLine.Length() < 3 ) continue;
438         TObjArray* array = currLine.Tokenize(" ");
439         resoCentrality.Add(new TObjString(array->At(0)->GetName()));
440         for ( Int_t ires=0; ires<3; ires++ ) {
441           resolution[ires].AddAt(((TObjString*)array->At(ires+1))->GetString().Atof(),ib);
442         }
443         delete array;
444         ib++;
445       } // ! EOF
446     } // file is open
447     for ( Int_t ires=0; ires<3; ires++ ) {
448       resolution[ires].Set(ib);
449     }
450   } // file exists
451   else {
452     // If no plane resolution file exist,
453     // use the centralities in terminateOption
454     // and set the resolution to 1 and errors to 0.
455
456     for ( Int_t ires=0; ires<3; ires++ ) {
457       resolution[ires].Set(centralityClasses->GetEntries());
458       Double_t resetVal = ( ires == 0 ) ? 1. : 0.;
459       resolution[ires].Reset(resetVal);
460     }
461     TArrayD currCos(3);
462     for ( Int_t icent=0; icent<centralityClasses->GetEntries(); icent++ ) {
463       resoCentrality.Add(new TObjString(centralityClasses->At(icent)->GetName()));
464       Int_t icos = 0;
465       Double_t sumRelErrSquare = 0.;
466       TString hResoTitle = "";
467       for ( Int_t isub=0; isub<resoEParray.GetEntries(); isub++ ) {
468         TString resoDet1 = resoEParray.At(isub)->GetName();
469         for ( Int_t jsub=isub+1; jsub<resoEParray.GetEntries(); jsub++ ) {
470           TString resoDet2 = resoEParray.At(jsub)->GetName();
471           TH1* histo = 0x0;
472           // Try both combinations: we want the first two values to contain the selectedEP!
473           TString resoDet = "";
474           for ( Int_t icomb=0; icomb<2; icomb++ ) {
475             TString hName = "hResoSub3";
476             resoDet = ( icomb == 0 ) ? resoDet1 + resoDet2 : resoDet2 + resoDet1;
477             hName += resoDet;
478             histo = static_cast<TH1*>(GetSum(physSel,trigClassName,resoCentrality.At(icent)->GetName(),hName));
479             if ( histo ) break;
480           }
481           if ( ! histo || histo->Integral() == 0 ) continue;
482           currCos[icos] = histo->GetMean();
483           if ( currCos[icos] <= 0. ) continue;
484           Double_t relErr = histo->GetMeanError() / currCos[icos];
485           sumRelErrSquare += relErr * relErr;
486           icos++;
487           hResoTitle += Form("%s ", resoDet.Data());
488         } // loop on sub events
489       } // loop on sub events
490       if ( icos != 3  ) {
491         printf("Warning: resolution cannot be estimated for %s\n", resoCentrality.At(icent)->GetName());
492         internalReso = kFALSE;
493         continue;
494       }
495       resolution[0][icent] = TMath::Sqrt((currCos[0]*currCos[1]/currCos[2]));
496       resolution[1][icent] = resolution[0][icent]/2.*TMath::Sqrt(sumRelErrSquare);
497       printf("Resolution (%s) %s  %g +- %g\n", hResoTitle.Data(), centralityClasses->At(icent)->GetName(), resolution[0][icent], resolution[1][icent]);
498     }
499   }
500
501   Bool_t hasResolution = ( externalReso || internalReso );
502
503
504   TArrayD resoCentrArray(resoCentrality.GetEntries()+1);
505   ib = 0;
506   for ( Int_t icent=0; icent<resoCentrality.GetEntries(); icent++ ) {
507     TArrayD range = GetCentralityRange(resoCentrality.At(icent)->GetName());
508     if ( icent == 0 ) resoCentrArray.AddAt(range[0],ib++);
509     resoCentrArray.AddAt(range[1], ib++);
510   }
511   
512   const Int_t kNresoCentr = resoCentrality.GetEntries();
513
514   //
515   // Create gridSparse array
516   //
517   AliCFGridSparse* gridSparseArray[kNsteps*kNresoCentr];
518   TString stepName[2];
519   for ( Int_t icent=0; icent<kNresoCentr; icent++ ) {
520     for ( Int_t istep=0; istep<kNsteps; istep++ ) {
521       gridSparseArray[istep*kNresoCentr+icent] = 0x0;
522     }
523
524     AliCFContainer* cfContainer = static_cast<AliCFContainer*> ( GetSum(physSel,trigClassName,resoCentrality.At(icent)->GetName(),Form("FlowSingleMuContainer%s", selectEP.Data())) );
525     if ( ! cfContainer ) continue;
526     for ( Int_t istep=0; istep<kNsteps; istep++ ) {
527       stepName[istep] = cfContainer->GetStepTitle(istep);
528       gridSparseArray[istep*kNresoCentr+icent] = cfContainer->GetGrid(istep);
529       if ( ! customBins ) customBins = gridSparseArray[istep*kNresoCentr+icent]->GetAxis(kHvarPt);
530     } // loop on step
531   } // loop on centrality
532   
533   
534   Bool_t isMC = furtherOpt.Contains("MC");
535   Int_t firstSrc = ( isMC ) ? 0 : kUnidentified;
536   Int_t lastSrc  = ( isMC ) ? kNtrackSources - 1 : kUnidentified;
537   Int_t srcColors[kNtrackSources] = {kBlack, kRed, kSpring, kTeal, kBlue, kViolet, kMagenta, kOrange};
538   if ( ! isMC ) srcColors[kUnidentified] = 1;
539   
540   
541   //
542   // Get histograms
543   //
544   const Int_t kNcentralities = centralityClasses->GetEntries();
545   const Int_t kNhistos = kNsteps*kNtrackSources*kNcentralities;
546   TObjArray ptVsDeltaPhiList(kNhistos);
547   TObjArray ptVsPhiList(kNhistos);
548   TArrayD wgtReso[3];
549   for ( Int_t ires=0; ires<3; ires++ ) {
550     wgtReso[ires].Set(kNhistos);
551     wgtReso[ires].Reset(0.);
552   }
553   for ( Int_t istep=0; istep<kNsteps; istep++ ) {
554     for ( Int_t isrc = firstSrc; isrc <= lastSrc; ++isrc ) {
555       TString baseName = Form("%s_%s", stepName[istep].Data(), fSrcKeys->At(isrc)->GetName());
556       for ( Int_t icent=0; icent<centralityClasses->GetEntries(); icent++ ) {
557         Int_t idx = kNcentralities * kNtrackSources * istep + kNcentralities * isrc + icent;
558         TH1* hPtVsDeltaPhi = 0x0;
559         TH1* hPtVsPhi = 0x0;
560         TString baseNameCent = Form("%s_%s", baseName.Data(), centralityClasses->At(icent)->GetName());
561         Double_t nMuons = 0.;
562         TString debugString = Form("\n%s", centralityClasses->At(icent)->GetName());
563         for ( Int_t rcent=0; rcent<kNresoCentr; rcent++ ) {
564           if ( resoCentrArray[rcent] < centralityArray[icent] ||
565               resoCentrArray[rcent+1] > centralityArray[icent+1] ) continue;
566           AliCFGridSparse* gridSparse = gridSparseArray[istep*kNresoCentr+rcent];
567           if ( ! gridSparse ||  gridSparse->GetEntries() == 0. ) continue;
568           
569           SetSparseRange(gridSparse, kHvarEta, "", -3.999, -2.501);
570           SetSparseRange(gridSparse, kHvarMotherType, "", isrc+1, isrc+1, "USEBIN");
571           TH1* auxHisto = gridSparse->Project(kHvarDeltaPhi,kHvarPt);
572           
573           Double_t currYield = auxHisto->Integral();
574           if ( currYield == 0. ) continue;
575           nMuons += currYield;
576           debugString += Form("\nAdding %s  yield %g  reso", resoCentrality.At(rcent)->GetName(), currYield);
577           for ( Int_t ires=0; ires<3; ires++ ) {
578             wgtReso[ires][idx] += resolution[ires][rcent] * currYield;
579             debugString += Form("  %g", resolution[ires][rcent]);
580           }
581           if ( ! hPtVsDeltaPhi ) hPtVsDeltaPhi = static_cast<TH1*>(auxHisto->Clone(Form("ptVsDeltaPhi_%s", baseNameCent.Data())));
582           else hPtVsDeltaPhi->Add(auxHisto);
583           delete auxHisto;
584           
585           auxHisto = gridSparse->Project(kHvarPhi,kHvarPt);
586           if ( ! hPtVsPhi ) hPtVsPhi = static_cast<TH1*>(auxHisto->Clone(Form("ptVsPhi_%s", baseNameCent.Data())));
587           else hPtVsPhi->Add(auxHisto);
588           delete auxHisto;
589           
590         } // loop on resolution centralities
591         
592         if ( nMuons == 0. ) continue;
593         debugString += Form("\nWeighted reso ");
594         for ( Int_t ires=0; ires<3; ires++ ) {
595           wgtReso[ires][idx] = wgtReso[ires][idx]/nMuons;
596           debugString += Form("  %g", wgtReso[ires][idx]);
597         }
598         AliDebug(1,debugString.Data());
599         
600         ptVsDeltaPhiList.AddAt(hPtVsDeltaPhi,idx);
601         ptVsPhiList.AddAt(hPtVsPhi,idx);
602       } // loop on centralities
603     } // loop on sources
604   } // loop on steps
605   
606   
607   
608   /////////////////////
609   // Plot resolution //
610   /////////////////////
611   if ( hasResolution ) {
612     currName = Form("reso_%s", selectEP.Data());
613     currName += externalReso ? Form("_external") : Form("_sub_%s_%s_%s", resoEParray.At(0)->GetName(), resoEParray.At(1)->GetName(), resoEParray.At(2)->GetName());
614     TGraphErrors* resoGraph = new TGraphErrors(kNresoCentr);
615     resoGraph->SetName(currName.Data());
616     currName.Append("_can");
617     can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
618     can->cd();
619     
620     for ( Int_t icent=0; icent<kNresoCentr; icent++ ) {
621       resoGraph->SetPoint(icent, 0.5*(resoCentrArray[icent+1] + resoCentrArray[icent]), resolution[0][icent]);
622       resoGraph->SetPointError(icent, 0.5*(resoCentrArray[icent+1] - resoCentrArray[icent]), resolution[1][icent]);
623     }
624     resoGraph->SetTitle(Form("Resolution %s", selectEP.Data()));
625     resoGraph->GetXaxis()->SetTitle("Centrality");
626     resoGraph->GetYaxis()->SetTitle("Resolution");
627     resoGraph->Draw("apz");
628     outList.Add(resoGraph);
629     
630     for ( Int_t istep=0; istep<kNsteps; istep++ ) {
631       for ( Int_t isrc = firstSrc; isrc <= lastSrc; ++isrc ) {
632         currName = Form("wgtReso_%s_%s", stepName[istep].Data(), fSrcKeys->At(isrc)->GetName());
633         resoGraph = new TGraphErrors();
634         resoGraph->SetName(currName.Data());
635         for ( Int_t icent=0; icent<kNcentralities; icent++ ) {
636           Int_t idx = kNcentralities * kNtrackSources * istep + kNcentralities * isrc + icent;
637           if ( wgtReso[0][idx] == 0. ) continue;
638           resoGraph->SetPoint(icent, 0.5*(centralityArray[icent+1] + centralityArray[icent]), wgtReso[0][idx]);
639           resoGraph->SetPointError(icent, 0.5*(centralityArray[icent+1] - centralityArray[icent]), wgtReso[1][idx]);
640         } // loop on centralities
641         if ( resoGraph->GetN() == 0 ) continue;
642         printf("New reso %s\n", resoGraph->GetName());
643         resoGraph->SetTitle(Form("Resolution %s", selectEP.Data()));
644         resoGraph->GetXaxis()->SetTitle("Centrality");
645         resoGraph->GetYaxis()->SetTitle("Resolution");
646         resoGraph->Draw("pz");
647         outList.Add(resoGraph);
648       } // loop on sources
649     } // loop on steps
650   }
651   
652   
653   
654   TString graphTypeName[3] = {"raw","correctStat","correctSyst"};
655   Int_t drawOrder[3] = {0, 2, 1};
656   TString drawOpt[3] = {"apz","pz","a2"};
657   Int_t nTypes = ( hasResolution ) ? 3 : 1;
658
659   TString histoName = "", histoPattern = "";
660   ///////////
661   // Flow  //
662   ///////////
663   
664   gStyle->SetOptFit(1111);
665   TString funcFormula = ( v2only ) ? "[0] * (1. + 2.*[1]*TMath::Cos(2*x))" : "[0] * (1. + 2.*([1]*TMath::Cos(x) + [2]*TMath::Cos(2*x) + [3]*TMath::Cos(3*x)))";
666   //
667   TF1* func = new TF1("funcFlow",funcFormula.Data(), 0., TMath::Pi());
668   if  ( v2only ) func->SetParNames("scale", "v2");
669   else func->SetParNames("scale","v1", "v2", "v3");
670   Int_t v2par = ( v2only ) ? 1 : 2;
671   
672   for ( Int_t istep=0; istep<kNsteps; istep++ ) {
673     for ( Int_t isrc = firstSrc; isrc <= lastSrc; ++isrc ) {
674       TString baseName = Form("%s_%s", stepName[istep].Data(), fSrcKeys->At(isrc)->GetName());
675       TGraphErrors* flowVsCentrality[3];
676       for ( Int_t itype=0; itype<nTypes; itype++ ) {
677         histoName = Form("v2VsCentrality_%s_%s", baseName.Data(), graphTypeName[itype].Data());
678         flowVsCentrality[itype] = new TGraphErrors();
679         flowVsCentrality[itype]->SetName(histoName.Data());
680         flowVsCentrality[itype]->SetTitle(histoName.Data());
681       }
682       for ( Int_t icent=0; icent<centralityClasses->GetEntries(); icent++ ) {
683         Int_t idx = kNcentralities * kNtrackSources * istep + kNcentralities * isrc + icent;
684         TH2* histo2D = static_cast<TH2*> (ptVsDeltaPhiList.At(idx));
685         if ( ! histo2D ) continue;
686         TString baseNameCent = Form("%s_%s", baseName.Data(), centralityClasses->At(icent)->GetName());
687
688         TAxis* ptAxis = histo2D->GetYaxis();
689         
690         Int_t ipad = 0;
691         TGraphErrors* flowVsPt[3];
692         for ( Int_t itype=0; itype<nTypes; itype++ ) {
693           histoName = Form("v2VsPt_%s_%s", baseNameCent.Data(), graphTypeName[itype].Data());
694           flowVsPt[itype] = new TGraphErrors();
695           flowVsPt[itype]->SetName(histoName.Data());
696           flowVsPt[itype]->SetTitle(histoName.Data());
697         }
698         
699         for ( Int_t ipt=0; ipt<=customBins->GetNbins(); ipt++ ) {
700           Int_t minCustomBin = ( ipt == 0 ) ? 1 : ipt;
701           Int_t maxCustomBin = ( ipt == 0 ) ? customBins->GetNbins() : ipt;
702           Int_t ptMinBin = ptAxis->FindBin(customBins->GetBinLowEdge(minCustomBin)+1.e-4);
703           Int_t ptMaxBin = ptAxis->FindBin(customBins->GetBinUpEdge(maxCustomBin)-1.e-4);
704           Double_t ptMin = ptAxis->GetBinLowEdge(ptMinBin);
705           Double_t ptMax = ptAxis->GetBinUpEdge(ptMaxBin);
706           TH1* projHisto = histo2D->ProjectionX(Form("%s_%s_ptBin%i", baseName.Data(), histo2D->GetName(), ipt), ptMinBin, ptMaxBin);
707           projHisto->SetTitle(Form("%g < p_{t} (GeV/c) < %g", ptMin, ptMax));
708           if ( projHisto->Integral() < 50 ) break;
709           if ( ipad%4 == 0 ) {
710             currName = projHisto->GetName();
711             currName.Append("_can");
712             showCan = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
713             showCan->Divide(2,2);
714             ipad = 0;
715           }
716           showCan->cd(++ipad);
717           if ( v2only ) func->SetParameters(projHisto->Integral(), 0.01);
718           else {
719             func->SetParameters(projHisto->Integral(), 0., 0.01, 0.);
720             func->FixParameter(1,0.);
721           }
722           projHisto->Fit(func,"R");
723           for ( Int_t itype=0; itype<nTypes; itype++ ) {
724             TGraphErrors* currGraph = ( ipt == 0 ) ? flowVsCentrality[itype] : flowVsPt[itype];
725             Double_t resoVal = ( itype == 0 ) ? 1. : wgtReso[0][idx];
726             Double_t resoErr = ( itype == 0 ) ? 0. : wgtReso[itype][idx];
727             Double_t rawVal = func->GetParameter(v2par);
728             Double_t rawErr = (itype == 2 ) ? 0. : func->GetParError(v2par);
729             Double_t yVal = rawVal/resoVal;
730             Double_t xVal = ( ipt == 0 ) ? 0.5*(centralityArray[icent+1] + centralityArray[icent]) : 0.5*(ptMax + ptMin);
731             Double_t xErr = ( ipt == 0 ) ? 0.5*(centralityArray[icent+1] - centralityArray[icent]) : 0.5*(ptMax - ptMin);
732             Int_t ipoint = currGraph->GetN();
733             currGraph->SetPoint(ipoint,xVal,yVal);
734             if ( itype > 0 && ipt != 0 ) {
735               // For the plot vs. pt the resolution error is fully correlated.
736               // Hence we do not use add it bin by bin, but rather write it
737               // on the title
738               resoErr = 0.;
739               TString graphTitle = currGraph->GetTitle();
740               if ( ! graphTitle.Contains("|") ) {
741                 graphTitle += Form(" | Rel. error on EP resolution: %.2g (stat.) %.2g (syst.)", wgtReso[1][idx]/wgtReso[0][idx], wgtReso[2][idx]/wgtReso[0][idx]);
742                 currGraph->SetTitle(graphTitle.Data());
743               }
744             }
745             Double_t rawErrRel = ( rawVal != 0. ) ? rawErr/rawVal : 0.;
746             Double_t resoErrRel = ( resoVal != 0. ) ? resoErr/resoVal : 0.;
747             Double_t yErr = TMath::Abs(yVal) * TMath::Sqrt(rawErrRel*rawErrRel+resoErrRel*resoErrRel);
748             currGraph->SetPointError(ipoint,xErr,yErr);
749           } // loop on type
750         } // loop on pt bins
751         for ( Int_t itype=0; itype<nTypes; itype++ ) {
752           Int_t currType = drawOrder[itype];
753           if ( itype < 2 ) {
754             currName = flowVsPt[currType]->GetName();
755             currName.Append("Can");
756             can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
757             igroup2++;
758             can->cd();
759           }
760           flowVsPt[currType]->GetXaxis()->SetTitle(ptAxis->GetTitle());
761           flowVsPt[currType]->GetYaxis()->SetTitle("v2");
762           flowVsPt[currType]->GetYaxis()->SetRangeUser(0., 0.25);
763           flowVsPt[currType]->SetFillStyle(0);
764           flowVsPt[currType]->Draw(drawOpt[currType].Data());
765           outList.Add(flowVsPt[currType]);
766         } // loop on types
767       } // loop on centrality
768       for ( Int_t itype=0; itype<nTypes; itype++ ) {
769         Int_t currType = drawOrder[itype];
770         if ( flowVsCentrality[currType]->GetN() == 0 ) continue;
771         if ( itype < 2 ) {
772           currName = flowVsCentrality[currType]->GetName();
773           currName.Append("Can");
774           can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
775           can->cd();
776         }
777         flowVsCentrality[currType]->GetXaxis()->SetTitle("Centrality");
778         flowVsCentrality[currType]->GetYaxis()->SetTitle("v2");
779         flowVsCentrality[currType]->SetFillStyle(0);
780         flowVsCentrality[currType]->GetYaxis()->SetRangeUser(0., 0.25);
781         flowVsCentrality[currType]->Draw(drawOpt[currType].Data());
782         outList.Add(flowVsCentrality[currType]);
783         igroup2++;
784       } // loop on type
785       igroup1++;
786       igroup2 = 0;
787     } // loop on track sources
788   } // loop on steps
789
790   ///////////////////////
791   // Event plane check //
792   ///////////////////////
793   igroup1++;
794   igroup2 = 0;
795   Int_t icolor=0;
796   currName ="eventPlaneDistrib";
797   can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
798   can->SetGridy();
799   igroup2++;
800
801   const Int_t kNfits = 2;
802   TString fitTypeName[kNfits] = {"Cos","Sin"};
803
804   Bool_t addOffsetToCentrality = kFALSE;
805   Double_t offsetStep = 0.1;
806
807   TLegend* leg = new TLegend(0.7,0.7,0.92,0.92);
808   leg->SetBorderSize(1);
809   
810   TGraphErrors* epCheck[kNfits];
811   for ( Int_t itype=0; itype<kNfits; itype++ ) {
812     histoName = Form("checkFourier_%s",fitTypeName[itype].Data());
813     epCheck[itype] = new TGraphErrors();
814     epCheck[itype]->SetName(histoName.Data());
815 //    epCheck[itype]->SetTitle(histoName.Data());
816   }
817   for ( Int_t icent=0; icent<centralityClasses->GetEntries(); icent++ ) {
818     TH1* histo = static_cast<TH1*> ( GetSum(physSel,trigClassName,centralityClasses->At(icent)->GetName(),Form("hEventPlane%s", selectEP.Data())) );
819     if ( ! histo ) continue;
820     histo->SetName(Form("%s_%s", histo->GetName(), centralityClasses->At(icent)->GetName()));
821     if ( histo->Integral() < 50. ) continue;
822     if ( histo->GetSumw2N() == 0 ) histo->Sumw2();
823     
824     for ( Int_t itype=0; itype<kNfits; itype++ ) {
825       TString checkFormula = Form("[0]*(1.+2.*[1]*TMath::%s(2.*x))",fitTypeName[itype].Data());
826       TF1* checkFunc = new TF1(Form("checkFunc_%s", fitTypeName[itype].Data()), checkFormula.Data(), histo->GetXaxis()->GetBinLowEdge(1), histo->GetXaxis()->GetBinUpEdge(histo->GetXaxis()->GetNbins()));
827       checkFunc->SetParameters(histo->Integral(), 0.);
828       histo->Fit(checkFunc,"R0");
829       TGraphErrors* currGraph = epCheck[itype];
830       Double_t yVal = checkFunc->GetParameter(1);
831       Double_t yErr = checkFunc->GetParError(1);
832       Double_t xVal = 0.5*(centralityArray[icent+1] + centralityArray[icent]);
833       Double_t xErr = 0.5*(centralityArray[icent+1] - centralityArray[icent]);
834       Int_t ipoint = currGraph->GetN();
835       currGraph->SetPoint(ipoint,xVal,yVal);
836       currGraph->SetPointError(ipoint,xErr,yErr);
837     } // loop on type
838     histo->Fit("pol0","R0");
839     TF1* pol0 = histo->GetFunction("pol0");
840     histo->Scale(1./pol0->GetParameter(0));
841     Double_t offset = ( addOffsetToCentrality ) ? offsetStep*(Double_t)icolor : 0.;
842     TGraphErrors* graph = new TGraphErrors();
843     graph->SetTitle(Form("%s %s", histo->GetTitle(), trigClassName.Data()));
844     for ( Int_t ipoint=0; ipoint<histo->GetXaxis()->GetNbins(); ipoint++ ) {
845       Int_t ibin = ipoint+1;
846       graph->SetPoint(ipoint, histo->GetXaxis()->GetBinCenter(ibin), histo->GetBinContent(ibin) + offset);
847       graph->SetPointError(ipoint, histo->GetXaxis()->GetBinWidth(ibin)/2., histo->GetBinError(ibin));
848     }
849     graph->GetXaxis()->SetTitle(histo->GetXaxis()->GetTitle());
850     graph->GetYaxis()->SetTitle("Data/Fit (pol0)");
851     TString legTitle = Form("%s %%", centralityClasses->At(icent)->GetName());
852     if ( addOffsetToCentrality ) {
853       graph->GetYaxis()->SetRangeUser(1.-offsetStep, 1.+1.5*offsetStep*(Double_t)centralityClasses->GetEntries());
854       graph->GetYaxis()->SetNdivisions(2*centralityClasses->GetEntries(),0,0);
855       legTitle += Form(" ( + %g)", offset);
856     }
857     Int_t currColor = ( icolor < kNtrackSources ) ? srcColors[icolor] : 20+icolor;
858     graph->SetLineColor(currColor);
859     graph->SetMarkerColor(currColor);
860     graph->SetMarkerStyle(20+icolor);
861     TString currDrawOpt = "pz";
862     if ( icolor == 0 ) currDrawOpt += "a";
863     graph->Draw(currDrawOpt.Data());
864     legTitle.ReplaceAll("_"," - ");
865     leg->AddEntry(graph,legTitle.Data(),"lp");
866     if ( addOffsetToCentrality ) {
867       TLatex* latex = new TLatex(histo->GetXaxis()->GetXmax()+0.5*histo->GetXaxis()->GetBinWidth(1), 1.+(Double_t)offset, Form("#chi^{2}/NDF = %.3g", pol0->GetChisquare()/(Double_t)pol0->GetNDF()));
868       latex->SetTextColor(currColor);
869       latex->SetTextSize(0.025);
870       latex->Draw("same");
871     }
872     icolor++;
873   }
874   leg->Draw("same");
875   
876   currName = "epCheckCan";
877   can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
878   can->SetGridy();
879   igroup2++;
880   leg = new TLegend(0.7,0.7,0.95,0.95);
881   leg->SetBorderSize(1);
882   for ( Int_t itype=0; itype<kNfits; itype++ ) {
883     epCheck[itype]->SetLineColor(srcColors[itype]);
884     epCheck[itype]->SetMarkerColor(srcColors[itype]);
885     epCheck[itype]->SetMarkerStyle(20+itype);
886     epCheck[itype]->GetXaxis()->SetTitle("Centrality (%)");
887     TString currDrawOpt = "pz";
888     if ( can->GetListOfPrimitives()->GetEntries() == 0 ) currDrawOpt.Prepend("a");
889     epCheck[itype]->Draw(currDrawOpt.Data());
890     leg->AddEntry(epCheck[itype],Form("<%s(2 #psi_{%s})>", fitTypeName[itype].Data(), selectEP.Data()),"lp");
891   }
892   leg->Draw("same");
893   
894   for ( Int_t istep=0; istep<kNsteps; istep++ ) {
895     for ( Int_t isrc = firstSrc; isrc <= lastSrc; ++isrc ) {
896       TString baseName = Form("%s_%s", stepName[istep].Data(), fSrcKeys->At(isrc)->GetName());
897       for ( Int_t icent=0; icent<centralityClasses->GetEntries(); icent++ ) {
898         Int_t idx = kNcentralities * kNtrackSources * istep + kNcentralities * isrc + icent;
899         TH2* histo2D = static_cast<TH2*> (ptVsPhiList.At(idx));
900         if ( ! histo2D ) continue;
901         TString baseNameCent = Form("%s_%s", baseName.Data(), centralityClasses->At(icent)->GetName());
902         
903         TAxis* ptAxis = histo2D->GetYaxis();
904         
905         Int_t ipad = 0;
906         TGraphErrors* detEffect[kNfits];
907         for ( Int_t itype=0; itype<kNfits; itype++ ) {
908           histoName = Form("checkResoVsPt_%s_%s", baseNameCent.Data(), fitTypeName[itype].Data());
909           detEffect[itype] = new TGraphErrors();
910           detEffect[itype]->SetName(histoName.Data());
911           detEffect[itype]->SetTitle(histoName.Data());
912         }
913         
914         for ( Int_t ipt=0; ipt<=customBins->GetNbins(); ipt++ ) {
915           Int_t minCustomBin = ( ipt == 0 ) ? 1 : ipt;
916           Int_t maxCustomBin = ( ipt == 0 ) ? customBins->GetNbins() : ipt;
917           Int_t ptMinBin = ptAxis->FindBin(customBins->GetBinLowEdge(minCustomBin)+1.e-4);
918           Int_t ptMaxBin = ptAxis->FindBin(customBins->GetBinUpEdge(maxCustomBin)-1.e-4);
919           Double_t ptMin = ptAxis->GetBinLowEdge(ptMinBin);
920           Double_t ptMax = ptAxis->GetBinUpEdge(ptMaxBin);
921           TH1* projHisto = histo2D->ProjectionX(Form("checkReso_%s_%s_ptBin%i", baseName.Data(), histo2D->GetName(), ipt), ptMinBin, ptMaxBin);
922           projHisto->SetTitle(Form("%g < p_{t} (GeV/c) < %g", ptMin, ptMax));
923           if ( projHisto->Integral() < 50 ) break;
924           if ( ipad%4 == 0 ) {
925             currName = projHisto->GetName();
926             currName.Append("_can");
927             showCan = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
928             showCan->Divide(2,2);
929             ipad = 0;
930           }
931           showCan->cd(++ipad);
932           for ( Int_t itype=0; itype<kNfits; itype++ ) {
933             TString checkFormula = Form("[0]*(1.+2.*[1]*TMath::%s(2.*x))",fitTypeName[itype].Data());
934             TF1* checkFunc = new TF1(Form("checkFunc_%s", fitTypeName[itype].Data()), checkFormula.Data(), projHisto->GetXaxis()->GetBinLowEdge(1), projHisto->GetXaxis()->GetBinUpEdge(projHisto->GetXaxis()->GetNbins()));
935             checkFunc->SetParameters(projHisto->Integral(), 0.);
936             projHisto->Fit(checkFunc,"R");
937             TGraphErrors* currGraph = detEffect[itype];
938             Double_t yVal = checkFunc->GetParameter(1);
939             Double_t yErr = checkFunc->GetParError(1);
940             Double_t xVal = 0.5*(ptMax + ptMin);
941             Double_t xErr = 0.5*(ptMax - ptMin);
942             Int_t ipoint = currGraph->GetN();
943             currGraph->SetPoint(ipoint,xVal,yVal);
944             currGraph->SetPointError(ipoint,xErr,yErr);
945           } // loop on type
946         } // loop on pt bins
947         for ( Int_t itype=0; itype<kNfits; itype++ ) {
948           currName = detEffect[itype]->GetName();
949           currName.Append("Can");
950           can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
951           igroup2++;
952           can->cd();
953           detEffect[itype]->GetXaxis()->SetTitle(ptAxis->GetTitle());
954           detEffect[itype]->GetYaxis()->SetTitle("v2");
955 //          detEffect[itype]->GetYaxis()->SetRangeUser(0., 0.25);
956           detEffect[itype]->SetFillStyle(0);
957           detEffect[itype]->Draw("ap");
958 //          outList.Add(detEffect[itype]);
959         } // loop on types
960       } // loop on centrality
961     } // loop on track sources
962   } // loop on steps
963
964   delete centralityClasses;
965
966   //////////////////////
967   // Event statistics //
968   //////////////////////
969   printf("\nTotal analyzed events:\n");
970   TString evtSel = Form("trigger:%s", trigClassName.Data());
971   fEventCounters->PrintSum(evtSel.Data());
972   printf("Physics selected analyzed events:\n");
973   evtSel = Form("trigger:%s/selected:yes", trigClassName.Data());
974   fEventCounters->PrintSum(evtSel.Data());
975   
976   TString countPhysSel = "any";
977   if ( physSel.Contains(fPhysSelKeys->At(kPhysSelPass)->GetName()) ) countPhysSel = "yes";
978   else if ( physSel.Contains(fPhysSelKeys->At(kPhysSelReject)->GetName()) ) countPhysSel="no";
979   countPhysSel.Prepend("selected:");
980   printf("Analyzed events vs. centrality:\n");
981   evtSel = Form("trigger:%s/%s", trigClassName.Data(), countPhysSel.Data());
982   fEventCounters->Print("centrality",evtSel.Data(),kTRUE);
983
984   if ( ! outFilename.IsNull() ) {
985     printf("\nWriting output file %s\n", outFilename.Data());
986     TFile* file = TFile::Open(outFilename.Data(), "RECREATE");
987     outList.Write();
988     file->Close();
989   }
990   
991   delete customBins;
992 }