Improve plotted histograms in Terminate (Diego)
[u/mrichter/AliRoot.git] / PWGPP / MUON / lite / AliAnalysisTaskMuonCuts.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: AliAnalysisTaskMuonCuts.cxx 47782 2011-02-24 18:37:31Z martinez $ */
17
18 /// Task for tuning cuts for single muons in the spectrometer.
19 /// The output is a list of histograms.
20 /// Task need to be run in steps:
21 /// Step 1: find <DCAx> and <DCAy> in MC
22 /// Step 2: find signal sigma_pDCA in MC
23 /// Step 3: find <DCAx> and <DCAy> in data
24 /// Step 4: apply sigma_pDCA cut on simulations and check the
25 ///         resulting kinematic variables
26 ///
27 /// \author Diego Stocco
28
29 #define AliAnalysisTaskMuonCuts_cxx
30
31 #include "AliAnalysisTaskMuonCuts.h"
32
33 // ROOT includes
34 #include "TROOT.h"
35 #include "TH1.h"
36 #include "TH2.h"
37 #include "TH3.h"
38 #include "TAxis.h"
39 #include "TCanvas.h"
40 #include "TLegend.h"
41 #include "TMath.h"
42 #include "TObjString.h"
43 #include "TObjArray.h"
44 #include "TF1.h"
45 #include "TStyle.h"
46 #include "TArrayI.h"
47 //#include "TMCProcess.h"
48
49 // STEER includes
50 #include "AliAODEvent.h"
51 #include "AliAODTrack.h"
52 #include "AliAODMCParticle.h"
53 #include "AliMCEvent.h"
54 #include "AliMCParticle.h"
55 //#include "AliStack.h"
56 #include "AliESDEvent.h"
57 #include "AliESDMuonTrack.h"
58 #include "AliInputEventHandler.h"
59
60 #include "AliAnalysisManager.h"
61
62 // PWG includes
63 #include "AliMergeableCollection.h"
64 #include "AliMuonEventCuts.h"
65 #include "AliMuonTrackCuts.h"
66 #include "AliAnalysisMuonUtility.h"
67
68
69 /// \cond CLASSIMP
70 ClassImp(AliAnalysisTaskMuonCuts) // Class implementation in ROOT context
71 /// \endcond
72
73
74 //________________________________________________________________________
75 AliAnalysisTaskMuonCuts::AliAnalysisTaskMuonCuts() :
76   AliVAnalysisMuon(),
77   fHistoTypeKeys(0x0),
78   fThetaAbsKeys(0x0),
79   fSigmaCuts(TArrayD())
80 {
81   /// Default ctor.
82 }
83
84 //________________________________________________________________________
85 AliAnalysisTaskMuonCuts::AliAnalysisTaskMuonCuts(const char *name, const AliMuonTrackCuts& cuts) :
86   AliVAnalysisMuon(name, cuts),
87   fHistoTypeKeys(0x0),
88   fThetaAbsKeys(0x0),
89   fSigmaCuts(TArrayD())
90 {
91   //
92   /// Constructor.
93   //
94   TString histoTypeKeys = "hDCAxVsP hDCAyVsP hPdcaVsP hPDCAVsPCheck hDCAVsPCheck hChiProbVsP hSigmaVsPt hSigmaVsEta";
95   fHistoTypeKeys = histoTypeKeys.Tokenize(" ");
96
97   TString thetaAbsKeys = "ThetaAbs23 ThetaAbs310";
98   fThetaAbsKeys = thetaAbsKeys.Tokenize(" ");
99
100   SetSigmaCuts();
101 }
102
103
104 //________________________________________________________________________
105 AliAnalysisTaskMuonCuts::~AliAnalysisTaskMuonCuts()
106 {
107   //
108   /// Destructor
109   //
110
111   delete fHistoTypeKeys;
112   delete fThetaAbsKeys;
113 }
114
115 //___________________________________________________________________________
116 void AliAnalysisTaskMuonCuts::MyUserCreateOutputObjects()
117 {
118   TH1* histo = 0x0;
119   TString histoName = "";
120   for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
121     for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
122       histoName = GetHistoName(kDCAxVsP, itheta, isrc);
123       histo = new TH2F(histoName.Data(), histoName.Data(), 50, 0., 200., 100, -50., 50.);
124       histo->SetXTitle("p (GeV/c)");
125       histo->SetYTitle("DCA_{x} (cm)");
126       AddObjectToCollection(histo);
127
128       histoName = GetHistoName(kDCAyVsP, itheta, isrc);
129       histo = new TH2F(histoName.Data(), histoName.Data(), 50, 0., 200., 100, -50., 50.);
130       histo->SetXTitle("p (GeV/c)");
131       histo->SetYTitle("DCA_{y} (cm)");
132       AddObjectToCollection(histo);
133
134       histoName = GetHistoName(kPdcaVsP, itheta, isrc);
135       histo = new TH2F(histoName.Data(), histoName.Data(), 100, 0., 400., 50, 0., 500.);
136       histo->SetXTitle("p (GeV/c)");
137       histo->SetYTitle("p#timesDCA (cm #times GeV/c)");
138       AddObjectToCollection(histo);
139
140       histoName = GetHistoName(kPDCAVsPCheck, itheta, isrc);
141       histo = new TH2F(histoName.Data(), histoName.Data(), 100, 0., 800., 100, 0., 5000.);
142       histo->SetXTitle("p (GeV/c)");
143       histo->SetYTitle("p#timesDCA (cm #times GeV/c)");
144       AddObjectToCollection(histo);
145
146       histoName = GetHistoName(kDCAVsPCheck, itheta, isrc);
147       histo = new TH2F(histoName.Data(), histoName.Data(), 100, 0., 800., 400, 0., 200.);
148       histo->SetXTitle("p (GeV/c)");
149       histo->SetYTitle("DCA (cm)");
150       AddObjectToCollection(histo);
151
152       histoName = GetHistoName(kChiProbVsP, itheta, isrc);
153       histo = new TH2F(histoName.Data(), histoName.Data(), 50, 0., 200., 50, 0., 1.);
154       histo->SetXTitle("p (GeV/c)");
155       histo->SetYTitle("Chisquare prob.");
156       AddObjectToCollection(histo);
157
158       histoName = GetHistoName(kSigmaVsPt, itheta, isrc);
159       histo = new TH2F(histoName.Data(), histoName.Data(), 100, 0., 100., fSigmaCuts.GetSize(), 0.5, 0.5+(Double_t)fSigmaCuts.GetSize());
160       histo->SetXTitle("p_{t} (GeV/c)");
161       histo->SetYTitle("#sigma_{p#timesDCA}");
162       for ( Int_t ibin=0; ibin<fSigmaCuts.GetSize(); ++ibin ) {
163         histo->GetYaxis()->SetBinLabel(ibin+1,Form("%g", fSigmaCuts[ibin]));
164       }
165       AddObjectToCollection(histo);
166
167       histoName = GetHistoName(kSigmaVsEta, itheta, isrc);
168       histo = new TH2F(histoName.Data(), histoName.Data(), 25, -4.5, -2., fSigmaCuts.GetSize(), 0.5, 0.5+(Double_t)fSigmaCuts.GetSize());
169       histo->SetXTitle("#eta");
170       histo->SetYTitle("#sigma_{p#timesDCA}");
171       for ( Int_t ibin=0; ibin<fSigmaCuts.GetSize(); ++ibin ) {
172         histo->GetYaxis()->SetBinLabel(ibin+1,Form("%g", fSigmaCuts[ibin]));
173       }
174       AddObjectToCollection(histo);
175     } // loop on track sources
176   } // loop on theta abs
177
178   fMuonTrackCuts->Print();
179
180 }
181
182 //________________________________________________________________________
183 TString AliAnalysisTaskMuonCuts::GetHistoName(Int_t histoTypeIndex, Int_t thetaAbsIndex, Int_t srcIndex)
184 {
185   /// Get local histogram name
186   TString histoName = Form("%s%s%s", fHistoTypeKeys->At(histoTypeIndex)->GetName(), fThetaAbsKeys->At(thetaAbsIndex)->GetName(), fSrcKeys->At(srcIndex)->GetName());
187
188   return histoName;
189 }
190
191 //________________________________________________________________________
192 void AliAnalysisTaskMuonCuts::ProcessEvent(TString physSel, const TObjArray& selectTrigClasses, TString centrality)
193 {
194   //
195   /// Fill histogram
196   //
197
198   TString histoName = "";
199   AliVParticle* track = 0x0;
200   Int_t nTracks = AliAnalysisMuonUtility::GetNTracks(InputEvent());
201   for (Int_t itrack = 0; itrack < nTracks; itrack++) {
202     track = AliAnalysisMuonUtility::GetTrack(itrack, InputEvent());
203     fMuonTrackCuts->CustomParam()->SetNSigmaPdca(1.e10);
204     if ( ! fMuonTrackCuts->IsSelected(track) ) continue;
205
206     Double_t thetaAbsEndDeg = AliAnalysisMuonUtility::GetThetaAbsDeg(track);
207     Int_t thetaAbsBin = ( thetaAbsEndDeg < 3. ) ? kThetaAbs23 : kThetaAbs310;
208
209     Int_t trackSrc = GetParticleType(track);
210
211     TVector3 dcaAtVz = fMuonTrackCuts->GetCorrectedDCA(track);
212     Double_t pTotMean = fMuonTrackCuts->GetAverageMomentum(track);
213     Double_t dca = dcaAtVz.Mag();
214     Double_t pDca = pTotMean * dca;
215
216     Double_t sigmaPdca = fMuonTrackCuts->IsThetaAbs23(track) ? fMuonTrackCuts->GetMuonTrackCutsParam().GetSigmaPdca23() : fMuonTrackCuts->GetMuonTrackCutsParam().GetSigmaPdca310();
217     Double_t chi2 = pDca / sigmaPdca;
218     chi2 *= chi2;
219     Double_t chiProb = TMath::Prob(chi2, 2);
220
221     Double_t pTot = track->P();
222     Double_t pt = track->Pt();
223     Double_t eta = track->Eta();
224
225     for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
226       TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
227       histoName = GetHistoName(kDCAxVsP, thetaAbsBin, trackSrc);
228       ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(pTot, dcaAtVz.X());
229
230       histoName = GetHistoName(kDCAyVsP, thetaAbsBin, trackSrc);
231       ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(pTot, dcaAtVz.Y());
232   
233       histoName = GetHistoName(kPdcaVsP, thetaAbsBin, trackSrc);
234       ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(pTot, pDca);
235
236       histoName = GetHistoName(kPDCAVsPCheck, thetaAbsBin, trackSrc);
237       ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(pTot, pDca);
238
239       histoName = GetHistoName(kDCAVsPCheck, thetaAbsBin, trackSrc);
240       ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(pTot, dca);
241
242       histoName = GetHistoName(kChiProbVsP, thetaAbsBin, trackSrc);
243       ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(pTot, chiProb);
244     } // loop on selected trigger classes
245
246     for ( Int_t isigma=0; isigma<fSigmaCuts.GetSize(); ++isigma) {
247       fMuonTrackCuts->CustomParam()->SetNSigmaPdca(fSigmaCuts[isigma]);
248       if ( ! fMuonTrackCuts->IsSelected(track) ) continue;
249       for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
250         TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
251         histoName = GetHistoName(kSigmaVsPt, thetaAbsBin, trackSrc);
252         ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(pt, isigma+1);
253         histoName = GetHistoName(kSigmaVsEta, thetaAbsBin, trackSrc);
254         ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(eta, isigma+1);
255       } // loop on selected trigger classes 
256     } // loop on sigmas
257   }
258 }
259
260 //________________________________________________________________________
261 void AliAnalysisTaskMuonCuts::SetSigmaCuts(Int_t nSigmaCuts, Double_t* sigmaCuts)
262 {
263   /// Set number of sigmas
264   if ( ! sigmaCuts || nSigmaCuts < 0 ) {
265 //     if ( defaultChiSquare ) {
266 //       Double_t cuts[] = {0.3, 0.1, 0.05, 0.03, 0.01, 0.005, 0.003, 0.001, 0.0005, 0.0003, 0.0001, 0.};
267 //       Int_t nCuts = sizeof(cuts)/sizeof(cuts[0]);
268 //       fSigmaCuts.Set(nCuts, cuts);
269 //     }
270 //     else {
271       Double_t cuts[] = {2., 3., 4., 5., 6., 7., 10., 15., 20., 25., 30., 1.e10};
272       Int_t nCuts = sizeof(cuts)/sizeof(cuts[0]);
273       fSigmaCuts.Set(nCuts, cuts);
274       //    }
275   }
276   else {
277     fSigmaCuts.Set(nSigmaCuts, sigmaCuts);
278   }
279 }
280
281 //________________________________________________________________________
282 void AliAnalysisTaskMuonCuts::Terminate(Option_t *) {
283   //
284   /// Draw some histogram at the end.
285   //
286
287   AliVAnalysisMuon::Terminate("");
288
289   if ( ! fMergeableCollection ) return;
290
291   TString physSel = fTerminateOptions->At(0)->GetName();
292   TString trigClassName = fTerminateOptions->At(1)->GetName();
293   TString centralityRange = fTerminateOptions->At(2)->GetName();
294   TString furtherOpt = fTerminateOptions->At(3)->GetName();
295   furtherOpt.ToUpper();
296
297   furtherOpt.ReplaceAll("  ", " ");
298   furtherOpt.ReplaceAll(" =", "=");
299   furtherOpt.ReplaceAll("= ", "=");
300   TObjArray* optArray = furtherOpt.Tokenize(" ");
301   Double_t refSigmaCut = 6.;
302   for ( Int_t iopt=0; iopt<optArray->GetEntries(); ++iopt ) {
303     TString currOpt = optArray->At(iopt)->GetName();
304     if ( currOpt.Contains("REFSIGMA") ) {
305       currOpt.Remove(0,currOpt.Index("=")+1);
306       refSigmaCut = currOpt.Atof();
307     }
308   }
309   delete optArray;
310   
311   fMuonTrackCuts->Print("param");
312
313   Int_t srcColors[kNtrackSources] = {kBlack, kRed, kSpring, kTeal, kBlue, kViolet, kMagenta, kOrange, kGray};
314
315   TCanvas* can = 0x0;
316   Int_t xshift = 100;
317   Int_t yshift = 100;
318   Int_t igroup1 = -1;
319   Int_t igroup2 = 0;
320
321   //////////////
322   // Reco DCA //
323   //////////////
324   igroup1++;
325   igroup2 = 0;
326   TString histoName = "", currName = "", histoPattern = "", drawOpt = "";
327   currName = Form("%s_recoDCA", GetName());
328   can = new TCanvas(currName.Data(),"Reco DCA",igroup1*xshift,igroup2*yshift,600,600);
329   can->Divide(2,2);
330   igroup2++;
331   Int_t recoDcaHisto[2] = {kDCAxVsP, kDCAyVsP};
332   TString dcaName[2] = {"DCAx", "DCAy"};
333   Double_t averageDca[4] = {0., 0.};
334   printf("\nAverage reconstructed DCA:\n");
335   TF1* fitFuncMeanDca = new TF1("fitFuncMeanDca","gausn",-20.,20.);
336   fitFuncMeanDca->SetParNames("Norm", "Mean", "Sigma");
337   for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
338     for ( Int_t ihisto=0; ihisto<2; ++ihisto ) {
339       TH2* histo = 0x0;
340       histoPattern = "";
341       histoPattern = Form("%s%s*", fHistoTypeKeys->At(recoDcaHisto[ihisto])->GetName(), fThetaAbsKeys->At(itheta)->GetName());
342       histo = (TH2*)GetSum(physSel, trigClassName, centralityRange, histoPattern);
343       if ( ! histo ) continue;
344
345       TH1* meanDcaVsP = histo->ProjectionX(Form("mean%sVsP_%s", dcaName[ihisto].Data(), fThetaAbsKeys->At(itheta)->GetName()));
346       meanDcaVsP->Reset();
347       meanDcaVsP->SetTitle(Form("Mean %s vs. p %s", dcaName[ihisto].Data(), fThetaAbsKeys->At(itheta)->GetName()));
348       meanDcaVsP->SetYTitle(Form("<%s> (cm)", dcaName[ihisto].Data()));
349       meanDcaVsP->SetStats(kFALSE);
350
351       Int_t nPbins = histo->GetXaxis()->GetNbins();
352       //Int_t nPadX = (Int_t)TMath::Sqrt(nPbins);
353       //Int_t nPadY = nPadX;
354       //if ( nPadX * nPadY < nPbins ) nPadX++;
355       TCanvas* meanDcaFitCan = 0x0;
356       Int_t nPadX = 5;
357       Int_t nPadY = 5;
358       Int_t ipad = 0;
359       Int_t ican = 0;
360
361       for ( Int_t ibin=2; ibin<=nPbins; ++ibin ) {
362         currName = Form("hMean%s_%s_%s", dcaName[ihisto].Data(), physSel.Data(), trigClassName.Data());
363         Int_t minBin = ( ibin == 0 ) ? 1 : ibin;
364         Int_t maxBin = ( ibin == 0 ) ? nPbins : ibin;
365         if ( ibin > 0 ) currName += Form("_pBin%i", ibin);
366         TH1* projHisto = histo->ProjectionY(currName.Data(), minBin, maxBin, "e");
367         projHisto->SetTitle(Form("%s %s %g < p < %g (GeV/c)", dcaName[ihisto].Data(), fThetaAbsKeys->At(itheta)->GetName(), meanDcaVsP->GetXaxis()->GetBinLowEdge(minBin), meanDcaVsP->GetXaxis()->GetBinUpEdge(maxBin)));
368
369         if ( projHisto->GetEntries() == 0 ) continue;
370         if ( ipad % (nPadX*nPadY) == 0 ) {
371           currName = histo->GetName();
372           currName += Form("Fit_can_%i", ican++);
373           meanDcaFitCan = new TCanvas(currName.Data(), currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
374           meanDcaFitCan->Divide(nPadX,nPadY);
375           ipad = 0;
376         }
377         meanDcaFitCan->cd(++ipad);
378         gPad->SetLogy();
379         if ( projHisto->Integral() > 50 ) {
380           fitFuncMeanDca->SetParameter(0, projHisto->Integral());
381           fitFuncMeanDca->SetParameter(1, projHisto->GetMean());
382           fitFuncMeanDca->SetParameter(2, projHisto->GetRMS());
383           Double_t fitDcaLim = ( ibin == 0 ) ?  40. :  40./((Double_t)ibin);
384           fitDcaLim = TMath::Max(5., fitDcaLim);
385           projHisto->Fit(fitFuncMeanDca, "RQ", "e", -fitDcaLim, fitDcaLim);
386           Double_t chi2 = fitFuncMeanDca->GetChisquare();
387           Double_t ndf = fitFuncMeanDca->GetNDF();
388           if ( ndf <= 0.) continue;
389           if ( chi2 / ndf > 100. ) continue;
390           Double_t meanDca = fitFuncMeanDca->GetParameter(1);
391           Double_t meanDcaErr = fitFuncMeanDca->GetParError(1);
392           meanDcaVsP->SetBinContent(ibin, meanDca);
393           meanDcaVsP->SetBinError(ibin, meanDcaErr);
394           //if ( ibin == 0 ) printf("%s  meanDca = %g +- %g\n", fThetaAbsKeys->At(itheta)->GetName(),  meanDca, meanDcaErr);
395         }
396         else projHisto->Draw("e");
397       } // loop on momentum bins
398       can->cd(2*itheta+ihisto+1);
399       //meanDcaVsP->SetLineColor(srcColors[isrc]);
400       //meanDcaVsP->SetMarkerColor(srcColors[isrc]);
401       //meanDcaVsP->SetMarkerStyle(20+isrc);
402       meanDcaVsP->Fit("pol0","Q");
403       TF1* trendFit = (TF1*)meanDcaVsP->GetListOfFunctions()->FindObject("pol0");
404       if ( trendFit ) {
405         averageDca[2*itheta+ihisto] = trendFit->GetParameter(0);
406         printf("  %s: mean %s = %g +- %g  (chi2/%i = %g)\n",  fThetaAbsKeys->At(itheta)->GetName(), dcaName[ihisto].Data(), trendFit->GetParameter(0), trendFit->GetParError(0), trendFit->GetNDF(), trendFit->GetChisquare()/((Double_t)trendFit->GetNDF()));
407       }
408         //drawOpt = "esame";
409       //leg->AddEntry(meanDcaVsP, fSrcKeys->At(isrc)->GetName(), "lp");
410       //} // loop on src
411     //can->cd(itheta+1);
412     //leg->Draw("same");
413
414       //can->cd(ipad++);
415       //histo->Draw();
416       //meanDca[ihisto] = histo->GetMean();
417     } // loop on histo type
418   } // loop on theta abs
419
420   for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {  
421     printf("muonCuts->SetMeanDCA(%g, %g, 0.); // %s\n", averageDca[2*itheta], averageDca[2*itheta+1], fThetaAbsKeys->At(itheta)->GetName());
422   }
423
424   //////////////
425   // Fit pDCA //
426   //////////////
427   Double_t nSigmaCut = fMuonTrackCuts->GetMuonTrackCutsParam().GetNSigmaPdca(); //6.;
428   Double_t sigmaMeasCut[2] = { fMuonTrackCuts->GetMuonTrackCutsParam().GetSigmaPdca23(), fMuonTrackCuts->GetMuonTrackCutsParam().GetSigmaPdca310()}; //{99., 54.}; //{120., 63.};
429   Double_t relPResolution = fMuonTrackCuts->GetMuonTrackCutsParam().GetRelPResolution(); //4.5e-4; //6.e-3;//8.e-4;
430   Double_t angleResolution = 535.*fMuonTrackCuts->GetMuonTrackCutsParam().GetSlopeResolution(); //6.e-4;
431   Double_t pMinCut = 0.1;
432   Double_t pMaxCut =  800.;
433   const Int_t kNCutFuncs = 2;
434   Int_t nShowFuncs = 1;
435   TString cutFormula = "[1]*TMath::Sqrt( ( [0] / ( 1. - [1]*[2]*x / ( 1.+[1]*[2]*x ) ) ) * ( [0] / ( 1. - [1]*[2]*x / ( 1.+[1]*[2]*x ) ) ) + [3]*[3]*x*x)";
436   Double_t cutParam[kNCutFuncs][4] = {{sigmaMeasCut[0], nSigmaCut, relPResolution, angleResolution}, {sigmaMeasCut[0], nSigmaCut, 0., 0.32}};
437   Int_t cutColor[kNCutFuncs] = {kBlack, kRed};
438   igroup1++;
439   igroup2 = 0;
440   can = new TCanvas("pdcaSigmaFit","Sigma vs P fit",igroup1*xshift,igroup2*yshift,600,600);
441   can->Divide(2,1);
442   igroup2++;
443   TF1* fitFunc = new TF1("fitFunc", "x*gausn", 0., 400.);
444   fitFunc->SetParNames("Norm", "Mean", "Sigma");
445   gStyle->SetOptFit(1111);
446   Double_t xMinFit[2] = {0., 0.};
447   Double_t xMaxFit[2] = {260., 150.}; //{320., 150.}; // {360., 180.};
448   printf("\nSigma p x DCA:\n");
449   Double_t averageSigmaPdca[kNtrackSources*kNthetaAbs] = {0.};
450   for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
451     TLegend* leg = new TLegend(0.7, 0.7, 0.9, 0.9);
452     leg->SetBorderSize(1);
453     for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
454       histoPattern = GetHistoName(kPdcaVsP, itheta, isrc);
455       TH2* histo = (TH2*)GetSum(physSel, trigClassName, centralityRange, histoPattern);
456       if ( ! histo ) continue;
457       if ( histo->Integral() < 200 ) {
458         delete histo;
459         continue;
460       }
461
462       TH1* sigmaVsP = histo->ProjectionX(Form("sigma%s_%s_%s", fHistoTypeKeys->At(kPdcaVsP)->GetName(), fThetaAbsKeys->At(itheta)->GetName(), fSrcKeys->At(isrc)->GetName()));
463       sigmaVsP->Reset();
464       sigmaVsP->SetTitle(Form("#sigma_{p#timesDCA} vs. p %s %s", fThetaAbsKeys->At(itheta)->GetName(), fSrcKeys->At(isrc)->GetName()));
465       sigmaVsP->SetYTitle("#sigma_{p#timesDCA} (cm #times GeV/c)");
466       sigmaVsP->SetStats(kFALSE);
467
468       Int_t nPbins = histo->GetXaxis()->GetNbins();
469       //Int_t nPadX = (Int_t)TMath::Sqrt(nPbins);
470       //Int_t nPadY = nPadX;
471       //if ( nPadX * nPadY < nPbins ) nPadX++;
472       TCanvas* pdcaFitCan = 0x0;
473       Int_t nPadX = 5;
474       Int_t nPadY = 5;
475       Int_t ipad = 0;
476       Int_t ican = 0;
477
478       for ( Int_t ibin=2; ibin<=nPbins; ++ibin ) {
479         currName = Form("hPDCA_%s_%s_%s_%s", fSrcKeys->At(isrc)->GetName(), physSel.Data(), trigClassName.Data(), fThetaAbsKeys->At(itheta)->GetName());
480         Int_t minBin = ( ibin == 0 ) ? 1 : ibin;
481         Int_t maxBin = ( ibin == 0 ) ? nPbins : ibin;
482         if ( ibin > 0 ) currName += Form("_pBin%i", ibin);
483         TH1* projHisto = histo->ProjectionY(currName.Data(), minBin, maxBin, "e");
484         if ( projHisto->GetEntries() == 0 ) continue;
485         projHisto->SetTitle(Form("P DCA %s %s %g < p < %g (GeV/c)", fThetaAbsKeys->At(itheta)->GetName(), fSrcKeys->At(isrc)->GetName(), sigmaVsP->GetXaxis()->GetBinLowEdge(minBin), sigmaVsP->GetXaxis()->GetBinUpEdge(maxBin)));
486         if ( ipad % (nPadX*nPadY) == 0 ) {
487           currName = histo->GetName();
488           currName += Form("Fit_can_%i", ican++);
489           pdcaFitCan = new TCanvas(currName.Data(), currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
490           pdcaFitCan->Divide(nPadX,nPadY);
491           ipad = 0;
492         }
493         pdcaFitCan->cd(++ipad);
494         if ( projHisto->Integral() > 0.) projHisto->Scale(1./projHisto->Integral());
495         gPad->SetLogy();
496         if ( projHisto->GetEntries() > 50 ) {
497           fitFunc->SetParameter(0, projHisto->Integral());
498           fitFunc->SetParameter(1, projHisto->GetMean());
499           fitFunc->SetParameter(2, projHisto->GetRMS());
500           projHisto->Fit(fitFunc, "RQ", "e", xMinFit[itheta], xMaxFit[itheta]);
501           Double_t chi2 = fitFunc->GetChisquare();
502           Double_t ndf = fitFunc->GetNDF();
503           if ( ndf <= 0.) continue;
504           if ( chi2 / ndf > 100. ) continue;
505           Double_t sigma = TMath::Abs(fitFunc->GetParameter(2));
506           Double_t sigmaErr = fitFunc->GetParError(2);
507           sigmaVsP->SetBinContent(ibin, sigma);
508           sigmaVsP->SetBinError(ibin, sigmaErr);
509           //if ( ibin == 0 ) printf("%s %s  sigma = %g +- %g\n", fThetaAbsKeys->At(itheta)->GetName(), fSrcKeys->At(isrc)->GetName(), sigma, sigmaErr);
510         }
511         else projHisto->Draw("e");
512       } // loop on momentum bins
513       can->cd(itheta+1);
514       sigmaVsP->SetLineColor(srcColors[isrc]);
515       sigmaVsP->SetMarkerColor(srcColors[isrc]);
516       sigmaVsP->SetMarkerStyle(20+isrc);
517       drawOpt = "e";
518       if ( gPad->GetListOfPrimitives()->GetEntries() != 0 ) drawOpt += "same";
519       sigmaVsP->Draw(drawOpt.Data());
520       sigmaVsP->Fit("pol0","Q",drawOpt.Data());
521       TF1* trendFit = (TF1*)sigmaVsP->GetListOfFunctions()->FindObject("pol0");
522       if ( trendFit ) {
523         averageSigmaPdca[kNtrackSources*itheta+isrc] = trendFit->GetParameter(0);
524         printf("  %s %s: Sigma = %g +- %g  (chi2/%i = %g)\n",  fThetaAbsKeys->At(itheta)->GetName(), fSrcKeys->At(isrc)->GetName(), trendFit->GetParameter(0), trendFit->GetParError(0), trendFit->GetNDF(), trendFit->GetChisquare()/((Double_t)trendFit->GetNDF()));
525       }
526       leg->AddEntry(sigmaVsP, fSrcKeys->At(isrc)->GetName(), "lp");
527     } // loop on src
528     can->cd(itheta+1);
529     leg->Draw("same");
530
531     for ( Int_t icut=0; icut<nShowFuncs; ++icut ) {
532       currName = Form("sigmaCut_%s_%i", fThetaAbsKeys->At(itheta)->GetName(), icut);
533       TF1* cutFunction = new TF1(currName.Data(), cutFormula.Data(), pMinCut, pMaxCut);
534       cutParam[icut][0] = sigmaMeasCut[itheta];
535       cutParam[icut][1] = 1.;
536       cutFunction->SetParameters(cutParam[icut]);
537       cutFunction->SetLineColor(cutColor[icut]);
538       cutFunction->SetLineWidth(2);
539       cutFunction->Draw("same");
540     }
541   } // loop on theta abs
542
543   for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
544     printf("muonCuts->SetSigmaPdca(%g, %g); // %s\n", averageSigmaPdca[isrc], averageSigmaPdca[kNtrackSources+isrc], fSrcKeys->At(isrc)->GetName());
545   }
546   printf("\n");
547
548   igroup2++;
549   for ( Int_t icheck=0; icheck<2; icheck++ ) {
550     Int_t hDraw = ( icheck == 0 ) ? kPDCAVsPCheck : kDCAVsPCheck;
551     for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
552       for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
553         histoPattern = GetHistoName(hDraw, itheta, isrc);
554         TH2* histoCheck = (TH2*)GetSum(physSel, trigClassName, centralityRange, histoPattern);
555         if ( ! histoCheck ) continue;
556         currName = histoCheck->GetName();
557         currName.Append("_plotCut");
558         TCanvas* pdcaVsPcan = new TCanvas(currName.Data(), currName.Data(), igroup1*xshift,(igroup2+icheck)*yshift,600,600);
559         pdcaVsPcan->SetRightMargin(0.12);
560         pdcaVsPcan->SetLogy();
561         pdcaVsPcan->SetLogz();
562         histoCheck->Draw("COLZ");
563
564         for ( Int_t icut=0; icut<nShowFuncs; ++icut ) {
565           currName = histoCheck->GetName();
566           currName.Append(Form("_cutFunc%i",icut));
567           TString currFormula = cutFormula;
568           if ( icheck == 1 ) currFormula.Append("/x");
569           TF1* cutFunction = new TF1(currName.Data(),currFormula.Data(), pMinCut, pMaxCut);
570           cutParam[icut][0] = sigmaMeasCut[itheta];
571           cutParam[icut][1] = nSigmaCut;
572           cutFunction->SetParameters(cutParam[icut]);
573           cutFunction->SetLineWidth(2);
574           cutFunction->SetLineColor(cutColor[icut]);
575           cutFunction->Draw("same");
576         } // loop on cut functions
577       } // loop on src
578     } //loop on theta
579   } // loop on check
580
581
582   ///////////////////////////
583   // Check Chi square prob //
584   ///////////////////////////
585   igroup1++;
586   igroup2 = 0;
587   for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
588     for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
589       histoPattern = GetHistoName(kChiProbVsP, itheta, isrc);
590       TH2* histo = (TH2*)GetSum(physSel, trigClassName, centralityRange, histoPattern);
591       if ( ! histo ) continue;
592
593       Int_t nPbins = histo->GetXaxis()->GetNbins();
594       //Int_t nPadX = (Int_t)TMath::Sqrt(nPbins);
595       //Int_t nPadY = nPadX;
596       //if ( nPadX * nPadY < nPbins ) nPadX++;
597       TCanvas* chi2probCan = 0x0;
598       Int_t nPadX = 5;
599       Int_t nPadY = 5;
600       Int_t ipad = 0;
601       Int_t ican = 0;
602
603       for ( Int_t ibin=0; ibin<=nPbins; ++ibin ) {
604         currName = Form("hChiProb_%s_%s_%s_%s", fSrcKeys->At(isrc)->GetName(), physSel.Data(), trigClassName.Data(), fThetaAbsKeys->At(itheta)->GetName());
605         Int_t minBin = ( ibin == 0 ) ? 1 : ibin;
606         Int_t maxBin = ( ibin == 0 ) ? nPbins : ibin;
607         if ( ibin > 0 ) currName += Form("_pBin%i", ibin);
608         TH1* projHisto = histo->ProjectionY(currName.Data(), minBin, maxBin);
609         projHisto->SetTitle(Form("Chisquare prob %s %s %g < p < %g (GeV/c)", fThetaAbsKeys->At(itheta)->GetName(), fSrcKeys->At(isrc)->GetName(), histo->GetXaxis()->GetBinLowEdge(minBin), histo->GetXaxis()->GetBinUpEdge(maxBin)));
610         if ( projHisto->GetEntries() == 0 ) continue;
611         if ( ipad % (nPadX*nPadY) == 0 ) {
612           currName = histo->GetName();
613           currName += Form("Fit_can_%i", ican++);
614           chi2probCan = new TCanvas(currName.Data(), currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
615           chi2probCan->Divide(nPadX,nPadY);
616           ipad = 0;
617         }
618         chi2probCan->cd(++ipad);
619         gPad->SetLogy();
620         projHisto->Draw("e");
621       } // loop on momentum bins
622     } // loop on src
623   } // loop on theta abs
624
625
626   //////////////////////
627   // Check sigma cuts //
628   //////////////////////
629   printf("\nReference sigma cut %g\n", refSigmaCut);
630
631   Float_t fracOfHeight = 0.35;
632   Float_t rightMargin = 0.03;
633   Int_t cutColors[14] = {kBlack, kRed, kBlue, kGreen, kCyan, kMagenta, kOrange, kViolet, kSpring, kGray, kSpring, kAzure, kPink, kYellow};
634   TArrayI orderCuts;
635   Int_t checkHistos[2] = {kSigmaVsPt, kSigmaVsEta};
636   Bool_t useCustomSigma = furtherOpt.Contains("CUSTOMSIGMA");
637
638   igroup1++;
639   igroup2 = 0;
640
641   for ( Int_t ihisto=0; ihisto<2; ++ihisto ) {
642     for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
643       can = 0x0;
644       for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
645         histoPattern = GetHistoName(checkHistos[ihisto], itheta, isrc);
646         TH2* histo = (TH2*)GetSum(physSel, trigClassName, centralityRange, histoPattern);
647         if ( ! histo ) continue;
648         if ( histo->Integral() == 0. ) {
649           delete histo;
650           continue;
651         }
652         if ( orderCuts.GetSize() == 0 ) {
653           // Re-order axis
654           TAxis* axis = histo->GetYaxis();
655           Int_t nSigmaCuts = ( useCustomSigma ) ? fSigmaCuts.GetSize() : axis->GetNbins();
656           orderCuts.Set(nSigmaCuts);
657           Int_t countBin = 0;
658           for ( Int_t isigma=0; isigma<axis->GetNbins(); ++isigma ) {
659             TString currLabel = axis->GetBinLabel(isigma+1);
660             Double_t currSigma = currLabel.Atof();
661             if ( useCustomSigma ) {
662               Bool_t foundMatch = kFALSE;
663               for ( Int_t jsigma=0; jsigma<fSigmaCuts.GetSize(); ++jsigma ) {
664                 if ( TMath::Abs(currSigma - fSigmaCuts[jsigma]) < 1e-4 ) {
665                   foundMatch = kTRUE;
666                   break;
667                 }
668               }
669               if ( ! foundMatch ) continue; 
670             }
671             Int_t currBin = ( TMath::Abs(currSigma - refSigmaCut) < 1e-4) ? 0 : ++countBin;
672             orderCuts[currBin] = isigma+1;
673           }
674         }
675         currName = histo->GetName();
676         currName.Append("_can");
677         can = new TCanvas(currName.Data(), currName.Data(), igroup1*xshift,igroup2*yshift, 600, 600);
678         TLegend* leg = new TLegend(0.6,0.6,0.9,0.9);
679         leg->SetBorderSize(1);
680         can->Divide(1,2,0,0);
681         can->cd(1);
682         gPad->SetPad(0., fracOfHeight, 0.99, 0.99);
683         //gPad->SetTopMargin(0.08);
684         gPad->SetTopMargin(0.03);
685         gPad->SetRightMargin(rightMargin);
686         gPad->SetLogy();
687
688         can->cd(2);
689         gPad->SetPad(0., 0., 0.99, fracOfHeight);
690         gPad->SetRightMargin(rightMargin);
691         gPad->SetBottomMargin(0.08/fracOfHeight);
692
693         TH1* refCutHisto = 0x0;
694         TString legTitle = "";
695         for ( Int_t isigma=0; isigma<orderCuts.GetSize(); ++isigma ) {
696           currName = histo->GetName();
697           currName.Append(Form("_sigma%i", isigma));
698           Int_t currBin = orderCuts[isigma];
699           TH1* projectHisto = histo->ProjectionX(currName.Data(), currBin, currBin);
700           projectHisto->SetLineColor(cutColors[isigma]);
701           projectHisto->SetMarkerColor(cutColors[isigma]);
702           projectHisto->SetMarkerStyle(20+isigma);
703           projectHisto->SetStats(kFALSE);
704           projectHisto->SetTitle("");
705           can->cd(1);
706           drawOpt = "e";
707           if ( gPad->GetListOfPrimitives()->GetEntries() != 0 ) drawOpt += "same";
708           projectHisto->Draw(drawOpt.Data());
709           TString currLabel = histo->GetYaxis()->GetBinLabel(currBin);
710           TString cutName = ( TMath::Abs(currLabel.Atof() ) > 1000. ) ? "Total" :  Form("%s #sigma cut", currLabel.Data());
711           leg->AddEntry(projectHisto, cutName.Data(), "lp");
712           if ( ! refCutHisto ) {
713             refCutHisto = projectHisto;
714             legTitle = Form("(pass n #sigma) / (pass %s #sigma)", currLabel.Data());
715             continue;
716           }
717           currName.Append("_ratio");
718           TH1* histoRatio = (TH1*)projectHisto->Clone(currName.Data());
719           histoRatio->Sumw2();
720           TH1* auxNum = projectHisto;
721           TH1* auxDen = refCutHisto;
722           Bool_t mustInvert = kFALSE;
723           if ( auxNum->Integral()>auxDen->Integral() ) {
724             auxNum = refCutHisto;
725             auxDen = projectHisto;
726             mustInvert = kTRUE;
727           }
728           histoRatio->Divide(auxNum,auxDen,1.,1.,"B");
729           if ( mustInvert ) {
730             TH1* auxHistoOne = static_cast<TH1*>(histoRatio->Clone("auxHistoOne"));
731             for ( Int_t ibin=1; ibin<=auxHistoOne->GetXaxis()->GetNbins(); ibin++ ) {
732               auxHistoOne->SetBinContent(ibin,1.);
733               auxHistoOne->SetBinError(ibin,0.);
734             }
735             histoRatio->Divide(auxHistoOne,histoRatio);
736             delete auxHistoOne;
737           }
738           histoRatio->SetTitle("");
739           histoRatio->GetYaxis()->SetTitle(legTitle.Data());
740           histoRatio->GetXaxis()->SetLabelSize(0.04/fracOfHeight);
741           histoRatio->GetXaxis()->SetTitleSize(0.035/fracOfHeight);
742           histoRatio->GetYaxis()->SetLabelSize(0.03/fracOfHeight);
743           histoRatio->GetYaxis()->SetTitleSize(0.03/fracOfHeight);
744           histoRatio->GetXaxis()->SetTitleOffset(1.);
745           histoRatio->GetYaxis()->SetTitleOffset(0.6);
746           histoRatio->GetYaxis()->SetRangeUser(0.5,1.5);
747
748           can->cd(2);
749           drawOpt = "e";
750           if ( gPad->GetListOfPrimitives()->GetEntries() != 0 ) drawOpt += "same";
751           histoRatio->Draw(drawOpt.Data());
752         }// loop on sigma cuts
753         can->cd(1);
754         leg->Draw("same");
755       } // loop on theta abs
756     } // loop on src
757   } // loop on histo type
758
759   igroup1++;
760   igroup2=0;
761   Double_t ptMin[] = {0., 2., 4., 15., 60.};
762   Int_t nPtMins = sizeof(ptMin)/sizeof(ptMin[0]);
763   for ( Int_t iptmin=0; iptmin<nPtMins; ++iptmin) {
764     for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
765       can = 0x0;
766       for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
767         TLegend* leg = 0x0;
768         histoName = GetHistoName(kSigmaVsPt, itheta, isrc);
769         for ( Int_t icent=1; icent<=GetCentralityClasses()->GetNbins(); ++icent ) {
770           TH2* histo = (TH2*)GetSum(physSel, trigClassName, GetCentralityClasses()->GetBinLabel(icent), histoName);
771           if ( ! histo ) continue;
772           Int_t ptMinBin = histo->GetXaxis()->FindBin(ptMin[iptmin]);
773           Int_t ptMaxBin = histo->GetXaxis()->GetNbins()+1;
774           currName = histo->GetName();
775           currName += Form("_%s_ptMin%i", GetCentralityClasses()->GetBinLabel(icent), TMath::Nint(ptMin[iptmin]));
776           TH1* projectHisto = histo->ProjectionY(currName.Data(), ptMinBin, ptMaxBin);
777           if ( projectHisto->GetMaximum() < 50. ) {
778             delete projectHisto;
779             continue;
780           }
781           if ( ! can ) {
782             currName = Form("CutEffect_%s_ptMin%i", fSrcKeys->At(isrc)->GetName(), TMath::Nint(ptMin[iptmin]));
783             can = new TCanvas(currName.Data(), currName.Data(), igroup1*xshift,igroup2*yshift, 600, 600);
784             can->Divide(2,1);
785           }
786           if ( ! leg ) {
787             leg = new TLegend(0.6,0.6,0.9,0.9);
788             leg->SetBorderSize(1);
789           }
790           can->cd(itheta+1);
791           projectHisto->SetTitle(Form("Cut effect %s %s %g < p_{t} < %g (GeV/c)", fSrcKeys->At(isrc)->GetName(), fThetaAbsKeys->At(itheta)->GetName(), histo->GetXaxis()->GetBinLowEdge(ptMinBin), histo->GetXaxis()->GetBinUpEdge(ptMaxBin)));
792           projectHisto->SetLineColor(cutColors[icent-1]);
793           projectHisto->SetMarkerColor(cutColors[icent-1]);
794           projectHisto->SetMarkerStyle(19+icent);
795           projectHisto->SetStats(0);
796           drawOpt = "e";
797           if ( gPad->GetListOfPrimitives()->GetEntries() != 0 ) drawOpt += "same";
798           projectHisto->Draw(drawOpt.Data());
799           leg->AddEntry(projectHisto, GetCentralityClasses()->GetBinLabel(icent), "lp");
800           Double_t keptEvents = projectHisto->GetBinContent(orderCuts[0]);
801           Double_t totalEvents = projectHisto->GetBinContent(orderCuts[orderCuts.GetSize()-1]);
802           Double_t accepted = ( totalEvents > 0. ) ? keptEvents / totalEvents : 1.;
803           Double_t acceptedErr = ( totalEvents > 0. ) ? TMath::Sqrt(accepted*(1.-accepted)/totalEvents) : 1.;
804           printf("%12s %11s %6s (pt>%g) rejected evts : %6.2f +- %6.3f %%\n", fSrcKeys->At(isrc)->GetName(), fThetaAbsKeys->At(itheta)->GetName(), GetCentralityClasses()->GetBinLabel(icent), ptMin[iptmin], (1.-accepted)*100., acceptedErr*100.);
805           //printf("  rejected %g  total %g   (%s vs %s)\n",totalEvents-keptEvents,totalEvents,projectHisto->GetXaxis()->GetBinLabel(orderCuts[0]),projectHisto->GetXaxis()->GetBinLabel(orderCuts[orderCuts.GetSize()-1]));
806         } // loop on centrality
807         if ( leg ) leg->Draw("same");
808       } // loop on theta abs
809     } // loop on sources
810   } // loop on pt min
811 }