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