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