Updated directory names on LDC in SDD calibration plotting macro (Mario Sitta)
[u/mrichter/AliRoot.git] / ITS / macrosSDD / PlotOutputQAtrainSDD.C
1 #if !defined(__CINT__) || defined(__MAKECINT__)
2 #include <TH1F.h>
3 #include <TH2F.h>
4 #include <TF1.h>
5 #include <TPad.h>
6 #include <TGraphErrors.h>
7 #include <TROOT.h>
8 #include <TFile.h>
9 #include <TTree.h>
10 #include <TGrid.h>
11 #include <TGridResult.h>
12 #include <TMath.h>
13 #include <TCanvas.h>
14 #include <TStyle.h>
15 #include <TLatex.h>
16 #include "AliITSgeomTGeo.h"
17 #endif
18
19 Double_t LangausFun(Double_t *x, Double_t *par);
20
21
22 void PlotOutputQAtrainSDD(TString option="local",
23                           Int_t nRun=0,
24                           TString period="LHC11a",
25                           TString recoPass="pass2_with_SDD",
26                           TString qaTrain="",
27                           TString fileName="QAresults.root"){
28
29   gStyle->SetOptStat(0);
30   //  gStyle->SetOptTitle(0);
31   TFile *f;
32   TString path;
33   Int_t year=2011;
34   if(period.Contains("LHC10")) year=2010;
35   else if(period.Contains("LHC09")) year=2009;
36
37   if(option.Contains("local")){
38     f=new TFile(fileName.Data());
39     printf("Opened file %s\n",f->GetName());
40   }else{
41     TGrid::Connect("alien:");
42     path=Form("/alice/data/%d/%s/%09d/ESDs/%s",year,period.Data(),nRun,recoPass.Data());    
43     printf("search in path %s\n",path.Data());
44     if(!gGrid||!gGrid->IsConnected()) {
45       printf("gGrid not found! exit macro\n");
46       return;
47     }
48     TGridResult *gr = gGrid->Query(path,fileName);
49     Int_t nFiles = gr->GetEntries();
50     if (nFiles < 1){
51       printf("QA file for run %d not found\n",nRun);
52       return;
53     }
54     printf("================>%d files found\n", nFiles);
55     if(qaTrain.Length()>0){
56       Int_t found=kFALSE;
57       for (Int_t iFil = 0; iFil <nFiles ; iFil++) { 
58         fileName=Form("%s",gr->GetKey(iFil,"turl"));
59         TString isMerged=fileName;
60         isMerged.Remove(isMerged.Sizeof()-16); 
61         isMerged.Remove(0,isMerged.Sizeof()-5);
62         if(!isMerged.Contains("QA")) continue;
63         if(fileName.Contains(qaTrain.Data())){
64           found=kTRUE;
65           break;
66         }
67       }
68       if(!found){
69         printf(" File from %s train not found\n",qaTrain.Data());
70         return;
71       }
72     }else{
73       Int_t theFile=0;
74       Int_t maxVer=0;
75       if (nFiles > 1){
76         for (Int_t iFil = 0; iFil <nFiles ; iFil++) { 
77           fileName=Form("%s",gr->GetKey(iFil,"turl"));
78           TString isMerged=fileName;
79           isMerged.Remove(isMerged.Sizeof()-16); 
80           isMerged.Remove(0,isMerged.Sizeof()-5);
81           if(!isMerged.Contains("QA")) continue;
82           TString cutFilename=fileName.ReplaceAll("/QAresults.root","");
83           Int_t last=cutFilename.Sizeof()-3;
84           cutFilename.Remove(0,last);
85           Int_t qaver=atoi(cutFilename.Data());
86           if(qaver>maxVer){
87             maxVer=qaver;
88             theFile=iFil;
89           }
90         }
91       }
92       fileName=Form("%s",gr->GetKey(theFile,"turl"));
93     }
94     f=TFile::Open(fileName.Data());
95   }
96
97   TDirectoryFile* df=(TDirectoryFile*)f->Get("SDD_Performance");
98   if(!df){
99     printf("SDD_Performance MISSING -> Exit\n");
100     return;
101   }
102   TList* l=(TList*)df->Get("coutputRP");
103   if(!df){
104     printf("coutputRP TList MISSING -> Exit\n");
105     return;
106   }
107
108   TH1F* hcllay=(TH1F*)l->FindObject("hCluInLay");
109   TH1F* hapmod=(TH1F*)l->FindObject("hAllPmod");
110   TH1F* hgpmod=(TH1F*)l->FindObject("hGoodPmod");
111   TH1F* hmpmod=(TH1F*)l->FindObject("hMissPmod");
112   TH1F* hbrmod=(TH1F*)l->FindObject("hBadRegmod");
113   TH1F* hskmod=(TH1F*)l->FindObject("hSkippedmod");
114   TH1F* hoamod=(TH1F*)l->FindObject("hOutAccmod");
115   TH1F* hnrmod=(TH1F*)l->FindObject("hNoRefitmod");
116
117   //  TH1F* hapxl=(TH1F*)l->FindObject("hAllPxloc");
118   TH1F* hgpxl=(TH1F*)l->FindObject("hGoodPxloc");
119   TH1F* hmpxl=(TH1F*)l->FindObject("hMissPxloc");
120   TH1F* hbrxl=(TH1F*)l->FindObject("hBadRegxloc");
121   //  TH1F* hapzl=(TH1F*)l->FindObject("hAllPzloc");
122   TH1F* hgpzl=(TH1F*)l->FindObject("hGoodPzloc");
123   TH1F* hmpzl=(TH1F*)l->FindObject("hMissPzloc");
124   TH1F* hbrzl=(TH1F*)l->FindObject("hBadRegzloc");
125
126   TH2F* hClSizAn=(TH2F*)l->FindObject("hCluSizAn");
127   TH2F* hClSizTb=(TH2F*)l->FindObject("hCluSizTb");
128
129   TH2F* hdedx3=(TH2F*)l->FindObject("hdEdxL3VsP");
130   TH2F* hdedx4=(TH2F*)l->FindObject("hdEdxL4VsP");
131   TH2F* hdedxmod=(TH2F*)l->FindObject("hdEdxVsMod");
132   
133
134   TH1F* hmodR=(TH1F*)l->FindObject("hRPMod");
135   TH1F* hmodT=(TH1F*)l->FindObject("hTPMod");
136   TH1F* hgamod=(TH1F*)l->FindObject("hGAMod");
137
138   TH2F* h2dmodR3=new TH2F("h2dmodR3","Rec Points, Layer 3",6,0.5,6.5,14,0.5,14.5);
139   TH2F* h2dmodR4=new TH2F("h2dmodR4","Rec Points, Layer 4",8,0.5,8.5,22,0.5,22.5);
140   TH2F* h2dmodT3=new TH2F("h2dmodT3","Track Points, Layer 3",6,0.5,6.5,14,0.5,14.5);
141   TH2F* h2dmodT4=new TH2F("h2dmodT4","Track Points, Layer 4",8,0.5,8.5,22,0.5,22.5);
142   TH2F* h2dmodR3N=new TH2F("h2dmodR3N","Rec Points/GoodAnode/Event, Layer 3",6,0.5,6.5,14,0.5,14.5);
143   TH2F* h2dmodR4N=new TH2F("h2dmodR4N","Rec Points/GoodAnode/Event, Layer 4",8,0.5,8.5,22,0.5,22.5);
144   TH2F* h2dmodT3N=new TH2F("h2dmodT3N","Track Points/GoodAnode/Event, Layer 3",6,0.5,6.5,14,0.5,14.5);
145   TH2F* h2dmodT4N=new TH2F("h2dmodT4N","Track Points/GoodAnode/Event, Layer 4",8,0.5,8.5,22,0.5,22.5);
146   TH1F* hmodRN=new TH1F("hmodRN","Normalized Rec Points per Module",260,239.5,499.5);
147   TH1F* hmodTN=new TH1F("hmodTN","Normalized Track Points per Module",260,239.5,499.5);
148
149   TH1F* hev=(TH1F*)l->FindObject("hNEvents");
150   Int_t nTotEvents=hev->GetBinContent(2);
151   Int_t nTrigEvents=hev->GetBinContent(3);
152   Int_t nEvents=nTotEvents;
153   printf("---- Statistics ----\n");
154   printf("Number of Events = %d\n",nTotEvents);
155   if(nTrigEvents>0){ 
156     printf("Number of Triggered Events = %d\n",nTrigEvents);
157     nEvents=nTrigEvents;
158   }else{
159     printf("No request on the trigger done when running the task\n");
160   }
161   Int_t bestMod=0;
162   for(Int_t iMod=0; iMod<260;iMod++){
163     Int_t gda=(Int_t)hgamod->GetBinContent(iMod+1);
164     if(gda>bestMod) bestMod=gda;
165   }
166   Int_t nChunks=1;
167   if(bestMod>512){
168     nChunks=(Int_t)(bestMod/512.+0.5);
169   }
170   printf("Chunks merged = %d\n",nChunks);
171   hgamod->Scale(1./nChunks);
172   TCanvas* cgan=new TCanvas("cgan","Good Anodes");
173   cgan->SetTickx();
174   cgan->SetTicky();
175   hgamod->SetMarkerStyle(20);
176   hgamod->SetMarkerSize(0.6);
177   hgamod->SetMinimum(-10.);
178   hgamod->Draw("P");
179   hgamod->GetXaxis()->SetTitle("SDD Module Id");
180   hgamod->GetYaxis()->SetTitle("Number of good anodes");
181   cgan->Update();
182
183   printf("---- Modules with > 2%% of bad anodes ----\n");
184   for(Int_t iMod=0; iMod<260; iMod++){
185     Int_t idMod=iMod+240;
186     Float_t rps=hmodR->GetBinContent(iMod+1);
187     Float_t tps=hmodT->GetBinContent(iMod+1);
188     Float_t ga=hgamod->GetBinContent(iMod+1);
189     if(ga<500){
190       printf("Module %d - Good Anodes = %d\n",idMod,(Int_t)ga);
191     }
192     Float_t rpsN=0.;
193     Float_t tpsN=0.;
194     Float_t erpsN=0.;
195     Float_t etpsN=0.;
196     if(ga>0 && nEvents>0){
197       rpsN=rps/ga/(Float_t)nEvents;
198       tpsN=tps/ga/(Float_t)nEvents;
199       erpsN=TMath::Sqrt(rps)/ga/(Float_t)nEvents;
200       etpsN=TMath::Sqrt(tps)/ga/(Float_t)nEvents;
201     }
202     hmodRN->SetBinContent(iMod+1,rpsN);
203     hmodTN->SetBinContent(iMod+1,tpsN);
204     hmodRN->SetBinError(iMod+1,erpsN);
205     hmodTN->SetBinError(iMod+1,etpsN);
206     Int_t iLay,iLad,iDet;
207     AliITSgeomTGeo::GetModuleId(idMod,iLay,iLad,iDet);
208     if(iLay==3){
209       h2dmodR3->SetBinContent(iDet,iLad,rps);
210       h2dmodT3->SetBinContent(iDet,iLad,tps);
211       h2dmodR3N->SetBinContent(iDet,iLad,rpsN);
212       h2dmodT3N->SetBinContent(iDet,iLad,tpsN);
213     }
214     else if(iLay==4){
215       h2dmodR4->SetBinContent(iDet,iLad,rps);
216       h2dmodT4->SetBinContent(iDet,iLad,tps);
217       h2dmodR4N->SetBinContent(iDet,iLad,rpsN);
218       h2dmodT4N->SetBinContent(iDet,iLad,tpsN);
219     }
220   }
221   if(nEvents<1) return;
222
223   gStyle->SetPalette(1);
224
225   if(hmodR->GetEntries()>0){
226     TCanvas* cmodR=new TCanvas("cmodR","RecPoint Occup",1200,1200);
227     cmodR->Divide(2,3);
228     cmodR->cd(1);
229     gPad->SetLeftMargin(0.14);
230     hmodR->Draw();
231     hmodR->GetXaxis()->SetTitle("SDD Module Id");
232     hmodR->GetYaxis()->SetTitle("RecPoints");
233     hmodR->GetYaxis()->SetTitleOffset(1.55);
234     cmodR->cd(2);
235     gPad->SetLeftMargin(0.14);
236     hmodRN->Draw("E");
237     hmodRN->GetXaxis()->SetTitle("SDD Module Id");
238     hmodRN->GetYaxis()->SetTitle("RecPoints/GoodAnode/Event");
239     hmodRN->GetYaxis()->SetTitleOffset(1.55);
240     cmodR->cd(3);
241     gPad->SetLeftMargin(0.14);
242     h2dmodR3->Draw("colz");
243     h2dmodR3->GetXaxis()->SetTitle("Detector");
244     h2dmodR3->GetYaxis()->SetTitle("Ladder");
245     cmodR->cd(4);
246     gPad->SetLeftMargin(0.14);
247     h2dmodR3N->Draw("colz");
248     h2dmodR3N->GetXaxis()->SetTitle("Detector");
249     h2dmodR3N->GetYaxis()->SetTitle("Ladder");
250     cmodR->cd(5);
251     gPad->SetLeftMargin(0.14);
252     h2dmodR4->Draw("colz");
253     h2dmodR4->GetXaxis()->SetTitle("Detector");
254     h2dmodR4->GetYaxis()->SetTitle("Ladder");
255     cmodR->cd(6);
256     gPad->SetLeftMargin(0.14);
257     gPad->SetLeftMargin(0.14);
258     h2dmodR4N->Draw("colz");
259     h2dmodR4N->GetXaxis()->SetTitle("Detector");
260     h2dmodR4N->GetYaxis()->SetTitle("Ladder");
261     cmodR->Update();
262   }
263
264
265   if(hmodT->GetEntries()>0){
266     TCanvas* cmodT=new TCanvas("cmodT","TrackPoint Occup",1200,1200);
267     cmodT->Divide(2,3);
268     cmodT->cd(1);
269     hmodT->Draw();
270     hmodT->GetXaxis()->SetTitle("SDD Module Id");
271     hmodT->GetYaxis()->SetTitle("TrackPoints");
272     hmodT->GetYaxis()->SetTitleOffset(1.4);
273     cmodT->cd(2);
274     gPad->SetLeftMargin(0.14);
275     hmodTN->Draw("E");
276     hmodTN->GetXaxis()->SetTitle("SDD Module Id");
277     hmodTN->GetYaxis()->SetTitle("TrackPoints");
278     hmodTN->GetYaxis()->SetTitleOffset(1.4);
279     cmodT->cd(3);
280     gPad->SetLeftMargin(0.14);
281     h2dmodT3->Draw("colz");
282     h2dmodT3->GetXaxis()->SetTitle("Detector");
283     h2dmodT3->GetYaxis()->SetTitle("Ladder");
284     cmodT->cd(4);
285     gPad->SetLeftMargin(0.14);
286     h2dmodT3N->Draw("colz");
287     h2dmodT3N->GetXaxis()->SetTitle("Detector");
288     h2dmodT3N->GetYaxis()->SetTitle("Ladder");  
289     cmodT->cd(5);
290     gPad->SetLeftMargin(0.14);
291     h2dmodT4->Draw("colz");
292     h2dmodT4->GetXaxis()->SetTitle("Detector");
293     h2dmodT4->GetYaxis()->SetTitle("Ladder");
294     cmodT->cd(6);
295     gPad->SetLeftMargin(0.14);
296     h2dmodT4N->Draw("colz");
297     h2dmodT4N->GetXaxis()->SetTitle("Detector");
298     h2dmodT4N->GetYaxis()->SetTitle("Ladder");
299     cmodT->Update();
300   }
301
302   TH1F* htplad3=(TH1F*)l->FindObject("hTPLad3");
303   TH1F* htplad4=(TH1F*)l->FindObject("hTPLad4");
304   TH1F* hgalad3=(TH1F*)l->FindObject("hGALad3");
305   TH1F* hgalad4=(TH1F*)l->FindObject("hGALad4");
306   TH1F* hnormOcc3=new TH1F("hnormOcc3","",14,-0.5,13.5);
307   TH1F* hnormOcc4=new TH1F("hnormOcc4","",22,-0.5,21.5);
308   Bool_t tpok=kFALSE;
309   for(Int_t ilad=0;ilad<14;ilad++){ 
310     Float_t occ=0.;
311     Float_t eocc=0.;
312     Int_t gd3=hgalad3->GetBinContent(ilad+1);
313     if(gd3>0){
314       occ=(Float_t)htplad3->GetBinContent(ilad+1)/(Float_t)gd3/(Float_t)nEvents;
315       eocc=TMath::Sqrt((Float_t)htplad3->GetBinContent(ilad+1))/(Float_t)gd3/(Float_t)nEvents;
316     }
317     hnormOcc3->SetBinContent(ilad+1,occ);
318     hnormOcc3->SetBinError(ilad+1,eocc);
319   }
320   for(Int_t ilad=0;ilad<22;ilad++){ 
321     Float_t occ=0.;
322     Float_t eocc=0.;
323     Int_t gd4=hgalad4->GetBinContent(ilad+1);
324     if(gd4>0){
325       occ=(Float_t)htplad4->GetBinContent(ilad+1)/(Float_t)gd4/(Float_t)nEvents;
326       eocc=TMath::Sqrt((Float_t)htplad4->GetBinContent(ilad+1))/(Float_t)gd4/(Float_t)nEvents;
327     }
328     hnormOcc4->SetBinContent(ilad+1,occ);
329     hnormOcc4->SetBinError(ilad+1,eocc);
330   }
331
332
333   if(tpok){
334     TCanvas* cn0=new TCanvas("cn0","Normalized Ladder Occupancy",1400,600);
335     cn0->Divide(2,1);
336     cn0->cd(1);
337     gPad->SetLeftMargin(0.14);
338     hnormOcc3->Draw();
339     hnormOcc3->GetXaxis()->SetTitle("Ladder number (layer 3)");
340     hnormOcc3->GetYaxis()->SetTitle("TrackPoints/GoodAnodes/Events");
341     hnormOcc3->GetYaxis()->SetTitleOffset(1.35);
342     cn0->cd(2);
343     gPad->SetLeftMargin(0.14);
344     hnormOcc4->Draw();
345     hnormOcc4->GetXaxis()->SetTitle("Ladder number (layer 4)");
346     hnormOcc4->GetYaxis()->SetTitle("TrackPoints/GoodAnode/Events");
347     hnormOcc4->GetYaxis()->SetTitleOffset(1.35);
348     cn0->Update();
349   }
350
351   if(hcllay){
352     Double_t norm=hcllay->GetBinContent(1);
353     if(norm>0.){
354       hcllay->Scale(1./norm);
355       hcllay->SetTitle("");
356       hcllay->GetXaxis()->SetRange(2,7);
357       hcllay->SetMinimum(0.);
358       hcllay->SetMaximum(1.1);
359       hcllay->SetMarkerStyle(23);
360       TCanvas* ceffL=new TCanvas("ceffL","PointPerLayer",1000,600);
361       ceffL->SetGridy();
362       hcllay->Draw();
363       hcllay->GetXaxis()->SetTitle("Layer");
364       hcllay->GetYaxis()->SetTitle("Fraction of tracks with point in layer");
365       ceffL->Update();
366     }
367   }
368
369   hgpmod->SetTitle("");
370   TCanvas* ceff0=new TCanvas("ceff0","ModuleIndexInfo",1000,600);
371   hgpmod->Draw();
372   hgpmod->GetXaxis()->SetTitle("SDD Module Id");
373   hgpmod->GetYaxis()->SetTitle("Number of tracks");
374   hmpmod->SetLineColor(2);
375   hmpmod->SetMarkerColor(2);
376   hmpmod->SetMarkerStyle(22);
377   hmpmod->SetMarkerSize(0.5);
378   hmpmod->Draw("psame");
379   hbrmod->SetLineColor(kGreen+1);
380   hbrmod->SetMarkerColor(kGreen+1);
381   hbrmod->SetMarkerStyle(20);
382   hbrmod->SetMarkerSize(0.5);
383   hbrmod->Draw("same");
384   hskmod->SetLineColor(kYellow);
385   hskmod->Draw("same");
386   hoamod->SetLineColor(4);
387   hoamod->Draw("same");
388   hnrmod->SetLineColor(6);
389   hnrmod->Draw("same");
390   TLatex* t1=new TLatex(0.7,0.85,"Good Point");
391   t1->SetNDC();
392   t1->SetTextColor(1);
393   t1->Draw();
394   TLatex* t2=new TLatex(0.7,0.8,"Missing Point");
395   t2->SetNDC();
396   t2->SetTextColor(2);
397   t2->Draw();
398   TLatex* t3=new TLatex(0.7,0.75,"Bad Region");
399   t3->SetNDC();
400   t3->SetTextColor(kGreen+1);
401   t3->Draw();
402   ceff0->Update();
403
404   TH1F* heff=new TH1F("heff","",260,239.5,499.5);
405   TH1F* hfracskip=new TH1F("hfracskip","",260,239.5,499.5);
406
407   for(Int_t imod=0; imod<260;imod++){
408     Float_t numer=hgpmod->GetBinContent(imod+1)+hbrmod->GetBinContent(imod+1)+hoamod->GetBinContent(imod+1)+hnrmod->GetBinContent(imod+1);
409     Float_t denom=hapmod->GetBinContent(imod+1)-hskmod->GetBinContent(imod+1);
410     Float_t eff=0.;
411     Float_t erreff=0.;
412     if(denom>0){
413       eff=numer/denom;
414       erreff=TMath::Sqrt(eff*(1-eff)/denom);
415     }
416     heff->SetBinContent(imod+1,eff);
417     heff->SetBinError(imod+1,erreff);
418     Float_t numer2=hskmod->GetBinContent(imod+1);
419     Float_t denom2=hapmod->GetBinContent(imod+1);
420     Float_t frac=0.;
421     Float_t efrac=0.;
422     if(denom2>0.){
423       frac=numer2/denom2;
424       efrac=TMath::Sqrt(frac*(1-frac)/denom2);
425     }
426     hfracskip->SetBinContent(imod+1,frac);
427     hfracskip->SetBinError(imod+1,efrac);
428
429   }
430
431   printf("---- Modules with efficiency < 90%% ----\n");
432   TCanvas* ceff1=new TCanvas("ceff1","Efficiency",1000,1000);
433   ceff1->Divide(1,2);
434   ceff1->cd(1);
435   heff->Draw();
436   heff->GetXaxis()->SetTitle("SDD Module Id");
437   heff->GetYaxis()->SetTitle("Fraction of tracks with point in good region");
438   for(Int_t ibin=1; ibin<=heff->GetNbinsX(); ibin++){
439     Float_t e=heff->GetBinContent(ibin);
440     if(e<0.9){
441       Int_t iMod=(Int_t)heff->GetBinCenter(ibin);
442       Int_t lay,lad,det;
443       AliITSgeomTGeo::GetModuleId(iMod,lay,lad,det);
444       printf("Module %d - Layer %d Ladder %2d Det %d  -   Eff. %.3f\n",iMod,lay,lad,det,heff->GetBinContent(ibin));
445     }
446   }
447   ceff1->cd(2);
448   hfracskip->Draw();
449   hfracskip->GetXaxis()->SetTitle("SDD Module Id");
450   hfracskip->GetYaxis()->SetTitle("Fraction of tracks with skipped SDD");
451   
452
453   if(hgpxl){
454     hgpxl->SetTitle("");
455     hgpzl->SetTitle("");
456     TCanvas* ceff2=new TCanvas("ceff2","LocalCoord",1000,600);
457     ceff2->Divide(2,1);
458     ceff2->cd(1);
459     hgpxl->Draw();
460     hgpxl->GetXaxis()->SetTitle("Xlocal (cm)");
461     hgpxl->GetYaxis()->SetTitle("Number of tracks");
462     hmpxl->SetLineColor(2);
463     hmpxl->SetMarkerColor(2);
464     hmpxl->SetMarkerStyle(22);
465     hmpxl->SetMarkerSize(0.5);
466     hmpxl->Draw("psame");
467     hbrxl->SetLineColor(kGreen+1);
468     hbrxl->SetMarkerColor(kGreen+1);
469     hbrxl->SetMarkerStyle(20);
470     hbrxl->SetMarkerSize(0.5);
471     hbrxl->Draw("same");
472     t1->Draw();
473     t2->Draw();
474     t3->Draw();
475     ceff2->cd(2);
476     hgpzl->Draw();
477     hgpzl->GetXaxis()->SetTitle("Zlocal (cm)");
478     hgpzl->GetYaxis()->SetTitle("Number of tracks");
479     hmpzl->SetLineColor(2);
480     hmpzl->SetMarkerColor(2);
481     hmpzl->SetMarkerStyle(22);
482     hmpzl->SetMarkerSize(0.5);
483     hmpzl->Draw("psame");
484     hbrzl->SetLineColor(kGreen+1);
485     hbrzl->SetMarkerColor(kGreen+1);
486     hbrzl->SetMarkerStyle(20);
487     hbrzl->SetMarkerSize(0.5);
488     hbrzl->Draw("same");
489     t1->Draw();
490     t2->Draw();
491     t3->Draw();
492     ceff2->Update();
493   }
494
495
496   if(hClSizAn && hClSizTb){
497     TCanvas* ccs=new TCanvas("ccs","Cluster Size",1200,600);
498     ccs->Divide(2,1);
499     ccs->cd(1);
500     gPad->SetLogz();
501     gPad->SetRightMargin(0.12);
502     hClSizAn->Draw("colz");
503     hClSizAn->GetXaxis()->SetTitle("Drift Time (ns)");
504     hClSizAn->GetYaxis()->SetTitle("Cluster Size - Anodes");
505     ccs->cd(2);
506     gPad->SetLogz();
507     gPad->SetRightMargin(0.12);
508     hClSizTb->Draw("colz");
509     hClSizTb->GetXaxis()->SetTitle("Drift Time (ns)");
510     hClSizTb->GetYaxis()->SetTitle("Cluster Size - Time Bins");
511     ccs->Update();
512   }
513
514   TH1F* htimR=(TH1F*)l->FindObject("hDrTimRP");
515   TH1F* htimT=(TH1F*)l->FindObject("hDrTimTPAll");
516   TH1F* htimTe=(TH1F*)l->FindObject("hDrTimTPExtra");
517   TH1F* htimTne=(TH1F*)l->FindObject("hDrTimTPNoExtra");
518   htimR->Rebin(4);
519   htimT->Rebin(4);
520   htimTe->Rebin(4);
521   htimTne->Rebin(4);
522   htimR->SetLineWidth(2);
523   htimT->SetLineWidth(2);
524   htimTe->SetLineWidth(2);
525   htimTne->SetLineWidth(2);
526
527   TCanvas* ctim=new TCanvas("ctim","DriftTime",1400,600);
528   ctim->Divide(2,1);
529   ctim->cd(1);
530   htimR->Draw(); 
531   htimR->GetYaxis()->SetTitleOffset(1.2);
532   htimR->GetXaxis()->SetTitle("Drift Time (ns)");
533   htimR->GetYaxis()->SetTitle("RecPoints");
534   ctim->cd(2);
535   htimT->Draw();
536   htimTe->SetLineColor(2);
537   htimTe->Draw("same");
538   htimTne->SetLineColor(4);
539   htimTne->Draw("same");
540   htimT->GetXaxis()->SetTitle("Drift Time (ns)");
541   htimT->GetYaxis()->SetTitle("TrackPoints");
542   htimT->GetYaxis()->SetTitleOffset(1.2);
543   TLatex* ta=new TLatex(0.5,0.85,"All Clusters");
544   ta->SetNDC();
545   ta->SetTextColor(1);
546   ta->Draw();
547   TLatex* te=new TLatex(0.5,0.8,"Extra Clusters");
548   te->SetNDC();
549   te->SetTextColor(2);
550   te->Draw();
551   TLatex* tn=new TLatex(0.5,0.75,"Non-Extra Clusters");
552   tn->SetNDC();
553   tn->SetTextColor(4);
554   tn->Draw();
555   ctim->Update();
556
557   TCanvas* cdedx=new TCanvas("cdedx","dedx",1400,600);
558   cdedx->Divide(3,1);
559   cdedx->cd(1);
560   gPad->SetLogz();
561   hdedx3->Draw("col");
562   hdedx3->GetXaxis()->SetTitle("P (GeV/c)");
563   hdedx3->GetYaxis()->SetTitle("dE/dx (keV/300 #mum) Layer 3");
564   hdedx3->GetYaxis()->SetTitleOffset(1.25);
565   cdedx->cd(2);
566   gPad->SetLogz();
567   hdedx4->Draw("col");
568   hdedx4->GetXaxis()->SetTitle("P (GeV/c)");
569   hdedx4->GetYaxis()->SetTitle("dE/dx (keV/300 #mum) Layer 4");
570   hdedx4->GetYaxis()->SetTitleOffset(1.25);
571   cdedx->cd(3);
572   gPad->SetLogz();
573   hdedxmod->Draw("col"); 
574   hdedxmod->GetXaxis()->SetTitle("SDD Module Id");
575   hdedxmod->GetYaxis()->SetTitle("dE/dx (keV/300 #mum)");
576   hdedxmod->GetYaxis()->SetTitleOffset(1.25);
577   cdedx->Update();
578
579   printf("---- dE/dx vs.DriftTime ----\n");
580   TCanvas* csig=new TCanvas("csig","dedx vs. DriftTime",1000,700);
581   csig->Divide(4,2);
582   TH1F* hSigTim[8];
583   TGraphErrors* gmpv=new TGraphErrors(0);
584   TGraphErrors* gsigg=new TGraphErrors(0);
585   TGraphErrors* gsigl=new TGraphErrors(0);
586   gmpv->SetTitle("");
587   gsigg->SetTitle("");
588   gsigl->SetTitle("");
589   Int_t iPoint=0;
590   TF1 *lfun = new TF1("LangausFun",LangausFun,50.,300.,4);
591   for(Int_t it=0; it<8; it++){
592     hSigTim[it]=(TH1F*)l->FindObject(Form("hSigTimeInt%d",it));
593     csig->cd(it+1);
594     hSigTim[it]->Draw();
595     if(hSigTim[it]->GetEntries()>200){
596       lfun->SetLineWidth(2);
597       lfun->SetParameter(0,5.);
598       lfun->SetParameter(1,80.);
599       lfun->SetParameter(2,hSigTim[it]->GetEntries()/10.);
600       lfun->SetParameter(3,10.);
601       lfun->SetParLimits(3,0.,20);
602
603       hSigTim[it]->Fit("LangausFun","QLR");
604       hSigTim[it]->GetXaxis()->SetTitle(Form("dE/dx, time interval %d",it+1));
605       hSigTim[it]->GetYaxis()->SetTitle("Events");
606       Float_t mpv=lfun->GetParameter(1);
607       Float_t maxf=lfun->GetMaximumX(70.,90.);
608       printf("mpv=%f   maxfunc=%f\n",mpv,maxf);
609       Float_t empv=lfun->GetParError(1);
610       Float_t sig=lfun->GetParameter(3);
611       Float_t esig=lfun->GetParError(3);
612       Float_t sigl=lfun->GetParameter(0);
613       Float_t esigl=lfun->GetParError(0);
614       gmpv->SetPoint(iPoint,(Float_t)it,mpv);
615       gmpv->SetPointError(iPoint,0.,empv);
616       gsigg->SetPoint(iPoint,(Float_t)it,sig);
617       gsigg->SetPointError(iPoint,0.,esig);
618       gsigl->SetPoint(iPoint,(Float_t)it,sigl);
619       gsigl->SetPointError(iPoint,0.,esigl);
620       ++iPoint;
621       gPad->Update();
622       printf("Bin %d - MPV=%.3f  \t SigmaLandau=%.3f  \t SigmaGaus=%.3f\n",it,mpv,sigl,sig);
623     }
624   }
625
626   TCanvas* cpars=new TCanvas("cpars","Params",800,900);
627   cpars->Divide(1,3,0.01,0.);
628   cpars->cd(1);
629   gPad->SetLeftMargin(0.14);
630   gPad->SetFrameLineWidth(2);
631   gPad->SetTickx();
632   gPad->SetTicky();
633   gmpv->SetMarkerStyle(20);
634   gmpv->SetMinimum(70);
635   gmpv->SetMaximum(90);
636   gmpv->GetXaxis()->SetLimits(-0.2,6.8);
637   gmpv->Draw("AP");
638   //  gmpv->GetXaxis()->SetTitle("Drift Time interval number");
639   gmpv->GetYaxis()->SetTitle("Landau MPV (keV)");
640   gmpv->GetXaxis()->SetTitleSize(0.05);
641   gmpv->GetYaxis()->SetTitleSize(0.05);
642   gmpv->GetYaxis()->SetTitleOffset(1.2);
643   cpars->cd(2);
644   gPad->SetLeftMargin(0.14);
645   gPad->SetFrameLineWidth(2);
646   gPad->SetTickx();
647   gPad->SetTicky();
648   gsigl->SetMarkerStyle(20);
649   gsigl->GetXaxis()->SetLimits(-0.2,6.8);
650   gsigl->SetMinimum(7.);
651   gsigl->SetMaximum(11.);
652   gsigl->Draw("AP");
653   //  gsigl->GetXaxis()->SetTitle("Drift Time interval number");
654   gsigl->GetYaxis()->SetTitle("#sigma_{Landau} (keV)");
655   gsigl->GetXaxis()->SetTitleSize(0.05);
656   gsigl->GetYaxis()->SetTitleSize(0.05);
657   gsigl->GetYaxis()->SetTitleOffset(1.2);
658   cpars->cd(3);
659   gPad->SetLeftMargin(0.14);
660   gPad->SetFrameLineWidth(2);
661   gPad->SetTickx();
662   gPad->SetTicky();
663   gsigg->SetMarkerStyle(20);
664   gsigg->GetXaxis()->SetLimits(-0.2,6.8);
665   gsigg->SetMinimum(3.);
666   gsigg->SetMaximum(7.);
667   gsigg->Draw("AP");
668   gsigg->GetXaxis()->SetTitle("Drift Time interval number");
669   gsigg->GetYaxis()->SetTitle("#sigma_{Gauss} (keV)");
670   gsigg->GetXaxis()->SetTitleSize(0.05);
671   gsigg->GetYaxis()->SetTitleSize(0.05);
672   gsigg->GetYaxis()->SetTitleOffset(1.2);
673   cpars->Update();
674
675 }
676
677
678 Double_t LangausFun(Double_t *x, Double_t *par) {
679
680   //Fit parameters:
681   //par[0]=Width (scale) parameter of Landau density
682   //par[1]=Most Probable (MP, location) parameter of Landau density
683   //par[2]=Total area (integral -inf to inf, normalization constant)
684   //par[3]=Width (sigma) of convoluted Gaussian function
685   //
686   //In the Landau distribution (represented by the CERNLIB approximation), 
687   //the maximum is located at x=-0.22278298 with the location parameter=0.
688   //This shift is corrected within this function, so that the actual
689   //maximum is identical to the MP parameter.
690
691   // Numeric constants
692   Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
693   Double_t mpshift  = -0.22278298;       // Landau maximum location
694
695   // Control constants
696   Double_t np = 100.0;      // number of convolution steps
697   Double_t sc =   5.0;      // convolution extends to +-sc Gaussian sigmas
698
699   // Variables
700   Double_t xx;
701   Double_t mpc;
702   Double_t fland;
703   Double_t sum = 0.0;
704   Double_t xlow,xupp;
705   Double_t step;
706   Double_t i;
707
708
709   // MP shift correction
710   mpc = par[1] - mpshift * par[0]; 
711
712   // Range of convolution integral
713   xlow = x[0] - sc * par[3];
714   xupp = x[0] + sc * par[3];
715
716   step = (xupp-xlow) / np;
717
718   // Convolution integral of Landau and Gaussian by sum
719   for(i=1.0; i<=np/2; i++) {
720     xx = xlow + (i-.5) * step;
721     fland = TMath::Landau(xx,mpc,par[0]) / par[0];
722     sum += fland * TMath::Gaus(x[0],xx,par[3]);
723
724     xx = xupp - (i-.5) * step;
725     fland = TMath::Landau(xx,mpc,par[0]) / par[0];
726     sum += fland * TMath::Gaus(x[0],xx,par[3]);
727   }
728
729   return (par[2] * step * sum * invsq2pi / par[3]);
730 }
731