]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGPP/TRD/AliTRDefficiency.cxx
Escape special characters for latex (Diego)
[u/mrichter/AliRoot.git] / PWGPP / TRD / AliTRDefficiency.cxx
CommitLineData
1ee39b3a 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/* $Id: AliTRDefficiency.cxx 27496 2008-07-22 08:35:45Z cblume $ */
17
18////////////////////////////////////////////////////////////////////////////
19// //
20// Reconstruction QA //
21// //
22// Authors: //
23// Markus Fasel <M.Fasel@gsi.de> //
24// //
25////////////////////////////////////////////////////////////////////////////
26
eb05d549 27#include <TROOT.h>
668a0654 28#include <TFile.h>
eb05d549 29#include <TStyle.h>
1ee39b3a 30#include <TClonesArray.h>
31#include <TObjArray.h>
32#include <TProfile.h>
e2e3cec2 33#include <TPad.h>
eb05d549 34#include <TCanvas.h>
a310e49b 35#include <TLegend.h>
eb05d549 36#include <THnSparse.h>
37#include <TH2.h>
38#include <TH3.h>
39#include <THStack.h>
1ee39b3a 40#include "TTreeStream.h"
95d47440 41#include "TPDGCode.h"
1ee39b3a 42
43#include "AliPID.h"
44#include "AliESDtrack.h"
45#include "AliTrackReference.h"
46#include "AliExternalTrackParam.h"
47#include "AliTracker.h"
48#include "AliAnalysisManager.h"
49
eb05d549 50#include "AliTRDgeometry.h"
51#include "AliTRDtrackV1.h"
1ee39b3a 52#include "Cal/AliTRDCalPID.h"
53#include "AliTRDefficiency.h"
54#include "info/AliTRDtrackInfo.h"
55
56ClassImp(AliTRDefficiency)
57
58//____________________________________________________________________
59AliTRDefficiency::AliTRDefficiency()
f2e89a4c 60 :AliTRDrecoTask()
db99a57a 61 ,fMissed(NULL)
eb05d549 62 ,fProj(NULL)
f8f46e4d 63{
64 //
65 // Default constructor
66 //
f2e89a4c 67 SetNameTitle("TRDefficiency", "TRD barrel tracking efficiency checker");
f8f46e4d 68}
69
f2e89a4c 70//____________________________________________________________________
f8f46e4d 71AliTRDefficiency::AliTRDefficiency(char* name)
72 :AliTRDrecoTask(name, "TRD barrel tracking efficiency checker")
db99a57a 73 ,fMissed(NULL)
eb05d549 74 ,fProj(NULL)
1ee39b3a 75{
76 //
77 // Default constructor
78 //
79}
80
81//____________________________________________________________________
82AliTRDefficiency::~AliTRDefficiency()
83{
84 // Destructor
85 if(fMissed){
86 fMissed->Delete();
87 delete fMissed;
88 }
89}
90
eb05d549 91// //____________________________________________________________________
92// void AliTRDefficiency::UserCreateOutputObjects()
93// {
94// //
95// // Create output objects
96// //
97//
98// const Int_t nbins = AliTRDCalPID::kNMom;
99// Float_t xbins[nbins+1] = {.5, .7, .9, 1.3, 1.7, 2.4, 3.5, 4.5, 5.5, 7., 9., 11.};
100//
101// TH1 *h = NULL;
102// fContainer = new TObjArray(); fContainer->SetOwner();
103// for(Int_t is=0; is<AliPID::kSPECIES; is++){
104// fContainer->Add(h = new TProfile(Form("h%s", AliTRDCalPID::GetPartSymb(is)), AliPID::ParticleShortName(is), nbins, xbins));
105// h->SetLineColor(AliTRDCalPID::GetPartColor(is));
106// h->SetMarkerColor(AliTRDCalPID::GetPartColor(is));
107// h->SetMarkerStyle(24);
108// }
109// fContainer->Add(h = new TProfile("h", "", nbins, xbins));
110// h->SetMarkerStyle(7);
111// PostData(1, fContainer);
112// }
113
1ee39b3a 114//____________________________________________________________________
eb05d549 115TH1* AliTRDefficiency::PlotBasicEff(const AliTRDtrackV1 *track)
1ee39b3a 116{
eb05d549 117// plot TRD efficiency based on ESD info
1ee39b3a 118
eb05d549 119 if(!fkESD){
120 AliDebug(4, "No ESD info.");
121 return NULL;
122 }
1ee39b3a 123
eb05d549 124 THnSparse *H(NULL);
125 if(!fContainer || !(H = ((THnSparse*)fContainer->FindObject("hEFF")))){
126 AliWarning("No output container defined.");
127 return NULL;
1ee39b3a 128 }
eb05d549 129 if(track) fkTrack = track;
130
95d47440 131 Double_t val[8]; memset(val, 0, 8*sizeof(Double_t));
eb05d549 132 ULong_t status(fkESD->GetStatus());
95d47440 133// printf(" %3d ITS[%c] TPC[%c] TRDin[%c] TRDout[%c] TRDStop[%c]\n", fkESD->GetId(),
134// (status&AliESDtrack::kITSout)?'y':'n',
135// (status&AliESDtrack::kTPCout)?'y':'n',
136// (status&AliESDtrack::kTRDin)?'y':'n',
137// (status&AliESDtrack::kTRDout)?'y':'n',
138// (status&AliESDtrack::kTRDStop)?'y':'n');
139
140
eb05d549 141 val[0] =((status&AliESDtrack::kTRDin)?1:0) +
142 ((status&AliESDtrack::kTRDStop)?1:0) +
143 ((status&AliESDtrack::kTRDout)?2:0);
ef09211d 144 val[1] = fkESD->Phi();
145 val[2] = fkESD->Eta();
146 val[3] = GetPtBin(fPt>0.?fPt:fkESD->Pt());
95d47440 147 val[4] = -2.; val[5] = fkMC?fkMC->GetPt():-1.;
148 if(fkMC && fkMC->GetNTrackRefs()>=2){ // TODO define trackable tracks more exactly
149 val[4] = -1.;
150 if(fkTrack && fkTrack->GetNumberOfTracklets()>0){ // TODO better define
151 if(fkMC->GetLabel() == fkMC->GetTRDlabel()) val[4] = 0.; // good tracks
152 else if(fkMC->GetLabel() == -fkMC->GetTRDlabel()) val[4] = 1.; // acceptable tracks
153 else val[4] = 2.; // wrongly matched tracks
154 }
155 val[5] = GetPtBin(fkMC->GetTrackRef()->Pt());
eb05d549 156 }
95d47440 157
158 val[6] = fkTrack?fkTrack->GetNumberOfTracklets():0;
668a0654 159 // down scale PID resolution
0f2c4c4f 160 Int_t spc(fSpecies); if(spc==-6) spc=0; if(spc==3) spc=2; if(spc==-3) spc=-2; if(spc>3) spc=3; if(spc<-3) spc=-3;
95d47440 161 val[7] = spc;
162 if(fkMC){
163 Int_t exactPID = 0;
164 switch(fkMC->GetPDG()){
165 case kElectron: exactPID = -1;break;
166 case kPositron: exactPID = 1;break;
167 case kPiPlus: exactPID = 2;break;
168 case kPiMinus: exactPID = -2;break;
169 case kProton: exactPID = 3;break;
170 case kProtonBar: exactPID = -3;break;
171 }
172 val[7] = exactPID;
173 }
174// printf("%3d label[%2d %+2d] PDG[%+4d] TR[%d] Trklts[%d] species[%d %d] pt%d[%6.4f] MC_pt%d[%6.4f] in[%c] out[%c] stop[%c] flag[%d]\n",
175// fkESD->GetId(), fkMC->GetLabel(), fkMC->GetTRDlabel(), fkMC->GetPDG(), fkMC->GetNTrackRefs(), Int_t(val[6]),
176// fSpecies, spc, Int_t(val[3]), fPt>0.?fPt:fkESD->Pt(), Int_t(val[5]), fkMC->GetTrackRef()?fkMC->GetTrackRef()->Pt():fkMC->GetPt(),
177// (status&AliESDtrack::kTRDin)?'x':'o', (status&AliESDtrack::kTRDout)?'x':'o', (status&AliESDtrack::kTRDStop)?'x':'o', Int_t(val[4]));
0f2c4c4f 178
95d47440 179/* if(val[0]<0.5 && TMath::Abs(val[4]+1)<0.1){
180
181 if(fkTrack) printf(" Seed[%s] tracklets[%d]\n", fkTrack->IsTPCseeded()?"TPC":"ITS", Int_t(val[5]));
182 }*/
0f2c4c4f 183 if(fkTrack) AliDebug(2, Form("%3d[%s] tracklets[%d] label[%2d %+2d] species[%d %d] in[%c] out[%c] stop[%c]",
184 fkESD->GetId(), fkTrack->IsTPCseeded()?"TPC":"ITS", Int_t(val[5]), fkMC->GetLabel(), fkMC->GetTRDlabel(), fSpecies, spc,
185 (status&AliESDtrack::kTRDin)?'x':'o', (status&AliESDtrack::kTRDout)?'x':'o', (status&AliESDtrack::kTRDStop)?'x':'o'));
eb05d549 186 H->Fill(val);
187 return NULL;
188}
189
190// //____________________________________________________________________
191// TH1* AliTRDefficiency::PlotMC(const AliTRDtrackV1 *track)
192// {
193// // plot TRD efficiency based on MC info
194//
195// if(!HasMC()) return NULL;
196// if(!fkESD){
197// AliDebug(4, "No ESD info.");
198// return NULL;
199// }
200// if(!fkMC){
201// AliDebug(4, "No MC info.");
202// return NULL;
203// }
204//
205// THnSparse *H(NULL);
206// if(!fContainer || !(H = ((THnSparse*)fContainer->FindObject("hMC")))){
207// AliWarning("No output container defined.");
208// return NULL;
209// }
210// if(track) fkTrack = track;
211// Double_t val[11]; memset(val, 0, 11*sizeof(Double_t));
212// ULong_t status(fkESD->GetStatus());
213// val[0] =((status&AliESDtrack::kTRDin)?1:0) +
214// ((status&AliESDtrack::kTRDStop)?1:0) +
215// ((status&AliESDtrack::kTRDout)?2:0);
216// val[1] = fkESD->Phi();
217// val[2] = fkESD->Eta();
218// val[3] = DebugLevel()>=1?GetPtBin(fkESD->Pt()):GetPtBinSignificant(fkESD->Pt());
219// if(fkTrack){ // read track status in debug mode with friends
220// val[4] = fkTrack->GetStatusTRD(-1);
221// for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++) val[5+ily]=fkTrack->GetStatusTRD(ily);
222// }
223// H->Fill(val);
224//
225// }
1ee39b3a 226
227//____________________________________________________________________
eb05d549 228void AliTRDefficiency::LocalUserExec(Option_t *)
1ee39b3a 229{
230 //
eb05d549 231 // Do it obsolete
1ee39b3a 232 //
233
eb05d549 234 Int_t labelsacc[10000];
1ee39b3a 235 memset(labelsacc, 0, sizeof(Int_t) * 10000);
236
5935a6da 237 fTracks = dynamic_cast<TObjArray *>(GetInputData(1));
238 if(!fTracks) return;
239 if(!fTracks->GetEntriesFast()) return;
240 else AliDebug(2, Form("Tracks[%d] for %s", fTracks->GetEntriesFast(), GetName()));
1ee39b3a 241 if(!fMissed){
242 fMissed = new TClonesArray("AliTRDtrackInfo", 10);
243 fMissed->SetOwner();
244 }
245
246 Float_t mom;
247 Int_t selection[10000], nselect = 0;
248 ULong_t status; Int_t pidx;
249 Int_t nTRD = 0, nTPC = 0, nMiss = 0;
db99a57a 250 AliTRDtrackInfo *track = NULL;
251 AliTrackReference *ref = NULL;
252 AliExternalTrackParam *esd = NULL;
1ee39b3a 253 for(Int_t itrk=0; itrk<fTracks->GetEntriesFast(); itrk++){
254 track = (AliTRDtrackInfo*)fTracks->UncheckedAt(itrk);
255
256 if(!track->HasESDtrack()) continue;
257 status = track->GetStatus();
258
259 // missing TPC propagation - interesting for SA
260 if(!(status&AliESDtrack::kTPCout)) continue;
261
262 // missing MC info.
263 if(HasMCdata() && track->GetNTrackRefs() <= 1) continue;
264
265 nTPC++;
266 selection[nselect++]=itrk;
267 ref = track->GetTrackRef(0);
268 esd = track->GetESDinfo()->GetOuterParam();
269 mom = ref ? ref->P(): esd->P();
270 pidx = AliTRDCalPID::GetPartIndex(track->GetPDG());
271 pidx = TMath::Max(pidx, 0);
8f7b1226 272 AliDebug(4, Form("PID: %d", pidx));
1ee39b3a 273
274 //Int_t n = track->GetNumberOfClusters();
275 // where are this tracklets ???
276 //if(ncls0 > ncls1) printf("%3d ESD[%3d] TRD[%3d|%3d]\n", itrk, ncls0, ncls1, n);
277 if(track->GetNumberOfClustersRefit()){
278 ((TProfile*)fContainer->At(pidx))->Fill(mom, 1.);
279 labelsacc[nTRD] = track->GetLabel();
280 nTRD++;
281 continue;
282 }
283
284
285
286 Float_t xmed, xleng;
287 Int_t iref = 1; Bool_t found = kFALSE;
288 while((ref = track->GetTrackRef(iref))){
289 xmed = .5*(ref->LocalX() + track->GetTrackRef(iref-1)->LocalX());
290 xleng= (ref->LocalX() - track->GetTrackRef(iref-1)->LocalX());
291 if(TMath::Abs(xmed - 298.5) < .5 &&
292 TMath::Abs(xleng - 3.7) < .1){
293 found = kTRUE;
294 break;
295 }
296 iref++;
297 }
298 if(!found){
299 nTPC--;
300 // track missing first layer. Maybe interesting for SA.
301 continue;
302 }
303 nselect--;
304 new ((*fMissed)[nMiss]) AliTRDtrackInfo(*track);
305 nMiss++;
306 }
307 AliDebug(2, Form("%3d Tracks: ESD[%3d] TPC[%3d] TRD[%3d | %5.2f%%] Off[%d]", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), fTracks->GetEntriesFast(), nTPC, nTRD, nTPC ? 1.E2*nTRD/float(nTPC) : 0., fMissed->GetEntriesFast()));
308
309
310 // Find double tracks
311 Float_t threshold = 10.;
db99a57a 312 AliTrackReference *refMiss = NULL;
313 AliExternalTrackParam *op = NULL;
314 AliTRDtrackInfo *tt = NULL;
1ee39b3a 315 for(Int_t imiss=0; imiss<nMiss; imiss++){
316 //printf("Searching missing %d ...\n", imiss);
317
318 // get outer param of missed
319 tt = (AliTRDtrackInfo*)fMissed->UncheckedAt(imiss);
320 op = tt->GetESDinfo()->GetOuterParam();
321 Double_t alpha = op->GetAlpha(), cosa = TMath::Cos(alpha), sina = TMath::Sin(alpha);
322
323 Double_t xyz[3], x0, y0, z0, x, y, z, dx, dy, dz, d;
324
325 Bool_t bFOUND = kFALSE;
326 for(Int_t iselect=0; iselect<nselect; iselect++){
327 track = (AliTRDtrackInfo*)fTracks->UncheckedAt(selection[iselect]);
328
329 // check first MC ... if available
330 d = 0;
331 for(Int_t iref=0; iref<track->GetNTrackRefs(); iref++){
332 if(!(ref = track->GetTrackRef(iref))) continue;
333 if((refMiss = tt->GetTrackRef(iref))){
334 dy = ref->LocalY() - refMiss->LocalY();
335 dz = ref->Z() - refMiss->Z();
336 } else {
337 // compare missOP with refTrackRef in LTC
338 x0 = ref->LocalX();
339 op->GetYAt(x0, AliTracker::GetBz(), y0);
340 op->GetZAt(x0, AliTracker::GetBz(), z0);
341 dy = y0 - ref->LocalY();
342 dz = z0 - ref->Z();
343 }
344 d += (dy*dy + dz*dz);
345 }
346 //printf("\td[%d] = %f N[%d]\n", selection[iselect], d, track->GetNTrackRefs());
347 if((track->GetNTrackRefs())){
348 d /= track->GetNTrackRefs();
349 if(d < threshold){
350 //printf("\t\tFound %2d in ref[%3d] : d[%f]\n", imiss, selection[iselect], d/track->GetNTrackRefs());
351 bFOUND = kTRUE; break;
352 }
353 }
354
355 // process outer param ... always available
356 // compare missOP with OP in GTC
357 esd = track->GetESDinfo()->GetOuterParam();
358 esd->GetXYZ(xyz);
359 x0 = esd->GetX();
360 op->GetYAt(x0, AliTracker::GetBz(), y0);
361 op->GetZAt(x0, AliTracker::GetBz(), z0);
362 x = x0*cosa - y0*sina;
363 y = x0*sina + y0*cosa;
364 z = z0;
365 dx=xyz[0]-x;
366 dy=xyz[1]-y;
367 dz=xyz[2]-z;
368 d = dx*dx+dy*dy+dz*dz;
369 //printf("\td[%d] = %f op\n", selection[iselect], d);
370 if(d < threshold){
371 //printf("\t\tFound %2d in op[%3d] : d[%f] dx[%5.2f] dy[%5.2f] dz[%5.2f]\n", imiss, selection[iselect], d, dx, dy, dz);
372 bFOUND = kTRUE; break;
373 }
374 }
375 if(bFOUND) nTPC--;
376 else{
377 ref = tt->GetTrackRef(0);
378 mom = ref ? ref->P(): op->P();
379 pidx = AliTRDCalPID::GetPartIndex(tt->GetPDG());
380 pidx = TMath::Max(pidx, 0);
381 ((TProfile*)fContainer->At(pidx))->Fill(mom, 0.);
382 AliDebug(2, Form(" NOT bFOUND Id[%d] Mom[%f]\n", tt->GetTrackId(), mom));
383 }
384 }
385
386 AliDebug(2, Form("%3d Tracks: ESD[%3d] TPC[%3d] TRD[%3d | %5.2f%%] Off[%d]", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), fTracks->GetEntriesFast(), nTPC, nTRD, nTPC ? 1.E2*nTRD/float(nTPC) : 0., fMissed->GetEntriesFast()));
387
388 //fMissed->Delete();
389 // check for double countings
390 Int_t indices[10000]; memset(indices, 0, sizeof(Int_t) * 10000);
391 TMath::Sort(nTRD, labelsacc, indices);
392 if(DebugLevel() > 2){
393 for(Int_t itk = 0; itk < nTRD - 1; itk++)
394 if(labelsacc[indices[itk]] ==labelsacc[indices[itk + 1]]) printf("Double counted MC track: %d\n", labelsacc[indices[itk]]);
395 }
1ee39b3a 396}
397
1ee39b3a 398//____________________________________________________________________
399Bool_t AliTRDefficiency::GetRefFigure(Int_t ifig)
400{
401// Steer reference figures
402
e2e3cec2 403 if(!gPad){
404 AliWarning("Please provide a canvas to draw results.");
405 return kFALSE;
406 }
407 gPad->SetLogx();
408
a310e49b 409 TLegend *leg(NULL);
410 Bool_t bFIRST(kTRUE);
411 TProfile *h(NULL);
1ee39b3a 412 switch(ifig){
413 case 0:
414 h = (TProfile*)fContainer->At(AliPID::kSPECIES);
415 for(Int_t is=0; is<AliPID::kSPECIES; is++){
416 h->Add((TProfile*)fContainer->At(is));
417 }
a310e49b 418 h->SetMarkerStyle(24);
1ee39b3a 419 h->SetMarkerColor(kBlack);
420 h->SetLineColor(kBlack);
e2e3cec2 421 h->SetTitle("TRD Efficiency integrated");
a310e49b 422 h->SetXTitle("p [GeV/c]");
e2e3cec2 423 h->GetXaxis()->SetMoreLogLabels();
a310e49b 424 h->SetYTitle("Efficiency");
425 h->GetYaxis()->CenterTitle();
1ee39b3a 426 h->Draw("e1");
427 break;
428 case 1:
429 bFIRST = kTRUE;
430 for(Int_t is=0; is<AliPID::kSPECIES; is++){
431 if(!(h = (TProfile*)fContainer->At(is))) continue;
a310e49b 432 h->SetMarkerStyle(24);
1ee39b3a 433 if(bFIRST){
434 h->Draw("e1");
a310e49b 435 h->SetXTitle("p [GeV/c]");
e2e3cec2 436 h->GetXaxis()->SetMoreLogLabels();
a310e49b 437 h->SetYTitle("Efficiency");
438 h->GetYaxis()->CenterTitle();
439 h->GetYaxis()->SetRangeUser(0.8, 1.05);
440 leg=new TLegend(.7, .2, .98, .6);
441 leg->SetHeader("Species");
442 leg->SetBorderSize(0);
443 leg->SetFillStyle(0);
444 leg->AddEntry(h, h->GetTitle(), "pl");
445 } else {
446 leg->AddEntry(h, h->GetTitle(), "pl");
447 h->Draw("same e1");
448 }
1ee39b3a 449 bFIRST = kFALSE;
450 }
7fe4e88b 451 if(leg) leg->Draw();
1ee39b3a 452 break;
453 }
454 return kTRUE;
455}
456
eb05d549 457//________________________________________________________
458TObjArray* AliTRDefficiency::Histos()
459{
460 //
461 // Define histograms
462 //
463
464 if(fContainer) return fContainer;
465
466 fContainer = new TObjArray(1); fContainer->SetOwner(kTRUE);
467 THnSparse *H(NULL);
468 TString st;
469
470 //++++++++++++++++++++++
471 // cluster to detector
472 if(!(H = (THnSparseI*)gROOT->FindObject("hEFF"))){
95d47440 473 const Int_t mdim(8);
474 Int_t nlabel(3);
475 const Char_t *eTitle[mdim] = {"status", "#phi [rad]", "eta", "p_{t} [bin]", "label", "p_{t}^{MC} [bin]", "N_{trklt}", "chg*spec*rc"};
476 const Int_t eNbins[mdim] = { 5, 144, 45, fNpt-1, nlabel, fNpt-1, AliTRDgeometry::kNlayer-2, 5};
477 const Double_t eMin[mdim] = { -0.5, -TMath::Pi(), -.9, -0.5, -1.5, -0.5, 1.5, -2.5},
478 eMax[mdim] = { 4.5, TMath::Pi(), .9, fNpt-1.5, 1.5, fNpt-1.5, 5.5, 2.5};
eb05d549 479 st = "basic efficiency;";
480 // define minimum info to be saved in non debug mode
95d47440 481 Int_t ndim=DebugLevel()>=1?mdim:(HasMCdata()?6:4);
eb05d549 482 for(Int_t idim(0); idim<ndim; idim++){ st += eTitle[idim]; st+=";";}
483 H = new THnSparseI("hEFF", st.Data(), ndim, eNbins, eMin, eMax);
eb05d549 484 } else H->Reset();
485 fContainer->AddAt(H, 0);
486
487 return fContainer;
488}
1ee39b3a 489
490//____________________________________________________________________
491Bool_t AliTRDefficiency::PostProcess()
492{
eb05d549 493// Fit, Project, Combine, Extract values from the containers filled during execution
494
495 if (!fContainer) {
496 AliError("ERROR: list not available");
497 return kFALSE;
498 }
499 if(!fProj){
500 AliInfo("Building array of projections ...");
95d47440 501 fProj = new TObjArray(2000); fProj->SetOwner(kTRUE);
eb05d549 502 }
0f2c4c4f 503 // set pt/p segmentation. guess from data
504 THnSparse *H(NULL);
505 if(!(H = (THnSparse*)fContainer->FindObject("hEFF"))){
506 AliError("Missing/Wrong data @ hEFF.");
507 return kFALSE;
508 }
509 fNpt=H->GetAxis(3)->GetNbins()+1;
510 if(!MakeMomSegmentation()) return kFALSE;
511
eb05d549 512 if(!MakeProjectionBasicEff()) return kFALSE;
513 return kTRUE;
514}
515
516//____________________________________________________________________
517Bool_t AliTRDefficiency::MakeProjectionBasicEff()
518{
519// Make basic efficiency plots
520
521 if(!fContainer || !fProj){
522 AliError("Missing data container.");
523 return kFALSE;
524 }
525 THnSparse *H(NULL);
526 if(!(H = (THnSparse*)fContainer->FindObject("hEFF"))){
527 AliError("Missing/Wrong data @ hEFF.");
528 return kFALSE;
529 }
530 Int_t ndim(H->GetNdimensions()); //Bool_t debug(ndim>Int_t(kNdimCl));
95d47440 531 TAxis *aa[11], *al(NULL), *apt(NULL); memset(aa, 0, sizeof(TAxis*) * 11);
eb05d549 532 for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
533 if(H->GetNdimensions() > 4) al = H->GetAxis(4);
95d47440 534 if(H->GetNdimensions() > 5) apt = H->GetAxis(5);
535 Int_t nlab=al?5:1;
eb05d549 536
537 // define rebinning strategy
538 //const Int_t nEtaPhi(4); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 5, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 1, 5};
95d47440 539 const Int_t nprojs(50);
540 AliTRDrecoProjection hp[nprojs]; TObjArray php(nprojs);
eb05d549 541 const Char_t *stat[] = {"!TRDin", "TRDin", "TRDin&TRDStop", "TRDin&TRDout", "TRDin&TRDout&TRDStop"};
95d47440 542 const Char_t *lab[] = {"MC-Miss", "MC-Trkble", "MC-Good", "MC-Accept", "MC-Wrong"};
eb05d549 543 Int_t ih(0);
544 for(Int_t ilab(0); ilab<nlab; ilab++){
545 for(Int_t istat(0); istat<5; istat++){
95d47440 546 hp[ih].Build(Form("HEff0%d%d", ilab, istat),
0f2c4c4f 547 Form("Efficiency :: Stat[#bf{%s}] %s", stat[istat], nlab>1?lab[ilab]:""), 2, 1, 3, aa);
eb05d549 548 //hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
95d47440 549 php.AddLast(&hp[ih++]);
550 if(!apt) continue;
551 hp[ih].Build(Form("HEff1%d%d", ilab, istat),
552 Form("Efficiency [MC] :: Stat[#bf{%s}] %s", stat[istat], lab[ilab]), 2, 1, 5, aa);
553 //hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
554 php.AddLast(&hp[ih++]);
eb05d549 555 }
556 }
557 AliInfo(Form("Build %3d 3D projections.", ih));
558
668a0654 559 AliTRDrecoProjection *pr0(NULL), *pr1(NULL);
eb05d549 560 Int_t istatus, ilab(0), coord[11]; memset(coord, 0, sizeof(Int_t) * 11); Double_t v = 0.;
561 for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
562 v = H->GetBinContent(ib, coord); if(v<1.) continue;
95d47440 563 istatus = coord[0]-1; if(istatus>4) continue;
564 ilab=al?coord[4]:0;
565 Int_t isel = (apt?2:1)*(ilab*5+istatus);
668a0654 566 if(isel>=ih){
95d47440 567 AliError(Form("Wrong selection %d [%3d] {stat[%d] lab[%d]}", isel, ih, istatus, ilab));
668a0654 568 return kFALSE;
569 }
570 if(!(pr0=(AliTRDrecoProjection*)php.At(isel))) {
571 AliError(Form("Missing projection %d", isel));
572 return kFALSE;
573 }
95d47440 574 if(strcmp(pr0->H()->GetName(), Form("HEff0%d%d", ilab, istatus))!=0){
575 AliError(Form("Projection mismatch :: request[HEff0%d%d] found[%s]", ilab, istatus, pr0->H()->GetName()));
668a0654 576 return kFALSE;
577 }
95d47440 578 //for(Int_t jh(0); jh<1/*np[isel]*/; jh++) ((AliTRDrecoProjection*)php.At(isel+jh))->Increment(coord, v);
579 AliDebug(2, Form("Found %s for selection stat[%d] lab[%d]", pr0->H()->GetName(), istatus, ilab));
580 ((AliTRDrecoProjection*)php.At(isel))->Increment(coord, v);
581 if(apt) ((AliTRDrecoProjection*)php.At(isel+1))->Increment(coord, v);
eb05d549 582 }
668a0654 583 if(HasDump3D()){
584 TDirectory *cwd = gDirectory;
95d47440 585 TFile::Open(Form("Dump%s_%s.root", GetName(), H->GetName()), "RECREATE");
668a0654 586 for(Int_t ip(0); ip<php.GetEntriesFast(); ip++){
587 if(!(pr0 = (AliTRDrecoProjection*)php.At(ip))) continue;
588 if(!pr0->H()) continue;
589 TH3 *h3=(TH3*)pr0->H()->Clone();
590 h3->Write();
591 }
592 gFile->Close();
593 cwd->cd();
594 }
595
0f2c4c4f 596 Int_t jh(0);
eb05d549 597 for(; ih--; ){
598 if(!hp[ih].H()) continue;
95d47440 599 if(hp[ih].H()->Integral()<5.) continue;
5047978d 600 for(Int_t ipt(0); ipt<=fNpt; ipt++) fProj->AddAt(hp[ih].Projection2Dbin(ipt), jh++);
eb05d549 601 }
602
95d47440 603 // PROCESS MC LABEL, CLONE
604 if(nlab>1){
605 // MISSED TRACKS
606 for(Int_t imc(0); imc<(apt?2:1); imc++){
607 AliTRDrecoProjection prMis, prMisOK, prS, prSs, prTRDf;
608 for(Int_t istat(0); istat<5; istat++) {
609 if(!(pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HEff%d0%d", imc, istat)))) continue;
610 if(!prMis.H()){
611 prMis = (*pr0);
612 prMis.SetNameTitle(Form("HEff%dM", imc), "Missed tracks");
613 prMis.H()->SetNameTitle(Form("HEff%dM", imc), Form("Efficiency [%s] :: Missed tracks", imc?"MC":"REC"));
614 } else prMis+=(*pr0);
615 if(istat==0) { // MC miss & !TRDin
616 prMisOK = (*pr0);
617 prMisOK.SetNameTitle(Form("HEff%dMok", imc), "Missed tracks/Not propagated");
618 prMisOK.H()->SetNameTitle(Form("HEff%dMok", imc), Form("Efficiency [%s] :: Missed tracks/Not propagated", imc?"MC":"REC"));
619 } else if(istat==1) { // MC miss & TRDin
620 prS = (*pr0);
621 prS.SetNameTitle(Form("HEff%dMS", imc), "Missed seeded tracks");
622 prS.H()->SetNameTitle(Form("HEff%dMS", imc), Form("Efficiency [%s] :: Missed tracks/Seeder propagated", imc?"MC":"REC"));
623 } else if(istat==2) { // MC miss & TRDstop
624 prSs = (*pr0);
625 prSs.SetNameTitle(Form("HEff%dMSs", imc), "Stop seeded tracks");
626 prSs.H()->SetNameTitle(Form("HEff%dMSs", imc), Form("Efficiency [%s] :: Missed tracks/Stop seeded tracks", imc?"MC":"REC"));
627 } else if(istat==3 || istat==4) { // MC miss & (TRDout || TRDout+TRDstop)
628 prTRDf = (*pr0);
629 prTRDf.SetNameTitle(Form("HEff%dMTtrd", imc), "Missed tracks/TRD fooled");
630 prTRDf.H()->SetNameTitle(Form("HEff%dMTtrd", imc), Form("Efficiency [%s] :: Missed tracks/TRD fooled", imc?"MC":"REC"));
631 } else {
632 if(Int_t(pr0->H()->Integral()) > 0) AliWarning(Form("Detected %3d MISSED entries for %s[%s]", Int_t(pr0->H()->Integral()), stat[istat], imc?"MC":"REC"));
633 }
634 }
635 for(Int_t ipt(-1); ipt<=fNpt; ipt++){
636 fProj->AddAt(prMis.Projection2Dbin(ipt, imc), jh++);
637 fProj->AddAt(prMisOK.Projection2Dbin(ipt, imc), jh++);
638 fProj->AddAt(prS.Projection2Dbin(ipt, imc), jh++);
639 fProj->AddAt(prSs.Projection2Dbin(ipt, imc), jh++);
640 fProj->AddAt(prTRDf.Projection2Dbin(ipt, imc), jh++);
641 }
642
643 // TRACKABLE TRACKS
644 AliTRDrecoProjection prAll, prTRDns, prTRDfail, prTRD, prTRDbad;//, prTRDstop;
645 for(ilab=1; ilab<=4; ilab++){
646 for(Int_t istat(0); istat<5; istat++) {
647 if(!(pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HEff%d%d%d", imc, ilab, istat)))){
648 AliError(Form("Missing prj. HEff%d%d%d", imc, ilab, istat));
649 continue;
650 }
651 if(!prAll.H()){
652 prAll = (*pr0);
653 prAll.SetNameTitle(Form("HEff%dAT", imc), "Trackable tracks");
654 prAll.H()->SetNameTitle(Form("HEff%dAT", imc), Form("Efficiency [%s] :: Trackable tracks", imc?"MC":"REC"));
655 } else prAll+=(*pr0);
656 if(ilab==1 && istat==0){ // MC trackable & !TRDin
657 prTRDns = (*pr0);
658 prTRDns.SetNameTitle(Form("HEff%dAnStrd", imc), "Trackable tracks/TRD not seeded");
659 prTRDns.H()->SetNameTitle(Form("HEff%dAnStrd", imc), Form("Efficiency [%s] :: Trackable tracks/TRD not seeded", imc?"MC":"REC"));
660 } else if(ilab==1 && istat==1){ // MC trackable & TRDin
661 prTRDfail = (*pr0);
662 prTRDfail.SetNameTitle(Form("HEff%dAFtrd", imc), "Trackable tracks/TRD failed");
663 prTRDfail.H()->SetNameTitle(Form("HEff%dAFtrd", imc), Form("Efficiency [%s] :: Trackable tracks/TRD failed", imc?"MC":"REC"));
664 } else if((ilab==2 && istat>=3) || // MC good & (TRDout || TRDout+TRDstop)
665 (ilab==3 && istat>=3) || // MC accept & (TRDout || TRDout+TRDstop)
666 (ilab==4 && istat==2)) { // MC bad & TRDin+TRDstop
667 if(!prTRD.H()){
668 prTRD = (*pr0);
669 prTRD.SetNameTitle(Form("HEff%dAtrd", imc), "Trackable tracks/TRD OK");
670 prTRD.H()->SetNameTitle(Form("HEff%dAtrd", imc), Form("Efficiency [%s] :: Trackable tracks/TRD OK", imc?"MC":"REC"));
671 } else prTRD+=(*pr0);
672 } else if(ilab==4 && istat>=3){ // MC trackable & (TRDout (+TRDstop))
673 if(!prTRDbad.H()){
674 prTRDbad = (*pr0);
675 prTRDbad.SetNameTitle(Form("HEff%dABtrd", imc), "Trackable tracks/TRD bad");
676 prTRDbad.H()->SetNameTitle(Form("HEff%dABtrd", imc), Form("Efficiency [%s] :: Trackable tracks/TRD bad", imc?"MC":"REC"));
677 } else prTRDbad+=(*pr0);
678 } else {
679 if(Int_t(pr0->H()->Integral()) > 0) AliWarning(Form("Detected %3d TRACKABLE entries for %s -> %s[%s]", Int_t(pr0->H()->Integral()), lab[ilab], stat[istat], imc?"MC":"REC"));
680 }
681 }
682 }
683 for(Int_t ipt(-1); ipt<=fNpt; ipt++){
684 fProj->AddAt(prAll.Projection2Dbin(ipt, imc), jh++);
685 fProj->AddAt(prTRDns.Projection2Dbin(ipt), jh++);
686 fProj->AddAt(prTRDfail.Projection2Dbin(ipt, imc), jh++);
687 fProj->AddAt(prTRD.Projection2Dbin(ipt, imc), jh++);
688 fProj->AddAt(prTRDbad.Projection2Dbin(ipt, imc), jh++);
689 }
690 } // END pt source LOOP
691 } // END PROCESS EFFICIENCY MC
eb05d549 692
95d47440 693 // PROCESS EFFICIENCY NOMC
694 const char suffix[] = {'A', 'S', 'T'};
695 const char *sname[] = {"All", "Stopped", "Tracked [Stop]"};
eb05d549 696 for(Int_t istat(0); istat<5; istat++) {
95d47440 697 if(!(pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HEff00%d", istat)))){
698 AliError(Form("Missing prj. HEff00%d", istat));
699 continue;
700 }
701 // sum over MC labels if available
702 for(ilab=1; ilab<nlab; ilab++){
703 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("HEff0%d%d", ilab, istat)))){
704 AliError(Form("Missing prj. HEff0%d%d", ilab, istat));
705 continue;
eb05d549 706 }
95d47440 707 (*pr0)+=(*pr1);
eb05d549 708 }
95d47440 709 pr0->H()->SetNameTitle(Form("HEff%d", istat), Form("Efficiency :: Stat[#bf{%s}]", stat[istat]));
710 for(Int_t ipt(0); ipt<=fNpt; ipt++) fProj->AddAt(pr0->Projection2Dbin(ipt), jh++);
711
712 if(istat>1 && (pr1 = (AliTRDrecoProjection*)php.FindObject("HEff001"))) (*pr1)+=(*pr0);
713 //if(istat>2 && (pr1 = (AliTRDrecoProjection*)php.FindObject("HEff02"))) (*pr1)+=(*pr0);
714 if(istat>3 && (pr1 = (AliTRDrecoProjection*)php.FindObject("HEff003"))) (*pr1)+=(*pr0);
eb05d549 715 }
0f2c4c4f 716 // Project 2D tracks
0f2c4c4f 717 for(Int_t istat(0); istat<3; istat++){
95d47440 718 if(!(pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HEff00%d", istat+1)))){
719 AliError(Form("Missing prj. HEff00%d", istat+1));
720 continue;
721 }
0f2c4c4f 722 pr0->H()->SetNameTitle(Form("HEff%c", suffix[istat]), Form("Efficiency :: %s Tracks", sname[istat]));
5047978d 723 for(Int_t ipt(-1); ipt<=fNpt; ipt++) fProj->AddAt(pr0->Projection2Dbin(ipt), jh++);
eb05d549 724 }
0f2c4c4f 725
726 // Efficiency
95d47440 727/* TH2F *h2T(NULL), *h2P(NULL);
0f2c4c4f 728 for(Int_t ipt(-1); ipt<=fNpt; ipt++){
729 if(!(h2T = (TH2F*)fProj->FindObject(Form("HEffT%d_2D", ipt)))) continue;
730 if(!(h2P = (TH2F*)fProj->FindObject(Form("HEffP%d_2D", ipt)))) continue;
95d47440 731 h2P->Divide(h2T); PutTrendValue(ipt<0?"pt":(Form("pt%d", ipt)), GetMeanStat(h2P, 0.01, ">"));
0f2c4c4f 732 }*/
eb05d549 733 AliInfo(Form("Done %3d 2D projections.", jh));
1ee39b3a 734 return kTRUE;
735}
eb05d549 736
737//____________________________________________________________________
738void AliTRDefficiency::MakeSummary()
739{
740// Build summary plots
741 if(!fProj){
742 AliError("Missing results");
743 return;
744 }
745 TVirtualPad *p(NULL); TCanvas *cOut(NULL);
746 TH2 *h2(NULL);
747 gStyle->SetPalette(1);
748
4b92200c 749 Int_t nbins(DebugLevel()==0?3:20);
eb05d549 750 //calculate true pt bin
751 Float_t ptBins[23]; ptBins[0] = 0.;
752 if(nbins==3){ // significant bins
753 ptBins[1] = 0.5;
754 ptBins[2] = 0.8;
755 ptBins[3] = 1.5;
756 ptBins[4] = 5.;
757 ptBins[5] = 10.;
758 } else if(nbins==20){ // logarithmic bins
759 ptBins[1] = 0.5;
760 Float_t dpt(0.002);
761 for(Int_t ipt(1); ipt<21; ipt++) ptBins[ipt+1] = ptBins[ipt]+(TMath::Exp(ipt*ipt*dpt)-1.);
762 ptBins[22] = 10.;
763 } else {
764 AliError(Form("Unknown no.[%d] of bins in the p_t spectrum", nbins));
765 return;// kFALSE;
766 }
767
33056e04 768 const Char_t cid[]={'T','P'};
eb05d549 769 cOut = new TCanvas(Form("%s_Eff", GetName()), "TRD Efficiency", 1536, 1536); cOut->Divide(2,2,1.e-5,1.e-5);
770 // tracking eff :: eta-phi dependence
771 for(Int_t it(0); it<2; it++){
668a0654 772 if(!(h2 = (TH2*)fProj->FindObject(Form("HEff%cEnN", cid[it])))) {
773 AliError(Form("Missing \"HEff%cEnN\".", cid[it]));
eb05d549 774 continue;
775 }
33056e04 776 h2->SetContour(10); h2->Scale(1.e2); SetRangeZ(h2, 80, 100, 5);
eb05d549 777 h2->GetZaxis()->SetTitle("Efficiency [%]"); h2->GetZaxis()->CenterTitle();
778 p=cOut->cd(it+1);p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712);
779 h2->Draw("colz");
33056e04 780 MakeDetectorPlot();
eb05d549 781 }
782 if(!(h2 = (TH2*)fProj->FindObject("HEff0En"))) {
783 AliError("Missing \"HEff0En\".");
784 return;
785 }
786 p=cOut->cd(3);p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712);
33056e04 787 h2->Draw("colz"); MakeDetectorPlot();
eb05d549 788 // tracking eff :: pt dependence
789 TH1 *h[2] = {0};
668a0654 790 if(!(h[0] = (TH1*)fProj->FindObject("HEffP_zN"))){
791 AliError("Missing \"HEffP_zN\".");
eb05d549 792 return;
793 }
668a0654 794 if(!(h[1] = (TH1*)fProj->FindObject("HEffT_zN"))){
795 AliError("Missing \"HEffT_zN\".");
eb05d549 796 return;
797 }
798 TH1 *h1[3] = {0};
799 Color_t color[] = {kGreen, kBlue, kRed};
800 for(Int_t il=0;il<3;il++){
801 h1[il]=new TH1F(Form("h1Eff%d", il), "", nbins+2, ptBins);
802 h1[il]->SetFillColor(color[il]);
803 h1[il]->SetFillStyle(il==2?3002:1001);
804 h1[il]->SetLineColor(color[il]);
668a0654 805 h1[il]->SetMarkerStyle(4);
806 h1[il]->SetMarkerColor(color[il]);
eb05d549 807 h1[il]->SetLineWidth(1);
808 }
809 for(Int_t ip(0);ip<=(nbins+1);ip++){
810 h1[0]->SetBinContent(ip+1, 1.e2*h[0]->GetBinContent(ip)); // propagated
811 h1[1]->SetBinContent(ip+1, 1.e2*(h[1]->GetBinContent(ip) - h[0]->GetBinContent(ip))); // stopped
812 h1[2]->SetBinContent(ip+1, 1.e2*(1 - h[1]->GetBinContent(ip))); // missed
813 }
33056e04 814 const Char_t *labEff[] = {"Propagated", "Stopped", "Missed"};
eb05d549 815 THStack *hs = new THStack("hEff","Tracking Efficiency;p_{t} [GeV/c];Efficiency [%]");
816 TLegend *leg = new TLegend(0.671371,0.1313559,0.9576613,0.2923729,NULL,"brNDC");
817 leg->SetBorderSize(0); leg->SetFillColor(kWhite); leg->SetFillStyle(1001);
818 for(Int_t ic(0); ic<3;ic++){ hs->Add(h1[ic]);leg->AddEntry(h1[ic], labEff[ic], "f");}
819 p=cOut->cd(4); p->SetLeftMargin(0.08266129); p->SetRightMargin(0.0141129);p->SetTopMargin(0.006355932);p->SetLogx();
668a0654 820 hs->Draw(); //hs->Draw("same nostack,e1p");
821 leg->Draw();
eb05d549 822 hs->GetXaxis()->SetRangeUser(0.6,10.);
823 hs->GetXaxis()->SetMoreLogLabels();
824 hs->GetXaxis()->SetTitleOffset(1.2);
825 hs->GetYaxis()->SetNdivisions(513);
668a0654 826 hs->SetMinimum(75.);
eb05d549 827 hs->GetYaxis()->CenterTitle();
828 cOut->SaveAs(Form("%s.gif", cOut->GetName()));
829
830 cOut = new TCanvas(Form("%s_MC", GetName()), "TRD Label", 1536, 1536); cOut->Divide(2,2,1.e-5,1.e-5);
831 for(Int_t ipad(0); ipad<3; ipad++){
832 p=cOut->cd(ipad+1);p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712);
668a0654 833 if(!(h2 = (TH2*)fProj->FindObject(Form("HEffLb%dEnN", ipad)))) continue;
eb05d549 834 h2->SetContour(10);
668a0654 835 h2->Scale(1.e2); SetRangeZ(h2, ipad==1?50:0., ipad==1?90.:50., ipad==1?0.01:0.01);
eb05d549 836 h2->GetZaxis()->SetTitle("Efficiency [%]"); h2->GetZaxis()->CenterTitle();
837 h2->Draw("colz");
33056e04 838 MakeDetectorPlot();
eb05d549 839 }
33056e04 840 color[0] = kRed; color[1] = kGreen; color[2] = kBlue;
eb05d549 841 for(Int_t il=0;il<3;il++){
668a0654 842 if(!(h[il] = (TH1D*)fProj->FindObject(Form("HEffLb%d_zN", il)))) continue;
eb05d549 843 h1[il]=new TH1F(Form("h1Lab%d", il), "", nbins+2, ptBins);
844 for(Int_t ip(0);ip<=(nbins+1);ip++) h1[il]->SetBinContent(ip+1, 1.e2*h[il]->GetBinContent(ip));
33056e04 845 h1[il]->SetFillColor(color[il]);
eb05d549 846 h1[il]->SetFillStyle(il==2?3002:1001);
33056e04 847 h1[il]->SetLineColor(color[il]);
eb05d549 848 h1[il]->SetLineWidth(1);
849 }
33056e04 850 const Char_t *labMC[] = {"TRD != ESD [bad]", "TRD == ESD [good]", "TRD == -ESD [accept]"};
eb05d549 851 leg = new TLegend(0.671371,0.1313559,0.9576613,0.2923729,NULL,"brNDC");
852 leg->SetBorderSize(0); leg->SetFillColor(kWhite); leg->SetFillStyle(1001);
853 hs = new THStack("hLab","TRD Label;p_{t} [GeV/c];Efficiency [%]");
854 hs->Add(h1[1]);leg->AddEntry(h1[1], labMC[1], "f"); // good
855 hs->Add(h1[2]);leg->AddEntry(h1[2], labMC[2], "f"); // accept
856 hs->Add(h1[0]);leg->AddEntry(h1[0], labMC[0], "f"); // bad
857 p=cOut->cd(4); p->SetLeftMargin(0.08266129); p->SetRightMargin(0.0141129);p->SetTopMargin(0.006355932); p->SetLogx();
668a0654 858 hs->Draw(/*"nostack,e1p"*/); leg->Draw();
eb05d549 859 cOut->Modified();cOut->Update();
860 hs->GetXaxis()->SetRangeUser(0.6,10.);
861 hs->GetXaxis()->SetMoreLogLabels();
862 hs->GetXaxis()->SetTitleOffset(1.2);
863 hs->GetYaxis()->SetNdivisions(513);
668a0654 864 hs->SetMinimum(50.);
eb05d549 865 hs->GetYaxis()->CenterTitle();
866 cOut->SaveAs(Form("%s.gif", cOut->GetName()));
867}