Updates by Ionut
[u/mrichter/AliRoot.git] / PWGPP / TRD / AliTRDcheckESD.cxx
1 /**************************************************************************
2 * Copyright(c) 1998-1999, 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 /////////////////////////////////////////////////////
17 //
18 // Check basic detector results at ESD level
19 //   - Geometrical efficiency  
20 //   - Tracking efficiency  
21 //   - PID efficiency  
22 //   - Refit efficiency  
23 //
24 // Author
25 //   Alex Bercuci <A.Bercuci@gsi.de>
26 //   Ionut Arsene <i.c.arsene@gsi.de>
27 //
28 //////////////////////////////////////////////////////
29
30 #include <TClonesArray.h>
31 #include <TCanvas.h>
32 #include <TObjArray.h>
33 #include <TPad.h>
34 #include <TLegend.h>
35 #include <TLatex.h>
36 #include <TLine.h>
37 #include <TF1.h>
38 #include <TH1D.h>
39 #include <TH2D.h>
40 #include <TH3D.h>
41 #include <TH2I.h>
42 #include <TH2F.h>
43 #include <TH3S.h>
44 #include <TH3F.h>
45 #include <TProfile2D.h>
46 #include <TProfile.h>
47 #include <TGraphErrors.h>
48 #include <TGraphAsymmErrors.h>
49 #include <TFile.h>
50 #include <TTree.h>
51 #include <TROOT.h>
52 #include <TChain.h>
53 #include <TParticle.h>
54 #include <TTimeStamp.h>
55 #include <TRandom.h>
56 #include <TString.h>
57
58 #include "AliLog.h"
59 #include "AliAnalysisManager.h"
60 #include "AliAnalysisCuts.h"
61 #include "AliPhysicsSelection.h"
62 #include "AliESDEvent.h"
63 #include "AliESDkink.h"
64 #include "AliMCEvent.h"
65 #include "AliESDInputHandler.h"
66 #include "AliMCEventHandler.h"
67 #include "AliESDpid.h"
68
69 #include "AliESDtrack.h"
70 #include "AliMCParticle.h"
71 #include "AliPID.h"
72 #include "AliStack.h"
73 #include "AliTrackReference.h"
74 //#include "AliESDCentrality.h"
75 #include "AliMultiplicity.h"
76 #include "AliCFContainer.h"
77
78 #include "AliTRDcheckESD.h"
79 #include <iostream>
80 using std::cout;
81 using std::endl;
82
83 ClassImp(AliTRDcheckESD)
84
85 const Float_t AliTRDcheckESD::fgkxTPC = 290.;
86 const Float_t AliTRDcheckESD::fgkxTOF = 365.;
87 const UChar_t AliTRDcheckESD::fgkNgraph[AliTRDcheckESD::kNrefs] ={
88 8, 4, 2, 20};
89 FILE* AliTRDcheckESD::fgFile = NULL;
90
91 const Float_t AliTRDcheckESD::fgkEvVertexZ = 15.;
92 const Int_t   AliTRDcheckESD::fgkEvVertexN = 1;
93 const Float_t AliTRDcheckESD::fgkTrkDCAxy  = 40.;
94 const Float_t AliTRDcheckESD::fgkTrkDCAz   = 15.;
95 const Int_t   AliTRDcheckESD::fgkNclTPC    = 100;
96 const Float_t AliTRDcheckESD::fgkPt        = 0.2;
97 const Float_t AliTRDcheckESD::fgkEta       = 0.9;
98 const Float_t AliTRDcheckESD::fgkQs        = 0.002;
99
100 //____________________________________________________________________
101 AliTRDcheckESD::AliTRDcheckESD():
102   AliAnalysisTaskSE()
103   ,fStatus(0)
104   ,fNRefFigures(0)
105   ,fESD(NULL)
106   ,fMC(NULL)
107   ,fESDpid(new AliESDpid)
108   ,fHistos(NULL)
109   ,fResults(NULL)
110   ,fExpertCF(NULL)
111   ,fMatchingPhiEtaCF(NULL)
112   ,fMatchingPtCF(NULL)
113   ,fBunchCrossingsCF(NULL)
114   ,fCentralityCF(NULL)
115   ,fQtotCF(NULL)
116   ,fPulseHeightCF(NULL)
117   ,fReferenceTrackFilter(NULL)
118   ,fPhysSelTriggersEnabled(kFALSE)
119   ,fUserEnabledTriggers("")
120   ,fNAssignedTriggers(0)
121 {
122   //
123   // Default constructor
124   //
125   for(Int_t i=0; i<kNTrdCfVariables; ++i) {
126     fExpertCFVars[i] = -1;
127     fMatchingPhiEtaCFVars[i] = -1;
128     fMatchingPtCFVars[i] = -1;
129     fBunchCrossingsCFVars[i] = -1;
130     fCentralityCFVars[i] = -1;
131     fQtotCFVars[i] = -1;
132     fPulseHeightCFVars[i] = -1;
133     fExpertCFVarsEnabled[i] = kFALSE;
134     fExpertCFVarNBins[i] = 0;
135     fExpertCFVarRanges[i][0] = -999.; fExpertCFVarRanges[i][1] = -999.;
136     fExpertCFVarBins[i] = "";
137   }
138   fExpertCFEnabledSteps[0] = kFALSE; fExpertCFEnabledSteps[1] = kFALSE; fExpertCFEnabledSteps[2] = kFALSE;
139   SetNameTitle("TRDcheckESD", "Check TRD @ ESD level");
140   SetMC(kTRUE);
141 }
142
143 //____________________________________________________________________
144 AliTRDcheckESD::AliTRDcheckESD(char* name):
145   AliAnalysisTaskSE(name)
146   ,fStatus(0)
147   ,fNRefFigures(0)
148   ,fESD(NULL)
149   ,fMC(NULL)
150   ,fESDpid(new AliESDpid)
151   ,fHistos(NULL)
152   ,fResults(NULL)
153   ,fExpertCF(NULL)
154   ,fMatchingPhiEtaCF(NULL)
155   ,fMatchingPtCF(NULL)
156   ,fBunchCrossingsCF(NULL)
157   ,fCentralityCF(NULL)
158   ,fQtotCF(NULL)
159   ,fPulseHeightCF(NULL)
160   ,fReferenceTrackFilter(NULL)
161   ,fPhysSelTriggersEnabled(kFALSE)
162   ,fUserEnabledTriggers("")
163   ,fNAssignedTriggers(0)
164 {
165   //
166   // Default constructor
167   //
168   for(Int_t i=0; i<kNTrdCfVariables; ++i) {
169     fExpertCFVars[i] = -1;
170     fMatchingPhiEtaCFVars[i] = -1;
171     fMatchingPtCFVars[i] = -1;
172     fBunchCrossingsCFVars[i] = -1;
173     fCentralityCFVars[i] = -1;
174     fQtotCFVars[i] = -1;
175     fPulseHeightCFVars[i] = -1;
176     fExpertCFVarsEnabled[i] = kFALSE;
177     fExpertCFVarNBins[i] = 0;
178     fExpertCFVarRanges[i][0] = -999.; fExpertCFVarRanges[i][1] = -999.;
179     fExpertCFVarBins[i] = "";
180   }
181   fExpertCFEnabledSteps[0] = kFALSE; fExpertCFEnabledSteps[1] = kFALSE; fExpertCFEnabledSteps[2] = kFALSE;
182   SetMC(kTRUE);
183   SetTitle("Check TRD @ ESD level");
184   DefineOutput(1, TObjArray::Class());
185 }
186
187 //____________________________________________________________________
188 AliTRDcheckESD::~AliTRDcheckESD()
189 {
190   // Destructor
191   if(fHistos && !(AliAnalysisManager::GetAnalysisManager() && AliAnalysisManager::GetAnalysisManager()->IsProofMode())){
192     if(fHistos->IsOwner()) fHistos->Delete();
193     delete fHistos;
194     fHistos = NULL;
195   }
196   
197   if(fResults){
198     fResults->Delete();
199     delete fResults;
200   }
201 }
202
203
204 //____________________________________________________________________
205 void AliTRDcheckESD::AddExpertCFVar(AliTRDcheckESD::ETrdCfVariables var, 
206                                     Int_t nbins, Double_t lowLim, Double_t highLim) {
207   //
208   // Configure variables for the expert CF container
209   //
210   fExpertCFVarsEnabled[var] = kTRUE;
211   fExpertCFVarNBins[var] = nbins;
212   fExpertCFVarRanges[var][0] = lowLim;
213   fExpertCFVarRanges[var][1] = highLim;
214 }
215
216
217 //____________________________________________________________________
218 void AliTRDcheckESD::AddExpertCFVar(AliTRDcheckESD::ETrdCfVariables var, 
219                                     const Char_t* bins) {
220   //
221   // Configure variables for the expert CF container
222   //
223   fExpertCFVarsEnabled[var] = kTRUE;
224   fExpertCFVarBins[var] = bins;
225 }
226
227
228 //____________________________________________________________________
229 void AliTRDcheckESD::UserCreateOutputObjects()
230 {       
231   //
232   // Create Output Containers (TObjectArray containing 1D histograms)
233   //
234   Histos();
235   PostData(1, fHistos);
236 }
237
238 //____________________________________________________________________
239 void AliTRDcheckESD::MakeSummaryFromCF(Double_t* trendValues, const Char_t* triggerName, Bool_t useIsolatedBC, Bool_t cutTOFbc){
240   //
241   // Draw summary plots for the ESDcheck task using the CF container
242   //
243
244   cout << "Make summary from CF" << endl;
245   TCanvas *cOut=0x0;
246   if(gROOT->FindObject("trackingSummary")) delete gROOT->FindObject("trackingSummary");
247   cOut = new TCanvas("trackingSummary", "Tracking summary for the ESD task", 1600, 1200);
248   cOut->cd();
249   PlotTrackingSummaryFromCF(trendValues, triggerName, useIsolatedBC, cutTOFbc);
250   cOut->SaveAs("trackingSummary.gif");
251   
252   if(gROOT->FindObject("pidSummary")) delete gROOT->FindObject("pidSummary");
253   cOut = new TCanvas("pidSummary", "PID summary for the ESD task", 1600, 1200);
254   cOut->cd();
255   //GetRefFigure(6);
256   PlotPidSummaryFromCF(trendValues, triggerName, useIsolatedBC, cutTOFbc);
257   cOut->SaveAs("pidSummary.gif");
258
259   if(gROOT->FindObject("centSummary")) delete gROOT->FindObject("centSummary");
260   cOut = new TCanvas("centSummary", "Centrality summary for the ESD task", 1600, 1200);
261   cOut->cd();
262   //GetRefFigure(7);
263   PlotCentSummaryFromCF(trendValues, triggerName, useIsolatedBC, cutTOFbc);
264   cOut->SaveAs("centSummary.gif");
265
266 }
267
268
269 //____________________________________________________________________
270 void AliTRDcheckESD::MakeSummary(Double_t* trendValues){
271   //
272   // Draw summary plots for the ESDcheck task
273   //
274   // Old method to draw summary pictures from histograms. Use the MakeSummaryFromCF() when CF container is present
275   
276   cout << "Make summary" << endl;
277   TCanvas *cTracking=0x0;
278   if(gROOT->FindObject("trackingSummary")) delete gROOT->FindObject("trackingSummary");
279   cTracking = new TCanvas("trackingSummary", "Tracking summary for the ESD task", 1600, 1200);
280   cTracking->cd();
281   //GetRefFigure(5);
282   if(PlotTrackingSummary(0, trendValues))
283     cTracking->SaveAs("trackingSummary.gif");
284   
285   TCanvas* cPid = 0x0;
286   if(gROOT->FindObject("pidSummary")) delete gROOT->FindObject("pidSummary");
287   cPid = new TCanvas("pidSummary", "PID summary for the ESD task", 1600, 1200);
288   cPid->cd();
289   //GetRefFigure(6);
290   if(PlotPidSummary(0, trendValues))
291     cPid->SaveAs("pidSummary.gif");
292
293   TCanvas* cCent=0x0;
294   if(gROOT->FindObject("centSummary")) delete gROOT->FindObject("centSummary");
295   cCent = new TCanvas("centSummary", "Centrality summary for the ESD task", 1600, 1200);
296   cCent->cd();
297   //GetRefFigure(7);
298   if(PlotCentSummary(trendValues))
299     cCent->SaveAs("centSummary.gif");
300 }
301
302 //____________________________________________________________________
303 Bool_t AliTRDcheckESD::GetRefFigure(Int_t ifig)
304 {
305   //
306   // Produce reference Plots during PostProcessing
307   //
308   if(ifig>=fNRefFigures){
309     AliWarning(Form("Ref plot %d not available. Valid only up to %d", ifig, fNRefFigures));
310     return kFALSE;
311   }
312   if(!gPad){
313     AliWarning("Please provide a canvas to draw results.");
314     return kFALSE;
315   } else {
316     gPad->SetLogx(0);gPad->SetLogy(0);
317     gPad->SetMargin(0.125, 0.015, 0.1, 0.015);
318   }
319
320   const Char_t *title[20];
321   TH1 *hF(NULL);
322   if((hF=(TH1S*)gROOT->FindObject("hFcheckESD"))) delete hF;
323   TLegend *leg(NULL);
324   TList *l(NULL); TVirtualPad *pad(NULL);
325   TGraphErrors *g(NULL);TGraphAsymmErrors *ga(NULL);
326   TObjArray *arr(NULL);
327   TLatex *lat=new TLatex();
328   lat->SetTextSize(0.07);
329   lat->SetTextColor(2);
330   TLine line;
331   TTimeStamp now;
332   switch(ifig){
333   case kNCl: // number of clusters/track
334     if(!(arr = (TObjArray*)fResults->At(kNCl))) return kFALSE;
335
336     leg = new TLegend(.83, .7, .99, .96);
337     leg->SetHeader("Species");
338     leg->SetBorderSize(0); leg->SetFillStyle(0);
339     for(Int_t ig(0); ig<fgkNgraph[kNCl-1]; ig++){
340       if(!(g = (TGraphErrors*)arr->At(ig))) return kFALSE;
341       if(!g->GetN()) continue;
342       g->Draw(ig?"pc":"apc"); leg->AddEntry(g, g->GetTitle(), "pl");
343       if(ig) continue;
344       hF=g->GetHistogram();
345       hF->SetXTitle("no of clusters");
346       hF->SetYTitle("entries"); 
347       hF->GetYaxis()->CenterTitle(1);
348       hF->GetYaxis()->SetTitleOffset(1.2);
349       hF->SetMinimum(5);
350     }
351     leg->Draw(); gPad->SetLogy();
352     break;
353   case kTRDstat: // Efficiency
354     if(!(arr = (TObjArray*)fResults->At(kTRDstat))) return kFALSE;
355     leg = new TLegend(.62, .77, .98, .98);
356     leg->SetHeader("TRD Efficiency");
357     leg->SetBorderSize(0); leg->SetFillStyle(0);
358     title[0] = "Geometrical (TRDin/TPCout)";
359     title[1] = "Tracking (TRDout/TRDin)";
360     title[2] = "PID (TRDpid/TRDin)";
361     title[3] = "Refit (TRDrefit/TRDin)";
362     hF = new TH1S("hFcheckESD", ";p [GeV/c];Efficiency", 10, 0.1, 10.);
363     hF->SetMaximum(1.4);
364     hF->GetXaxis()->SetMoreLogLabels();
365     hF->GetYaxis()->CenterTitle(1);
366     hF->Draw("p");
367     for(Int_t ig(0); ig<fgkNgraph[kTRDstat-1]; ig++){
368       if(!(g = (TGraphErrors*)arr->At(ig))) return kFALSE;
369       g->Draw("pl"); leg->AddEntry(g, title[ig], "pl");
370       //PutTrendValue(name[id], g->GetMean(2));
371       //PutTrendValue(Form("%sRMS", name[id]), g->GetRMS(2));
372     }
373     leg->Draw(); gPad->SetLogx();
374     break;
375   case kTRDmom: // Energy loss
376     if(!(arr = (TObjArray*)fResults->At(kTRDmom))) return kFALSE;
377     leg = new TLegend(.65, .7, .95, .99);
378     leg->SetHeader("Energy Loss");
379     leg->SetBorderSize(1); leg->SetFillColor(0);
380     title[0] = "Max & 90% quantile";
381     title[1] = "Mean & 60% quantile";
382     hF = new TH1S("hFcheckESD", ";layer;#Delta E", 6, -0.5, 5.5);
383     hF->SetMaximum(1.3);hF->SetMinimum(-.3);
384     hF->Draw("p");
385     for(Int_t ig(0); ig<fgkNgraph[kTRDmom-1]; ig++){
386       if(!(ga = (TGraphAsymmErrors*)arr->At(ig))) return kFALSE;
387       ga->Draw("pl"); leg->AddEntry(ga, title[ig], "pl");
388       //PutTrendValue(name[id], g->GetMean(2));
389       //PutTrendValue(Form("%sRMS", name[id]), g->GetRMS(2));
390     }
391     leg->Draw();gPad->SetLogx(kFALSE);
392     break;
393   case kPtRes: // Pt resolution @ vertex
394     if(!(arr = (TObjArray*)fResults->At(kPtRes))) return kFALSE;
395     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
396     pad = ((TVirtualPad*)l->At(0)); pad->cd(); pad->SetLogx();
397     pad->SetMargin(0.1, 0.022, 0.1, 0.023);
398     hF = new TH1S("hFcheckESD", "ITS+TPC+TRD;p_{t} [GeV/c];#Delta p_{t} / p_{t} [%]", 10, 0.2, 10.);
399     hF->SetMaximum(10.);hF->SetMinimum(-3.);
400     hF->GetXaxis()->SetMoreLogLabels();
401     hF->GetXaxis()->SetTitleOffset(1.2);
402     hF->GetYaxis()->CenterTitle();
403     hF->Draw("p");
404     //for(Int_t ig(0); ig<fgkNgraph[kPtRes-1]/2; ig++){
405     for(Int_t ig(2); ig<6; ig++){
406       if(!(g = (TGraphErrors*)arr->At(ig))) continue;
407       if(!g->GetN()) continue;
408       g->Draw("pl");
409       //PutTrendValue(name[id], g->GetMean(2));
410       //PutTrendValue(Form("%sRMS", name[id]), g->GetRMS(2));
411     }
412     pad = ((TVirtualPad*)l->At(1)); pad->cd(); pad->SetLogx();
413     pad->SetMargin(0.1, 0.22, 0.1, 0.023);
414     hF = (TH1*)hF->Clone("hFcheckESD1");
415     hF->SetTitle("ITS+TPC");
416     hF->SetMaximum(10.);hF->SetMinimum(-3.);
417     hF->Draw("p");
418     leg = new TLegend(.78, .1, .99, .98);
419     leg->SetHeader("P_{t} @ DCA");
420     leg->SetBorderSize(1); leg->SetFillColor(0);
421     leg->SetTextAlign(22);
422     leg->SetTextFont(12);
423     leg->SetTextSize(0.03813559);
424     {
425       Int_t nPlots(0);
426       //for(Int_t ig(fgkNgraph[kPtRes-1]/2); ig<fgkNgraph[kPtRes-1]; ig++){
427       for(Int_t ig(12); ig<16; ig++){
428         if(!(g = (TGraphErrors*)arr->At(ig))) continue;
429         if(!g->GetN()) continue;
430         nPlots++;
431         g->Draw("pl"); leg->AddEntry(g, g->GetTitle(), "pl");
432         //PutTrendValue(name[id], g->GetMean(2));
433         //PutTrendValue(Form("%sRMS", name[id]), g->GetRMS(2));
434       }
435       if(nPlots) leg->Draw();
436     }
437     break;
438   case 5: // plot a 3x3 canvas with tracking related histograms
439     PlotTrackingSummary(0);
440     break;
441     
442   case 6: // plot a 3x3 canvas with PID related histograms
443     PlotPidSummary(0);
444     break;
445
446   case 7: // plot a 3x3 canvas with centrality dependence histograms
447     PlotCentSummary();
448     break;
449
450   }
451   return kTRUE;
452 }
453
454 //____________________________________________________________________
455 void AliTRDcheckESD::UserExec(Option_t *){
456   //
457   // Run the Analysis
458   //
459   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
460   fMC = MCEvent();
461
462   if(!fESD){
463     AliError("ESD event missing.");
464     return;
465   }
466   
467   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
468   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
469   if(!inputHandler) return;
470   
471   if(!fPhysSelTriggersEnabled) {
472     InitializeCFContainers();
473     fPhysSelTriggersEnabled = kTRUE;
474   }
475     
476   UInt_t isSelected = AliVEvent::kAny;
477   if(inputHandler){
478     if(inputHandler->GetEventSelection()) {
479       isSelected = inputHandler->IsEventSelected();
480     }
481   }
482   if(!isSelected) return;
483
484   TString triggerClasses = fESD->GetFiredTriggerClasses();
485   //  cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++triggers fired:  " << triggerClasses.Data() << endl;
486   TObjArray* triggers = triggerClasses.Tokenize(" ");
487   TObjArray* userTriggers = fUserEnabledTriggers.Tokenize(";");
488   if(triggers->GetEntries()<1) return;
489   Bool_t hasGoodTriggers = kFALSE;
490   Int_t triggerIndices[kNMaxAssignedTriggers] = {0};
491   Int_t nTrigFired=0;
492   for(Int_t i=0; i<triggers->GetEntries(); ++i) {
493     //    cout << "check trigger " << triggers->At(i)->GetName() << endl;
494     TString trigStr=triggers->At(i)->GetName();
495     if(!trigStr.Contains("NOTRD") && !trigStr.Contains("MUON")) hasGoodTriggers = kTRUE;    // check wheter TRD was read out in this event
496     if(i>=kNMaxAssignedTriggers) continue;
497     //    triggerIndices[i] = GetTriggerIndex(triggers->At(i)->GetName(), kFALSE);
498     for(Int_t j=0;j<userTriggers->GetEntries();++j) {
499       TString userTrigStr=userTriggers->At(j)->GetName();
500       if(trigStr.Contains(userTrigStr.Data())) {
501         triggerIndices[nTrigFired] = GetTriggerIndex(userTrigStr.Data(), kFALSE);
502         if(triggerIndices[nTrigFired]==-1) triggerIndices[nTrigFired]=1;  // 0-assigned to all other triggers
503         ++nTrigFired;
504       }
505     }
506     triggerIndices[nTrigFired] = GetTriggerIndex(trigStr.Data(), kFALSE);
507     if(triggerIndices[nTrigFired]==-1) triggerIndices[nTrigFired]=1;  // 0-assigned to all other triggers
508     ++nTrigFired;
509   } 
510   //  Int_t nTRDtracks = fESD->GetNumberOfTrdTracks();
511   //  Int_t nGlobalTracks = fESD->GetNumberOfTracks();
512   //cout << "TRD/All tracks: " << nTRDtracks << "/" << nGlobalTracks << endl;
513   for(Int_t i=0; i<nTrigFired; ++i) 
514     ((TH1F*)fHistos->At(kTriggerDefs))->Fill(triggerIndices[i]);
515
516   if(!hasGoodTriggers) {
517     PostData(1, fHistos);
518     return;
519   }
520   
521   // Get MC information if available
522   AliStack * fStack = NULL;
523   if(HasMC()){
524     if(!fMC){ 
525       AliWarning("MC event missing");
526       SetMC(kFALSE);
527     } else {
528       if(!(fStack = fMC->Stack())){
529         AliWarning("MC stack missing");
530         SetMC(kFALSE);
531       }
532     }
533   }
534   TH1 *h(NULL);
535   
536   Double_t values[kNTrdCfVariables];      // array where the CF container variables are stored
537   values[kEventVtxZ] = fESD->GetPrimaryVertex()->GetZv();
538   values[kEventBC] = fESD->GetBunchCrossNumber();
539   
540   const AliMultiplicity* mult=fESD->GetMultiplicity();
541   Double_t itsNTracklets = mult->GetNumberOfTracklets();
542   if(itsNTracklets<1) return;
543   Int_t multLimits[6] = {0, 700, 1400, 2100, 2800, 3500};
544   Int_t centralityClass = 0;
545   for(Int_t iCent=0; iCent<5; ++iCent) {
546     if(itsNTracklets>=multLimits[iCent] && itsNTracklets<multLimits[iCent+1])
547       centralityClass=iCent+1;
548   }
549   values[kEventMult] = itsNTracklets;
550   if(centralityClass == 0) return;
551   
552   Double_t* valuesMatchingPhiEtaCF = new Double_t[fMatchingPhiEtaCF->GetNVar()];
553   Double_t* valuesMatchingPtCF = new Double_t[fMatchingPtCF->GetNVar()];
554   Double_t* valuesBCCF = new Double_t[fBunchCrossingsCF->GetNVar()];
555   Double_t* valuesCentCF = new Double_t[fCentralityCF->GetNVar()];
556   Double_t* valuesQtotCF = new Double_t[fQtotCF->GetNVar()];
557   Double_t* valuesPHCF = new Double_t[fPulseHeightCF->GetNVar()];
558   Double_t* valuesExpertCF = (fExpertCF ? new Double_t[fExpertCF->GetNVar()] : 0x0);
559   
560   AliESDtrack *esdTrack(NULL);
561   for(Int_t itrk = 0; itrk < fESD->GetNumberOfTracks(); itrk++){
562     esdTrack = fESD->GetTrack(itrk);
563     //    cout << "track pt/eta: " << esdTrack->Pt() << "/" << esdTrack->Eta() << endl;
564     Float_t dcaxy,dcaz;
565     esdTrack->GetImpactParameters(dcaxy,dcaz);
566     //    cout << "dca xy/z: " << dcaxy << "/" << dcaz << endl;
567     //    cout << "TPC ncls: " << esdTrack->GetTPCNcls() << endl;
568     //    cout << "ITS hit map: ";
569     //    for (Int_t iC=0; iC<6; iC++) {
570     //      cout << ((esdTrack->GetITSClusterMap())&(1<<(iC)) ? "1" : "0") << flush;
571     //    }
572     //    cout << endl;
573     if(!fReferenceTrackFilter->IsSelected(esdTrack)) continue;
574     //    cout << "track passed" << endl;    
575
576     ULong_t status = esdTrack->GetStatus(); //PrintStatus(status);
577             
578     // pid quality
579     Bool_t kBarrel = Bool_t(status & AliESDtrack::kTRDin);
580
581     // find position and momentum of the track at entrance in TRD
582     Double_t rTRD[6] = {298.0, 311.0, 324.0, 337.0, 350.0, 363.0};
583     Double_t localCoord[6][3] = {{0.0}};
584     Bool_t localCoordGood[6];
585     for(Int_t il=0;il<6;++il) localCoordGood[il] = esdTrack->GetXYZAt(rTRD[il], fESD->GetMagneticField(), localCoord[il]);
586     Double_t localMom[6][3] = {{0.0}};
587     Bool_t localMomGood[6];
588     for(Int_t il=0; il<6; ++il) localMomGood[il] = esdTrack->GetPxPyPzAt(rTRD[il], fESD->GetMagneticField(), localMom[il]);
589     //Double_t localPhi = (localMomGood ? TMath::ATan2(localMom[1], localMom[0]) : 0.0);
590     Double_t localSagitaPhi[6] = {-999.};
591     for(Int_t il=0; il<6; ++il) localSagitaPhi[il] = (localCoordGood[il] ? TMath::ATan2(localCoord[il][1], localCoord[il][0]) : -999.);
592
593     values[kTrackTOFBC]      = esdTrack->GetTOFBunchCrossing(fESD->GetMagneticField());
594     Float_t dcaXY=0.0; Float_t dcaZ=0.0;
595     esdTrack->GetImpactParameters(dcaXY, dcaZ);
596     values[kTrackDCAxy]  = dcaXY;
597     values[kTrackDCAz]   = dcaZ;
598     values[kTrackCharge] = esdTrack->Charge();
599     values[kTrackPhi]    = localSagitaPhi[0];
600     values[kTrackEta]    = esdTrack->Eta();
601     values[kTrackPt]     = esdTrack->Pt();
602     values[kTrackP]      = esdTrack->P();
603     values[kTrackTrdTracklets] = esdTrack->GetTRDntracklets();
604     values[kTrackTrdClusters] = esdTrack->GetTRDncls();
605     for(Int_t i=0; i<6; ++i) values[kTrackQtot+i] = 0.0;
606         
607     if(localCoordGood[0] && localMomGood[0]) {
608       for(Int_t itrig=0; itrig<nTrigFired; ++itrig) {
609         values[kEventTrigger] = triggerIndices[itrig];
610         if((fMatchingPhiEtaCF->GetVar("trigger")<0 && itrig==0) || (fMatchingPhiEtaCF->GetVar("trigger")>=0)) {
611           for(Int_t iv=0; iv<fMatchingPhiEtaCF->GetNVar(); ++iv) valuesMatchingPhiEtaCF[iv] = values[fMatchingPhiEtaCFVars[iv]];
612           fMatchingPhiEtaCF->Fill(valuesMatchingPhiEtaCF, 0);
613         }
614         if((fMatchingPtCF->GetVar("trigger")<0 && itrig==0) || (fMatchingPtCF->GetVar("trigger")>=0)) {
615           for(Int_t iv=0; iv<fMatchingPtCF->GetNVar(); ++iv) valuesMatchingPtCF[iv] = values[fMatchingPtCFVars[iv]];
616           fMatchingPtCF->Fill(valuesMatchingPtCF, 0);
617         }
618         if(values[kTrackPt]>1.0 && values[kTrackPt]<3.0)
619           if((fBunchCrossingsCF->GetVar("trigger")<0 && itrig==0) || (fBunchCrossingsCF->GetVar("trigger")>=0)) {
620             for(Int_t iv=0; iv<fBunchCrossingsCF->GetNVar(); ++iv) valuesBCCF[iv] = values[fBunchCrossingsCFVars[iv]];
621             fBunchCrossingsCF->Fill(valuesBCCF, 0);
622           }
623         if(fExpertCF) {
624           if((fExpertCF->GetVar("trigger")<0 && itrig==0) || (fExpertCF->GetVar("trigger")>=0))
625             if(fExpertCF->GetStep("TPC")>=0 && fExpertCF->GetStep("TPC")<3) {
626               for(Int_t iv=0; iv<fExpertCF->GetNVar(); ++iv) valuesExpertCF[iv] = values[fExpertCFVars[iv]];
627               fExpertCF->Fill(valuesExpertCF, fExpertCF->GetStep("TPC"));
628             }
629         }
630       }
631     }
632             
633     // TRD reference tracks
634     if(values[kTrackTrdTracklets]>=1) {
635       // (slicePH,sliceNo) distribution and Qtot from slices
636       for(Int_t iPlane=0; iPlane<6; iPlane++) {
637         values[kTrackQtot+iPlane] = fgkQs*esdTrack->GetTRDslice(iPlane, 0);
638         values[kTrackPhi] = localSagitaPhi[iPlane];
639         for(Int_t itrig=0; itrig<nTrigFired; ++itrig) {
640           values[kEventTrigger] = triggerIndices[itrig];
641           if((fCentralityCF->GetVar("trigger")<0 && itrig==0) || (fCentralityCF->GetVar("trigger")>=0)) {
642             for(Int_t iv=0; iv<fCentralityCF->GetNVar(); ++iv) valuesCentCF[iv] = values[fCentralityCFVars[iv]];
643             valuesCentCF[fCentralityCF->GetNVar()-1] = values[kTrackQtot+iPlane];
644             fCentralityCF->Fill(valuesCentCF, 0);
645           }
646           if(values[kTrackTrdTracklets]>=4)
647             if((fQtotCF->GetVar("trigger")<0 && itrig==0) || (fQtotCF->GetVar("trigger")>=0)) {
648               for(Int_t iv=0; iv<fQtotCF->GetNVar()-2; ++iv) valuesQtotCF[iv] = values[fQtotCFVars[iv]];
649               valuesQtotCF[fQtotCF->GetNVar()-2] = values[kTrackQtot+iPlane];
650               valuesQtotCF[fQtotCF->GetNVar()-1] = iPlane;
651               fQtotCF->Fill(valuesQtotCF, 0);
652             }
653         }
654         for(Int_t iSlice=0; iSlice<8; iSlice++) {
655           if(esdTrack->GetTRDslice(iPlane, iSlice)>20.) {
656             values[kTrackPHslice+iSlice] =  fgkQs*esdTrack->GetTRDslice(iPlane, iSlice);
657             h = (TH2F*)fHistos->At(kPHSlice); h->Fill(iSlice, values[kTrackPHslice+iSlice]);
658             h = (TH2F*)fHistos->At(kPHSlice+centralityClass); h->Fill(iSlice, values[kTrackPHslice+iSlice]);
659             for(Int_t itrig=0; itrig<triggers->GetEntries(); ++itrig) {
660               values[kEventTrigger] = triggerIndices[itrig];
661               if((fPulseHeightCF->GetVar("trigger")<0 && itrig==0) || (fPulseHeightCF->GetVar("trigger")>=0)) {
662                 for(Int_t iv=0; iv<fPulseHeightCF->GetNVar()-2; ++iv) valuesPHCF[iv] = values[fPulseHeightCFVars[iv]];
663                 valuesPHCF[fPulseHeightCF->GetNVar()-2] = values[kTrackPHslice+iSlice];
664                 valuesPHCF[fPulseHeightCF->GetNVar()-1] = iSlice;
665                 fPulseHeightCF->Fill(valuesPHCF, 0);
666               }
667             }
668           }
669         }
670       }
671       values[kTrackPhi] = localSagitaPhi[0];
672             
673       if(localCoordGood[0] && localMomGood[0]) {
674         for(Int_t itrig=0; itrig<nTrigFired; ++itrig) {
675           values[kEventTrigger] = triggerIndices[itrig];
676           if((fMatchingPhiEtaCF->GetVar("trigger")<0 && itrig==0) || (fMatchingPhiEtaCF->GetVar("trigger")>=0)) {
677             for(Int_t iv=0; iv<fMatchingPhiEtaCF->GetNVar(); ++iv) valuesMatchingPhiEtaCF[iv] = values[fMatchingPhiEtaCFVars[iv]];
678             fMatchingPhiEtaCF->Fill(valuesMatchingPhiEtaCF, 1);
679           }
680           if((fMatchingPtCF->GetVar("trigger")<0 && itrig==0) || (fMatchingPtCF->GetVar("trigger")>=0)) {
681             for(Int_t iv=0; iv<fMatchingPtCF->GetNVar(); ++iv) valuesMatchingPtCF[iv] = values[fMatchingPtCFVars[iv]];
682             fMatchingPtCF->Fill(valuesMatchingPtCF, 1);
683           }
684           if(values[kTrackPt]>1.0 && values[kTrackPt]<3.0)
685             if((fBunchCrossingsCF->GetVar("trigger")<0 && itrig==0) || (fBunchCrossingsCF->GetVar("trigger")>=0)) {
686               for(Int_t iv=0; iv<fBunchCrossingsCF->GetNVar(); ++iv) valuesBCCF[iv] = values[fBunchCrossingsCFVars[iv]];
687               fBunchCrossingsCF->Fill(valuesBCCF, 1);
688             }
689           if(fExpertCF) {
690             if((fExpertCF->GetVar("trigger")<0 && itrig==0) || (fExpertCF->GetVar("trigger")>=0)) {
691               if(fExpertCF->GetStep("TRD")>=0 && fExpertCF->GetStep("TRD")<3) {
692                 for(Int_t iv=0; iv<fExpertCF->GetNVar(); ++iv) valuesExpertCF[iv] = values[fExpertCFVars[iv]];
693                 fExpertCF->Fill(valuesExpertCF, fExpertCF->GetStep("TRD"));
694               }
695             }
696           } 
697         }
698         if(Bool_t(status & AliESDtrack::kTOFpid)) {
699           for(Int_t itrig=0; itrig<nTrigFired; ++itrig) {
700             values[kEventTrigger] = triggerIndices[itrig];
701             if((fMatchingPhiEtaCF->GetVar("trigger")<0 && itrig==0) || (fMatchingPhiEtaCF->GetVar("trigger")>=0)) {
702               for(Int_t iv=0; iv<fMatchingPhiEtaCF->GetNVar(); ++iv) valuesMatchingPhiEtaCF[iv] = values[fMatchingPhiEtaCFVars[iv]];
703               fMatchingPhiEtaCF->Fill(valuesMatchingPhiEtaCF, 2);
704             }
705             if((fMatchingPtCF->GetVar("trigger")<0 && itrig==0) || (fMatchingPtCF->GetVar("trigger")>=0)) {
706               for(Int_t iv=0; iv<fMatchingPtCF->GetNVar(); ++iv) valuesMatchingPtCF[iv] = values[fMatchingPtCFVars[iv]];
707               fMatchingPtCF->Fill(valuesMatchingPtCF, 2);
708             }
709             if(values[kTrackPt]>1.0 && values[kTrackPt]<3.0)
710               if((fBunchCrossingsCF->GetVar("trigger")<0 && itrig==0) || (fBunchCrossingsCF->GetVar("trigger")>=0)) {
711                 for(Int_t iv=0; iv<fBunchCrossingsCF->GetNVar(); ++iv) valuesBCCF[iv] = values[fBunchCrossingsCFVars[iv]];
712                 fBunchCrossingsCF->Fill(valuesBCCF, 2);
713               }
714             if(fExpertCF) {
715               if((fExpertCF->GetVar("trigger")<0 && itrig==0) || (fExpertCF->GetVar("trigger")>=0)) {
716                 if(fExpertCF->GetStep("TOF")>=0 && fExpertCF->GetStep("TOF")<3) {
717                   for(Int_t iv=0; iv<fExpertCF->GetNVar(); ++iv) valuesExpertCF[iv] = values[fExpertCFVars[iv]];
718                   fExpertCF->Fill(valuesExpertCF, fExpertCF->GetStep("TOF"));
719                 }
720               }
721             }
722           }
723         }
724       }
725     }  // end if nTRDtrkl>=1
726     
727     // look at external track param
728     const AliExternalTrackParam *op = esdTrack->GetOuterParam();
729     const AliExternalTrackParam *ip = esdTrack->GetInnerParam();
730
731     Double_t pt(0.), pt0(0.), ptTRD(0.); 
732     // read MC info if available
733     Bool_t kFOUND(kFALSE), kPhysPrim(kFALSE);
734     AliMCParticle *mcParticle(NULL);
735     if(HasMC()){
736       AliTrackReference *ref(NULL); 
737       Int_t fLabel(esdTrack->GetLabel());
738       Int_t fIdx(TMath::Abs(fLabel));
739       if(!fStack || fIdx > fStack->GetNtrack()) continue; 
740       
741       // read MC particle 
742       if(!(mcParticle = (AliMCParticle*) fMC->GetTrack(fIdx))) {
743         AliWarning(Form("MC particle missing. Label[ %d].", fLabel));
744         continue;
745       }
746   
747       pt   = esdTrack->Pt();
748       pt0  = mcParticle->Pt();
749       //Double_t eta0 = mcParticle->Eta();
750       //Double_t phi0 = mcParticle->Phi();
751       kPhysPrim = fMC->IsPhysicalPrimary(fIdx);
752
753       // read track references
754       Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
755       if(!nRefs){
756         AliWarning(Form("No TR found for track @ Label[%d].", fLabel));
757         continue;
758       }
759       Int_t iref = 0;
760       while(iref<nRefs){
761         ref = mcParticle->GetTrackReference(iref);
762         if(ref->LocalX() > fgkxTPC) break;
763         ref=NULL; iref++;
764       }
765       if(ref){ 
766         if(ref->LocalX() > fgkxTOF){ // track skipping TRD fiducial volume
767           ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
768         }
769       } else { // track stopped in TPC 
770         ref = mcParticle->GetTrackReference(TMath::Max(iref-1, 0));
771       }
772       ptTRD = ref->Pt();kFOUND=kTRUE;
773     } else { // use reconstructed values
774       if(op){
775         Double_t x(op->GetX());
776         if(x<fgkxTOF && x>fgkxTPC){
777           ptTRD=op->Pt();
778           kFOUND=kTRUE;
779         }
780       }
781
782       if(!kFOUND && ip){
783         ptTRD=ip->Pt();
784         kFOUND=kTRUE;
785       }
786     }     // end if(HasMC())
787
788     if(kFOUND){
789       h = (TH2I*)fHistos->At(kTRDstat);
790       if(status & AliESDtrack::kTPCout) h->Fill(ptTRD, kTPCout);
791       if(status & AliESDtrack::kTRDin) h->Fill(ptTRD, kTRDin);
792       if(kBarrel && (status & AliESDtrack::kTRDout)) h->Fill(ptTRD, kTRDout);
793       if(kBarrel && (status & AliESDtrack::kTRDpid)) h->Fill(ptTRD, kTRDpid);
794       if(kBarrel && (status & AliESDtrack::kTRDrefit)) h->Fill(ptTRD, kTRDref);
795     }
796     Int_t idx(HasMC() ? Pdg2Idx(TMath::Abs(mcParticle->PdgCode())): 0)
797          ,sgn(esdTrack->Charge()<0?0:1);
798     if(kBarrel && kPhysPrim) {
799       TH3 *h3 = (TH3S*)fHistos->At(kPtRes);
800       Int_t offset = (status & AliESDtrack::kTRDrefit) ? 0 : 10; 
801       h3->Fill(pt0, 1.e2*(pt/pt0-1.), 
802         offset + 2*idx + sgn);
803     }
804     ((TH1*)fHistos->At(kNCl))->Fill(esdTrack->GetTRDncls(), 2*idx + sgn);
805     if(ip){
806       h = (TH2I*)fHistos->At(kTRDmom);
807       Float_t pTRD(0.);
808       for(Int_t ily=6; ily--;){
809         if((pTRD=esdTrack->GetTRDmomentum(ily))<0.) continue;
810         h->Fill(ip->GetP()-pTRD, ily);
811       }
812     }
813   }  // end loop over tracks
814   
815   
816   
817   delete [] valuesMatchingPhiEtaCF;
818   delete [] valuesMatchingPtCF;
819   delete [] valuesBCCF;
820   delete [] valuesCentCF;
821   delete [] valuesQtotCF;
822   delete [] valuesPHCF;
823   if(valuesExpertCF) delete [] valuesExpertCF;
824   
825   PostData(1, fHistos);
826 }
827
828 //____________________________________________________________________
829 TObjArray* AliTRDcheckESD::Histos()
830 {
831 // Retrieve histograms array if already build or build it
832
833   if(fHistos) return fHistos;
834
835   fHistos = new TObjArray(kNhistos+1);
836   fHistos->SetOwner(kTRUE);
837
838   TH1 *h = NULL;
839
840   // clusters per track
841   const Int_t kNpt(30);
842   Float_t pt(0.2);
843   Float_t binsPt[kNpt+1];
844   for(Int_t i=0;i<kNpt+1; i++,pt+=(TMath::Exp(i*i*.001)-1.)) binsPt[i]=pt;
845   if(!(h = (TH2I*)gROOT->FindObject("hNCl"))){
846     h = new TH2I("hNCl", "Clusters per TRD track;N_{cl}^{TRD};SPECIES;entries", 60, 0., 180., 10, -0.5, 9.5);
847     TAxis *ay(h->GetYaxis());
848     ay->SetLabelOffset(0.015);
849     for(Int_t i(0); i<AliPID::kSPECIES; i++){
850       ay->SetBinLabel(2*i+1, Form("%s^{-}", AliPID::ParticleLatexName(i)));
851       ay->SetBinLabel(2*i+2, Form("%s^{+}", AliPID::ParticleLatexName(i)));
852     }
853   } else h->Reset();
854   fHistos->AddAt(h, kNCl); fNRefFigures++;
855
856   // status bits histogram
857   const Int_t kNbits(5);
858   Float_t bits(.5);
859   Float_t binsBits[kNbits+1];
860   for(Int_t i=0; i<kNbits+1; i++,bits+=1.) binsBits[i]=bits;
861   if(!(h = (TH2I*)gROOT->FindObject("hTRDstat"))){
862     h = new TH2I("hTRDstat", "TRD status bits;p_{t} @ TRD [GeV/c];status;entries", kNpt, binsPt, kNbits, binsBits);
863     TAxis *ay(h->GetYaxis());
864     ay->SetBinLabel(1, "kTPCout");
865     ay->SetBinLabel(2, "kTRDin");
866     ay->SetBinLabel(3, "kTRDout");
867     ay->SetBinLabel(4, "kTRDpid");
868     ay->SetBinLabel(5, "kTRDrefit");
869   } else h->Reset();
870   fHistos->AddAt(h, kTRDstat);
871
872   // energy loss
873   if(!(h = (TH2I*)gROOT->FindObject("hTRDmom"))){
874     h = new TH2I("hTRDmom", "TRD energy loss;p_{inner} - p_{ly} [GeV/c];ly;entries", 100, -1., 2., 6, -0.5, 5.5);
875   } else h->Reset();
876   fHistos->AddAt(h, kTRDmom);
877   //if(!HasMC()) return fHistos;
878
879   // pt resolution
880   const Int_t kNdpt(100), kNspec(4*AliPID::kSPECIES);
881   Float_t dpt(-3.), spec(-0.5);
882   Float_t binsDPt[kNdpt+1], binsSpec[kNspec+1];
883   for(Int_t i=0; i<kNdpt+1; i++,dpt+=6.e-2) binsDPt[i]=dpt;
884   for(Int_t i=0; i<kNspec+1; i++,spec+=1.) binsSpec[i]=spec;
885   if(!(h = (TH3S*)gROOT->FindObject("hPtRes"))){
886     h = new TH3S("hPtRes", "P_{t} resolution @ DCA;p_{t}^{MC} [GeV/c];#Delta p_{t}/p_{t}^{MC} [%];SPECIES", kNpt, binsPt, kNdpt, binsDPt, kNspec, binsSpec);
887     TAxis *az(h->GetZaxis());
888     az->SetLabelOffset(0.015);
889     for(Int_t i(0); i<AliPID::kSPECIES; i++){
890       az->SetBinLabel(2*i+1, Form("%s^{-}", AliPID::ParticleLatexName(i)));
891       az->SetBinLabel(2*i+2, Form("%s^{+}", AliPID::ParticleLatexName(i)));
892       az->SetBinLabel(10+2*i+1, Form("%s^{-}", AliPID::ParticleLatexName(i)));
893       az->SetBinLabel(10+2*i+2, Form("%s^{+}", AliPID::ParticleLatexName(i)));
894     }
895   } else h->Reset();
896   fHistos->AddAt(h, kPtRes);
897
898   // TPC event vertex distribution
899   if(!(h = (TH1F*)gROOT->FindObject("hTPCVertex"))){
900     h = new TH1F("hTPCVertex", "Event vertex Z coord. from TPC tracks", 100, -25., 25.);
901   } else h->Reset();
902   fHistos->AddAt(h, kTPCVertex);
903   
904   // Event vertex
905   if(!(h = (TH1F*)gROOT->FindObject("hEventVertex"))){
906     h = new TH1F("hEventVertex", "Event vertex Z coord.", 100, -25., 25.);
907   } else h->Reset();
908   fHistos->AddAt(h, kEventVertex);
909   
910   // Number of all tracks
911   if(!(h = (TH1I*)gROOT->FindObject("hNTracksAll"))){
912     h = new TH1I("hNTracksAll", "Number of tracks per event, event vertex cuts", 5000, 0, 5000);
913   } else h->Reset();
914   fHistos->AddAt(h, kNTracksAll);
915   
916   // Number of tracks in acceptance and DCA cut
917   if(!(h = (TH1I*)gROOT->FindObject("hNTracksAcc"))){
918     h = new TH1I("hNTracksAcc", Form("Number of tracks per event, |#eta|<%.1f, |DCAxy|<%.1f, |DCAz|<%.1f",
919                                      fgkEta, fgkTrkDCAxy, fgkTrkDCAz), 5000, 0, 5000);
920   } else h->Reset();
921   fHistos->AddAt(h, kNTracksAcc);
922   
923   // Number of tracks in TPC (Ncls>10)
924   if(!(h = (TH1I*)gROOT->FindObject("hNTracksTPC"))){
925     h = new TH1I("hNTracksTPC", Form("Number of tracks per event, |#eta|<%.1f, pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",
926                                      fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 5000, 0, 5000);
927   } else h->Reset();
928   fHistos->AddAt(h, kNTracksTPC);
929   
930   // Distribution of DCA-xy
931   if(!(h = (TH1F*)gROOT->FindObject("hDCAxy"))){
932     h = new TH1F("hDCAxy", "Distribution of transverse DCA", 100, -100., 100.);
933   } else h->Reset();
934   fHistos->AddAt(h, kDCAxy);
935   
936   // Distribution of DCA-z
937   if(!(h = (TH1F*)gROOT->FindObject("hDCAz"))){
938     h = new TH1F("hDCAz", "Distribution of longitudinal DCA", 100, -100., 100.);
939   } else h->Reset();
940   fHistos->AddAt(h, kDCAz);
941   
942   Double_t binPtLimits[33] = {0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9,
943                              1.0, 1.1, 1.2, 1.3, 1.4, 
944                              1.6, 1.8, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0,
945                              3.4, 3.8, 4.2, 4.6, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0};
946   // Pt distributions
947   if(!(h = (TH1F*)gROOT->FindObject("hPt1"))){
948     h = new TH1F("hPt1", Form("dN/dpt, |#eta|<%.1f and pt>%.1f", fgkEta, fgkPt), 32, binPtLimits);
949   } else h->Reset();
950   fHistos->AddAt(h, kPt1);
951   
952   if(!(h = (TH1F*)gROOT->FindObject("hPt2"))){
953     h = new TH1F("hPt2", Form("dN/dpt, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f",
954                               fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz), 32, binPtLimits);
955   } else h->Reset();
956   fHistos->AddAt(h, kPt2);
957   
958   if(!(h = (TH1F*)gROOT->FindObject("hPt3pos"))){
959     h = new TH1F("hPt3pos", Form("dN/dpt (positives), |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",
960                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 32, binPtLimits);
961   } else h->Reset();
962   fHistos->AddAt(h, kPt3pos);
963   
964   if(!(h = (TH1F*)gROOT->FindObject("hPt3neg"))){
965     h = new TH1F("hPt3neg", Form("dN/dpt (negatives), |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",
966                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 32, binPtLimits);
967   } else h->Reset();
968   fHistos->AddAt(h, kPt3neg);
969   
970   if(!(h = (TH1F*)gROOT->FindObject("hPt4pos"))){
971     h = new TH1F("hPt4pos", Form("dN/dpt (positives), |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d, nTRDtracklets#geq 1",
972                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 32, binPtLimits);
973   } else h->Reset();
974   fHistos->AddAt(h, kPt4pos);
975   
976   if(!(h = (TH1F*)gROOT->FindObject("hPt4neg"))){
977     h = new TH1F("hPt4pos", Form("dN/dpt (negatives), |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d, nTRDtracklets#geq 1",
978                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 32, binPtLimits);
979   } else h->Reset();
980   fHistos->AddAt(h, kPt4neg);
981   
982   // theta distribution of TRD tracks
983   if(!(h = (TH1F*)gROOT->FindObject("hTheta"))){
984     h = new TH1F("hTheta", Form("dN/d#theta, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d, nTRDtracklets#geq 1",
985                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 220,.5,2.7);
986   } else h->Reset();
987   fHistos->AddAt(h, kTheta);
988   
989   // phi distribution of TRD tracks
990   if(!(h = (TH1F*)gROOT->FindObject("hPhi"))){
991     h = new TH1F("hPhi", Form("dN/d#varphi, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d, nTRDtracklets#geq 1",
992                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 157,0,6.28);
993   } else h->Reset();
994   fHistos->AddAt(h, kPhi);
995   
996   // TPC cluster distribution
997   if(!(h = (TH1F*)gROOT->FindObject("hNTPCCl"))){
998     h = new TH1I("hNTPCCl", Form("Number of TPC clusters/track, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f",
999                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz), 160, 0, 160);
1000   } else h->Reset();
1001   fHistos->AddAt(h, kNTPCCl);
1002   
1003   if(!(h = (TH1I*)gROOT->FindObject("hNTPCCl2"))){
1004     h = new TH1F("hNTPCCl2", Form("Number of TPC clusters/track, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, pt>1.0 GeV/c",
1005                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz), 160, 0, 160);
1006   } else h->Reset();
1007   fHistos->AddAt(h, kNTPCCl2);
1008   
1009   // dE/dx vs P for TPC reference tracks
1010   if(!(h = (TH2F*)gROOT->FindObject("hTPCDedx"))){
1011     h = new TH2F("hTPCDedx", Form("TPC dE/dx vs P, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",
1012                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100, 0.1,10.1, 150, 0, 150.);
1013   } else h->Reset();
1014   fHistos->AddAt(h, kTPCDedx);
1015   
1016   // eta,phi distribution of TPC reference tracks
1017   if(!(h = (TH2F*)gROOT->FindObject("hEtaPhi"))){
1018     h = new TH2F("hEtaPhi", Form("TPC (#eta,#varphi), |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",
1019                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 50, -1, 1, 157, 0, 6.28);
1020   } else h->Reset();
1021   fHistos->AddAt(h, kEtaPhi);
1022   
1023   // Nclusters vs eta distribution for TPC tracks
1024   if(!(h = (TH2F*)gROOT->FindObject("hEtaNclsTPC"))){
1025     h = new TH2F("hEtaNclsTPC", Form("TPC Nclusters vs. #eta, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f",
1026                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz), 50, -1, 1, 160, 0, 160.);
1027   } else h->Reset();
1028   fHistos->AddAt(h, kEtaNclsTPC);
1029   
1030   // Nclusters vs phi distribution for TPC reference tracks
1031   if(!(h = (TH2F*)gROOT->FindObject("hPhiNclsTPC"))){
1032     h = new TH2F("hPhiNclsTPC", Form("TPC Nclusters vs. #varphi, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f",
1033                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz), 157, 0, 6.28, 160, 0, 160.);
1034   } else h->Reset();
1035   fHistos->AddAt(h, kPhiNclsTPC);
1036
1037   // SPD multiplicity distribution
1038   if(!(h = (TH1F*)gROOT->FindObject("hSPDMult"))){
1039     h = new TH1F("hSPDMult", "SPD multiplicity", 10000, -0.5, 9999.5);
1040   } else h->Reset();
1041   fHistos->AddAt(h, kSPDMult);
1042
1043   // Ntracklets/track vs P for TRD reference tracks
1044   Double_t binsP[19] = {0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.7, 2.0,
1045                         2.5, 3.0, 3.5, 4.0, 5.0, 6.0, 7.0, 9.0, 12.0};
1046   for(Int_t iCent=0; iCent<=5; ++iCent) {
1047     if(!(h = (TH2F*)gROOT->FindObject(Form("hNTrackletsTRD_cent%d",iCent+1)))){
1048       h = new TH2F(Form("hNTrackletsTRD_cent%d",iCent+1), Form("TRD Ntracklets/track vs. P, centrality %d, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",
1049                                                                iCent+1, fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 18, binsP, 7, -0.5, 6.5);
1050     } else h->Reset();
1051     fHistos->AddAt(h, kNTrackletsTRD+iCent);
1052   }
1053   
1054   // Nclusters/track vs P for TRD reference tracks
1055   for(Int_t iCent=0; iCent<=5; ++iCent) {
1056     if(!(h = (TH2F*)gROOT->FindObject(Form("hNClsTrackTRD_cent%d",iCent+1)))){
1057       h = new TH2F(Form("hNClsTrackTRD_cent%d",iCent+1), Form("TRD Nclusters/track vs. P, centrality %d, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",
1058                                                               iCent+1, fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 18, binsP, 300, 0., 300.);
1059     } else h->Reset();
1060     fHistos->AddAt(h, kNClsTrackTRD+iCent);
1061   }  
1062
1063   // <PH> vs slice number for TRD reference tracklets
1064   for(Int_t iCent=0; iCent<=5; ++iCent) {
1065     if(!(h = (TH2F*)gROOT->FindObject(Form("hPHSlice_cent%d",iCent+1)))){
1066       h = new TH2F(Form("hPHSlice_cent%d",iCent+1), Form("<PH> vs sliceNo, centrality %d, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",
1067                                     iCent+1, fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 8, -0.5, 7.5, 200, 0., 2000.);
1068     } else h->Reset();
1069     fHistos->AddAt(h, kPHSlice+iCent);
1070   }
1071
1072   // <PH> vs slice number for TRD reference tracklets, from TPC pions
1073   for(Int_t iCent=0; iCent<=5; ++iCent) {
1074     if(!(h = (TH2F*)gROOT->FindObject(Form("hPHSliceTPCpions_cent%d",iCent+1)))){
1075       h = new TH2F(Form("hPHSliceTPCpions_cent%d",iCent+1), Form("<PH> vs sliceNo, centrality %d, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d, TPC pions",
1076                                                                  iCent+1, fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 8, -0.5, 7.5, 100, 0., 2000.);
1077     } else h->Reset();
1078     fHistos->AddAt(h, kPHSliceTPCpions+iCent);
1079   }
1080
1081   // TPC dE/dx vs P for TRD reference tracks, pions
1082   for(Int_t iCent=0; iCent<=5; ++iCent) {
1083     if(!(h = (TH2F*)gROOT->FindObject(Form("hTPCdedxPions_cent%d",iCent+1)))){
1084       h = new TH2F(Form("hTPCdedxPions_cent%d",iCent+1), Form("TPC dE/dx vs P, TPC pions, centrality %d, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",
1085                                          iCent+1, fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100, 0.1,10.1, 100, 0,100.);
1086     } else h->Reset();
1087     fHistos->AddAt(h, kTPCdedxPions+iCent);
1088   }
1089
1090   // <PH> vs slice number for TRD reference tracklets, from TPC electrons
1091   for(Int_t iCent=0; iCent<=5; ++iCent) {
1092     if(!(h = (TH2F*)gROOT->FindObject(Form("hPHSliceTPCelectrons_cent%d",iCent+1)))){
1093       h = new TH2F(Form("hPHSliceTPCelectrons_cent%d",iCent+1), Form("<PH> vs sliceNo, centrality %d,|#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d, TPC electrons",
1094                                                 iCent+1, fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 8, -0.5, 7.5, 100, 0., 2000.);
1095     } else h->Reset();
1096     fHistos->AddAt(h, kPHSliceTPCelectrons+iCent);
1097   }
1098
1099   // TPC dE/dx vs P for TRD reference tracks, electrons
1100   for(Int_t iCent=0; iCent<=5; ++iCent) {
1101     if(!(h = (TH2F*)gROOT->FindObject(Form("hTPCdedxElectrons_cent%d",iCent+1)))){
1102       h = new TH2F(Form("hTPCdedxElectrons_cent%d",iCent+1), Form("TPC dE/dx vs P, TPC electrons, centrality %d, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",
1103                                              iCent+1, fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100, 0.1,10.1, 100, 0,100.);
1104     } else h->Reset();
1105     fHistos->AddAt(h, kTPCdedxElectrons+iCent);
1106   }
1107
1108   // Qtot vs P for TRD reference tracklets
1109   for(Int_t iCent=0; iCent<=5; ++iCent) {
1110     if(!(h = (TH2F*)gROOT->FindObject(Form("hQtotP_cent%d",iCent+1)))){
1111       h = new TH2F(Form("hQtotP_cent%d",iCent+1), Form("Qtot(from slices) vs P, centrality %d, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",
1112                                   iCent+1, fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 18, binsP, 400, 0., 20);
1113     } else h->Reset();
1114     fHistos->AddAt(h, kQtotP+iCent);
1115   }
1116   
1117   // (X,Y,momentum) distribution after AliESDtrack::PropagateTo(r=300.)
1118   if(!(h = (TH3F*)gROOT->FindObject("hPropagXYvsP"))){
1119     h = new TH3F("hPropagXYvsP", Form("(x,y) vs P after AliESDtrack::PropagateTo(r=300.), |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",
1120                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100,-500,500, 100,-500,500, 10, 0.,10.);
1121   } else h->Reset();
1122   fHistos->AddAt(h, kPropagXYvsP);
1123   
1124   // (R,Z,momentum) distribution after AliESDtrack::PropagateTo(r=300.)
1125   if(!(h = (TH3F*)gROOT->FindObject("hPropagRZvsP"))){
1126     h = new TH3F("hPropagRZvsP", Form("(r,z) vs P after AliESDtrack::PropagateTo(r=300.), |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",
1127                  fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100,-350., 350., 100,0.,500., 10, 0.,10.);
1128   } else h->Reset();
1129   fHistos->AddAt(h, kPropagRZvsP);
1130   
1131   Double_t etaBinLimits[101];   
1132   for(Int_t i=0; i<101; i++) etaBinLimits[i] = -1.0 + i*2.0/100.;
1133   Double_t phiBinLimits[151];
1134   for(Int_t i=0; i<151; i++) phiBinLimits[i] = -1.1*TMath::Pi() + i*2.2*TMath::Pi()/150.;
1135   // (eta,detector phi,P) distribution of reference TPC positive tracks
1136   for(Int_t iCent=0; iCent<=5; ++iCent) {
1137     if(!(h = (TH3F*)gROOT->FindObject(Form("hTPCRefTracksPos_cent%d",iCent+1)))){
1138       h = new TH3F(Form("hTPCRefTracksPos_cent%d",iCent+1), Form("(#eta,detector #varphi,p) for TPC positive reference tracks, centrality %d, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",
1139                                             iCent+1, fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100, etaBinLimits, 150, phiBinLimits, 32, binPtLimits);
1140     } else h->Reset();
1141     fHistos->AddAt(h, kTPCRefTracksPos+iCent);
1142   }
1143   
1144   // (eta,detector phi,P) distribution of reference TPC negative tracks
1145   for(Int_t iCent=0; iCent<=5; ++iCent) {
1146     if(!(h = (TH3F*)gROOT->FindObject(Form("hTPCRefTracksNeg_cent%d",iCent+1)))){
1147       h = new TH3F(Form("hTPCRefTracksNeg_cent%d",iCent+1), Form("(#eta,detector #varphi,p) for TPC negative reference tracks, centrality %d, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d",
1148                                             iCent+1, fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100, etaBinLimits, 150, phiBinLimits, 32, binPtLimits);
1149     } else h->Reset();
1150     fHistos->AddAt(h, kTPCRefTracksNeg+iCent);
1151   }
1152   
1153   // (eta,detector phi,P) distribution of reference TRD positive tracks
1154   for(Int_t iCent=0; iCent<=5; ++iCent) {
1155     if(!(h = (TH3F*)gROOT->FindObject(Form("hTRDRefTracksPos_cent%d",iCent+1)))){
1156       h = new TH3F(Form("hTRDRefTracksPos_cent%d",iCent+1), Form("(#eta,detector #varphi,p) for TRD positive reference tracks, centrality %d, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d, nTRDtracklets#geq1",
1157                                             iCent+1, fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100, etaBinLimits, 150, phiBinLimits, 32, binPtLimits);
1158     } else h->Reset();
1159     fHistos->AddAt(h, kTRDRefTracksPos+iCent);
1160   }
1161   
1162   // (eta,detector phi,P) distribution of reference TRD negative tracks
1163   for(Int_t iCent=0; iCent<=5; ++iCent) {
1164     if(!(h = (TH3F*)gROOT->FindObject(Form("hTRDRefTracksNeg_cent%d",iCent+1)))){
1165       h = new TH3F(Form("hTRDRefTracksNeg_cent%d",iCent+1), Form("(#eta,detector #varphi,p) for TRD negative reference tracks, centrality %d, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d, nTRDtracklets#geq1",
1166                                             iCent+1, fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100, etaBinLimits, 150, phiBinLimits, 32, binPtLimits);
1167     } else h->Reset();
1168     fHistos->AddAt(h, kTRDRefTracksNeg+iCent);
1169   }
1170
1171   // (eta,detector phi,P) distribution of reference TRD positive tracks with 4 tracklets
1172   for(Int_t iCent=0; iCent<=5; ++iCent) {
1173     if(!(h = (TH3F*)gROOT->FindObject(Form("hTRDRefTracksPos4_cent%d",iCent+1)))){
1174       h = new TH3F(Form("hTRDRefTracksPos4_cent%d",iCent+1), Form("(#eta,detector #varphi,p) for TRD positive reference tracks, centrality %d, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d, nTRDtracklets#geq1",
1175                                                                  iCent+1, fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100, etaBinLimits, 150, phiBinLimits, 32, binPtLimits);
1176     } else h->Reset();
1177     fHistos->AddAt(h, kTRDRefTracksPos4+iCent);
1178   }
1179
1180   // (eta,detector phi,P) distribution of reference TRD positive tracks with 5 tracklets
1181   for(Int_t iCent=0; iCent<=5; ++iCent) {
1182     if(!(h = (TH3F*)gROOT->FindObject(Form("hTRDRefTracksPos5_cent%d",iCent+1)))){
1183       h = new TH3F(Form("hTRDRefTracksPos5_cent%d",iCent+1), Form("(#eta,detector #varphi,p) for TRD positive reference tracks, centrality %d, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d, nTRDtracklets#geq1",
1184                                                                   iCent+1, fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100, etaBinLimits, 150, phiBinLimits, 32, binPtLimits);
1185     } else h->Reset();
1186     fHistos->AddAt(h, kTRDRefTracksPos5+iCent);
1187   }
1188
1189   // (eta,detector phi,P) distribution of reference TRD positive tracks with 6 tracklets
1190   for(Int_t iCent=0; iCent<=5; ++iCent) {
1191     if(!(h = (TH3F*)gROOT->FindObject(Form("hTRDRefTracksPos6_cent%d",iCent+1)))){
1192       h = new TH3F(Form("hTRDRefTracksPos6_cent%d",iCent+1), Form("(#eta,detector #varphi,p) for TRD positive reference tracks, centrality %d, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d, nTRDtracklets#geq1",
1193                                                                   iCent+1, fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100, etaBinLimits, 150, phiBinLimits, 32, binPtLimits);
1194     } else h->Reset();
1195     fHistos->AddAt(h, kTRDRefTracksPos6+iCent);
1196   }
1197   
1198   // (eta,detector phi,P) distribution of reference TRD negative tracks with 4 tracklets
1199   for(Int_t iCent=0; iCent<=5; ++iCent) {
1200     if(!(h = (TH3F*)gROOT->FindObject(Form("hTRDRefTracksNeg4_cent%d",iCent+1)))){
1201       h = new TH3F(Form("hTRDRefTracksNeg4_cent%d",iCent+1), Form("(#eta,detector #varphi,p) for TRD negative reference tracks, centrality %d, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d, nTRDtracklets#geq1",
1202                                                                  iCent+1, fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100, etaBinLimits, 150, phiBinLimits, 32, binPtLimits);
1203     } else h->Reset();
1204     fHistos->AddAt(h, kTRDRefTracksNeg4+iCent);
1205   }
1206
1207   // (eta,detector phi,P) distribution of reference TRD negative tracks with 5 tracklets
1208   for(Int_t iCent=0; iCent<=5; ++iCent) {
1209     if(!(h = (TH3F*)gROOT->FindObject(Form("hTRDRefTracksNeg5_cent%d",iCent+1)))){
1210       h = new TH3F(Form("hTRDRefTracksNeg5_cent%d",iCent+1), Form("(#eta,detector #varphi,p) for TRD negative reference tracks, centrality %d, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d, nTRDtracklets#geq1",
1211                                                                   iCent+1, fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100, etaBinLimits, 150, phiBinLimits, 32, binPtLimits);
1212     } else h->Reset();
1213     fHistos->AddAt(h, kTRDRefTracksNeg5+iCent);
1214   }
1215
1216   // (eta,detector phi,P) distribution of reference TRD negative tracks with 6 tracklets
1217   for(Int_t iCent=0; iCent<=5; ++iCent) {
1218     if(!(h = (TH3F*)gROOT->FindObject(Form("hTRDRefTracksNeg6_cent%d",iCent+1)))){
1219       h = new TH3F(Form("hTRDRefTracksNeg6_cent%d",iCent+1), Form("(#eta,detector #varphi,p) for TRD negative reference tracks, centrality %d, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d, nTRDtracklets#geq1",
1220                                                                   iCent+1, fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100, etaBinLimits, 150, phiBinLimits, 32, binPtLimits);
1221     } else h->Reset();
1222     fHistos->AddAt(h, kTRDRefTracksNeg6+iCent);
1223   }
1224
1225
1226   // (eta,detector phi) profile of average number of TRD tracklets/track
1227   for(Int_t iCent=0; iCent<=5; ++iCent) {
1228     if(!(h = (TProfile2D*)gROOT->FindObject(Form("hTRDEtaPhiAvNtrkl_cent%d",iCent+1)))){
1229       h = new TProfile2D(Form("hTRDEtaPhiAvNtrkl_cent%d",iCent+1), Form("<Ntracklets/track> vs (#eta,detector #varphi) for TRD reference tracks, centrality %d, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d, nTRDtracklets#geq1",
1230                                                    iCent+1, fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100, -1.0, 1.0, 150, -1.1*TMath::Pi(), 1.1*TMath::Pi());
1231     } else h->Reset();
1232     fHistos->AddAt(h, kTRDEtaPhiAvNtrkl+iCent);
1233   }
1234
1235   // (eta,delta phi) profile of average number of TRD tracklets/track
1236   for(Int_t iCent=0; iCent<=5; ++iCent) {
1237     if(!(h = (TProfile2D*)gROOT->FindObject(Form("hTRDEtaDeltaPhiAvNtrkl_cent%d",iCent+1)))){
1238       h = new TProfile2D(Form("hTRDEtaDeltaPhiAvNtrkl_cent%d",iCent+1), Form("<Ntracklets/track> vs (#eta, #Delta#varphi) for TRD reference tracks, centrality %d, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d, nTRDtracklets#geq1",
1239                                                         iCent+1, fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100, -1.0, 1.0, 50, -0.4*TMath::Pi(), 0.4*TMath::Pi());
1240     } else h->Reset();
1241     fHistos->AddAt(h, kTRDEtaDeltaPhiAvNtrkl+iCent);
1242   }  
1243
1244   // (eta, detector phi) profile of average tracklet Qtot from slices
1245   for(Int_t iCent=0; iCent<=5; ++iCent) {
1246     for(Int_t iLayer=0;iLayer<6;iLayer++) {
1247       if(!(h = (TProfile2D*)gROOT->FindObject(Form("hTRDEtaPhiAvQtot_Layer%d_cent%d",iLayer,iCent+1)))) {
1248         h = new TProfile2D(Form("hTRDEtaPhiAvQtot_Layer%d_cent%d",iLayer,iCent+1),
1249                            Form("<Q_{tot}> vs (#eta, detector #varphi) for TRD reference tracks (layer %d), centrality %d, |#eta|<%.1f and pt>%.1f, |DCAxy|<%.1f, |DCAz|<%.1f, TPC nclusters>%d, nTRDtracklets#geq1",
1250                                 iLayer, iCent+1, fgkEta, fgkPt, fgkTrkDCAxy, fgkTrkDCAz, fgkNclTPC), 100, -1.0, 1.0, 150, -1.1*TMath::Pi(), 1.1*TMath::Pi());
1251       } else h->Reset();
1252       fHistos->AddAt(h, kTRDEtaPhiAvQtot+iCent*6+iLayer);
1253     }
1254   }
1255   
1256   // Trigger definitions
1257   if(!(h=(TH1F*)gROOT->FindObject("hTriggerDefs"))) {
1258     h = new TH1F("hTriggerDefs", "Trigger definitions", kNMaxAssignedTriggers, 0.5, 0.5+Float_t(kNMaxAssignedTriggers));
1259   }
1260   else h->Reset();
1261   fHistos->AddAt(h, kTriggerDefs);
1262
1263   // dummy histo
1264   if(!(h=(TH1F*)gROOT->FindObject("hDummy"))) {
1265     h = new TH1F("hDummy", "Dummy hist", 10, 0., 1.);
1266   }
1267   else h->Reset();
1268   fHistos->AddAt(h, 0);
1269   
1270   fMatchingPhiEtaCF = CreateCFContainer("MatchingPhiEta", "CF container with TRD-TPC matching data");
1271   fHistos->AddAt(fMatchingPhiEtaCF, kMatchingPhiEtaCF);
1272   fMatchingPtCF = CreateCFContainer("MatchingPt", "CF container with TRD-TPC matching data");
1273   fHistos->AddAt(fMatchingPtCF, kMatchingPtCF);
1274   fBunchCrossingsCF = CreateCFContainer("BunchCrossingsCF", "CF container with bunch crossings dependent data");
1275   fHistos->AddAt(fBunchCrossingsCF, kBunchCrossingsCF);
1276   fCentralityCF = CreateCFContainer("CentralityCF", "CF container with TRD-TPC matching data");
1277   fHistos->AddAt(fCentralityCF, kCentralityCF);
1278   fQtotCF = CreateCFContainer("QtotCF", "CF container with TRD tracklet charge data");
1279   fHistos->AddAt(fQtotCF, kQtotCF);
1280   fPulseHeightCF = CreateCFContainer("PulseHeightCF", "CF container with TRD tracklet PH data");
1281   fHistos->AddAt(fPulseHeightCF, kPulseHeightCF);
1282   fExpertCF = CreateCFContainer("ExpertCF", "CF container with customized information");
1283   if(fExpertCF) fHistos->AddAt(fExpertCF, kExpertCF);
1284   
1285   return fHistos;
1286 }
1287
1288
1289 //__________________________________________________________________________________________________________
1290 void AliTRDcheckESD::InitializeCFContainers() {
1291   //
1292   //  Initialize the CF container
1293   //
1294   AliAnalysisManager* man=AliAnalysisManager::GetAnalysisManager();
1295   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
1296   if(!inputHandler) return;
1297
1298   GetTriggerIndex("All triggers", kTRUE);
1299
1300   AliPhysicsSelection* physSel = (AliPhysicsSelection*)inputHandler->GetEventSelection();
1301   const TList* trigList = (physSel ? physSel->GetCollisionTriggerClasses() : 0x0);
1302   const TList* bgTrigList = (physSel ? physSel->GetBGTriggerClasses() : 0x0);
1303   
1304   // Add collision triggers from PhysicsSelection
1305   if(trigList) {
1306     for(Int_t it=0; it<trigList->GetEntries(); ++it) {
1307       TString trigName = trigList->At(it)->GetName();
1308       TObjArray* arr = trigName.Tokenize(" ");
1309       trigName = arr->At(0)->GetName();
1310       trigName.Remove(0,1);
1311       TObjArray* arr2 = trigName.Tokenize(",");
1312       for(Int_t jt=0; jt<arr2->GetEntries(); ++jt) {
1313         // Assign an index into the trigger histogram and the CF container for this trigger
1314         GetTriggerIndex(arr2->At(jt)->GetName(), kTRUE);
1315       }
1316     }
1317   }
1318   // Add background triggers from PhysicsSelection
1319   if(bgTrigList) {
1320     for(Int_t it=0; it<bgTrigList->GetEntries(); ++it) {
1321       TString trigName = bgTrigList->At(it)->GetName();
1322       TObjArray* arr = trigName.Tokenize(" ");
1323       trigName = arr->At(0)->GetName();
1324       trigName.Remove(0,1);
1325       TObjArray* arr2 = trigName.Tokenize(",");
1326       for(Int_t jt=0; jt<arr2->GetEntries(); ++jt) {
1327         // Assign an index into the trigger histogram and the CF container for this trigger
1328         GetTriggerIndex(arr2->At(jt)->GetName(), kTRUE);
1329       }
1330     }
1331   }
1332     
1333   // Add user enabled triggers
1334   TObjArray* arr = fUserEnabledTriggers.Tokenize(";");
1335   for(Int_t it=0; it<arr->GetEntries(); ++it) {
1336     GetTriggerIndex(arr->At(it)->GetName(), kTRUE);
1337   }
1338 }
1339
1340
1341 //__________________________________________________________________________________________________________
1342 AliCFContainer* AliTRDcheckESD::CreateCFContainer(const Char_t* name, const Char_t* title) {
1343   //
1344   //  make a CF container
1345   //
1346   // create a CF container and add it to the list of histograms
1347   Int_t nbinsCf[kNTrdCfVariables];
1348   for(Int_t i=0;i<kNTrdCfVariables;++i) nbinsCf[i]=0;
1349   nbinsCf[kEventVtxZ]         =    5;
1350   nbinsCf[kEventMult]         =    5;
1351   nbinsCf[kEventTrigger]      =    kNMaxAssignedTriggers;
1352   nbinsCf[kEventBC]           = 3500;
1353   nbinsCf[kTrackTOFBC]        =    2;
1354   nbinsCf[kTrackDCAxy]        =    9;
1355   nbinsCf[kTrackDCAz]         =    5;
1356   nbinsCf[kTrackCharge]       =    2;
1357   nbinsCf[kTrackPhi]          =  180;
1358   nbinsCf[kTrackEta]          =   90;
1359   nbinsCf[kTrackPt]           =   18;
1360   nbinsCf[kTrackP]            =   17;
1361   nbinsCf[kTrackTrdTracklets] =    7;
1362   nbinsCf[kTrackTrdClusters]  =  200;
1363   for(Int_t i=0;i<6;++i) nbinsCf[kTrackQtot+i] = 100;
1364   for(Int_t i=0;i<8;++i) nbinsCf[kTrackPHslice+i] = 400;
1365   //Double_t evVtxLims[2]      = {-10.,+10.};
1366   Double_t evMultLims[6]     = {0.0, 700., 1400., 2100., 2800., 3500.};
1367   Double_t evTriggerLims[2]  = {0.5, 0.5+Float_t(kNMaxAssignedTriggers)};
1368   Double_t evBCLims[2]       = {-0.5, +3499.5};
1369   //Double_t trkTOFBClims[3]   = {-0.5, 0.5, 5.5};
1370   //Double_t trkDCAxyLims[10]  = {-10.0,  -6.0, -3.0,  -2.0,  -1.0,   
1371   //                               +1.0, +2.0,  +3.0,  +6.0, +10.0};    
1372   //Double_t trkDCAzLims[2]    = {-15.0, +15.0};    
1373   Double_t trkChargeLims[2]  = {-1.5, +1.5};
1374   Double_t trkPhiLims[2]     = {-1.001*TMath::Pi(), +1.001*TMath::Pi()};
1375   Double_t trkEtaLims[2]     = {-0.9, +0.9};
1376   Double_t trkPtLims[19]     = {0.0, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75, 2.0, 2.5, 3.0, 
1377                                 3.5, 4.0, 4.5,  5.0, 6.0,  7.0, 8.0,  9.0, 10.0};
1378   Double_t trkPLims[18]      = {0.0, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.7, 2.0, 2.5, 
1379                                 3.0, 3.5, 4.0, 5.0, 6.0, 7.0, 9.0, 12.0};
1380   Double_t trkTrdNLims[2]    = {-0.5, 6.5};
1381   Double_t trkTrdNclsLims[2] = {-0.5, 199.5};
1382   Double_t trkQtotLims[2]    = {0.0, 20.};
1383   /*Double_t trkQtotLims[20]   = {0.0, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0,
1384                                 5.5, 6.0, 6.5, 7.0, 8.0, 9.0,10.0,12.0,15.0,20.0};*/
1385   const Char_t* varNames[kNTrdCfVariables] = {"vtxZ", "multiplicity", "trigger", "BC", "TOFBC", "DCAxy", "DCAz",
1386     "charge", "phi", "eta", "pt", "P", "tracklets", "clusters", 
1387     "PH0", "PH1", "PH2", "PH3", "PH4", "PH5", "PH6", "PH7",
1388     "Qtot0", "Qtot1", "Qtot2", "Qtot3", "Qtot4", "Qtot5"};
1389   
1390   AliCFContainer* cf;
1391   TString nameStr=name;
1392   if(nameStr.Contains("MatchingPhiEta")) {
1393     fMatchingPhiEtaCFVars[0] = kTrackCharge;       fMatchingPhiEtaCFVars[1] = kTrackPhi; fMatchingPhiEtaCFVars[2] = kTrackEta;
1394     fMatchingPhiEtaCFVars[3] = kTrackTrdTracklets; fMatchingPhiEtaCFVars[4] = kEventTrigger;
1395     const Int_t nVars = 5;
1396     Int_t nBins[nVars]; for(Int_t i=0; i<nVars; ++i) nBins[i] = nbinsCf[fMatchingPhiEtaCFVars[i]];
1397     cf = new AliCFContainer(name, title, 3, nVars, nBins);
1398     cf->SetBinLimits(0, trkChargeLims[0], trkChargeLims[1]);
1399     cf->SetBinLimits(1, trkPhiLims[0], trkPhiLims[1]);
1400     cf->SetBinLimits(2, trkEtaLims[0], trkEtaLims[1]);
1401     cf->SetBinLimits(3, trkTrdNLims[0], trkTrdNLims[1]);
1402     cf->SetBinLimits(4, evTriggerLims[0], evTriggerLims[1]);
1403     for(Int_t i=0; i<nVars; ++i) cf->SetVarTitle(i, varNames[fMatchingPhiEtaCFVars[i]]);
1404     cf->SetStepTitle(0, "TPC");
1405     cf->SetStepTitle(1, "TRD");
1406     cf->SetStepTitle(2, "TOF");
1407     return cf;
1408   }
1409   if(nameStr.Contains("MatchingPt")) {
1410     fMatchingPtCFVars[0] = kEventMult; fMatchingPtCFVars[1] = kTrackCharge;        fMatchingPtCFVars[2] = kTrackPhi; 
1411     fMatchingPtCFVars[3] = kTrackPt;   fMatchingPtCFVars[4] = kTrackTrdTracklets;  fMatchingPtCFVars[5] = kEventTrigger;
1412     const Int_t nVars = 6;
1413     Int_t nBins[nVars]; for(Int_t i=0; i<nVars; ++i) nBins[i] = nbinsCf[fMatchingPtCFVars[i]];
1414     cf = new AliCFContainer(name, title, 3, nVars, nBins);
1415     cf->SetBinLimits(0, evMultLims);
1416     cf->SetBinLimits(1, trkChargeLims[0], trkChargeLims[1]);
1417     cf->SetBinLimits(2, trkPhiLims[0], trkPhiLims[1]);
1418     cf->SetBinLimits(3, trkPtLims);
1419     cf->SetBinLimits(4, trkTrdNLims[0], trkTrdNLims[1]);
1420     cf->SetBinLimits(5, evTriggerLims[0], evTriggerLims[1]);
1421     for(Int_t i=0; i<nVars; ++i) cf->SetVarTitle(i, varNames[fMatchingPtCFVars[i]]);
1422     cf->SetStepTitle(0, "TPC");
1423     cf->SetStepTitle(1, "TRD");
1424     cf->SetStepTitle(2, "TOF");
1425     return cf;
1426   }
1427   if(nameStr.Contains("BunchCrossings")) {
1428     fBunchCrossingsCFVars[0] = kEventBC; fBunchCrossingsCFVars[1] = kTrackPhi;
1429     const Int_t nVars = 2;
1430     Int_t nBins[nVars]; for(Int_t i=0; i<nVars; ++i) nBins[i] = nbinsCf[fBunchCrossingsCFVars[i]];
1431     cf = new AliCFContainer(name, title, 3, nVars, nBins);
1432     cf->SetBinLimits(0, evBCLims[0], evBCLims[1]);
1433     cf->SetBinLimits(1, trkPhiLims[0], trkPhiLims[1]);
1434     for(Int_t i=0; i<nVars; ++i) cf->SetVarTitle(i, varNames[fBunchCrossingsCFVars[i]]);
1435     cf->SetStepTitle(0, "TPC");
1436     cf->SetStepTitle(1, "TRD");
1437     cf->SetStepTitle(2, "TOF");
1438     return cf;
1439   }
1440   if(nameStr.Contains("Centrality")) {
1441     fCentralityCFVars[0] = kEventMult; fCentralityCFVars[1] = kTrackP; fCentralityCFVars[2] = kTrackTrdClusters; 
1442     fCentralityCFVars[3] = kTrackQtot; fCentralityCFVars[4] = kEventTrigger;
1443     const Int_t nVars = 5;
1444     Int_t nBins[nVars]; for(Int_t i=0; i<nVars; ++i) nBins[i] = nbinsCf[fCentralityCFVars[i]];
1445     cf = new AliCFContainer(name, title, 1, nVars, nBins);
1446     cf->SetBinLimits(0, evMultLims);
1447     cf->SetBinLimits(1, trkPLims);
1448     cf->SetBinLimits(2, trkTrdNclsLims[0], trkTrdNclsLims[1]);
1449     cf->SetBinLimits(3, trkQtotLims[0], trkQtotLims[1]);
1450     cf->SetBinLimits(4, evTriggerLims[0], evTriggerLims[1]);
1451     for(Int_t i=0; i<nVars; ++i) cf->SetVarTitle(i, varNames[fCentralityCFVars[i]]);
1452     cf->SetStepTitle(0, "TRD");
1453     return cf;
1454   }
1455   if(nameStr.Contains("Qtot")) {
1456     fQtotCFVars[0] = kTrackPhi; fQtotCFVars[1] = kTrackEta; fQtotCFVars[2] = kTrackQtot; 
1457     fQtotCFVars[3] = kEventTrigger;
1458     const Int_t nVars = 5;
1459     Int_t nBins[nVars]; for(Int_t i=0; i<nVars-1; ++i) nBins[i] = nbinsCf[fQtotCFVars[i]];
1460     nBins[2] = 50;
1461     nBins[nVars-1] = 6;
1462     cf = new AliCFContainer(name, title, 1, nVars, nBins);
1463     cf->SetBinLimits(0, trkPhiLims[0], trkPhiLims[1]);
1464     cf->SetBinLimits(1, trkEtaLims[0], trkEtaLims[1]);
1465     cf->SetBinLimits(2, trkQtotLims[0], trkQtotLims[1]);
1466     cf->SetBinLimits(3, evTriggerLims[0], evTriggerLims[1]);
1467     cf->SetBinLimits(4, -0.5, 5.5);
1468     for(Int_t i=0; i<nVars-1; ++i) cf->SetVarTitle(i, varNames[fQtotCFVars[i]]);
1469     cf->SetVarTitle(nVars-1, "layer");
1470     cf->SetStepTitle(0, "TRD");
1471     return cf;
1472   }
1473   if(nameStr.Contains("PulseHeight")) {
1474     fPulseHeightCFVars[0] = kTrackP; fPulseHeightCFVars[1] = kTrackPHslice; fPulseHeightCFVars[2] = kEventTrigger;
1475     const Int_t nVars = 4;
1476     Int_t nBins[nVars]; for(Int_t i=0; i<nVars-1; ++i) nBins[i] = nbinsCf[fPulseHeightCFVars[i]];
1477     nBins[nVars-1] = 8;
1478     cf = new AliCFContainer(name, title, 1, nVars, nBins);
1479     //cf->SetBinLimits(0, evTriggerLims[0], evTriggerLims[1]);
1480     cf->SetBinLimits(0, trkPLims);
1481     cf->SetBinLimits(1, trkQtotLims[0], trkQtotLims[1]);
1482     cf->SetBinLimits(2, evTriggerLims[0], evTriggerLims[1]);
1483     cf->SetBinLimits(3, -0.5, 7.5);
1484     for(Int_t i=0; i<nVars-1; ++i) cf->SetVarTitle(i, varNames[fPulseHeightCFVars[i]]);
1485     cf->SetVarTitle(nVars-1, "slice");
1486     cf->SetStepTitle(0, "TRD");
1487     return cf;
1488   }
1489   if(nameStr.Contains("Expert")) {
1490     Int_t nVars = 0;
1491     Int_t nBins[kNTrdCfVariables];
1492     for(Int_t ivar=0; ivar<kNTrdCfVariables; ++ivar) {
1493       if(!fExpertCFVarsEnabled[ivar]) continue;
1494       if(fExpertCFVarBins[ivar][0]=='\0') {
1495         nBins[nVars] = fExpertCFVarNBins[ivar];
1496         nVars++;
1497       }
1498       else {
1499         TObjArray* arr = fExpertCFVarBins[ivar].Tokenize(";");
1500         nBins[nVars] = arr->GetEntries()-1;
1501         if(nBins[nVars]>0) nVars++;
1502       }
1503     }
1504     if(nVars<1) return 0x0;
1505     Int_t nSteps = 0; for(Int_t i=0; i<3; ++i) if(fExpertCFEnabledSteps[i]) nSteps++;
1506     if(nSteps<1) return 0x0;
1507     cf = new AliCFContainer(name, title, nSteps, nVars, nBins);
1508     Int_t iUsedVar = 0;
1509     for(Int_t ivar=0; ivar<kNTrdCfVariables; ++ivar) {
1510       if(!fExpertCFVarsEnabled[ivar]) continue;
1511       if(fExpertCFVarBins[ivar][0]=='\0')
1512         cf->SetBinLimits(iUsedVar, fExpertCFVarRanges[ivar][0], fExpertCFVarRanges[ivar][1]);
1513       else {
1514         TObjArray* arr = fExpertCFVarBins[ivar].Tokenize(";");
1515         if(arr->GetEntries()-1>0) {
1516           Double_t* binLims = new Double_t[arr->GetEntries()];
1517           for(Int_t ib=0;ib<arr->GetEntries();++ib) {
1518             TString binStr = arr->At(ib)->GetName();
1519             binLims[ib] = binStr.Atof();
1520           }
1521           cf->SetBinLimits(iUsedVar++, binLims);
1522         }
1523       }
1524       cf->SetVarTitle(iUsedVar, varNames[ivar]);
1525     }
1526     const Char_t* stepNames[3] = {"TPC","TRD","TOF"};
1527     Int_t iUsedStep = 0;
1528     for(Int_t istep=0; istep<3; ++istep) {
1529       if(fExpertCFEnabledSteps[istep]) cf->SetStepTitle(iUsedStep++, stepNames[istep]);
1530     }
1531     return cf;
1532   }  
1533   return 0x0;
1534 }
1535
1536
1537 //____________________________________________________________________
1538 Bool_t AliTRDcheckESD::Load(const Char_t *file, const Char_t *dir, const Char_t *name)
1539 {
1540 // Load data from performance file
1541
1542   if(!TFile::Open(file)){
1543     AliWarning(Form("Couldn't open file %s.", file));
1544     return kFALSE;
1545   }
1546   if(dir){
1547     if(!gFile->cd(dir)){
1548       AliWarning(Form("Couldn't cd to %s in %s.", dir, file));
1549       return kFALSE;
1550     }
1551   }
1552   TObjArray *o(NULL);
1553   const Char_t *tn=(name ? name : GetName());
1554   if(!(o = (TObjArray*)gDirectory->Get(tn))){
1555     AliWarning(Form("Missing histogram container %s.", tn));
1556     return kFALSE;
1557   }
1558   fHistos = (TObjArray*)o->Clone(GetName());
1559   fMatchingPhiEtaCF = (AliCFContainer*)fHistos->At(kMatchingPhiEtaCF);
1560   fMatchingPtCF = (AliCFContainer*)fHistos->At(kMatchingPtCF);
1561   fBunchCrossingsCF = (AliCFContainer*)fHistos->At(kBunchCrossingsCF);
1562   fCentralityCF = (AliCFContainer*)fHistos->At(kCentralityCF);
1563   fQtotCF = (AliCFContainer*)fHistos->At(kQtotCF);
1564   fPulseHeightCF = (AliCFContainer*)fHistos->At(kPulseHeightCF);
1565   fExpertCF = (AliCFContainer*)fHistos->At(kExpertCF);
1566   
1567   /*
1568   TObjArray *cfs(NULL);
1569   if(!(cfs = (TObjArray*)gDirectory->Get(Form("%s_CF", tn)))){
1570     AliWarning(Form("Missing CFs container %s_CF.", tn));
1571     fCfList = NULL;
1572     //return kFALSE;
1573   } 
1574   else
1575     fCfList = (TObjArray*)cfs->Clone(Form("%s_CF_clone", GetName()));
1576   */
1577   
1578   gFile->Close();
1579   return kTRUE;
1580 }
1581
1582 //_______________________________________________________
1583 Bool_t AliTRDcheckESD::PutTrendValue(const Char_t *name, Double_t val)
1584 {
1585 // Dump trending value to default file
1586
1587   if(!fgFile){
1588     fgFile = fopen("TRD.Performance.txt", "at");
1589   }
1590   fprintf(fgFile, "%s_%s %f\n", GetName(), name, val);
1591   return kTRUE;
1592 }
1593
1594 //____________________________________________________________________
1595 void AliTRDcheckESD::Terminate(Option_t *)
1596 {
1597   // Steer post-processing 
1598   if(!fHistos){
1599     fHistos = dynamic_cast<TObjArray *>(GetOutputData(1));
1600     if(!fHistos){
1601       AliError("Histogram container not found in output");
1602       return;
1603     }
1604   }
1605
1606   //  fNRefFigures = 15;
1607   //  return;
1608
1609   const Char_t *name[kNrefs] = {
1610     "Ncl", "Eff", "Eloss", "PtResDCA"
1611   };
1612
1613   TObjArray *arr(NULL); TGraph *g(NULL);
1614   if(!fResults){
1615     fResults = new TObjArray(kNrefs);
1616     fResults->SetOwner();
1617     fResults->SetName("results");
1618     for(Int_t iref(0); iref<kNrefs; iref++){
1619       fResults->AddAt(arr = new TObjArray(fgkNgraph[iref]), iref);
1620       arr->SetName(name[iref]);  arr->SetOwner();
1621       switch(iref+1){
1622       case kNCl:
1623         for(Int_t ig(0); ig<fgkNgraph[iref]; ig++){
1624           arr->AddAt(g = new TGraphErrors(), ig);
1625           g->SetLineColor(ig+1); 
1626           g->SetMarkerColor(ig+1); 
1627           g->SetMarkerStyle(ig+20); 
1628           g->SetName(Form("s%d", ig));
1629           switch(ig){
1630           case 0: g->SetTitle("ALL"); break;
1631           case 1: g->SetTitle("NEG"); break;
1632           case 2: g->SetTitle("POS"); break;
1633           default: g->SetTitle(AliPID::ParticleLatexName(ig-3)); break;
1634           };
1635         }
1636         break;
1637       case kTRDmom:
1638         for(Int_t ig(0); ig<fgkNgraph[iref]; ig++){
1639           arr->AddAt(g = new TGraphAsymmErrors(), ig);
1640           g->SetLineColor(ig+1); 
1641           g->SetMarkerColor(ig+1); 
1642           g->SetMarkerStyle(ig+20); 
1643         }
1644         break;
1645       case kPtRes:
1646         for(Int_t idx(0); idx<AliPID::kSPECIES; idx++){
1647           Int_t ig(2*idx);
1648           arr->AddAt(g = new TGraphErrors(), ig);
1649           g->SetLineColor(kRed-idx); 
1650           g->SetMarkerColor(kRed-idx); 
1651           g->SetMarkerStyle(20+idx); 
1652           g->SetNameTitle(Form("s%d", ig), Form("res %s", AliPID::ParticleLatexName(idx)));
1653           arr->AddAt(g = new TGraphErrors(), ig+1);
1654           g->SetLineColor(kBlue-idx); 
1655           g->SetMarkerColor(kBlue-idx); 
1656           g->SetMarkerStyle(20+idx); 
1657           g->SetNameTitle(Form("m%d", ig+1), Form("sys %s", AliPID::ParticleLatexName(idx)));
1658
1659           ig+=10;
1660           arr->AddAt(g = new TGraphErrors(), ig);
1661           g->SetLineColor(kRed-idx); 
1662           g->SetMarkerColor(kRed-idx); 
1663           g->SetMarkerStyle(20+idx); 
1664           g->SetNameTitle(Form("s%d", ig), Form("sigma %s", AliPID::ParticleLatexName(idx)));
1665           arr->AddAt(g = new TGraphErrors(), ig+1);
1666           g->SetLineColor(kBlue-idx); 
1667           g->SetMarkerColor(kBlue-idx); 
1668           g->SetMarkerStyle(20+idx); 
1669           g->SetNameTitle(Form("m%d", ig+1), Form("mean %s", AliPID::ParticleLatexName(idx)));
1670         }
1671         break;
1672       default:
1673         for(Int_t ig(0); ig<fgkNgraph[iref]; ig++){
1674           arr->AddAt(g = new TGraphErrors(), ig);
1675           g->SetLineColor(ig+1); 
1676           g->SetMarkerColor(ig+1); 
1677           g->SetMarkerStyle(ig+20); 
1678         }
1679         break;
1680       }
1681     }
1682   }
1683   TH1 *h1[2] = {NULL, NULL};
1684   TH2I *h2(NULL);
1685   TAxis *ax(NULL);
1686
1687   // No of clusters
1688   if(!(h2 = (TH2I*)fHistos->At(kNCl))) return;
1689   ax = h2->GetXaxis();
1690   arr = (TObjArray*)fResults->At(kNCl);
1691   
1692   // All tracks
1693   h1[0] = h2->ProjectionX("Ncl_px");
1694   TGraphErrors *ge=(TGraphErrors*)arr->At(0);
1695   for(Int_t ib=2; ib<=ax->GetNbins(); ib++){
1696     ge->SetPoint(ib-2, ax->GetBinCenter(ib), h1[0]->GetBinContent(ib));
1697   }
1698   
1699   // All charged tracks
1700   TH1 *hNclCh[2] = {(TH1D*)h1[0]->Clone("NEG"), (TH1D*)h1[0]->Clone("POS")};
1701   hNclCh[0]->Reset();hNclCh[1]->Reset();
1702   for(Int_t is(1); is<=AliPID::kSPECIES; is++){
1703     hNclCh[0]->Add(h2->ProjectionX("Ncl_px", 2*is-1, 2*is-1)); // neg
1704     hNclCh[1]->Add(h2->ProjectionX("Ncl_px", 2*is, 2*is));     // pos
1705   }
1706   if(Int_t(hNclCh[0]->GetEntries())){
1707     ge=(TGraphErrors*)arr->At(1);
1708     for(Int_t ib=2; ib<=ax->GetNbins(); ib++){
1709       ge->SetPoint(ib-2, ax->GetBinCenter(ib), hNclCh[0]->GetBinContent(ib));
1710     }
1711   }
1712   
1713   if(Int_t(hNclCh[1]->GetEntries())){
1714     ge=(TGraphErrors*)arr->At(2);
1715     for(Int_t ib=2; ib<=ax->GetNbins(); ib++){
1716       ge->SetPoint(ib-2, ax->GetBinCenter(ib), hNclCh[1]->GetBinContent(ib));
1717     }
1718   }
1719   // Species wise
1720   for(Int_t is(1); is<=AliPID::kSPECIES; is++){
1721     h1[0] = h2->ProjectionX("Ncl_px", 2*is-1, 2*is);
1722     if(!Int_t(h1[0]->GetEntries())) continue;
1723     ge=(TGraphErrors*)arr->At(2+is);
1724     for(Int_t ib=2; ib<=ax->GetNbins(); ib++){
1725       ge->SetPoint(ib-2, ax->GetBinCenter(ib), h1[0]->GetBinContent(ib));
1726     }
1727   }
1728   
1729   fNRefFigures = 1;
1730
1731   // EFFICIENCY
1732   // geometrical efficiency
1733   if(!(h2 = (TH2I*)fHistos->At(kTRDstat))) return;
1734   arr = (TObjArray*)fResults->At(kTRDstat-1);
1735   h1[0] = h2->ProjectionX("checkESDx0", kTPCout, kTPCout);
1736   h1[1] = h2->ProjectionX("checkESDx1", kTRDin, kTRDin);
1737   Process(h1, (TGraphErrors*)arr->At(0));
1738   delete h1[0];delete h1[1];
1739   // tracking efficiency
1740   h1[0] = h2->ProjectionX("checkESDx0", kTRDin, kTRDin);
1741   h1[1] = h2->ProjectionX("checkESDx1", kTRDout, kTRDout);
1742   Process(h1, (TGraphErrors*)arr->At(1));
1743   delete h1[1];
1744   // PID efficiency
1745   h1[1] = h2->ProjectionX("checkESDx1", kTRDpid, kTRDpid);
1746   Process(h1, (TGraphErrors*)arr->At(2));
1747   delete h1[1];
1748   // Refit efficiency
1749   h1[1] = h2->ProjectionX("checkESDx1", kTRDref, kTRDref);
1750   Process(h1, (TGraphErrors*)arr->At(3));
1751   delete h1[1];
1752   fNRefFigures++;
1753
1754   // ENERGY LOSS
1755   if(!(h2 = dynamic_cast<TH2I*>(fHistos->At(kTRDmom)))) return;
1756   arr = (TObjArray*)fResults->At(kTRDmom-1);
1757   TGraphAsymmErrors *g06 = (TGraphAsymmErrors*)arr->At(0), *g09 = (TGraphAsymmErrors*)arr->At(1);
1758   ax=h2->GetXaxis();
1759   const Int_t nq(4);
1760   const Double_t xq[nq] = {0.05, 0.2, 0.8, 0.95};
1761   Double_t yq[nq];
1762   for(Int_t ily=6; ily--;){
1763     h1[0] = h2->ProjectionX("checkESDp0", ily+1, ily+1);
1764     h1[0]->GetQuantiles(nq,yq,xq);
1765     g06->SetPoint(ily, Float_t(ily), ax->GetBinCenter(h1[0]->GetMaximumBin()));
1766     g06->SetPointError(ily, 0., 0., TMath::Abs(yq[0]), yq[3]);
1767     g09->SetPoint(ily, Float_t(ily), h1[0]->GetMean());
1768     g09->SetPointError(ily, 0., 0., TMath::Abs(yq[1]), yq[2]);
1769
1770     //printf(" max[%f] mean[%f] q[%f %f %f %f]\n", ax->GetBinCenter(h1[0]->GetMaximumBin()), h1[0]->GetMean(), yq[0], yq[1], yq[2], yq[3]);
1771     delete h1[0];
1772   }
1773   fNRefFigures++;
1774 //  if(!HasMC()) return;
1775
1776   // Pt RESOLUTION @ DCA
1777   TH3S* h3(NULL); TGraphErrors *gg[2] = {NULL,NULL};
1778   if(!(h3 = dynamic_cast<TH3S*>(fHistos->At(kPtRes)))) return;
1779   arr = (TObjArray*)fResults->At(kPtRes-1);
1780   TAxis *az(h3->GetZaxis());
1781   for(Int_t i(0); i<AliPID::kSPECIES; i++){
1782     Int_t idx(2*i);
1783     az->SetRange(idx+1, idx+2); 
1784     gg[1] = (TGraphErrors*)arr->At(idx);
1785     gg[0] = (TGraphErrors*)arr->At(idx+1);
1786     Process2D((TH2*)h3->Project3D("yx"), gg);
1787
1788     idx+=10;
1789     az->SetRange(idx+1, idx+2); 
1790     gg[1] = (TGraphErrors*)arr->At(idx);
1791     gg[0] = (TGraphErrors*)arr->At(idx+1);
1792     Process2D((TH2*)h3->Project3D("yx"), gg);
1793   }
1794   fNRefFigures++;
1795   
1796   fNRefFigures++;
1797   // 3x3 tracking summary canvases for every centrality class
1798   fNRefFigures++;
1799   // 3x3 PID summary canvases for every centrality class
1800   fNRefFigures++;
1801   // 3x3 for centrality dependent pictures
1802   fNRefFigures++;
1803   
1804   //DoTrending();  
1805 }
1806
1807 //____________________________________________________________________
1808 Int_t AliTRDcheckESD::Pdg2Idx(Int_t pdg) const
1809 {
1810   //
1811   // Helper function converting PDG code into AliPID index
1812   //
1813   switch(pdg){
1814   case kElectron: 
1815   case kPositron: return AliPID::kElectron;  
1816   case kMuonPlus:
1817   case kMuonMinus: return AliPID::kMuon;  
1818   case kPiPlus: 
1819   case kPiMinus: return AliPID::kPion;  
1820   case kKPlus: 
1821   case kKMinus: return AliPID::kKaon;
1822   case kProton: 
1823   case kProtonBar: return AliPID::kProton;
1824   } 
1825   return -1;
1826 }
1827
1828 //____________________________________________________________________
1829 void AliTRDcheckESD::Process(TH1 **h1, TGraphErrors *g)
1830 {
1831 // Generic function to process one reference plot
1832
1833   Int_t n1 = 0, n2 = 0, ip=0;
1834   Double_t eff = 0.;
1835
1836   TAxis *ax = h1[0]->GetXaxis();
1837   for(Int_t ib=1; ib<=ax->GetNbins(); ib++){
1838     if(!(n1 = (Int_t)h1[0]->GetBinContent(ib))) continue;
1839     n2 = (Int_t)h1[1]->GetBinContent(ib);
1840     eff = n2/Float_t(n1);
1841
1842     ip=g->GetN();
1843     g->SetPoint(ip, ax->GetBinCenter(ib), eff);
1844     g->SetPointError(ip, 0., n2 ? eff*TMath::Sqrt(1./n1+1./n2) : 0.);
1845   }
1846 }  
1847 //________________________________________________________
1848 void AliTRDcheckESD::Process2D(TH2 * const h2, TGraphErrors **g)
1849 {
1850   //
1851   // Do the processing
1852   //
1853
1854   Int_t n = 0;
1855   if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);
1856   if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);
1857   TF1 f("fg", "gaus", -3.,3.);
1858   for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){
1859     Double_t x = h2->GetXaxis()->GetBinCenter(ibin);
1860     TH1D *h = h2->ProjectionY("py", ibin, ibin);
1861     if(h->GetEntries()<100) continue;
1862     //AdjustF1(h, f);
1863
1864     h->Fit(&f, "QN");
1865     Int_t ip = g[0]->GetN();
1866     g[0]->SetPoint(ip, x, f.GetParameter(1));
1867     g[0]->SetPointError(ip, 0., f.GetParError(1));
1868     g[1]->SetPoint(ip, x, f.GetParameter(2));
1869     g[1]->SetPointError(ip, 0., f.GetParError(2));
1870   }
1871   return;
1872 }
1873 //____________________________________________________________________
1874 void AliTRDcheckESD::PrintStatus(ULong_t status)
1875 {
1876 // Dump track status to stdout
1877
1878   printf("ITS[i(%d) o(%d) r(%d)] TPC[i(%d) o(%d) r(%d) p(%d)] TRD[i(%d) o(%d) r(%d) p(%d) s(%d)] HMPID[o(%d) p(%d)]\n"
1879     ,Bool_t(status & AliESDtrack::kITSin)
1880     ,Bool_t(status & AliESDtrack::kITSout)
1881     ,Bool_t(status & AliESDtrack::kITSrefit)
1882     ,Bool_t(status & AliESDtrack::kTPCin)
1883     ,Bool_t(status & AliESDtrack::kTPCout)
1884     ,Bool_t(status & AliESDtrack::kTPCrefit)
1885     ,Bool_t(status & AliESDtrack::kTPCpid)
1886     ,Bool_t(status & AliESDtrack::kTRDin)
1887     ,Bool_t(status & AliESDtrack::kTRDout)
1888     ,Bool_t(status & AliESDtrack::kTRDrefit)
1889     ,Bool_t(status & AliESDtrack::kTRDpid)
1890     ,Bool_t(status & AliESDtrack::kTRDStop)
1891     ,Bool_t(status & AliESDtrack::kHMPIDout)
1892     ,Bool_t(status & AliESDtrack::kHMPIDpid)
1893   );
1894 }
1895
1896 //____________________________________________________________________
1897 TH1D* AliTRDcheckESD::Proj2D(TH2* hist, TH1* fitErr) {
1898   //
1899   // project the PH vs Slice 2D-histo into a 1D histo
1900   //
1901   
1902   TH1D* hProjection = (TH1D*)hist->ProjectionX(Form("hProjection_%f", gRandom->Rndm()));
1903   hProjection->Reset();
1904   
1905   TF1* fitLandau = new TF1("landauFunc","landau",0.,2000.);
1906   TH1D *hD;
1907   for(Int_t iBin=1;iBin<=hist->GetXaxis()->GetNbins();iBin++) {
1908     if(gROOT->FindObject("projection"))
1909       delete gROOT->FindObject("projection");
1910     hD = (TH1D*)hist->ProjectionY("projection",iBin,iBin);
1911     hD->Rebin(4);
1912     if(hD->Integral()>10) {
1913       fitLandau->SetParameter(1, hD->GetBinCenter(hD->GetMaximumBin()));
1914       fitLandau->SetParLimits(1, 0.2*hD->GetBinCenter(hD->GetMaximumBin()), 3.0*hD->GetBinCenter(hD->GetMaximumBin()));
1915       fitLandau->SetParameter(0, 1000.);
1916       fitLandau->SetParLimits(0, 1., 10000000.);
1917       fitLandau->SetParameter(2, 0.5*hD->GetBinCenter(hD->GetMaximumBin()));
1918       fitLandau->SetParLimits(2, 0.01*hD->GetBinCenter(hD->GetMaximumBin()), 1.0*hD->GetBinCenter(hD->GetMaximumBin()));
1919       hD->Fit(fitLandau, "Q0", "", hD->GetXaxis()->GetXmin(), hD->GetXaxis()->GetXmax());
1920       hD->Fit(fitLandau, "Q0", "", hD->GetXaxis()->GetXmin(), hD->GetXaxis()->GetXmax());
1921       hProjection->SetBinContent(iBin, fitLandau->GetParameter(1));
1922       hProjection->SetBinError(iBin, fitLandau->GetParameter(2));
1923       if(fitErr) {
1924         fitErr->SetBinContent(iBin, fitLandau->GetParameter(1));
1925         fitErr->SetBinError(iBin, fitLandau->GetParError(1));
1926       }
1927     }
1928     else{
1929       hProjection->SetBinContent(iBin, 0);
1930       hProjection->SetBinError(iBin, 0);
1931     }
1932   }
1933   return hProjection;
1934 }
1935
1936 //____________________________________________________________________
1937 TH2F* AliTRDcheckESD::Proj3D(TH3* hist, TH2* accMap, Int_t zbinLow, Int_t zbinHigh, Float_t &entries) {
1938   //
1939   //  Project a 3D histogram to a 2D histogram in the Z axis interval [zbinLow,zbinHigh] 
1940   //  Return the 2D histogram and also the number of entries into this projection (entries)
1941
1942   Int_t nBinsX = hist->GetXaxis()->GetNbins();   // X and Y axis bins are assumed to be all equal
1943   Float_t minX = hist->GetXaxis()->GetXmin();
1944   Float_t maxX = hist->GetXaxis()->GetXmax();
1945   Int_t nBinsY = hist->GetYaxis()->GetNbins();
1946   Float_t minY = hist->GetYaxis()->GetXmin();
1947   Float_t maxY = hist->GetYaxis()->GetXmax();
1948   Int_t nBinsZ = hist->GetZaxis()->GetNbins();  // Z axis bins (pt) might have different widths
1949
1950   TH2F* projHisto = (TH2F*)gROOT->FindObject("projHisto");
1951   if(projHisto) 
1952     projHisto->Reset();
1953   else
1954     projHisto = new TH2F("projHisto", "projection", nBinsX, minX, maxX, nBinsY, minY, maxY);
1955
1956   for(Int_t iZ=1; iZ<=nBinsZ; iZ++) {
1957     if(iZ<zbinLow) continue;
1958     if(iZ>zbinHigh) continue;
1959     for(Int_t iX=1; iX<=nBinsX; iX++) {
1960       for(Int_t iY=1; iY<=nBinsY; iY++) {
1961         if(accMap) {
1962           if(accMap->GetBinContent(iX,iY)>0.1)
1963             projHisto->SetBinContent(iX, iY, projHisto->GetBinContent(iX, iY)+hist->GetBinContent(iX,iY,iZ));
1964         }
1965         else    // no acc. cut 
1966           projHisto->SetBinContent(iX, iY, projHisto->GetBinContent(iX, iY)+hist->GetBinContent(iX,iY,iZ));
1967         // count only the entries which are inside the acceptance map
1968         if(accMap) {
1969           if(accMap->GetBinContent(iX,iY)>0.1)
1970             entries+=hist->GetBinContent(iX,iY,iZ);
1971         }
1972         else    // no acc. cut
1973           entries+=hist->GetBinContent(iX,iY,iZ);
1974       }
1975     }
1976   }
1977   return projHisto;
1978 }
1979
1980 //____________________________________________________________________
1981 void AliTRDcheckESD::CheckActiveSM(TH1D* phiProj, Bool_t activeSM[18]) {
1982   //
1983   // Check the active super-modules
1984   //
1985   Double_t entries[18] = {0.0};
1986   Double_t smPhiLimits[19];
1987   for(Int_t ism=0; ism<=18; ++ism) smPhiLimits[ism] = -TMath::Pi() + (2.0*TMath::Pi()/18.0)*ism;
1988   for(Int_t phiBin=1; phiBin<=phiProj->GetXaxis()->GetNbins(); ++phiBin) {
1989     Double_t phi = phiProj->GetBinCenter(phiBin);
1990     Int_t sm = -1;
1991     for(Int_t ism=0; ism<18; ++ism) 
1992       if(phi>=smPhiLimits[ism] && phi<smPhiLimits[ism+1]) sm = ism;
1993     if(sm==-1) continue;
1994     entries[sm] += phiProj->GetBinContent(phiBin);
1995   }
1996   Double_t avEntries = Double_t(phiProj->Integral())/18.0;
1997   for(Int_t ism=0; ism<18; ++ism) 
1998     if(entries[ism]>0.5*avEntries) activeSM[ism] = kTRUE;
1999 }
2000
2001
2002 //__________________________________________________________________________________________________
2003 TH1F* AliTRDcheckESD::EfficiencyFromPhiPt(AliCFContainer* cf, Int_t minNtrkl, Int_t maxNtrkl, 
2004                                           Int_t stepNom, Int_t stepDenom, const Char_t* varStr/*="pt"*/, const Char_t* type/*="TPCTRD"*/) {
2005   //
2006   // Use the CF container to extract the efficiency vs pt
2007   //
2008   Int_t varTrackPhi = cf->GetVar("phi");
2009   Int_t var = cf->GetVar(varStr);
2010     
2011   TH1D* phiProj = (TH1D*)cf->Project(1, varTrackPhi);
2012   Bool_t activeSM[18] = {kFALSE};
2013   CheckActiveSM(phiProj, activeSM); delete phiProj;
2014   Double_t smPhiLimits[19];
2015   for(Int_t ism=0; ism<=18; ++ism) smPhiLimits[ism] = -TMath::Pi() + (2.0*TMath::Pi()/18.0)*ism;
2016   
2017   TString effTypeStr = type;
2018   if(effTypeStr.Contains("TRDTOF")) {
2019     if(minNtrkl>-1 && minNtrkl<7 && maxNtrkl>-1 && maxNtrkl<7)
2020       cf->SetRangeUser(cf->GetVar("tracklets"), Double_t(minNtrkl), Double_t(maxNtrkl));
2021   }
2022   TH2D* hDenomPhiVar = (TH2D*)cf->Project(stepDenom, var, varTrackPhi);
2023   if(effTypeStr.Contains("TPCTRD")) {
2024     if(minNtrkl>-1 && minNtrkl<7 && maxNtrkl>-1 && maxNtrkl<7)
2025       cf->SetRangeUser(cf->GetVar("tracklets"), Double_t(minNtrkl), Double_t(maxNtrkl));
2026   }
2027   TH2D* hNomPhiVar = (TH2D*)cf->Project(stepNom, var, varTrackPhi);
2028   cf->SetRangeUser(cf->GetVar("tracklets"), 0.0, 6.0);
2029   
2030   TH1F* hEff = new TH1F(Form("hEff%s_%d_%d_%f", varStr, stepNom, stepDenom, gRandom->Rndm()), "", 
2031                         hNomPhiVar->GetXaxis()->GetNbins(), hNomPhiVar->GetXaxis()->GetXbins()->GetArray());
2032   for(Int_t ivar=1; ivar<=hEff->GetXaxis()->GetNbins(); ++ivar) {
2033     Double_t nom = 0.0; Double_t denom = 0.0;
2034     Double_t eff = 0.0; Double_t err = 0.0;
2035     for(Int_t iphi=1; iphi<=hNomPhiVar->GetYaxis()->GetNbins(); ++iphi) {
2036       Double_t phi = hNomPhiVar->GetYaxis()->GetBinCenter(iphi);
2037       Bool_t isActive = kFALSE;
2038       for(Int_t ism=0; ism<18; ++ism) 
2039         if(phi>=smPhiLimits[ism] && phi<smPhiLimits[ism+1] && activeSM[ism]) 
2040           isActive = kTRUE;
2041       if(!isActive) continue;
2042       nom += hNomPhiVar->GetBinContent(ivar, iphi);
2043       denom += hDenomPhiVar->GetBinContent(ivar, iphi);
2044     }
2045     eff = (denom>0.001 ? nom/denom : 0.0);
2046     err = (denom>0.001 && (denom-nom)>0.001 && nom>0.001 ? (TMath::Sqrt(nom*(denom-nom)/denom/denom/denom)) : 0.0);
2047     hEff->SetBinContent(ivar, eff);
2048     hEff->SetBinError(ivar, err);
2049   }   // end loop over pt bins
2050   delete hNomPhiVar; delete hDenomPhiVar;
2051   return hEff;
2052 }
2053
2054
2055 //____________________________________________________________________
2056 TH1F* AliTRDcheckESD::EfficiencyTRD(TH3* tpc3D, TH3* trd3D, Bool_t useAcceptance) {
2057   //
2058   // Calculate the TRD-TPC matching efficiency as function of pt
2059   //
2060   
2061   if(!tpc3D || !trd3D) return NULL;
2062   Int_t nBinsZ = trd3D->GetZaxis()->GetNbins();
2063   // project everything on the eta-phi map to obtain an acceptance map
2064   Float_t nada = 0.;
2065   TH2F *trdAcc = (useAcceptance ? (TH2F*)Proj3D(trd3D, 0x0, 1, nBinsZ, nada)->Clone(Form("trdAcc%f", gRandom->Rndm())) : 0x0);
2066   TH1D *phiProj = (trdAcc ? trdAcc->ProjectionY(Form("phiProj%f", gRandom->Rndm())) : 0x0);
2067   
2068   // prepare the acceptance map
2069   Bool_t activeSM[18] = {kFALSE};
2070   Double_t smPhiLimits[19];
2071   for(Int_t ism=0; ism<=18; ++ism) smPhiLimits[ism] = -TMath::Pi() + (2.0*TMath::Pi()/18.0)*ism;
2072   if(phiProj) {
2073     CheckActiveSM(phiProj, activeSM);   // get the active SMs
2074     trdAcc->Reset();
2075     // Put 1 entry in every bin which belongs to an active SM
2076     for(Int_t iY=1; iY<=trdAcc->GetYaxis()->GetNbins(); ++iY) {
2077       Double_t phi = trdAcc->GetYaxis()->GetBinCenter(iY);
2078       Bool_t isActive = kFALSE;
2079       for(Int_t ism=0; ism<18; ++ism) {
2080         if(phi>=smPhiLimits[ism] && phi<smPhiLimits[ism+1] && activeSM[ism]) {
2081           isActive = kTRUE;
2082         }
2083       }
2084       if(!isActive) continue;
2085       for(Int_t iX=1; iX<=trdAcc->GetXaxis()->GetNbins(); ++iX) 
2086         if(trdAcc->GetXaxis()->GetBinCenter(iX)>=-0.85 && trdAcc->GetXaxis()->GetBinCenter(iX)<=0.85) trdAcc->SetBinContent(iX, iY, 1.0);
2087     }  // end for over Y(phi) bins
2088   }  // end if phiProj
2089     
2090   // get the bin limits from the Z axis of 3D histos
2091   Float_t *ptBinLimits = new Float_t[nBinsZ+1];
2092   for(Int_t i=1; i<=nBinsZ; i++) {
2093     ptBinLimits[i-1] = trd3D->GetZaxis()->GetBinLowEdge(i);
2094   }
2095   ptBinLimits[nBinsZ] = trd3D->GetZaxis()->GetBinUpEdge(nBinsZ);
2096   
2097   TH1F *efficiency = new TH1F(Form("eff%d", Int_t(1000000.0*gRandom->Rndm())), "TRD-TPC matching efficiency", nBinsZ, ptBinLimits);
2098   
2099   // loop over Z bins
2100   Bool_t effGood = kFALSE;
2101   for(Int_t i=1; i<=nBinsZ; i++) {
2102     Float_t tpcEntries = 0.0; Float_t trdEntries = 0.0;
2103     Proj3D(tpc3D, trdAcc, i, i, tpcEntries);
2104     Proj3D(trd3D, trdAcc, i, i, trdEntries);
2105     Float_t ratio = 0;
2106     if(tpcEntries>0) ratio = trdEntries/tpcEntries;
2107     Float_t error = 0;
2108     if(tpcEntries>0 && trdEntries>0 && (tpcEntries-trdEntries)>=0.0) 
2109       error = TMath::Sqrt(trdEntries*(tpcEntries-trdEntries)/tpcEntries/tpcEntries/tpcEntries);
2110     if(ratio>0.001) {
2111       efficiency->SetBinContent(i,ratio);
2112       efficiency->SetBinError(i,error);
2113       effGood = kTRUE;
2114     }
2115   }     // end loop over Z bins
2116   if(!effGood) return 0x0;
2117   
2118   return efficiency;
2119 }
2120
2121 //__________________________________________________________________________________________________
2122 void AliTRDcheckESD::PlotCentSummaryFromCF(Double_t* /*trendValues*/, const Char_t* /*triggerName*/, Bool_t /*useIsolatedBC*/, Bool_t /*cutTOFbc*/) {
2123   //
2124   // Make the centrality summary figure from the CF container 
2125   // 
2126   if(!fMatchingPtCF) return;
2127   AliCFContainer* cf = 0x0;    
2128   
2129   TLatex* lat=new TLatex();
2130   lat->SetTextSize(0.06);
2131   lat->SetTextColor(2);
2132
2133   gPad->SetTopMargin(0.05); gPad->SetBottomMargin(0.001); gPad->SetLeftMargin(0.001); gPad->SetRightMargin(0.001);
2134   gPad->Divide(3,3,0.,0.);
2135   TList* l=gPad->GetListOfPrimitives();
2136   TVirtualPad* pad=0x0;
2137   
2138   //if(cutTOFbc) cf->SetRangeUser(stepTOFBC, 0.0, 0.0);
2139   
2140   if(gROOT->FindObject("rangeEffPt")) delete gROOT->FindObject("rangeEffPt");
2141   TH2F* rangeEffPt=new TH2F("rangeEffPt", "",10,0.,10.,10,0.,1.3);
2142   rangeEffPt->SetStats(kFALSE);
2143   SetStyle(rangeEffPt->GetXaxis(), "p_{T} [GeV/c]", 0.07, 0.8, kTRUE, 0.05);
2144   SetStyle(rangeEffPt->GetYaxis(), "efficiency", 0.07, 0.8, kTRUE, 0.05);
2145   
2146   Int_t padsForEffs[5] = {0,3,6,1,4};
2147   for(Int_t iCent=1; iCent<6; ++iCent) {
2148     pad = ((TVirtualPad*)l->At(padsForEffs[iCent-1])); pad->cd();
2149     pad->SetLeftMargin(0.15); pad->SetRightMargin(0.02); pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);
2150     pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);
2151     rangeEffPt->Draw();
2152     TLine line;
2153     line.SetLineStyle(2);
2154     line.SetLineWidth(2);
2155     line.DrawLine(rangeEffPt->GetXaxis()->GetXmin(), 0.7, rangeEffPt->GetXaxis()->GetXmax(), 0.7);
2156     line.DrawLine(rangeEffPt->GetXaxis()->GetXmin(), 0.9, rangeEffPt->GetXaxis()->GetXmax(), 0.9);
2157     
2158     cf = fMatchingPtCF;
2159     cf->SetRangeUser(cf->GetVar("multiplicity"), Double_t(iCent), Double_t(iCent), kTRUE);
2160        
2161     cf->SetRangeUser(cf->GetVar("charge"), +1.0, +1.0);  
2162     TH1F* hEffPosAll = EfficiencyFromPhiPt(cf,  0,  6, 1, 0);
2163     TH1F* hEffPosTrk4 = EfficiencyFromPhiPt(cf, 4,  4, 1, 0);  
2164     TH1F* hEffPosTrk5 = EfficiencyFromPhiPt(cf, 5,  5, 1, 0);
2165     TH1F* hEffPosTrk6 = EfficiencyFromPhiPt(cf, 6,  6, 1, 0);
2166      
2167     cf->SetRangeUser(cf->GetVar("charge"), -1.0, -1.0);  
2168     TH1F* hEffNegAll = EfficiencyFromPhiPt(cf, 0, 6, 1, 0);  
2169     TH1F* hEffNegTrk4 = EfficiencyFromPhiPt(cf, 4, 4, 1, 0);  
2170     TH1F* hEffNegTrk5 = EfficiencyFromPhiPt(cf, 5, 5, 1, 0);  
2171     TH1F* hEffNegTrk6 = EfficiencyFromPhiPt(cf, 6, 6, 1, 0);
2172     cf->SetRangeUser(cf->GetVar("tracklets"), 0.0, 6.0);  
2173     cf->SetRangeUser(cf->GetVar("charge"), -1.0, +1.0);  
2174     
2175     SetStyle(hEffPosAll,  1, kRed, 1, 24, kRed, 1);
2176     SetStyle(hEffPosTrk4, 1, kRed, 1, 25, kRed, 1);
2177     SetStyle(hEffPosTrk5, 1, kRed, 1, 26, kRed, 1);
2178     SetStyle(hEffPosTrk6, 1, kRed, 1, 27, kRed, 1);
2179     SetStyle(hEffNegAll,  1, kBlue, 1, 24, kBlue, 1);
2180     SetStyle(hEffNegTrk4, 1, kBlue, 1, 25, kBlue, 1);
2181     SetStyle(hEffNegTrk5, 1, kBlue, 1, 26, kBlue, 1);
2182     SetStyle(hEffNegTrk6, 1, kBlue, 1, 27, kBlue, 1);
2183     hEffPosAll->Draw("same");
2184     hEffNegAll->Draw("same");
2185     hEffPosTrk4->Draw("same");
2186     hEffNegTrk4->Draw("same");
2187     hEffPosTrk5->Draw("same");
2188     hEffNegTrk5->Draw("same");
2189     hEffPosTrk6->Draw("same");
2190     hEffNegTrk6->Draw("same");    
2191         
2192     TLegend* leg=new TLegend(0.18, 0.7, 0.77, 0.89);
2193     if(iCent==1) {
2194       leg->SetFillColor(0);
2195       leg->SetNColumns(2);
2196       leg->SetMargin(0.1);
2197       leg->SetBorderSize(0);
2198       leg->AddEntry(hEffPosAll,  "pos. (#geq 1 tracklet)", "p");
2199       leg->AddEntry(hEffNegAll,  "neg. (#geq 1 tracklet)", "p");
2200       leg->AddEntry(hEffPosTrk4, "pos. (4 tracklets)", "p");
2201       leg->AddEntry(hEffNegTrk4, "neg. (4 tracklets)", "p");
2202       leg->AddEntry(hEffPosTrk5, "pos. (5 tracklets)", "p");
2203       leg->AddEntry(hEffNegTrk5, "neg. (5 tracklets)", "p");
2204       leg->AddEntry(hEffPosTrk6, "pos. (6 tracklets)", "p");     
2205       leg->AddEntry(hEffNegTrk6, "neg. (6 tracklets)", "p");
2206       leg->Draw();
2207     }
2208     lat->DrawLatex(0.2, 1.32, Form("%.0f < SPD tracklets < %.0f", cf->GetAxis(cf->GetVar("multiplicity"),0)->GetBinLowEdge(iCent), cf->GetAxis(cf->GetVar("multiplicity"),0)->GetBinUpEdge(iCent)));
2209   }   // end for loop over multiplicity classes
2210   
2211   // Reset the modified user ranges of the CF container
2212   cf->SetRangeUser(cf->GetVar("multiplicity"), 0., 3500.);
2213      
2214   // Cluster distributions in all multiplicity classes
2215   pad = ((TVirtualPad*)l->At(2)); pad->cd();
2216   pad->SetLeftMargin(0.15); pad->SetRightMargin(0.02);
2217   pad->SetTopMargin(0.02); pad->SetBottomMargin(0.15);
2218   pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);
2219   if(gROOT->FindObject("rangeNcls")) delete gROOT->FindObject("rangeNcls");
2220   TH2F* rangeNcls = new TH2F("rangeNcls", "", 10, 0.0, 199.9, 10, 0.0, 1.199);
2221   SetStyle(rangeNcls->GetXaxis(), "# TRD clusters", 0.07, 0.8, kTRUE, 0.05);
2222   SetStyle(rangeNcls->GetYaxis(), "entries (a.u.)", 0.07, 0.8, kTRUE, 0.05);
2223   rangeNcls->SetStats(kFALSE);
2224   rangeNcls->Draw();
2225     
2226   TH1D* hNcls[6]={0x0};
2227   TLegend* legCls=new TLegend(0.7, 0.75, 0.97, 0.97);
2228   legCls->SetBorderSize(0);
2229   legCls->SetFillColor(0);
2230   legCls->SetMargin(0.15);
2231   cf = fCentralityCF;
2232   for(Int_t iCent=0; iCent<6; ++iCent) {
2233     if(iCent>0)
2234       cf->SetRangeUser(cf->GetVar("multiplicity"), Double_t(iCent), Double_t(iCent), kTRUE);
2235     hNcls[iCent] = (TH1D*)cf->Project(0, cf->GetVar("clusters"));
2236     if(!hNcls[iCent]) continue;
2237     
2238     hNcls[iCent]->SetLineColor(iCent<4 ? iCent+1 : iCent+2);
2239     Double_t maximum = hNcls[iCent]->GetMaximum();
2240     if(maximum>1.0)
2241       hNcls[iCent]->Scale(1.0/maximum);
2242     hNcls[iCent]->SetStats(kFALSE);
2243     hNcls[iCent]->SetTitle("");
2244     hNcls[iCent]->SetLineWidth(2);
2245     
2246     if(hNcls[iCent]->Integral()>0.01) {
2247       hNcls[iCent]->Draw("same");
2248       legCls->AddEntry(hNcls[iCent], (iCent==0 ? "all centralities" : Form("%.0f < SPD tracklets < %.0f", cf->GetAxis(cf->GetVar("multiplicity"),0)->GetBinLowEdge(iCent), 
2249                                                                            cf->GetAxis(cf->GetVar("multiplicity"),0)->GetBinUpEdge(iCent))), "l");
2250     }
2251   }
2252   legCls->Draw();
2253   cf->SetRangeUser(cf->GetVar("multiplicity"), 0.0, 6.0, kTRUE);
2254   
2255   // Qtot vs P
2256   pad = ((TVirtualPad*)l->At(5)); pad->cd();
2257   pad->SetLeftMargin(0.15); pad->SetRightMargin(0.02);
2258   pad->SetTopMargin(0.02); pad->SetBottomMargin(0.15);
2259   pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);
2260   if(gROOT->FindObject("rangeQtot")) delete gROOT->FindObject("rangeQtot");
2261   TH2F* rangeQtot = new TH2F("rangeQtot", "", 10, 0.0, 9.999, 10, 0.0, 1.199);
2262   SetStyle(rangeQtot->GetXaxis(), "Q_{tot} (a.u.)", 0.07, 0.8, kTRUE, 0.05);
2263   SetStyle(rangeQtot->GetYaxis(), "entries (a.u.)", 0.07, 0.8, kTRUE, 0.05);
2264   rangeQtot->SetStats(kFALSE);
2265   rangeQtot->Draw();
2266   
2267   TH1D* hQtot[6]={0x0};
2268   TLegend* leg2=new TLegend(0.6, 0.7, 0.9, 0.97);
2269   leg2->SetFillColor(0);
2270   leg2->SetBorderSize(0);
2271   
2272   for(Int_t iCent=0; iCent<6; ++iCent) {
2273     if(iCent>0)
2274       cf->SetRangeUser(cf->GetVar("multiplicity"), Double_t(iCent), Double_t(iCent), kTRUE);
2275     
2276     hQtot[iCent] = (TH1D*)cf->Project(0, cf->GetVar("Qtot0"));
2277     if(!hQtot[iCent]) continue;
2278     hQtot[iCent]->SetBinContent(1, 0);
2279     
2280     Double_t maximum = hQtot[iCent]->GetMaximum();
2281     if(maximum>1.0)
2282       hQtot[iCent]->Scale(1.0/maximum);
2283     hQtot[iCent]->SetLineColor(iCent<4 ? iCent+1 : iCent+2);
2284     hQtot[iCent]->SetStats(kFALSE);
2285     hQtot[iCent]->SetTitle("");
2286     hQtot[iCent]->SetLineWidth(2);
2287     if(hQtot[iCent]->Integral()>0.01) {
2288       hQtot[iCent]->Draw(iCent==0 ? "" : "same");
2289       leg2->AddEntry(hQtot[iCent], (iCent==0 ? "all centralities" : Form("%.0f < SPD tracklets < %.0f", cf->GetAxis(cf->GetVar("multiplicity"),0)->GetBinLowEdge(iCent), 
2290                                                                            cf->GetAxis(cf->GetVar("multiplicity"),0)->GetBinUpEdge(iCent))), "l");
2291     }
2292   }
2293   leg2->Draw();
2294   cf->SetRangeUser(cf->GetVar("multiplicity"), 0.0, 5.0, kTRUE);
2295   //if(cutTOFbc) cf->SetRangeUser(stepTOFBC, -1000.0, +1000.0);  // reset the cut on TOFbc
2296 }
2297
2298
2299 //_________________________________________________________________
2300 void AliTRDcheckESD::PlotTrackingSummaryFromCF(Double_t* trendValues, const Char_t* /*triggerName*/, Bool_t /*useIsolatedBC*/, Bool_t /*cutTOFbc*/) {
2301   //
2302   //  Plot tracking summary
2303   //
2304   if(!fMatchingPhiEtaCF || !fMatchingPtCF || !fCentralityCF || !fBunchCrossingsCF) return;
2305   AliCFContainer* cf = 0x0;  
2306   
2307   TLatex *lat=new TLatex();
2308   lat->SetTextSize(0.06);
2309   lat->SetTextColor(2);
2310   
2311   gPad->SetTopMargin(0.05); gPad->SetBottomMargin(0.001);
2312   gPad->SetLeftMargin(0.001); gPad->SetRightMargin(0.001);
2313   gPad->Divide(3,3,0.,0.);
2314   TList* l=gPad->GetListOfPrimitives();
2315   
2316   // eta-phi distr. for positive TPC tracks
2317   TVirtualPad* pad = ((TVirtualPad*)l->At(0)); pad->cd();
2318   pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);
2319   pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);
2320   pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);
2321  
2322   //cf->SetRangeUser(stepDCAxy, -0.999, +0.999);
2323   //cf->SetRangeUser(stepDCAz, -3.0, +3.0);
2324   //if(cutTOFbc) cf->SetRangeUser(stepTOFBC, 0.0, 0.0);
2325   
2326   // find all the isolated bunch crossings with entries
2327   TH2D* hTPCrefPos = 0x0; TH2D* hTRDrefPos = 0x0; TH2D* hTOFrefPos = 0x0;
2328   TH2D* hTPCrefNeg = 0x0; TH2D* hTRDrefNeg = 0x0; TH2D* hTOFrefNeg = 0x0;
2329   cf = fMatchingPhiEtaCF;
2330   cf->SetRangeUser(cf->GetVar("tracklets"), 0.0, 6.0);
2331   cf->SetRangeUser(cf->GetVar("charge"), +1.0, +1.0);      // positive charges
2332   hTPCrefPos = (TH2D*)cf->Project(0, cf->GetVar("eta"), cf->GetVar("phi"));
2333   hTRDrefPos = (TH2D*)cf->Project(1, cf->GetVar("eta"), cf->GetVar("phi"));
2334   hTOFrefPos = (TH2D*)cf->Project(2, cf->GetVar("eta"), cf->GetVar("phi"));
2335   cf->SetRangeUser(cf->GetVar("charge"), -1.0, -1.0);      // negative charges
2336   hTPCrefNeg = (TH2D*)cf->Project(0, cf->GetVar("eta"), cf->GetVar("phi"));
2337   hTRDrefNeg = (TH2D*)cf->Project(1, cf->GetVar("eta"), cf->GetVar("phi"));
2338   hTOFrefNeg = (TH2D*)cf->Project(2, cf->GetVar("eta"), cf->GetVar("phi"));
2339   cf->SetRangeUser(cf->GetVar("charge"), -1.0, +1.0);      // reset charge cut
2340     
2341   if(gROOT->FindObject("rangeEtaPhi")) delete gROOT->FindObject("rangeEtaPhi");
2342   TH2F* rangeEtaPhi = new TH2F("rangeEtaPhi", "", 10, -0.99, +0.99, 10, -3.4, +3.4);
2343   SetStyle(rangeEtaPhi->GetXaxis(), "#eta", 0.07, 0.8, kTRUE, 0.05);
2344   SetStyle(rangeEtaPhi->GetYaxis(), "detector #varphi", 0.07, 0.8, kTRUE, 0.05);
2345   rangeEtaPhi->SetStats(kFALSE);  
2346   
2347   //----------------------------------------------
2348   // eta-phi efficiency for positive TRD tracks
2349   pad = ((TVirtualPad*)l->At(0)); pad->cd();
2350   pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);
2351   pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);
2352   pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);
2353   rangeEtaPhi->Draw();
2354   
2355   TH2D* hTRDeffPos = (hTRDrefPos ? (TH2D*)hTRDrefPos->Clone("hTRDeffPos") : 0x0);
2356   if(hTRDeffPos) {
2357     hTRDeffPos->Reset();
2358     hTRDeffPos->SetStats(kFALSE);
2359     hTRDeffPos->Divide(hTRDrefPos, hTPCrefPos);
2360     hTRDeffPos->SetMaximum(1.0);
2361     hTRDeffPos->Draw("samecolz");
2362     lat->DrawLatex(-0.9, 3.6, "TPC-TRD matching for positive tracks");
2363     DrawTRDGrid();
2364   }
2365   
2366   //----------------------------------------------
2367   // eta-phi efficiency for negative TRD tracks
2368   pad = ((TVirtualPad*)l->At(3)); pad->cd();
2369   pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);
2370   pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);
2371   pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);
2372   rangeEtaPhi->Draw();
2373   
2374   TH2D* hTRDeffNeg = (hTRDrefNeg ? (TH2D*)hTRDrefNeg->Clone("hTRDeffNeg") : 0x0);
2375   if(hTRDeffNeg) {
2376     hTRDeffNeg->Reset();
2377     hTRDeffNeg->SetStats(kFALSE);
2378     hTRDeffNeg->Divide(hTRDrefNeg, hTPCrefNeg);
2379     hTRDeffNeg->SetMaximum(1.0);
2380     hTRDeffNeg->Draw("samecolz");
2381     lat->DrawLatex(-0.9, 3.6, "TPC-TRD matching for negative tracks");
2382     DrawTRDGrid();  
2383   }
2384   
2385   //----------------------------------------------
2386   // eta-phi TRD-TOF matching efficiency for positive tracks
2387   pad = ((TVirtualPad*)l->At(1)); pad->cd();
2388   pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);
2389   pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);
2390   pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);
2391   rangeEtaPhi->Draw();
2392   
2393   TH2D* hTOFeffPos = (hTOFrefPos ? (TH2D*)hTOFrefPos->Clone("hTOFeffPos") : 0x0);
2394   if(hTOFeffPos) {
2395     hTOFeffPos->Reset();
2396     hTOFeffPos->SetStats(kFALSE);
2397     hTOFeffPos->Divide(hTOFrefPos, hTRDrefPos);
2398     hTOFeffPos->SetMaximum(1.0);
2399     hTOFeffPos->Draw("samecolz");
2400     lat->DrawLatex(-0.9, 3.6, "TRD-TOF matching for positive tracks");
2401     DrawTRDGrid();
2402   }
2403   
2404   //----------------------------------------------
2405   // eta-phi TRD-TOF matching efficiency for negative tracks
2406   pad = ((TVirtualPad*)l->At(4)); pad->cd();
2407   pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);
2408   pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);
2409   pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);
2410   rangeEtaPhi->Draw();
2411   
2412   TH2D* hTOFeffNeg = (hTOFrefNeg ? (TH2D*)hTOFrefNeg->Clone("hTOFeffNeg") : 0x0);
2413   if(hTOFeffNeg) {
2414     hTOFeffNeg->Reset();
2415     hTOFeffNeg->SetStats(kFALSE);
2416     hTOFeffNeg->Divide(hTOFrefNeg, hTRDrefNeg);
2417     hTOFeffNeg->SetMaximum(1.0);
2418     hTOFeffNeg->Draw("samecolz");
2419     lat->DrawLatex(-0.9, 3.6, "TRD-TOF matching for negative tracks");
2420     DrawTRDGrid();
2421   }
2422   
2423   if(hTRDrefPos) delete hTRDrefPos; if(hTPCrefPos) delete hTPCrefPos; if(hTOFrefPos) delete hTOFrefPos;
2424   if(hTRDrefNeg) delete hTRDrefNeg; if(hTPCrefNeg) delete hTPCrefNeg; if(hTOFrefNeg) delete hTOFrefNeg;
2425   
2426   // switch to the Pt cf container
2427   cf = fMatchingPtCF;
2428   cf->SetRangeUser(cf->GetVar("charge"), +1.0, +1.0);  
2429   TH1F* hTRDEffPtPosAll = EfficiencyFromPhiPt(cf, 0, 6, 1, 0);
2430   TH1F* hTOFEffPtPosAll = EfficiencyFromPhiPt(cf, 0, 6, 2, 1, "pt", "TRDTOF");
2431   TH1F* hTRDEffPtPosTrk4 = EfficiencyFromPhiPt(cf, 4, 4, 1, 0);
2432   TH1F* hTOFEffPtPosTrk4 = EfficiencyFromPhiPt(cf, 4, 4, 2, 1, "pt", "TRDTOF");
2433   TH1F* hTRDEffPtPosTrk5 = EfficiencyFromPhiPt(cf, 5, 5, 1, 0);
2434   TH1F* hTOFEffPtPosTrk5 = EfficiencyFromPhiPt(cf, 5, 5, 2, 1, "pt", "TRDTOF");
2435   TH1F* hTRDEffPtPosTrk6 = EfficiencyFromPhiPt(cf, 6, 6, 1, 0);
2436   TH1F* hTOFEffPtPosTrk6 = EfficiencyFromPhiPt(cf, 6, 6, 2, 1, "pt", "TRDTOF");
2437   
2438   cf->SetRangeUser(cf->GetVar("charge"), -1.0, -1.0);  
2439   TH1F* hTRDEffPtNegAll = EfficiencyFromPhiPt(cf, 0, 6, 1, 0);
2440   TH1F* hTOFEffPtNegAll = EfficiencyFromPhiPt(cf, 0, 6, 2, 1, "pt", "TRDTOF");
2441   TH1F* hTRDEffPtNegTrk4 = EfficiencyFromPhiPt(cf, 4, 4, 1, 0);
2442   TH1F* hTOFEffPtNegTrk4 = EfficiencyFromPhiPt(cf, 4, 4, 2, 1, "pt", "TRDTOF");
2443   TH1F* hTRDEffPtNegTrk5 = EfficiencyFromPhiPt(cf, 5, 5, 1, 0);
2444   TH1F* hTOFEffPtNegTrk5 = EfficiencyFromPhiPt(cf, 5, 5, 2, 1, "pt", "TRDTOF");
2445   TH1F* hTRDEffPtNegTrk6 = EfficiencyFromPhiPt(cf, 6, 6, 1, 0);
2446   TH1F* hTOFEffPtNegTrk6 = EfficiencyFromPhiPt(cf, 6, 6, 2, 1, "pt", "TRDTOF");
2447   cf->SetRangeUser(cf->GetVar("charge"), -1.0, +1.0);  
2448   
2449   TF1* funcConst = new TF1("constFunc", "[0]", 1.0, 3.0);
2450   if(trendValues) {
2451     if(hTRDEffPtPosAll && hTRDEffPtPosAll->Integral()>0.1) {
2452       hTRDEffPtPosAll->Fit(funcConst, "Q0ME", "goff", 1.0, 3.0);
2453       trendValues[0] = funcConst->GetParameter(0);
2454       trendValues[1] = funcConst->GetParError(0);
2455     }
2456   }
2457   if(trendValues) { 
2458     if(hTRDEffPtNegAll && hTRDEffPtNegAll->Integral()>0.1) {
2459       hTRDEffPtNegAll->Fit(funcConst, "Q0ME", "goff", 1.0, 3.0);
2460       trendValues[2] = funcConst->GetParameter(0);
2461       trendValues[3] = funcConst->GetParError(0);
2462     }
2463   }
2464   if(trendValues) { 
2465     if(hTOFEffPtPosAll && hTOFEffPtPosAll->Integral()>0.1) {
2466       hTOFEffPtPosAll->Fit(funcConst, "Q0ME", "goff", 1.0, 3.0);
2467       trendValues[4] = funcConst->GetParameter(0);
2468       trendValues[5] = funcConst->GetParError(0);
2469     }
2470   }
2471   if(trendValues) { 
2472     if(hTOFEffPtNegAll && hTOFEffPtNegAll->Integral()>0.1) {
2473       hTOFEffPtNegAll->Fit(funcConst, "Q0ME", "goff", 1.0, 3.0);
2474       trendValues[6] = funcConst->GetParameter(0);
2475       trendValues[7] = funcConst->GetParError(0);
2476     }
2477   }
2478   
2479   //---------------------------------------------------------
2480   // TPC-TRD matching efficiency vs pt
2481   pad = ((TVirtualPad*)l->At(6)); pad->cd();
2482   pad->SetLeftMargin(0.15); pad->SetRightMargin(0.02);
2483   pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);
2484   pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);
2485   
2486   if(gROOT->FindObject("rangeEffPt2")) delete gROOT->FindObject("rangeEffPt2");
2487   TH2F* rangeEffPt=new TH2F("rangeEffPt2", "",10,0.,10.,10,0.,1.4);
2488   rangeEffPt->SetStats(kFALSE);
2489   SetStyle(rangeEffPt->GetXaxis(), "p_{T} [GeV/c]", 0.07, 0.8, kTRUE, 0.05);
2490   SetStyle(rangeEffPt->GetYaxis(), "efficiency", 0.07, 0.8, kTRUE, 0.05);
2491   rangeEffPt->Draw();
2492   lat->DrawLatex(0.2, 1.42, "TPC-TRD matching efficiency");
2493   //++++++++++++++++++
2494   TLine line;
2495   line.SetLineStyle(2);
2496   line.SetLineWidth(2);
2497   line.DrawLine(rangeEffPt->GetXaxis()->GetXmin(), 0.7, rangeEffPt->GetXaxis()->GetXmax(), 0.7);
2498   line.DrawLine(rangeEffPt->GetXaxis()->GetXmin(), 0.9, rangeEffPt->GetXaxis()->GetXmax(), 0.9);
2499   TLegend* leg=new TLegend(0.2, 0.7, 0.7, 0.89);
2500   leg->SetNColumns(2);
2501   leg->SetMargin(0.15);
2502   leg->SetBorderSize(0);
2503   leg->SetFillColor(0);
2504
2505   SetStyle(hTRDEffPtPosAll, 1, kRed, 1, 24, kRed, 1);
2506   SetStyle(hTRDEffPtNegAll, 1, kBlue, 1, 24, kBlue, 1);
2507   SetStyle(hTRDEffPtPosTrk4, 1, kRed, 1, 25, kRed, 1);
2508   SetStyle(hTRDEffPtNegTrk4, 1, kBlue, 1, 25, kBlue, 1);
2509   SetStyle(hTRDEffPtPosTrk5, 1, kRed, 1, 26, kRed, 1);
2510   SetStyle(hTRDEffPtNegTrk5, 1, kBlue, 1, 26, kBlue, 1);
2511   SetStyle(hTRDEffPtPosTrk6, 1, kRed, 1, 27, kRed, 1);
2512   SetStyle(hTRDEffPtNegTrk6, 1, kBlue, 1, 27, kBlue, 1);
2513   if(hTRDEffPtPosAll) {hTRDEffPtPosAll->Draw("same"); leg->AddEntry(hTRDEffPtPosAll, "pos. (#geq 1 tracklet)", "p");}
2514   if(hTRDEffPtNegAll) {hTRDEffPtNegAll->Draw("same"); leg->AddEntry(hTRDEffPtNegAll, "neg. (#geq 1 tracklet)", "p");}
2515   hTRDEffPtPosTrk4->Draw("same"); leg->AddEntry(hTRDEffPtPosTrk4, "pos. (4 tracklets)", "p");
2516   hTRDEffPtNegTrk4->Draw("same"); leg->AddEntry(hTRDEffPtNegTrk4, "neg. (4 tracklets)", "p");
2517   hTRDEffPtPosTrk5->Draw("same"); leg->AddEntry(hTRDEffPtPosTrk5, "pos. (5 tracklets)", "p");
2518   hTRDEffPtNegTrk5->Draw("same"); leg->AddEntry(hTRDEffPtNegTrk5, "neg. (5 tracklets)", "p");
2519   hTRDEffPtPosTrk6->Draw("same"); leg->AddEntry(hTRDEffPtPosTrk6, "pos. (6 tracklets)", "p");
2520   hTRDEffPtNegTrk6->Draw("same"); leg->AddEntry(hTRDEffPtNegTrk6, "neg. (6 tracklets)", "p");
2521   
2522   leg->Draw();
2523   
2524   //---------------------------------------------------------
2525   // TRD-TOF matching efficiency vs pt
2526   pad = ((TVirtualPad*)l->At(7)); pad->cd();
2527   pad->SetLeftMargin(0.15); pad->SetRightMargin(0.02);
2528   pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);
2529   pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);
2530   
2531   rangeEffPt->Draw();
2532   lat->DrawLatex(0.2, 1.42, "TRD-TOF matching efficiency");
2533   SetStyle(hTOFEffPtPosAll, 1, kRed, 1, 24, kRed, 1);
2534   SetStyle(hTOFEffPtPosTrk4, 1, kRed, 1, 25, kRed, 1);
2535   SetStyle(hTOFEffPtPosTrk5, 1, kRed, 1, 26, kRed, 1);
2536   SetStyle(hTOFEffPtPosTrk6, 1, kRed, 1, 27, kRed, 1);
2537   SetStyle(hTOFEffPtNegAll, 1, kBlue, 1, 24, kBlue, 1);
2538   SetStyle(hTOFEffPtNegTrk4, 1, kBlue, 1, 25, kBlue, 1);
2539   SetStyle(hTOFEffPtNegTrk5, 1, kBlue, 1, 26, kBlue, 1);
2540   SetStyle(hTOFEffPtNegTrk6, 1, kBlue, 1, 27, kBlue, 1);
2541   if(hTOFEffPtPosAll) hTOFEffPtPosAll->Draw("same"); 
2542   hTOFEffPtPosTrk4->Draw("same"); 
2543   hTOFEffPtPosTrk5->Draw("same"); 
2544   hTOFEffPtPosTrk6->Draw("same"); 
2545   if(hTOFEffPtNegAll) hTOFEffPtNegAll->Draw("same"); 
2546   hTOFEffPtNegTrk4->Draw("same"); 
2547   hTOFEffPtNegTrk5->Draw("same"); 
2548   hTOFEffPtNegTrk6->Draw("same");  
2549     
2550   //-----------------------------------------------------
2551   // <ntracklets> vs (phi,eta)
2552   pad = ((TVirtualPad*)l->At(2)); pad->cd();
2553   pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);
2554   pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);
2555   pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);
2556   
2557   rangeEtaPhi->Draw();
2558   lat->DrawLatex(-0.9, 3.6, "TRD <N_{tracklets}>");
2559   
2560   cf = fMatchingPhiEtaCF;
2561   TH3D* hNtracklets = (TH3D*)cf->Project(1, cf->GetVar("phi"), cf->GetVar("eta"), cf->GetVar("tracklets"));
2562   
2563   TProfile2D* hNtrackletsProf = hNtracklets->Project3DProfile();
2564   delete hNtracklets;
2565   if(hNtrackletsProf) {
2566     hNtrackletsProf->SetStats(kFALSE);
2567     hNtrackletsProf->SetMinimum(0.);
2568     hNtrackletsProf->SetMaximum(6.);
2569     hNtrackletsProf->Draw("samecolz");
2570     DrawTRDGrid();
2571   }
2572   
2573   // calculate the trend value for tracklets/track
2574   cf = fMatchingPtCF;
2575   TH2D* hNtrackletsVsP = (TH2D*)cf->Project(1, cf->GetVar("pt"), cf->GetVar("tracklets"));
2576   if(trendValues &&  hNtrackletsVsP && hNtrackletsVsP->GetEntries()>0.1) {
2577     TProfile* hNtrackletsVsPprof = hNtrackletsVsP->ProfileX("hNtrackletsVsPprof");
2578     hNtrackletsVsPprof->Fit(funcConst, "QME0", "goff", 1.0, 3.0);
2579     trendValues[8] = funcConst->GetParameter(0);
2580     trendValues[9] = funcConst->GetParError(0);
2581     delete hNtrackletsVsP;
2582   }
2583       
2584   //--------------------------------------------------------------
2585   // Nclusters per TRD track vs momentum
2586   pad = ((TVirtualPad*)l->At(5)); pad->cd();
2587   pad->SetLeftMargin(0.15); pad->SetRightMargin(0.12);
2588   pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);
2589   pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);
2590   pad->SetLogz();
2591   
2592   if(gROOT->FindObject("rangeNclsP")) delete gROOT->FindObject("rangeNclsP");
2593   TH2F* rangeNclsP = new TH2F("rangeNclsP", "", 10, 0.0, 11.99, 10, 0.0, 199.0);
2594   SetStyle(rangeNclsP->GetXaxis(), "p [GeV/c]", 0.07, 0.8, kTRUE, 0.05);
2595   SetStyle(rangeNclsP->GetYaxis(), "#clusters", 0.07, 0.8, kTRUE, 0.05);
2596   rangeNclsP->SetStats(kFALSE);
2597   rangeNclsP->Draw();
2598   lat->DrawLatex(1.0, 205., "TRD Clusters / track");
2599   
2600   cf = fCentralityCF;
2601   TH2D* hNclsVsP = (TH2D*)cf->Project(0, cf->GetVar("P"), cf->GetVar("clusters"));
2602   if(hNclsVsP) {
2603     hNclsVsP->SetStats(kFALSE);
2604     hNclsVsP->Draw("samecolz");
2605   }    
2606   
2607   if(trendValues && hNclsVsP && hNclsVsP->GetEntries()>10) {
2608     TProfile* hNclsVsPprof = hNclsVsP->ProfileX("hNclsVsPprof");
2609     hNclsVsPprof->Fit(funcConst, "QME0", "goff", 1.0, 3.0);
2610     trendValues[10] = funcConst->GetParameter(0);
2611     trendValues[11] = funcConst->GetParError(0);
2612   }
2613     
2614   //--------------------------------------------------------------
2615   // TRD-TPC and TOF-TRD matching efficiency vs bunch crossing
2616   pad = ((TVirtualPad*)l->At(8)); pad->cd();
2617   pad->SetLeftMargin(0.15); pad->SetRightMargin(0.02);
2618   pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);
2619   pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);
2620
2621   cf = fBunchCrossingsCF;
2622   TH1F* hTRDEffBC = EfficiencyFromPhiPt(cf, 0, 6, 1, 0, "BC", "TPCTRD");
2623   TH1F* hTOFEffBC = EfficiencyFromPhiPt(cf, 0, 6, 2, 1, "BC", "TRDTOF");
2624   
2625   if(gROOT->FindObject("rangeBC")) delete gROOT->FindObject("rangeBC");
2626   TH2F* rangeBC = new TH2F("rangeBC", "", 10, -0.5, 3499.5, 10, 0.0, 1.4);
2627   rangeBC->SetStats(kFALSE);
2628   SetStyle(rangeBC->GetXaxis(), "Bunch crossing", 0.07, 0.8, kTRUE, 0.05);
2629   SetStyle(rangeBC->GetYaxis(), "efficiency", 0.07, 0.8, kTRUE, 0.05);
2630   rangeBC->Draw();
2631   
2632   TLegend* legBC=new TLegend(0.8, 0.7, 0.95, 0.89);
2633   legBC->SetBorderSize(0);
2634   legBC->SetMargin(0.15);
2635   legBC->SetFillColor(0);
2636   if(hTRDEffBC) {
2637     hTRDEffBC->SetStats(kFALSE);
2638     SetStyle(hTRDEffBC, 1, kRed, 2, 24, kRed, 1); legBC->AddEntry(hTRDEffBC, "TPC-TRD", "p");
2639     SetStyle(hTOFEffBC, 1, kBlue, 2, 24, kBlue, 1); legBC->AddEntry(hTOFEffBC, "TRD-TOF", "p");
2640     hTRDEffBC->Draw("same");
2641     hTOFEffBC->Draw("same");
2642     legBC->Draw();
2643     lat->DrawLatex(200., 1.42, "Matching efficiency at 1<p_{T}<3 GeV/c");
2644   }
2645     
2646   // reset the user range on the event multiplicity
2647   //if(cutTOFbc) cf->SetRangeUser(stepTOFBC, -1000.0, +1000.0);  // reset the cut on TOFbc
2648   
2649   delete funcConst;
2650 }
2651
2652
2653 //_________________________________________________________________
2654 void AliTRDcheckESD::PlotPidSummaryFromCF(Double_t* trendValues, const Char_t* /*triggerName*/, Bool_t /*useIsolatedBC*/, Bool_t /*cutTOFbc*/) {
2655   //
2656   // Centrality summary
2657   //
2658   if(!fQtotCF || !fPulseHeightCF || !fCentralityCF) return;
2659   
2660   AliCFContainer* cf = 0x0;
2661   
2662   TLatex *lat=new TLatex();
2663   lat->SetTextSize(0.07);
2664   lat->SetTextColor(2);
2665   gPad->SetTopMargin(0.05); gPad->SetBottomMargin(0.001);
2666   gPad->SetLeftMargin(0.001); gPad->SetRightMargin(0.001);
2667   gPad->Divide(3,3,0.,0.);
2668   TList* l=gPad->GetListOfPrimitives();
2669   
2670   //cf->SetRangeUser(stepDCAxy, -0.999, +0.999);
2671   //cf->SetRangeUser(stepDCAz, -3.0, +3.0);
2672   //if(cutTOFbc) cf->SetRangeUser(stepTOFBC, 0.0, 0.0);
2673   
2674   if(gROOT->FindObject("rangeEtaPhi2")) delete gROOT->FindObject("rangeEtaPhi2");
2675   TH2F* rangeEtaPhi = new TH2F("rangeEtaPhi2", "", 10, -0.99, +0.99, 10, -3.4, +3.4);
2676   SetStyle(rangeEtaPhi->GetXaxis(), "#eta", 0.07, 0.8, kTRUE, 0.05);
2677   SetStyle(rangeEtaPhi->GetYaxis(), "detector #varphi", 0.07, 0.8, kTRUE, 0.05);
2678   rangeEtaPhi->SetStats(kFALSE);  
2679   
2680   // eta-phi distr. for <Qtot> in layer 0
2681   TVirtualPad* pad;
2682   TProfile2D* hProf2D;
2683   cf = fQtotCF;
2684   for(Int_t iLayer=0; iLayer<6; ++iLayer) {
2685     pad = ((TVirtualPad*)l->At((iLayer<3 ? iLayer*3 : (iLayer-3)*3+1))); pad->cd();
2686     pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);
2687     pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);
2688     pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);
2689     rangeEtaPhi->Draw();
2690     
2691     cf->SetRangeUser(cf->GetVar("layer"), Double_t(iLayer), Double_t(iLayer));
2692     TH3D* hQtotEtaPhi = (TH3D*)cf->Project(0, cf->GetVar("phi"), cf->GetVar("eta"), cf->GetVar("Qtot0"));
2693     hProf2D = (hQtotEtaPhi ? hQtotEtaPhi->Project3DProfile() : 0x0);
2694     if(hQtotEtaPhi) delete hQtotEtaPhi;
2695     
2696     if(hProf2D) {
2697       hProf2D->SetName(Form("Qtot_layer%d",iLayer));
2698       hProf2D->SetStats(kFALSE);
2699       hProf2D->SetMinimum(0.);
2700       hProf2D->SetMaximum(4.);
2701       hProf2D->Draw("samecolz");
2702     }
2703     lat->DrawLatex(-0.9, 3.6, Form("TRD <Q_{tot}> Layer %d", iLayer));
2704     DrawTRDGrid();
2705   }
2706   cf->SetRangeUser(cf->GetVar("layer"), 0.0, 5.0);
2707     
2708   // PH versus slice number
2709   pad = ((TVirtualPad*)l->At(2)); pad->cd();
2710   pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);
2711   pad->SetTopMargin(0.03); pad->SetBottomMargin(0.15);
2712   pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);
2713   
2714   if(gROOT->FindObject("rangePHslice")) delete gROOT->FindObject("rangePHslice");
2715   TH2F* rangePHslice=new TH2F("rangePHslice", "", 8, -0.5, 7.5, 10, 0.0, 6.);
2716   rangePHslice->SetStats(kFALSE);
2717   SetStyle(rangePHslice->GetXaxis(), "slice", 0.07, 0.8, kTRUE, 0.05);
2718   SetStyle(rangePHslice->GetYaxis(), "PH", 0.07, 0.8, kTRUE, 0.05);
2719   rangePHslice->Draw();
2720   
2721   TF1* funcPol1 = new TF1("funcPol1", "[0]+[1]*x", 2.9, 6.4);
2722   
2723   cf = fPulseHeightCF;
2724   TH2D* hPH = (TH2D*)cf->Project(0, cf->GetVar("slice"), cf->GetVar("PH0"));
2725   TH1D* hSliceErr = new TH1D(Form("hSliceErr%f", gRandom->Rndm()), "", hPH->GetXaxis()->GetNbins(), hPH->GetXaxis()->GetXbins()->GetArray());
2726   TH1D* hLandauFit = Proj2D(hPH, hSliceErr);
2727   hPH->SetStats(kFALSE);
2728   hPH->Draw("samecolz");
2729   if(trendValues) {
2730     hSliceErr->Fit(funcPol1, "QME0", "goff", 2.9, 6.4);
2731     trendValues[12] = funcPol1->GetParameter(0);  // PH plateau
2732     trendValues[13] = funcPol1->GetParError(0);   // PH plateau
2733     trendValues[14] = funcPol1->GetParameter(1);  // PH slope
2734     trendValues[15] = funcPol1->GetParError(1);   // PH slope
2735   }
2736   hLandauFit->SetLineWidth(2);
2737   hLandauFit->SetLineStyle(2);
2738   hLandauFit->Draw("same");
2739   
2740   delete funcPol1; delete hSliceErr;
2741   
2742   // Qtot vs P
2743   pad = ((TVirtualPad*)l->At(5)); pad->cd();
2744   pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);
2745   pad->SetTopMargin(0.03); pad->SetBottomMargin(0.15);
2746   pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);
2747   pad->SetLogz();
2748   
2749   if(gROOT->FindObject("rangeQtotP")) delete gROOT->FindObject("rangeQtotP");
2750   TH2F* rangeQtotP = new TH2F("rangeQtotP", "", 10, 0.0, 11.99, 10, 0.0, 11.99);
2751   SetStyle(rangeQtotP->GetXaxis(), "P [GeV/c]", 0.07, 0.8, kTRUE, 0.05);
2752   SetStyle(rangeQtotP->GetYaxis(), "Q_{tot}", 0.07, 0.8, kTRUE, 0.05);
2753   rangeQtotP->SetStats(kFALSE);
2754   rangeQtotP->Draw();
2755   
2756   cf = fCentralityCF;
2757   TH2D* hQtotP = (TH2D*)cf->Project(0, cf->GetVar("P"), cf->GetVar("Qtot0"));
2758     
2759   if(hQtotP)
2760     for(Int_t i=1; i<=hQtotP->GetXaxis()->GetNbins(); ++i) 
2761       hQtotP->SetBinContent(i, 1, 0.0);  
2762   TH1D* hQtotProj = (hQtotP ? Proj2D(hQtotP) : 0x0);
2763   if(hQtotProj) SetStyle(hQtotProj, 2, kBlue, 2, 1, kBlue, 1);
2764   if(trendValues && hQtotProj && hQtotProj->GetEntries()>2) {
2765     trendValues[16] = hQtotProj->GetBinContent(hQtotProj->FindBin(1.0));   // Landau MPV at 1GeV/c
2766     trendValues[17] = hQtotProj->GetBinError(hQtotProj->FindBin(1.0));     // Landau width at 1 GeV/c
2767   }
2768   if(hQtotP) {
2769     hQtotP->SetStats(kFALSE);
2770     hQtotP->Draw("samecolz");
2771     hQtotProj->Draw("same");
2772   }
2773   //if(cutTOFbc) cf->SetRangeUser(stepTOFBC, -1000.0, +1000.0);  // reset the cut on TOFbc
2774 }
2775
2776
2777 //_________________________________________________________________
2778 Bool_t AliTRDcheckESD::PlotCentSummary(Double_t* /*trendValues*/) {
2779
2780   Bool_t isGoodForSaving=kFALSE;
2781   
2782   TLatex* lat=new TLatex();
2783   lat->SetTextSize(0.06);
2784   lat->SetTextColor(2);
2785
2786   gPad->SetTopMargin(0.05); gPad->SetBottomMargin(0.001);
2787   gPad->SetLeftMargin(0.001); gPad->SetRightMargin(0.001);
2788   gPad->Divide(3,3,0.,0.);
2789   TList* l=gPad->GetListOfPrimitives();
2790   
2791   TPad* pad=0x0;  
2792
2793   if(gROOT->FindObject("rangeEffPt")) delete gROOT->FindObject("rangeEffPt");
2794   TH2F* rangeEffPt=new TH2F("rangeEffPt", "",10,0.,10.,10,0.,1.4);
2795   rangeEffPt->SetStats(kFALSE);
2796   SetStyle(rangeEffPt->GetXaxis(), "p_{T} [GeV/c]", 0.07, 0.8, kTRUE, 0.05);
2797   SetStyle(rangeEffPt->GetYaxis(), "efficiency", 0.07, 0.8, kTRUE, 0.05);
2798   
2799   TH3F *h3(NULL), *h3p(NULL), *h3n(NULL);
2800   Int_t padsForEffs[5] = {0,3,6,1,4};
2801   for(Int_t iCent=1; iCent<6; ++iCent) {
2802     // TPC-TRD matching efficiencies
2803     pad = ((TPad*)l->At(padsForEffs[iCent-1])); pad->cd();
2804     pad->SetLeftMargin(0.15); pad->SetRightMargin(0.02); pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);
2805     pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);
2806     
2807     if(!(h3p = dynamic_cast<TH3F*>(fHistos->At(kTPCRefTracksPos+iCent)))) continue;
2808     if(!(h3n = dynamic_cast<TH3F*>(fHistos->At(kTPCRefTracksNeg+iCent)))) continue;
2809     // =============================================
2810     if(!(h3 = dynamic_cast<TH3F*>(fHistos->At(kTRDRefTracksPos+iCent)))) continue;
2811     TH1F* hFeffP = EfficiencyTRD(h3p, h3, kTRUE);
2812     //
2813     if(!(h3 = dynamic_cast<TH3F*>(fHistos->At(kTRDRefTracksNeg+iCent)))) continue;
2814     TH1F* hFeffN = EfficiencyTRD(h3n, h3, kTRUE);
2815     // =============================================
2816     if(!(h3 = dynamic_cast<TH3F*>(fHistos->At(kTRDRefTracksPos4+iCent)))) continue;
2817     TH1F* hFeffP4 = EfficiencyTRD(h3p, h3, kTRUE);
2818     //
2819     if(!(h3 = dynamic_cast<TH3F*>(fHistos->At(kTRDRefTracksNeg4+iCent)))) continue;
2820     TH1F* hFeffN4 = EfficiencyTRD(h3n, h3, kTRUE);
2821     // =============================================
2822     if(!(h3 = dynamic_cast<TH3F*>(fHistos->At(kTRDRefTracksPos5+iCent)))) continue;
2823     TH1F* hFeffP5 = EfficiencyTRD(h3p, h3, kTRUE);
2824     //
2825     if(!(h3 = dynamic_cast<TH3F*>(fHistos->At(kTRDRefTracksNeg5+iCent)))) continue;
2826     TH1F* hFeffN5 = EfficiencyTRD(h3n, h3, kTRUE);
2827     // =============================================
2828     if(!(h3 = dynamic_cast<TH3F*>(fHistos->At(kTRDRefTracksPos6+iCent)))) continue;
2829     TH1F* hFeffP6 = EfficiencyTRD(h3p, h3, kTRUE);
2830     //
2831     if(!(h3 = dynamic_cast<TH3F*>(fHistos->At(kTRDRefTracksNeg6+iCent)))) continue;
2832     TH1F* hFeffN6 = EfficiencyTRD(h3n, h3, kTRUE);
2833   
2834     rangeEffPt->Draw();
2835     
2836     TLine line;
2837     line.SetLineStyle(2);
2838     line.SetLineWidth(2);
2839     line.DrawLine(rangeEffPt->GetXaxis()->GetXmin(), 0.7, rangeEffPt->GetXaxis()->GetXmax(), 0.7);
2840     line.DrawLine(rangeEffPt->GetXaxis()->GetXmin(), 0.9, rangeEffPt->GetXaxis()->GetXmax(), 0.9);
2841     line.SetLineStyle(1);
2842     line.SetLineWidth(1);
2843     line.DrawLine(rangeEffPt->GetXaxis()->GetXmin(), 1.0, rangeEffPt->GetXaxis()->GetXmax(), 1.0);
2844     if(hFeffP) SetStyle(hFeffP, 1, kRed, 1, 24, kRed, 1);
2845     if(hFeffP4) SetStyle(hFeffP4, 1, kRed, 1, 25, kRed, 1);
2846     if(hFeffP5) SetStyle(hFeffP5, 1, kRed, 1, 26, kRed, 1);
2847     if(hFeffP6) SetStyle(hFeffP6, 1, kRed, 1, 27, kRed, 1);
2848     if(hFeffN) SetStyle(hFeffN, 1, kBlue, 1, 24, kBlue, 1);
2849     if(hFeffN4) SetStyle(hFeffN4, 1, kBlue, 1, 25, kBlue, 1);
2850     if(hFeffN5) SetStyle(hFeffN5, 1, kBlue, 1, 26, kBlue, 1);
2851     if(hFeffN6) SetStyle(hFeffN6, 1, kBlue, 1, 27, kBlue, 1);
2852     
2853     TLegend* leg=new TLegend(0.16, 0.7, 0.61, 0.89);
2854     leg->SetFillColor(0);
2855     leg->SetNColumns(2);
2856     leg->SetTextSize(0.039);
2857     leg->SetMargin(0.1);
2858     if(hFeffP && hFeffP->Integral()>0.001) {
2859       isGoodForSaving = kTRUE;
2860       hFeffP->Draw("same");
2861       leg->AddEntry(hFeffP, "pos. (#geq 1 trcklt)", "p");
2862     }
2863     if(hFeffN && hFeffN->Integral()>0.001) {
2864       isGoodForSaving = kTRUE;
2865       hFeffN->Draw("same");
2866       leg->AddEntry(hFeffN, "neg. (#geq 1 trcklt)", "p");
2867     }
2868     if(hFeffP4 && hFeffP4->Integral()>0.001) {
2869       isGoodForSaving = kTRUE;
2870       hFeffP4->Draw("same");
2871       leg->AddEntry(hFeffP4, "pos. (4 trcklts)", "p");
2872     }
2873     if(hFeffN4 && hFeffN4->Integral()>0.001) {
2874       isGoodForSaving = kTRUE;
2875       hFeffN4->Draw("same");
2876       leg->AddEntry(hFeffN4, "neg. (4 trcklts)", "p");
2877     }
2878     if(hFeffP5 && hFeffP5->Integral()>0.001) {
2879       isGoodForSaving = kTRUE;
2880       hFeffP5->Draw("same");
2881       leg->AddEntry(hFeffP5, "pos. (5 trcklts)", "p");
2882     }
2883     if(hFeffN5 && hFeffN5->Integral()>0.001) {
2884       isGoodForSaving = kTRUE;
2885       hFeffN5->Draw("same");
2886       leg->AddEntry(hFeffN5, "neg. (5 trcklts)", "p");
2887     }
2888     if(hFeffP6 && hFeffP6->Integral()>0.001) {
2889       isGoodForSaving = kTRUE;
2890       hFeffP6->Draw("same");
2891       leg->AddEntry(hFeffP6, "pos. (6 trcklts)", "p");
2892     }
2893     if(hFeffN6 && hFeffN6->Integral()>0.001) {
2894       isGoodForSaving = kTRUE;
2895       hFeffN6->Draw("same");
2896       leg->AddEntry(hFeffN6, "neg. (6 trklts)", "p");
2897     }
2898     
2899     if(isGoodForSaving) {
2900       if(iCent==1) leg->Draw();
2901       lat->DrawLatex(5.6, 1.3, Form("Centrality class %d", iCent));
2902       lat->DrawLatex(0.5, 1.42, "TPC-TRD matching efficiency");
2903     }
2904   }   // end loop over multiplicity intervals
2905
2906   // Number of clusters per TRD track
2907   pad = ((TPad*)l->At(2)); pad->cd();
2908   pad->SetLeftMargin(0.15); pad->SetRightMargin(0.02);
2909   pad->SetTopMargin(0.02); pad->SetBottomMargin(0.15);
2910   pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);
2911   if(gROOT->FindObject("rangeNcls")) delete gROOT->FindObject("rangeNcls");
2912   TH2F* rangeNcls = new TH2F("rangeNcls", "", 10, 0.0, 199.9, 10, 0.0, 1.199);
2913   SetStyle(rangeNcls->GetXaxis(), "# TRD clusters", 0.07, 0.8, kTRUE, 0.05);
2914   SetStyle(rangeNcls->GetYaxis(), "entries (a.u.)", 0.07, 0.8, kTRUE, 0.05);
2915   rangeNcls->SetStats(kFALSE);
2916   rangeNcls->Draw();
2917   
2918   TH2F* h2F[6]; TH1D* proj[6];
2919   TLegend* leg=new TLegend(0.2, 0.7, 0.5, 0.95);
2920   leg->SetFillColor(0);
2921   Bool_t isGood=kFALSE;
2922   for(Int_t iCent=0; iCent<6; ++iCent) {
2923     h2F[iCent] = dynamic_cast<TH2F*>(fHistos->At(kNClsTrackTRD+iCent));
2924     proj[iCent] = (h2F[iCent] && h2F[iCent]->GetEntries()>10 ? h2F[iCent]->ProjectionY(Form("projCent%d",iCent)) : 0x0);
2925     if(proj[iCent]) {
2926       proj[iCent]->SetLineColor(iCent<4 ? iCent+1 : iCent+2);
2927       Double_t maximum = proj[iCent]->GetMaximum();
2928       if(maximum>1.0)
2929         proj[iCent]->Scale(1.0/maximum);
2930       proj[iCent]->SetStats(kFALSE);
2931       proj[iCent]->Draw("same");
2932       leg->AddEntry(proj[iCent], (iCent==0 ? "all centralities" : Form("centrality class %d", iCent)), "l");
2933       isGood = kTRUE;
2934     }
2935   }
2936   if(isGood) leg->Draw();
2937   isGoodForSaving = isGoodForSaving || isGood;
2938   
2939   // Qtot vs P
2940   pad = ((TPad*)l->At(5)); pad->cd();
2941   pad->SetLeftMargin(0.15); pad->SetRightMargin(0.02);
2942   pad->SetTopMargin(0.02); pad->SetBottomMargin(0.15);
2943   pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);
2944   if(gROOT->FindObject("rangeQtot")) delete gROOT->FindObject("rangeQtot");
2945   TH2F* rangeQtot = new TH2F("rangeQtot", "", 10, 0.0, 9.999, 10, 0.0, 1.199);
2946   SetStyle(rangeQtot->GetXaxis(), "Q_{tot} (a.u.)", 0.07, 0.8, kTRUE, 0.05);
2947   SetStyle(rangeQtot->GetYaxis(), "entries (a.u.)", 0.07, 0.8, kTRUE, 0.05);
2948   rangeQtot->SetStats(kFALSE);
2949   rangeQtot->Draw();
2950   
2951   TH1D* projQ[6];
2952   TLegend* leg2=new TLegend(0.6, 0.7, 0.9, 0.95);
2953   leg2->SetFillColor(0);
2954   isGood = kFALSE;
2955   for(Int_t iCent=0; iCent<6; ++iCent) {  
2956     h2F[iCent] = dynamic_cast<TH2F*>(fHistos->At(kQtotP+iCent));
2957     projQ[iCent] = (h2F[iCent] && h2F[iCent]->GetEntries()>10 ? h2F[iCent]->ProjectionY(Form("projQCent%d",iCent)) : 0x0);
2958     if(projQ[iCent]) {
2959       projQ[iCent]->SetLineColor(iCent<4 ? iCent+1 : iCent+2);
2960       Double_t maximum = projQ[iCent]->GetMaximum();
2961       if(maximum>1.0)
2962         projQ[iCent]->Scale(1.0/maximum);
2963       projQ[iCent]->SetStats(kFALSE);
2964       projQ[iCent]->Draw("same");
2965       leg2->AddEntry(projQ[iCent], (iCent==0 ? "all centralities" : Form("centrality class %d", iCent)), "l");
2966       isGood = kTRUE;
2967     }
2968   }
2969   if(isGood) leg2->Draw();
2970   isGoodForSaving = isGoodForSaving || isGood;
2971   return isGoodForSaving;
2972 }
2973
2974
2975 //_________________________________________________________________
2976 Bool_t AliTRDcheckESD::PlotTrackingSummary(Int_t centralityClass, Double_t* trendValues) {
2977
2978   Bool_t isGoodForSaving=kFALSE;
2979   
2980   TLatex *lat=new TLatex();
2981   lat->SetTextSize(0.07);
2982   lat->SetTextColor(2);
2983   gPad->SetTopMargin(0.05); gPad->SetBottomMargin(0.001);
2984   gPad->SetLeftMargin(0.001); gPad->SetRightMargin(0.001);
2985   gPad->Divide(3,3,0.,0.);
2986   TList* l=gPad->GetListOfPrimitives();
2987   // eta-phi distr. for positive TPC tracks
2988   TVirtualPad* pad = ((TVirtualPad*)l->At(0)); pad->cd();
2989   pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);
2990   pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);
2991   pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);
2992   if(gROOT->FindObject("rangeEtaPhi")) delete gROOT->FindObject("rangeEtaPhi");
2993   TH2F* rangeEtaPhi = new TH2F("rangeEtaPhi", "", 10, -0.99, +0.99, 10, -3.4, +3.4);
2994   SetStyle(rangeEtaPhi->GetXaxis(), "#eta", 0.07, 0.8, kTRUE, 0.05);
2995   SetStyle(rangeEtaPhi->GetYaxis(), "detector #varphi", 0.07, 0.8, kTRUE, 0.05);
2996   rangeEtaPhi->SetStats(kFALSE);
2997   rangeEtaPhi->Draw();
2998   lat->DrawLatex(-0.9, 3.6, "TPC positive ref. tracks");
2999   
3000   TH3F* h3F = dynamic_cast<TH3F*>(fHistos->At(kTPCRefTracksPos+centralityClass));
3001   TH2F* h2FtpcP = 0x0;
3002   Float_t nada=0.0;
3003   if(h3F && h3F->GetEntries()>10) {
3004     h2FtpcP = (TH2F*)Proj3D((TH3F*)h3F, 0x0, 1, h3F->GetZaxis()->GetNbins(), nada)->Clone();
3005     h2FtpcP->SetStats(kFALSE);
3006     h2FtpcP->Draw("samecolz"); isGoodForSaving = kTRUE;
3007     isGoodForSaving = kTRUE;
3008   }
3009   //-----------------
3010   // eta-phi distr. for negative TPC tracks
3011   pad = ((TVirtualPad*)l->At(1)); pad->cd();
3012   pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);
3013   pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);
3014   pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);
3015   rangeEtaPhi->Draw();
3016   lat->DrawLatex(-0.9, 3.6, "TPC negative ref. tracks");
3017   
3018   h3F = dynamic_cast<TH3F*>(fHistos->At(kTPCRefTracksNeg+centralityClass));
3019   TH2F* h2FtpcN = 0x0;
3020   if(h3F && h3F->GetEntries()>10) {
3021     h2FtpcN = (TH2F*)Proj3D((TH3F*)h3F, 0x0, 1, h3F->GetZaxis()->GetNbins(), nada)->Clone();
3022     h2FtpcN->SetStats(kFALSE);
3023     h2FtpcN->Draw("samecolz"); 
3024     isGoodForSaving = kTRUE;
3025   }
3026   //----------------------------------------------
3027   // eta-phi distr. for positive TRD tracks
3028   pad = ((TVirtualPad*)l->At(3)); pad->cd();
3029   pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);
3030   pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);
3031   pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);
3032   rangeEtaPhi->Draw();
3033   lat->DrawLatex(-0.9, 3.6, "TRD positive ref. tracks");
3034   
3035   h3F = dynamic_cast<TH3F*>(fHistos->At(kTRDRefTracksPos+centralityClass));
3036   TH2F* h2FtrdP = 0x0;
3037   if(h3F && h3F->GetEntries()>10) {
3038     h2FtrdP = (TH2F*)Proj3D((TH3F*)h3F, 0x0, 1, h3F->GetZaxis()->GetNbins(), nada)->Clone();
3039     h2FtrdP->SetStats(kFALSE);
3040     h2FtrdP->SetMaximum((h2FtpcP ? h2FtpcP->GetMaximum() : h2FtrdP->GetMaximum()));
3041     h2FtrdP->Draw("samecolz"); 
3042     isGoodForSaving=kTRUE;
3043   }
3044   //--------------------------------------------
3045   // eta-phi distr. for negative TRD tracks
3046   pad = ((TVirtualPad*)l->At(4)); pad->cd();
3047   pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);
3048   pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);
3049   pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);
3050   rangeEtaPhi->Draw();
3051   lat->DrawLatex(-0.9, 3.6, "TRD negative ref. tracks");
3052   
3053   h3F = dynamic_cast<TH3F*>(fHistos->At(kTRDRefTracksNeg+centralityClass));
3054   TH2F* h2FtrdN = 0x0;
3055   if(h3F && h3F->GetEntries()>10) {
3056     h2FtrdN = (TH2F*)Proj3D((TH3F*)h3F, 0x0, 1, h3F->GetZaxis()->GetNbins(), nada)->Clone();
3057     h2FtrdN->SetStats(kFALSE);
3058     h2FtrdN->SetMaximum(h2FtpcN ? h2FtpcN->GetMaximum() : h2FtrdN->GetMaximum());
3059     h2FtrdN->Draw("samecolz"); 
3060     isGoodForSaving=kTRUE;
3061   }
3062   //----------------------------------------------
3063   // eta-phi efficiency for positive TRD tracks
3064   pad = ((TVirtualPad*)l->At(6)); pad->cd();
3065   pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);
3066   pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);
3067   pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);
3068   rangeEtaPhi->Draw();
3069   lat->DrawLatex(-0.9, 3.6, "Efficiency positive tracks");
3070   
3071   TH2F* h2Feff = (h2FtrdP ? (TH2F*)h2FtrdP->Clone("h2FeffPos") : 0x0);
3072   if(h2Feff) {
3073     h2Feff->Reset();
3074     h2Feff->SetStats(kFALSE);
3075     h2Feff->Divide(h2FtrdP, h2FtpcP);
3076     h2Feff->SetMaximum(1.0);
3077     if(h2Feff->GetEntries()>1) { h2Feff->Draw("samecolz"); isGoodForSaving=kTRUE; }
3078   }
3079   //-------------------------------------------------
3080   // eta-phi efficiency for negative TRD tracks
3081   pad = ((TVirtualPad*)l->At(7)); pad->cd();
3082   pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);
3083   pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);
3084   pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);
3085   rangeEtaPhi->Draw();
3086   lat->DrawLatex(-0.9, 3.6, "Efficiency negative tracks");
3087   
3088   h2Feff = (h2FtrdN ? (TH2F*)h2FtrdN->Clone("h2FeffNeg") : 0x0);
3089   if(h2Feff) {
3090     h2Feff->Reset();
3091     h2Feff->SetStats(kFALSE);
3092     h2Feff->Divide(h2FtrdN, h2FtpcN);
3093     h2Feff->SetMaximum(1.0);
3094     if(h2Feff->GetEntries()>0.1) { h2Feff->Draw("samecolz"); isGoodForSaving=kTRUE; }
3095   }
3096   //-----------------------------------------------------
3097   // <ntracklets> vs (phi,eta)
3098   pad = ((TVirtualPad*)l->At(2)); pad->cd();
3099   pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);
3100   pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);
3101   pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);
3102   rangeEtaPhi->Draw();
3103   lat->DrawLatex(-0.9, 3.6, "TRD <N_{tracklets}>");
3104   
3105   TProfile2D* hProf2D;
3106   if((hProf2D = dynamic_cast<TProfile2D*>(fHistos->At(kTRDEtaPhiAvNtrkl+centralityClass)))) {
3107     if(hProf2D->GetEntries()>10) {
3108       hProf2D->SetStats(kFALSE);
3109       hProf2D->SetMinimum(0.);
3110       hProf2D->SetMaximum(6.);
3111       if(hProf2D->GetEntries()>1) { hProf2D->Draw("samecolz"); isGoodForSaving = kTRUE; }
3112     }
3113   }
3114   //---------------------------------------------------------
3115   // TPC-TRD matching efficiency vs pt
3116   pad = ((TVirtualPad*)l->At(5)); pad->cd();
3117   pad->SetLeftMargin(0.15); pad->SetRightMargin(0.02);
3118   pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);
3119   pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);
3120   TH1F* hFeffP = EfficiencyTRD(dynamic_cast<TH3F*>(fHistos->At(kTPCRefTracksPos+centralityClass)),
3121               dynamic_cast<TH3F*>(fHistos->At(kTRDRefTracksPos+centralityClass)), kTRUE);
3122   TH1F* hFeffN = EfficiencyTRD(dynamic_cast<TH3F*>(fHistos->At(kTPCRefTracksNeg+centralityClass)),
3123               dynamic_cast<TH3F*>(fHistos->At(kTRDRefTracksNeg+centralityClass)), kTRUE);
3124   TH1F* hFeffP4 = EfficiencyTRD(dynamic_cast<TH3F*>(fHistos->At(kTPCRefTracksPos+centralityClass)),
3125         dynamic_cast<TH3F*>(fHistos->At(kTRDRefTracksPos4+centralityClass)), kTRUE);
3126   TH1F* hFeffN4 = EfficiencyTRD(dynamic_cast<TH3F*>(fHistos->At(kTPCRefTracksNeg+centralityClass)),
3127         dynamic_cast<TH3F*>(fHistos->At(kTRDRefTracksNeg4+centralityClass)), kTRUE);
3128   TH1F* hFeffP5 = EfficiencyTRD(dynamic_cast<TH3F*>(fHistos->At(kTPCRefTracksPos+centralityClass)),
3129         dynamic_cast<TH3F*>(fHistos->At(kTRDRefTracksPos5+centralityClass)), kTRUE);
3130   TH1F* hFeffN5 = EfficiencyTRD(dynamic_cast<TH3F*>(fHistos->At(kTPCRefTracksNeg+centralityClass)),
3131         dynamic_cast<TH3F*>(fHistos->At(kTRDRefTracksNeg5+centralityClass)), kTRUE);
3132   TH1F* hFeffP6 = EfficiencyTRD(dynamic_cast<TH3F*>(fHistos->At(kTPCRefTracksPos+centralityClass)),
3133         dynamic_cast<TH3F*>(fHistos->At(kTRDRefTracksPos6+centralityClass)), kTRUE);
3134   TH1F* hFeffN6 = EfficiencyTRD(dynamic_cast<TH3F*>(fHistos->At(kTPCRefTracksNeg+centralityClass)),
3135         dynamic_cast<TH3F*>(fHistos->At(kTRDRefTracksNeg6+centralityClass)), kTRUE);
3136   
3137   TF1* funcConst = new TF1("funcConst", "[0]", 1.0, 3.0);
3138   
3139   if(gROOT->FindObject("rangeEffPt2")) delete gROOT->FindObject("rangeEffPt2");
3140   TH2F* rangeEffPt2=new TH2F("rangeEffPt2", "",10,0.,10.,10,0.,1.4);
3141   rangeEffPt2->SetStats(kFALSE);
3142   SetStyle(rangeEffPt2->GetXaxis(), "p_{T} [GeV/c]", 0.07, 0.8, kTRUE, 0.05);
3143   SetStyle(rangeEffPt2->GetYaxis(), "efficiency", 0.07, 0.8, kTRUE, 0.05);
3144   rangeEffPt2->Draw();
3145   lat->DrawLatex(0.5, 1.42, "TRD-TPC matching efficiency");
3146   //++++++++++++++++++
3147   TLine line;
3148   line.SetLineStyle(2);
3149   line.SetLineWidth(2);
3150   line.DrawLine(rangeEffPt2->GetXaxis()->GetXmin(), 0.7, rangeEffPt2->GetXaxis()->GetXmax(), 0.7);
3151   line.DrawLine(rangeEffPt2->GetXaxis()->GetXmin(), 0.9, rangeEffPt2->GetXaxis()->GetXmax(), 0.9);
3152   line.SetLineStyle(1);
3153   line.SetLineWidth(1);
3154   line.DrawLine(rangeEffPt2->GetXaxis()->GetXmin(), 1.0, rangeEffPt2->GetXaxis()->GetXmax(), 1.0);
3155   TLegend* leg=new TLegend(0.2, 0.7, 0.6, 0.89);
3156   leg->SetNColumns(2);
3157   leg->SetFillColor(0);
3158   if(hFeffP){
3159     hFeffP->SetMarkerStyle(24);
3160     hFeffP->SetMarkerColor(2);
3161     hFeffP->SetLineColor(2);
3162     if(trendValues && hFeffP->GetEntries()>1) {
3163       hFeffP->Fit(funcConst, "QME0", "goff", 1.0, 3.0);
3164       trendValues[0] = funcConst->GetParameter(0);
3165       trendValues[1] = funcConst->GetParError(0);
3166     }
3167     if(hFeffP->Integral()>0.001) {
3168       hFeffP->Draw("same"); 
3169       leg->AddEntry(hFeffP, "positives (#geq 1 tracklet)", "p");
3170     }
3171   }
3172   if(hFeffN){
3173     hFeffN->SetMarkerStyle(24);
3174     hFeffN->SetMarkerColor(4);
3175     hFeffN->SetLineColor(4);
3176     if(trendValues && hFeffN->GetEntries()>1) {
3177       hFeffN->Fit(funcConst, "QME0", "goff", 1.0, 3.0);
3178       trendValues[2] = funcConst->GetParameter(0);
3179       trendValues[3] = funcConst->GetParError(0);
3180     }
3181     if(hFeffN->Integral()>0.001) {
3182       hFeffN->Draw("same"); 
3183       leg->AddEntry(hFeffN, "negatives (#geq 1 tracklet)", "p");
3184     }
3185   }
3186   if(hFeffP4){
3187     hFeffP4->SetMarkerStyle(25);
3188     hFeffP4->SetMarkerColor(2);
3189     hFeffP4->SetLineColor(2);
3190     if(hFeffP4->Integral()>0.001) {
3191       hFeffP4->Draw("same"); 
3192       leg->AddEntry(hFeffP4, "positives (4 tracklets)", "p");
3193     }
3194   }
3195   if(hFeffN4){
3196     hFeffN4->SetMarkerStyle(25);
3197     hFeffN4->SetMarkerColor(4);
3198     hFeffN4->SetLineColor(4);
3199     if(hFeffN4->Integral()>0.001) {
3200       hFeffN4->Draw("same"); 
3201       leg->AddEntry(hFeffN4, "negatives (4 tracklets)", "p");
3202     }
3203   }
3204   if(hFeffP5){
3205     hFeffP5->SetMarkerStyle(26);
3206     hFeffP5->SetMarkerColor(2);
3207     hFeffP5->SetLineColor(2);
3208     if(hFeffP5->Integral()>0.001) {
3209       hFeffP5->Draw("same"); 
3210       leg->AddEntry(hFeffP5, "positives (5 tracklets)", "p");
3211     }
3212   }
3213   if(hFeffN5){
3214     hFeffN5->SetMarkerStyle(26);
3215     hFeffN5->SetMarkerColor(4);
3216     hFeffN5->SetLineColor(4);
3217     if(hFeffN5->Integral()>0.001) {
3218       hFeffN5->Draw("same"); 
3219       leg->AddEntry(hFeffN5, "negatives (5 tracklets)", "p");
3220     }
3221   }
3222   if(hFeffP6){
3223     hFeffP6->SetMarkerStyle(27);
3224     hFeffP6->SetMarkerColor(2);
3225     hFeffP6->SetLineColor(2);
3226     if(hFeffP6->Integral()>0.001) {
3227       hFeffP6->Draw("same"); 
3228       leg->AddEntry(hFeffP6, "positives (6 tracklets)", "p");
3229     }
3230   }
3231   if(hFeffN6){
3232     hFeffN6->SetMarkerStyle(27);
3233     hFeffN6->SetMarkerColor(4);
3234     hFeffN6->SetLineColor(4);
3235     if(hFeffN6->Integral()>0.001) {
3236       hFeffN6->Draw("same"); 
3237       leg->AddEntry(hFeffN6, "negatives (6 tracklets)", "p");
3238     }
3239   }
3240   leg->Draw();
3241   
3242   //--------------------------------------------------------------
3243   // Nclusters per TRD track
3244   pad = ((TVirtualPad*)l->At(8)); pad->cd();
3245   pad->SetLeftMargin(0.15); pad->SetRightMargin(0.12);
3246   pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);
3247   pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);
3248   pad->SetLogz();
3249   if(gROOT->FindObject("rangeNclsP")) delete gROOT->FindObject("rangeNclsP");
3250   TH2F* rangeNclsP = new TH2F("rangeNclsP", "", 10, 0.0, 11.99, 10, 0.0, 199.0);
3251   SetStyle(rangeNclsP->GetXaxis(), "p [GeV/c]", 0.07, 0.8, kTRUE, 0.05);
3252   SetStyle(rangeNclsP->GetYaxis(), "#clusters", 0.07, 0.8, kTRUE, 0.05);
3253   rangeNclsP->SetStats(kFALSE);
3254   rangeNclsP->Draw();
3255   lat->DrawLatex(1.0, 205., "TRD Clusters / track");
3256   
3257   TH2F* hNclsVsP=0x0;
3258   if((hNclsVsP = dynamic_cast<TH2F*>(fHistos->At(kNClsTrackTRD+centralityClass)))) {
3259     hNclsVsP->SetStats(kFALSE);
3260     if(hNclsVsP->GetEntries()>10) {
3261       hNclsVsP->Draw("samecolz"); isGoodForSaving=kTRUE;
3262       if(trendValues) {
3263         TProfile* h2FProf = hNclsVsP->ProfileX("nclsVsPprof");
3264         h2FProf->Fit(funcConst, "QME0", "goff", 1.0, 3.0);
3265         trendValues[4] = funcConst->GetParameter(0);
3266         trendValues[5] = funcConst->GetParError(0);
3267       }
3268     }
3269   }
3270   
3271   delete funcConst;
3272   return isGoodForSaving;
3273 }
3274
3275
3276 //_________________________________________________________________
3277 Bool_t AliTRDcheckESD::PlotPidSummary(Int_t centralityClass, Double_t* trendValues) {
3278
3279   Bool_t isGoodForSaving=kFALSE;
3280   
3281   TLatex *lat=new TLatex();
3282   lat->SetTextSize(0.07);
3283   lat->SetTextColor(2);
3284   gPad->SetTopMargin(0.05); gPad->SetBottomMargin(0.001);
3285   gPad->SetLeftMargin(0.001); gPad->SetRightMargin(0.001);
3286   gPad->Divide(3,3,0.,0.);
3287   TList* l=gPad->GetListOfPrimitives();
3288   // eta-phi distr. for <Qtot> in layer 0
3289   TVirtualPad* pad;
3290   TProfile2D* hProf2D;
3291   if(gROOT->FindObject("rangeEtaPhi2")) delete gROOT->FindObject("rangeEtaPhi2");
3292   TH2F* rangeEtaPhi = new TH2F("rangeEtaPhi2", "", 10, -0.99, +0.99, 10, -3.4, +3.4);
3293   SetStyle(rangeEtaPhi->GetXaxis(), "#eta", 0.07, 0.8, kTRUE, 0.05);
3294   SetStyle(rangeEtaPhi->GetYaxis(), "detector #varphi", 0.07, 0.8, kTRUE, 0.05);
3295   rangeEtaPhi->SetStats(kFALSE);
3296   
3297   for(Int_t iLayer=0; iLayer<6; ++iLayer) {
3298     pad = ((TVirtualPad*)l->At((iLayer<3 ? iLayer*3 : (iLayer-3)*3+1))); pad->cd();
3299     pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);
3300     pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);
3301     pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);
3302     rangeEtaPhi->Draw();
3303     lat->DrawLatex(-0.9, 3.6, Form("TRD <Q_{tot}> Layer %d", iLayer));
3304     
3305     if(!(hProf2D = dynamic_cast<TProfile2D*>(fHistos->At(kTRDEtaPhiAvQtot+6*centralityClass+iLayer)))) continue;
3306     if(hProf2D && hProf2D->GetEntries()>10) {
3307       hProf2D->SetStats(kFALSE);
3308       hProf2D->SetMinimum(0.);
3309       hProf2D->SetMaximum(4.);
3310       if(hProf2D->GetEntries()>10) { hProf2D->Draw("samecolz"); isGoodForSaving=kTRUE; }
3311     }
3312   }
3313     
3314   // PH versus slice number
3315   pad = ((TVirtualPad*)l->At(2)); pad->cd();
3316   pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);
3317   pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);
3318   pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);
3319   if(gROOT->FindObject("rangePHslice")) delete gROOT->FindObject("rangePHslice");
3320   TH2F* rangePHslice = new TH2F("rangePHslice", "", 10, -0.5, 7.5, 10, 0.0, 2000.0);
3321   SetStyle(rangePHslice->GetXaxis(), "slice", 0.07, 0.8, kTRUE, 0.05);
3322   SetStyle(rangePHslice->GetYaxis(), "PH", 0.07, 0.8, kTRUE, 0.05);
3323   rangePHslice->SetStats(kFALSE);
3324   rangePHslice->Draw();
3325   
3326   TF1* funcPol1 = new TF1("funcPol1", "[0]+[1]*x", 2.9, 6.4);
3327   
3328   TH2F* h2F;
3329   TH1D* hF;
3330   if((h2F = dynamic_cast<TH2F*>(fHistos->At(kPHSlice+centralityClass)))) {
3331     if(h2F && h2F->GetEntries()>10) {
3332       hF = Proj2D(h2F);
3333       h2F->SetStats(kFALSE);
3334       h2F->Draw("samecolz");
3335       isGoodForSaving=kTRUE;
3336       if(trendValues) {
3337         hF->Fit(funcPol1, "QME0", "goff", 2.9, 6.4);
3338         trendValues[6] = funcPol1->GetParameter(0);
3339         trendValues[7] = funcPol1->GetParError(0);
3340         trendValues[8] = funcPol1->GetParameter(1);
3341         trendValues[9] = funcPol1->GetParError(1);
3342       }
3343       hF->SetLineWidth(2);
3344       hF->Draw("same");
3345     }
3346   }
3347   delete funcPol1;
3348
3349   // Qtot vs P
3350   pad = ((TVirtualPad*)l->At(5)); pad->cd();
3351   pad->SetLeftMargin(0.15); pad->SetRightMargin(0.1);
3352   pad->SetTopMargin(0.1); pad->SetBottomMargin(0.15);
3353   pad->SetGridx(kFALSE); pad->SetGridy(kFALSE);
3354   pad->SetLogz();
3355   if(gROOT->FindObject("rangeQtotP")) delete gROOT->FindObject("rangeQtotP");
3356   TH2F* rangeQtotP = new TH2F("rangeQtotP", "", 10, 0.0, 11.99, 10, 0.0, 11.99);
3357   SetStyle(rangeQtotP->GetXaxis(), "P [GeV/c]", 0.07, 0.8, kTRUE, 0.05);
3358   SetStyle(rangeQtotP->GetYaxis(), "Q_{tot}", 0.07, 0.8, kTRUE, 0.05);
3359   rangeQtotP->SetStats(kFALSE);
3360   rangeQtotP->Draw();
3361   
3362   if((h2F = dynamic_cast<TH2F*>(fHistos->At(kQtotP+centralityClass)))) {
3363     if(h2F && h2F->GetEntries()>10) {
3364       h2F->SetStats(kFALSE);
3365       h2F->Draw("samecolz");
3366       isGoodForSaving=kTRUE;