Fixes for std:: need with the trunk of Root (Jochen, Yves)
[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 "AliMuonTrackCuts.h"
80
81
82 /// \cond CLASSIMP
83 ClassImp(AliAnalysisTaskFlowSingleMu) // Class implementation in ROOT context
84 /// \endcond
85
86 using std::ifstream;
87
88 //________________________________________________________________________
89 AliAnalysisTaskFlowSingleMu::AliAnalysisTaskFlowSingleMu() :
90   AliVAnalysisMuon(),
91   fEPKeys(0x0),
92   fRandom(0x0)
93 {
94   /// Default ctor.
95 }
96
97 //________________________________________________________________________
98 AliAnalysisTaskFlowSingleMu::AliAnalysisTaskFlowSingleMu(const char *name, const AliMuonTrackCuts& cuts) :
99   AliVAnalysisMuon(name, cuts),
100   fEPKeys(0x0),
101   fRandom(0x0)
102 {
103   //
104   /// Constructor.
105   //
106   TString epMethods = "V0A V0C TPC V0Cr1 V0Cr4";
107   fEPKeys = epMethods.Tokenize(" ");
108 }
109
110
111 //________________________________________________________________________
112 AliAnalysisTaskFlowSingleMu::~AliAnalysisTaskFlowSingleMu()
113 {
114   //
115   /// Destructor
116   //
117   delete fEPKeys;
118   delete fRandom;
119 }
120
121 //___________________________________________________________________________
122 void AliAnalysisTaskFlowSingleMu::MyUserCreateOutputObjects()
123 {
124   //
125   /// Create outputs
126   //
127
128   if ( ! fRandom ) fRandom = new TRandom3();
129   // Set seed here, to avoid having the same seed for all jobs in grid
130   fRandom->SetSeed(0);
131
132   Int_t nPtBins = 160;
133   Double_t ptMin = 0., ptMax = 80.;
134   TString ptName("Pt"), ptTitle("p_{t}"), ptUnits("GeV/c");
135   
136   Int_t nEtaBins = 25;
137   Double_t etaMin = -4.5, etaMax = -2.;
138   TString etaName("Eta"), etaTitle("#eta"), etaUnits("");
139   
140   Int_t nPhiBins = 36;
141   Double_t phiMin = 0.; Double_t phiMax = 2.*TMath::Pi();
142   TString phiName("Phi"), phiTitle("#phi"), phiUnits("rad");
143
144   Int_t nDeltaPhiBins = 18;
145   Double_t deltaPhiMin = 0.; Double_t deltaPhiMax = TMath::Pi();
146   TString deltaPhiName("DeltaPhi"), deltaPhiTitle("#phi_{#mu} - #psi_{plane}"), deltaPhiUnits("rad");
147     
148   Int_t nChargeBins = 2;
149   Double_t chargeMin = -2, chargeMax = 2.;
150   TString chargeName("Charge"), chargeTitle("charge"), chargeUnits("e");
151   
152   Int_t nMotherTypeBins = kNtrackSources;
153   Double_t motherTypeMin = -0.5, motherTypeMax = (Double_t)kNtrackSources - 0.5;
154   TString motherType("MotherType"), motherTypeTitle("motherType"), motherTypeUnits("");
155
156   Int_t nbins[kNvars] = {nPtBins, nEtaBins, nPhiBins, nDeltaPhiBins, nChargeBins, nMotherTypeBins};
157   Double_t xmin[kNvars] = {ptMin, etaMin, phiMin, deltaPhiMin, chargeMin, motherTypeMin};
158   Double_t xmax[kNvars] = {ptMax, etaMax, phiMax, deltaPhiMax, chargeMax, motherTypeMax};
159   TString axisTitle[kNvars] = {ptTitle, etaTitle, phiTitle, deltaPhiTitle, chargeTitle, motherTypeTitle};
160   TString axisUnits[kNvars] = {ptUnits, etaUnits, phiUnits, deltaPhiUnits, chargeUnits, motherTypeUnits};
161
162   Int_t iobj=0;
163   TH1* histo = 0x0;
164   TString currTitle = "";
165   for ( Int_t iep=0; iep<fEPKeys->GetEntries(); iep++ ) {
166     TString currEP = fEPKeys->At(iep)->GetName();
167     AliCFContainer* cfContainer = new AliCFContainer(Form("FlowSingleMuContainer%s", currEP.Data()),"Container for tracks",kNsteps,kNvars,nbins);
168   
169     TString currName = "";
170     for ( Int_t idim = 0; idim<kNvars; idim++){
171       currName = Form("%s (%s)", axisTitle[idim].Data(), axisUnits[idim].Data());
172       currName.ReplaceAll("()","");
173
174       cfContainer->SetVarTitle(idim, currName.Data());
175       cfContainer->SetBinLimits(idim, xmin[idim], xmax[idim]);
176     }
177     
178     TString stepTitle[kNsteps] = {"reconstructed", "generated"};
179
180     TAxis* currAxis = 0x0;
181     for (Int_t istep=0; istep<kNsteps; istep++){
182       cfContainer->SetStepTitle(istep, stepTitle[istep].Data());
183       AliCFGridSparse* gridSparse = cfContainer->GetGrid(istep);
184
185       currAxis = gridSparse->GetAxis(kHvarMotherType);
186       for ( Int_t ibin=0; ibin<fSrcKeys->GetEntries(); ibin++ ) {
187         currAxis->SetBinLabel(ibin+1, fSrcKeys->At(ibin)->GetName());
188       }
189     }
190
191     AddObjectToCollection(cfContainer, iobj++);
192
193     Bool_t isTPC = currEP.Contains("TPC");
194     Double_t minPsiPlane = ( isTPC ) ? 0. : -TMath::Pi()/2.;
195     Double_t maxPsiPlane = ( isTPC ) ? TMath::Pi() : TMath::Pi()/2.;
196     
197     histo = new TH1D(Form("hEventPlane%s", currEP.Data()), Form("Event plane distribution: %s", currEP.Data()), 18, minPsiPlane, maxPsiPlane);
198     histo->GetXaxis()->SetTitle(Form("#psi_{%s} (rad)", currEP.Data()));
199     AddObjectToCollection(histo, iobj++);
200
201     for ( Int_t jep=iep+1; jep<fEPKeys->GetEntries(); jep++ ) {
202       TString auxEP = fEPKeys->At(jep)->GetName();
203       currTitle = Form("cos(#psi_{%s} - #psi_{%s})",currEP.Data(),auxEP.Data());
204       histo = new TH1D(Form("hResoSub3%s%s",currEP.Data(),auxEP.Data()), currTitle.Data(), 100, -1.,1.);
205       histo->GetXaxis()->SetTitle(currTitle.Data());
206       AddObjectToCollection(histo, iobj++);
207     }
208
209     currTitle = Form("cos(#psi_{1} - #psi_{2})");
210     histo = new TH1D(Form("hResoSub2%s",currEP.Data()), Form("%s: %s", currEP.Data(), currTitle.Data()), 100, -1.,1.);
211     histo->GetXaxis()->SetTitle(currTitle.Data());
212
213     AddObjectToCollection(histo, iobj++);
214   } // loop on EPs
215
216   fMuonTrackCuts->Print("mask");
217 }
218
219 //________________________________________________________________________
220 void AliAnalysisTaskFlowSingleMu::ProcessEvent(TString physSel, const TObjArray& selectTrigClasses, TString centrality)
221 {
222   //
223   /// Fill output objects
224   //
225
226   //
227   // Base event cuts
228   //
229   if ( GetVertexSPD()->GetNContributors() < fMinNvtxContirbutors ) return;
230   //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
231   //if ( isPileupFromSPD ) return;
232   if ( TMath::Abs(GetVertexSPD()->GetZ()) > 10. ) 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     } // loop on trigger
279   } // loop on EP
280
281   //
282   // Fill track info
283   //
284   Double_t containerInput[kNvars];
285   AliVParticle* track = 0x0;
286   for ( Int_t istep = 0; istep<2; ++istep ) {
287     Int_t nTracks = ( istep == kStepReconstructed ) ? GetNTracks() : GetNMCTracks();
288     for (Int_t itrack = 0; itrack < nTracks; itrack++) {
289       track = ( istep == kStepReconstructed ) ? GetTrack(itrack) : GetMCTrack(itrack);
290       
291       Bool_t isSelected = ( istep == kStepReconstructed ) ? fMuonTrackCuts->IsSelected(track) : ( TMath::Abs(track->PdgCode()) == 13 );
292       if ( ! isSelected ) continue;
293       
294       Int_t trackSrc = ( istep == kStepReconstructed ) ? GetParticleType(track) : RecoTrackMother(track);
295  
296
297       containerInput[kHvarPt]         = track->Pt();
298       containerInput[kHvarEta]        = track->Eta();
299       containerInput[kHvarPhi]        = track->Phi();
300       containerInput[kHvarCharge]     = track->Charge()/3.;
301       containerInput[kHvarMotherType] = (Double_t)trackSrc;
302
303       for ( Int_t iep=0; iep<fEPKeys->GetEntriesFast(); iep++ ) {
304         if ( psiPlane[iep] < kLowestVal ) continue; 
305         TString currEP = fEPKeys->At(iep)->GetName();
306         Double_t deltaPhi = containerInput[kHvarPhi] - psiPlane[iep];
307
308         // The difference between phi - psiPlane must be between (0, pi)
309         // Since the reaction plain is symmetric in pi,
310         // we can do in such a way that:
311         // -pi < delta_phi + N*pi < pi
312         // We choose N accrodingly
313         Double_t nPi = TMath::Floor( 1. - deltaPhi / TMath::Pi() );
314         deltaPhi += nPi * TMath::Pi();
315         containerInput[kHvarDeltaPhi] = deltaPhi;
316
317         for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
318           TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
319           if ( istep == kStepReconstructed && ! TrackPtCutMatchTrigClass(track, trigClassName) ) continue;
320           ((AliCFContainer*)GetMergeableObject(physSel, trigClassName, centrality, Form("FlowSingleMuContainer%s",currEP.Data())))->Fill(containerInput,istep);
321         } // loop on selected trigger classes
322       } // loop on event plane
323     } // loop on tracks
324   } // loop on container steps
325 }
326
327
328 //________________________________________________________________________
329 TArrayD AliAnalysisTaskFlowSingleMu::GetCentralityRange(TString sRange)
330 {
331   //
332   /// Get the range from string
333   //
334   TArrayD centralityRange(2);
335   centralityRange.Reset(-999.);
336   TObjArray* array = sRange.Tokenize("_");
337   if ( array->GetEntries() == 2 ) {
338     for ( Int_t ientry=0; ientry<2; ientry++ ) {
339       centralityRange[ientry] = ((TObjString*)array->At(ientry))->GetString().Atof();
340     }
341   }
342   delete array;
343   return centralityRange;
344 }
345
346 //________________________________________________________________________
347 void AliAnalysisTaskFlowSingleMu::Terminate(Option_t *) {
348   //
349   /// Draw some histograms at the end.
350   //
351
352   AliVAnalysisMuon::Terminate("");
353
354   if ( ! fMergeableCollection ) return;
355   
356   TString physSel = fTerminateOptions->At(0)->GetName();
357   TString trigClassName = fTerminateOptions->At(1)->GetName();
358   TString centralityRange = fTerminateOptions->At(2)->GetName();
359   TString furtherOpt = fTerminateOptions->At(3)->GetName();
360
361   TObjArray* centralityClasses = centralityRange.Tokenize(",");
362   TArrayD centralityArray(centralityClasses->GetEntries()+1);
363   Int_t ib = 0;
364   for ( Int_t icent=0; icent<centralityClasses->GetEntries(); icent++ ) {
365     TArrayD range = GetCentralityRange(centralityClasses->At(icent)->GetName());
366     if ( icent == 0 ) centralityArray.AddAt(range[0],ib++);
367     centralityArray.AddAt(range[1],ib++);
368   }
369   
370   TString selectEP = "";
371   TObjArray resoEParray;
372   resoEParray.SetOwner();
373   
374   TString outFilename = "";
375   TObjArray* optArr = furtherOpt.Tokenize(" ");
376   TString currName = "";
377   TString resoFilename = "";
378   for ( Int_t iopt=0; iopt<optArr->GetEntries(); iopt++ ) {
379     currName = optArr->At(iopt)->GetName();
380     if ( currName.Contains(".root") ) outFilename = currName;
381     else if ( currName.Contains(".txt") ) resoFilename = currName;
382     else {
383       for ( Int_t iep=0; iep<fEPKeys->GetEntries(); iep++ ) {
384         if ( ! currName.CompareTo(fEPKeys->At(iep)->GetName()) ) resoEParray.Add(new TObjString(currName.Data()));
385       }
386     }
387   }
388   delete optArr;
389
390   if ( resoEParray.At(0) ) selectEP = resoEParray.At(0)->GetName();
391
392   furtherOpt.ToUpper();
393
394   Bool_t v2only = furtherOpt.Contains("V2ONLY");
395
396   TAxis* customBins = 0x0;
397   
398   //Double_t myBins[] = {0., 0.5, 1., 1.5, 2., 3., 4., 5., 6., 8., 10.};
399   Double_t myBins[] = {4., 5., 6., 8., 10.}; // REMEMBER TO CHOOSE
400   //Double_t myBins[] = {2., 3., 4., 5., 6., 8., 10.};
401   //Double_t myBins[] = {0., 0.5, 1., 1.5, 2., 3., 4., 10.};
402   //Double_t myBins[] = {2., 3., 4., 5., 6., 8., 10., 80.};
403   Int_t nMyBins = sizeof(myBins)/sizeof(myBins[0])-1;
404   customBins = new TAxis(nMyBins, myBins); // REMEMBER TO CHOOSE
405
406   TCanvas *can = 0x0, *showCan = 0x0;
407   Int_t xshift = 100;
408   Int_t yshift = 20;
409   Int_t igroup1 = 0;
410   Int_t igroup2 = 0;
411
412   //
413   /// Read plane resolution file (if present)
414   //
415   TObjArray resoCentrality;
416   resoCentrality.SetOwner();
417   TArrayD resolution[3];
418   for ( Int_t ires=0; ires<3; ires++ ) {
419     resolution[ires].Set(50);
420   }
421
422   Bool_t externalReso = ( ! resoFilename.IsNull() && ! gSystem->AccessPathName(resoFilename.Data()) );
423   Bool_t internalReso = kTRUE;
424   if ( externalReso ) {
425     // If plane resolution file exists, fill resolution
426     TString currLine = "";
427     ifstream inFile(resoFilename.Data());
428     if (inFile.is_open()) {
429       ib = 0;
430       while (! inFile.eof() ) {
431         currLine.ReadLine(inFile,kFALSE);
432         if ( currLine.BeginsWith("#") || currLine.Length() < 3 ) continue;
433         TObjArray* array = currLine.Tokenize(" ");
434         resoCentrality.Add(new TObjString(array->At(0)->GetName()));
435         for ( Int_t ires=0; ires<3; ires++ ) {
436           resolution[ires].AddAt(((TObjString*)array->At(ires+1))->GetString().Atof(),ib);
437         }
438         delete array;
439         ib++;
440       } // ! EOF
441     } // file is open
442     for ( Int_t ires=0; ires<3; ires++ ) {
443       resolution[ires].Set(ib);
444     }
445   } // file exists
446   else {
447     // If no plane resolution file exist,
448     // use the centralities in terminateOption
449     // and set the resolution to 1 and errors to 0.
450
451     for ( Int_t ires=0; ires<3; ires++ ) {
452       resolution[ires].Set(centralityClasses->GetEntries());
453       Double_t resetVal = ( ires == 0 ) ? 1. : 0.;
454       resolution[ires].Reset(resetVal);
455     }
456     TArrayD currCos(3);
457     for ( Int_t icent=0; icent<centralityClasses->GetEntries(); icent++ ) {
458       resoCentrality.Add(new TObjString(centralityClasses->At(icent)->GetName()));
459       Int_t icos = 0;
460       Double_t sumRelErrSquare = 0.;
461       for ( Int_t isub=0; isub<resoEParray.GetEntries(); isub++ ) {
462         for ( Int_t jsub=0; jsub<resoEParray.GetEntries(); jsub++ ) {
463           TH1* histo = static_cast<TH1*>(GetSum(physSel,trigClassName,resoCentrality.At(icent)->GetName(),Form("hResoSub3%s%s", resoEParray.At(isub)->GetName(),resoEParray.At(jsub)->GetName())));
464           if ( ! histo || histo->Integral() == 0 ) continue;
465           currCos[icos] = histo->GetMean();
466           if ( currCos[icos] <= 0. ) continue;
467           Double_t relErr = histo->GetMeanError() / currCos[icos];
468           sumRelErrSquare += relErr * relErr;
469           icos++;
470         } // loop on sub events
471       } // loop on sub events
472       if ( icos != 3  ) {
473         printf("Warning: resolution cannot be estimated for %s\n", resoCentrality.At(icent)->GetName());
474         internalReso = kFALSE;
475         continue;
476       }
477       resolution[0][icent] = TMath::Sqrt((currCos[0]*currCos[1]/currCos[2]));
478       resolution[1][icent] = resolution[0][icent]/2.*TMath::Sqrt(sumRelErrSquare);
479       printf("Resolution %s  %g +- %g\n", centralityClasses->At(icent)->GetName(), resolution[0][icent], resolution[1][icent]);
480     }
481   }
482
483   Bool_t hasResolution = ( externalReso || internalReso );
484
485
486   TArrayD resoCentrArray(resoCentrality.GetEntries()+1);
487   ib = 0;
488   for ( Int_t icent=0; icent<resoCentrality.GetEntries(); icent++ ) {
489     TArrayD range = GetCentralityRange(resoCentrality.At(icent)->GetName());
490     if ( icent == 0 ) resoCentrArray.AddAt(range[0],ib++);
491     resoCentrArray.AddAt(range[1], ib++);
492   }
493
494   if ( hasResolution ) {
495     currName = externalReso ? Form("externalReso") : Form("reso_%s_%s_%s", resoEParray.At(0)->GetName(), resoEParray.At(1)->GetName(), resoEParray.At(2)->GetName());
496     can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
497     TGraphErrors* resoGraph = new TGraphErrors(centralityClasses->GetEntries());
498     for ( Int_t icent=0; icent<resoCentrality.GetEntries(); icent++ ) {
499
500       resoGraph->SetPoint(icent, 0.5*(resoCentrArray[icent+1] + resoCentrArray[icent]), resolution[0][icent]);
501       resoGraph->SetPointError(icent, 0.5*(resoCentrArray[icent+1] - resoCentrArray[icent]), resolution[1][icent]);
502     }
503     resoGraph->SetTitle(Form("Resolution %s", selectEP.Data()));
504     resoGraph->GetXaxis()->SetTitle("Centrality");
505     resoGraph->GetYaxis()->SetTitle("Resolution");
506     resoGraph->Draw("apz");
507   }
508
509   const Int_t kNresoCentr = resoCentrality.GetEntries();
510
511   //
512   // Create gridSparse array
513   //
514   AliCFGridSparse* gridSparseArray[kNsteps*kNresoCentr];
515   TString stepName[2];
516   for ( Int_t icent=0; icent<kNresoCentr; icent++ ) {
517     for ( Int_t istep=0; istep<kNsteps; istep++ ) {
518       gridSparseArray[istep*kNresoCentr+icent] = 0x0;
519     }
520
521     AliCFContainer* cfContainer = static_cast<AliCFContainer*> ( GetSum(physSel,trigClassName,resoCentrality.At(icent)->GetName(),Form("FlowSingleMuContainer%s", selectEP.Data())) );
522     if ( ! cfContainer ) continue;
523     for ( Int_t istep=0; istep<kNsteps; istep++ ) {
524       stepName[istep] = cfContainer->GetStepTitle(istep);
525       gridSparseArray[istep*kNresoCentr+icent] = cfContainer->GetGrid(istep);
526     } // loop on step
527   } // loop on centrality
528
529   Int_t srcColors[kNtrackSources] = {kBlack, kRed, kSpring, kTeal, kBlue, kViolet, kMagenta, kOrange};
530   
531   Bool_t isMC = furtherOpt.Contains("MC");
532   Int_t firstSrc = ( isMC ) ? 0 : kUnidentified;
533   Int_t lastSrc  = ( isMC ) ? kNtrackSources - 1 : kUnidentified;
534   if ( ! isMC ) srcColors[kUnidentified] = 1;
535
536   TList outList;
537
538   TString graphTypeName[3] = {"raw","correctStat","correctSyst"};
539   Int_t drawOrder[3] = {0, 2, 1};
540   TString drawOpt[3] = {"apz","pz","a2"};
541   Int_t nTypes = ( hasResolution ) ? 3 : 1;
542
543   TString histoName = "", histoPattern = "";
544   ///////////
545   // Flow  //
546   ///////////
547   gStyle->SetOptFit(1111);
548   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)))";
549   //
550   TF1* func = new TF1("funcFlow",funcFormula.Data(), 0., TMath::Pi());
551   if  ( v2only ) func->SetParNames("scale", "v2");
552   else func->SetParNames("scale","v1", "v2", "v3");
553   Int_t v2par = ( v2only ) ? 1 : 2;
554
555   for ( Int_t istep=0; istep<kNsteps; istep++ ) {
556     for ( Int_t isrc = firstSrc; isrc <= lastSrc; ++isrc ) {
557       TString baseName = Form("%s_%s", stepName[istep].Data(), fSrcKeys->At(isrc)->GetName());
558       TGraphErrors* flowVsCentrality[3];
559       for ( Int_t itype=0; itype<nTypes; itype++ ) {
560         histoName = Form("v2VsCentrality_%s_%s", baseName.Data(), graphTypeName[itype].Data());
561         flowVsCentrality[itype] = new TGraphErrors();
562         flowVsCentrality[itype]->SetName(histoName.Data());
563         flowVsCentrality[itype]->SetTitle(histoName.Data());
564       }
565       for ( Int_t icent=0; icent<centralityClasses->GetEntries(); icent++ ) {
566         TH2* hDeltaPhi2D = 0x0;
567         TString baseNameCent = Form("%s_%s", baseName.Data(), centralityClasses->At(icent)->GetName());
568         TArrayD weightedReso(3);
569         weightedReso.Reset(0.);
570         Double_t nMuons = 0.;
571         TString debugString = Form("\n%s", centralityClasses->At(icent)->GetName());
572         for ( Int_t rcent=0; rcent<kNresoCentr; rcent++ ) {
573           if ( resoCentrArray[rcent] < centralityArray[icent] ||
574                resoCentrArray[rcent+1] > centralityArray[icent+1] ) continue;
575           AliCFGridSparse* gridSparse = gridSparseArray[istep*kNresoCentr+rcent];
576           if ( ! gridSparse ||  gridSparse->GetEntries() == 0. ) continue;
577
578           SetSparseRange(gridSparse, kHvarEta, "", -3.999, -2.501);
579           SetSparseRange(gridSparse, kHvarMotherType, "", isrc+1, isrc+1, "USEBIN");
580           TH1* auxHisto = gridSparse->Project(kHvarDeltaPhi,kHvarPt);
581
582           Double_t currYield = auxHisto->Integral();
583           if ( currYield == 0. ) continue;
584           nMuons += currYield;
585           debugString += Form("\nAdding %s  yield %g  reso", resoCentrality.At(rcent)->GetName(), currYield);
586           for ( Int_t ires=0; ires<3; ires++ ) {
587             weightedReso[ires] += resolution[ires][rcent] * currYield;
588             debugString += Form("  %g", resolution[ires][rcent]);
589           }
590           histoName = Form("ptVsDeltaPhi_%s", baseNameCent.Data());
591           if ( ! hDeltaPhi2D ) hDeltaPhi2D = static_cast<TH2*>(auxHisto->Clone(histoName.Data()));
592           else hDeltaPhi2D->Add(auxHisto);
593           delete auxHisto;
594         }
595
596         if ( ! hDeltaPhi2D ) continue;
597
598         debugString += Form("\nWeighted reso ");
599         for ( Int_t ires=0; ires<3; ires++ ) {
600           weightedReso[ires] = weightedReso[ires]/nMuons;
601           debugString += Form("  %g", weightedReso[ires]);
602         }
603         AliDebug(1,debugString.Data());
604
605         TAxis* ptAxis = hDeltaPhi2D->GetYaxis();
606
607         Int_t ipad = 0;
608         TGraphErrors* flowVsPt[3];
609         for ( Int_t itype=0; itype<nTypes; itype++ ) {
610           histoName = Form("v2VsPt_%s_%s", baseNameCent.Data(), graphTypeName[itype].Data());
611           flowVsPt[itype] = new TGraphErrors();
612           flowVsPt[itype]->SetName(histoName.Data());
613           flowVsPt[itype]->SetTitle(histoName.Data());
614         }
615
616         if ( ! customBins ) customBins = static_cast<TAxis*>(ptAxis->Clone());
617
618
619         for ( Int_t ipt=0; ipt<=customBins->GetNbins(); ipt++ ) {
620           Int_t minCustomBin = ( ipt == 0 ) ? 1 : ipt;
621           Int_t maxCustomBin = ( ipt == 0 ) ? customBins->GetNbins() : ipt;
622           Int_t ptMinBin = ptAxis->FindBin(customBins->GetBinLowEdge(minCustomBin)+1.e-4);
623           Int_t ptMaxBin = ptAxis->FindBin(customBins->GetBinUpEdge(maxCustomBin)-1.e-4);
624           Double_t ptMin = ptAxis->GetBinLowEdge(ptMinBin);
625           Double_t ptMax = ptAxis->GetBinUpEdge(ptMaxBin);
626           TH1* projHisto = hDeltaPhi2D->ProjectionX(Form("checkHisto_%s_ptBin%i", baseNameCent.Data(), ipt), ptMinBin, ptMaxBin);
627           projHisto->SetTitle(Form("%g < p_{t} (GeV/c) < %g", ptMin, ptMax));
628           if ( projHisto->Integral() < 50 ) break;
629           if ( ipad%4 == 0 ) {
630             currName = projHisto->GetName();
631             currName.Append("_can");
632             showCan = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
633             showCan->Divide(2,2);
634             ipad = 0;
635           }
636           showCan->cd(++ipad);
637           if ( v2only ) func->SetParameters(projHisto->Integral(), 0.01);
638           else {
639             func->SetParameters(projHisto->Integral(), 0., 0.01, 0.);
640             func->FixParameter(1,0.);
641           }
642           projHisto->Fit(func,"R");
643           for ( Int_t itype=0; itype<nTypes; itype++ ) {
644             TGraphErrors* currGraph = ( ipt == 0 ) ? flowVsCentrality[itype] : flowVsPt[itype];
645             Double_t resoVal = ( itype == 0 ) ? 1. : weightedReso[0];
646             Double_t resoErr = ( itype == 0 ) ? 0. : weightedReso[itype];
647             Double_t rawVal = func->GetParameter(v2par);
648             Double_t rawErr = (itype == 2 ) ? 0. : func->GetParError(v2par);
649             Double_t yVal = rawVal/resoVal;
650             Double_t xVal = ( ipt == 0 ) ? 0.5*(centralityArray[icent+1] + centralityArray[icent]) : 0.5*(ptMax + ptMin);
651             Double_t xErr = ( ipt == 0 ) ? 0.5*(centralityArray[icent+1] - centralityArray[icent]) : 0.5*(ptMax - ptMin);
652             Int_t ipoint = currGraph->GetN();
653             currGraph->SetPoint(ipoint,xVal,yVal);
654             if ( itype > 0 && ipt != 0 ) {
655               // For the plot vs. pt the resolution error is fully correlated.
656               // Hence we do not use add it bin by bin, but rather write it
657               // on the title
658               resoErr = 0.;
659               TString graphTitle = currGraph->GetTitle();
660               if ( ! graphTitle.Contains("|") ) {
661                 graphTitle += Form(" | Rel. error on EP resolution: %.2g (stat.) %.2g (syst.)", weightedReso[1]/weightedReso[0], weightedReso[2]/weightedReso[0]);
662                 currGraph->SetTitle(graphTitle.Data());
663               }
664             }
665             Double_t rawErrRel = ( rawVal != 0. ) ? rawErr/rawVal : 0.;
666             Double_t resoErrRel = ( resoVal != 0. ) ? resoErr/resoVal : 0.;
667             Double_t yErr = TMath::Abs(yVal) * TMath::Sqrt(rawErrRel*rawErrRel+resoErrRel*resoErrRel);
668             currGraph->SetPointError(ipoint,xErr,yErr);
669           } // loop on type
670         } // loop on pt bins
671         for ( Int_t itype=0; itype<nTypes; itype++ ) {
672           Int_t currType = drawOrder[itype];
673           if ( itype < 2 ) {
674             currName = flowVsPt[currType]->GetName();
675             currName.Append("Can");
676             can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
677             igroup2++;
678             can->cd();
679           }
680           flowVsPt[currType]->GetXaxis()->SetTitle(ptAxis->GetTitle());
681           flowVsPt[currType]->GetYaxis()->SetTitle("v2");
682           flowVsPt[currType]->GetYaxis()->SetRangeUser(0., 0.25);
683           flowVsPt[currType]->SetFillStyle(0);
684           flowVsPt[currType]->Draw(drawOpt[currType].Data());
685           outList.Add(flowVsPt[currType]);
686         } // loop on types
687       } // loop on centrality
688       for ( Int_t itype=0; itype<nTypes; itype++ ) {
689         Int_t currType = drawOrder[itype];
690         if ( flowVsCentrality[currType]->GetN() == 0 ) continue;
691         if ( itype < 2 ) {
692           currName = flowVsCentrality[currType]->GetName();
693           currName.Append("Can");
694           can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
695           can->cd();
696         }
697         flowVsCentrality[currType]->GetXaxis()->SetTitle("Centrality");
698         flowVsCentrality[currType]->GetYaxis()->SetTitle("v2");
699         flowVsCentrality[currType]->SetFillStyle(0);
700         flowVsCentrality[currType]->GetYaxis()->SetRangeUser(0., 0.25);
701         flowVsCentrality[currType]->Draw(drawOpt[currType].Data());
702         outList.Add(flowVsCentrality[currType]);
703         igroup2++;
704       } // loop on type
705       igroup1++;
706       igroup2 = 0;
707     } // loop on track sources      
708   } // loop on steps
709
710   delete customBins;
711
712
713   ///////////////////////
714   // Event plane check //
715   ///////////////////////
716   igroup1++;
717   igroup2 = 0;
718   Int_t icolor=0;
719   currName ="eventPlaneDistrib";
720   can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
721   can->SetGridy();
722   igroup2++;
723
724   Bool_t addOffsetToCentrality = kFALSE;
725
726   TLegend* leg = new TLegend(0.7,0.7,0.92,0.92);
727   leg->SetBorderSize(1);
728   for ( Int_t icent=0; icent<centralityClasses->GetEntries(); icent++ ) {
729     TH1* histo = static_cast<TH1*> ( GetSum(physSel,trigClassName,centralityClasses->At(icent)->GetName(),Form("hEventPlane%s", selectEP.Data())) );
730     if ( ! histo ) continue;
731     histo->SetName(Form("%s_%s", histo->GetName(), centralityClasses->At(icent)->GetName()));
732     if ( histo->Integral() < 50. ) continue;
733     if ( histo->GetSumw2N() == 0 ) histo->Sumw2();
734     histo->Fit("pol0","R0");
735     TF1* pol0 = histo->GetFunction("pol0");
736     histo->Scale(1./pol0->GetParameter(0));
737     Int_t offset = ( addOffsetToCentrality ) ? icolor : 0;
738     TGraphErrors* graph = new TGraphErrors();
739     graph->SetTitle(Form("%s %s", histo->GetTitle(), trigClassName.Data()));
740     for ( Int_t ipoint=0; ipoint<histo->GetXaxis()->GetNbins(); ipoint++ ) {
741       Int_t ibin = ipoint+1;
742       graph->SetPoint(ipoint, histo->GetXaxis()->GetBinCenter(ibin), histo->GetBinContent(ibin) + offset);
743       graph->SetPointError(ipoint, histo->GetXaxis()->GetBinWidth(ibin)/2., histo->GetBinError(ibin));
744     }
745     graph->GetXaxis()->SetTitle(histo->GetXaxis()->GetTitle());
746     graph->GetYaxis()->SetTitle("Data/Fit (pol0)");
747     TString legTitle = Form("%s %%", centralityClasses->At(icent)->GetName());
748     if ( addOffsetToCentrality ) {
749       graph->GetYaxis()->SetRangeUser(0, 1.5*(Double_t)centralityClasses->GetEntries());
750       graph->GetYaxis()->SetNdivisions(2*centralityClasses->GetEntries(),0,0);
751       legTitle += Form(" ( + %i)", offset);
752     }
753     Int_t currColor = ( icolor < kNtrackSources ) ? srcColors[icolor] : 20+icolor;
754     graph->SetLineColor(currColor);
755     graph->SetMarkerColor(currColor);
756     graph->SetMarkerStyle(20+icolor);
757     TString currDrawOpt = "pz";
758     if ( icolor == 0 ) currDrawOpt += "a";
759     graph->Draw(currDrawOpt.Data());
760     legTitle.ReplaceAll("_"," - ");
761     leg->AddEntry(graph,legTitle.Data(),"lp");
762     if ( addOffsetToCentrality ) {
763       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()));
764       latex->SetTextColor(currColor);
765       latex->SetTextSize(0.025);
766       latex->Draw("same");
767     }
768     icolor++;
769   }
770   leg->Draw("same");
771
772
773   //////////////////////
774   // Event statistics //
775   //////////////////////
776   printf("\nTotal analyzed events:\n");
777   TString evtSel = Form("trigger:%s", trigClassName.Data());
778   fEventCounters->PrintSum(evtSel.Data());
779   printf("Physics selected analyzed events:\n");
780   evtSel = Form("trigger:%s/selected:yes", trigClassName.Data());
781   fEventCounters->PrintSum(evtSel.Data());
782   
783   TString countPhysSel = "any";
784   if ( physSel.Contains(fPhysSelKeys->At(kPhysSelPass)->GetName()) ) countPhysSel = "yes";
785   else if ( physSel.Contains(fPhysSelKeys->At(kPhysSelReject)->GetName()) ) countPhysSel="no";
786   countPhysSel.Prepend("selected:");
787   printf("Analyzed events vs. centrality:\n");
788   evtSel = Form("trigger:%s/%s", trigClassName.Data(), countPhysSel.Data());
789   fEventCounters->Print("centrality",evtSel.Data(),kTRUE);
790
791   if ( ! outFilename.IsNull() ) {
792     printf("\nWriting output file %s\n", outFilename.Data());
793     TFile* file = TFile::Open(outFilename.Data(), "RECREATE");
794     outList.Write();
795     file->Close();
796   }
797 }