]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/macrosSDD/PlotOutputQAtrainSDD.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Macro to display the output of SDD task in QA-train
[u/mrichter/AliRoot.git] / ITS / macrosSDD / PlotOutputQAtrainSDD.C
CommitLineData
7541e187 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
19Double_t LangausFun(Double_t *x, Double_t *par);
20
21
22void 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
193 gStyle->SetPalette(1);
194
195 if(hmodR->GetEntries()>0){
196 TCanvas* cmodR=new TCanvas("cmodR","RecPoint Occup",1200,1200);
197 cmodR->Divide(2,3);
198 cmodR->cd(1);
199 gPad->SetLeftMargin(0.14);
200 hmodR->Draw();
201 hmodR->GetXaxis()->SetTitle("SDD Module Id");
202 hmodR->GetYaxis()->SetTitle("RecPoints");
203 hmodR->GetYaxis()->SetTitleOffset(1.55);
204 cmodR->cd(2);
205 gPad->SetLeftMargin(0.14);
206 hmodRN->Draw("E");
207 hmodRN->GetXaxis()->SetTitle("SDD Module Id");
208 hmodRN->GetYaxis()->SetTitle("RecPoints/GoodAnode/Event");
209 hmodRN->GetYaxis()->SetTitleOffset(1.55);
210 cmodR->cd(3);
211 gPad->SetLeftMargin(0.14);
212 h2dmodR3->Draw("colz");
213 h2dmodR3->GetXaxis()->SetTitle("Detector");
214 h2dmodR3->GetYaxis()->SetTitle("Ladder");
215 cmodR->cd(4);
216 gPad->SetLeftMargin(0.14);
217 h2dmodR3N->Draw("colz");
218 h2dmodR3N->GetXaxis()->SetTitle("Detector");
219 h2dmodR3N->GetYaxis()->SetTitle("Ladder");
220 cmodR->cd(5);
221 gPad->SetLeftMargin(0.14);
222 h2dmodR4->Draw("colz");
223 h2dmodR4->GetXaxis()->SetTitle("Detector");
224 h2dmodR4->GetYaxis()->SetTitle("Ladder");
225 cmodR->cd(6);
226 gPad->SetLeftMargin(0.14);
227 gPad->SetLeftMargin(0.14);
228 h2dmodR4N->Draw("colz");
229 h2dmodR4N->GetXaxis()->SetTitle("Detector");
230 h2dmodR4N->GetYaxis()->SetTitle("Ladder");
231 cmodR->Update();
232 }
233
234 TCanvas* cmodT=new TCanvas("cmodT","TrackPoint Occup",1200,1200);
235 cmodT->Divide(2,3);
236 cmodT->cd(1);
237 hmodT->Draw();
238 hmodT->GetXaxis()->SetTitle("SDD Module Id");
239 hmodT->GetYaxis()->SetTitle("TrackPoints");
240 hmodT->GetYaxis()->SetTitleOffset(1.4);
241 cmodT->cd(2);
242 gPad->SetLeftMargin(0.14);
243 hmodTN->Draw("E");
244 hmodTN->GetXaxis()->SetTitle("SDD Module Id");
245 hmodTN->GetYaxis()->SetTitle("TrackPoints");
246 hmodTN->GetYaxis()->SetTitleOffset(1.4);
247 cmodT->cd(3);
248 gPad->SetLeftMargin(0.14);
249 h2dmodT3->Draw("colz");
250 h2dmodT3->GetXaxis()->SetTitle("Detector");
251 h2dmodT3->GetYaxis()->SetTitle("Ladder");
252 cmodT->cd(4);
253 gPad->SetLeftMargin(0.14);
254 h2dmodT3N->Draw("colz");
255 h2dmodT3N->GetXaxis()->SetTitle("Detector");
256 h2dmodT3N->GetYaxis()->SetTitle("Ladder");
257 cmodT->cd(5);
258 gPad->SetLeftMargin(0.14);
259 h2dmodT4->Draw("colz");
260 h2dmodT4->GetXaxis()->SetTitle("Detector");
261 h2dmodT4->GetYaxis()->SetTitle("Ladder");
262 cmodT->cd(6);
263 gPad->SetLeftMargin(0.14);
264 h2dmodT4N->Draw("colz");
265 h2dmodT4N->GetXaxis()->SetTitle("Detector");
266 h2dmodT4N->GetYaxis()->SetTitle("Ladder");
267 cmodT->Update();
268
269 TH1F* htplad3=(TH1F*)l->FindObject("hTPLad3");
270 TH1F* htplad4=(TH1F*)l->FindObject("hTPLad4");
271 TH1F* hgalad3=(TH1F*)l->FindObject("hGALad3");
272 TH1F* hgalad4=(TH1F*)l->FindObject("hGALad4");
273 TH1F* hnormOcc3=new TH1F("hnormOcc3","",14,-0.5,13.5);
274 TH1F* hnormOcc4=new TH1F("hnormOcc4","",22,-0.5,21.5);
275 Bool_t tpok=kFALSE;
276 for(Int_t ilad=0;ilad<14;ilad++){
277 Float_t occ=0.;
278 Float_t eocc=0.;
279 Int_t gd3=hgalad3->GetBinContent(ilad+1);
280 if(gd3>0){
281 occ=(Float_t)htplad3->GetBinContent(ilad+1)/(Float_t)gd3/(Float_t)nEvents;
282 eocc=TMath::Sqrt((Float_t)htplad3->GetBinContent(ilad+1))/(Float_t)gd3/(Float_t)nEvents;
283 }
284 hnormOcc3->SetBinContent(ilad+1,occ);
285 hnormOcc3->SetBinError(ilad+1,eocc);
286 }
287 for(Int_t ilad=0;ilad<22;ilad++){
288 Float_t occ=0.;
289 Float_t eocc=0.;
290 Int_t gd4=hgalad4->GetBinContent(ilad+1);
291 if(gd4>0){
292 occ=(Float_t)htplad4->GetBinContent(ilad+1)/(Float_t)gd4/(Float_t)nEvents;
293 eocc=TMath::Sqrt((Float_t)htplad4->GetBinContent(ilad+1))/(Float_t)gd4/(Float_t)nEvents;
294 }
295 hnormOcc4->SetBinContent(ilad+1,occ);
296 hnormOcc4->SetBinError(ilad+1,eocc);
297 }
298
299
300 if(tpok){
301 TCanvas* cn0=new TCanvas("cn0","Normalized Ladder Occupancy",1400,600);
302 cn0->Divide(2,1);
303 cn0->cd(1);
304 gPad->SetLeftMargin(0.14);
305 hnormOcc3->Draw();
306 hnormOcc3->GetXaxis()->SetTitle("Ladder number (layer 3)");
307 hnormOcc3->GetYaxis()->SetTitle("TrackPoints/GoodAnodes/Events");
308 hnormOcc3->GetYaxis()->SetTitleOffset(1.35);
309 cn0->cd(2);
310 gPad->SetLeftMargin(0.14);
311 hnormOcc4->Draw();
312 hnormOcc4->GetXaxis()->SetTitle("Ladder number (layer 4)");
313 hnormOcc4->GetYaxis()->SetTitle("TrackPoints/GoodAnode/Events");
314 hnormOcc4->GetYaxis()->SetTitleOffset(1.35);
315 cn0->Update();
316 }
317 if(hcllay){
318 Double_t norm=hcllay->GetBinContent(1);
319 hcllay->Scale(1./norm);
320 hcllay->SetTitle("");
321 hcllay->GetXaxis()->SetRange(2,7);
322 hcllay->SetMinimum(0.);
323 hcllay->SetMaximum(1.1);
324 hcllay->SetMarkerStyle(23);
325 TCanvas* ceffL=new TCanvas("ceffL","PointPerLayer",1000,600);
326 ceffL->SetGridy();
327 hcllay->Draw();
328 hcllay->GetXaxis()->SetTitle("Layer");
329 hcllay->GetYaxis()->SetTitle("Fraction of tracks with point in layer");
330 ceffL->Update();
331 }
332
333 hgpmod->SetTitle("");
334 TCanvas* ceff0=new TCanvas("ceff0","ModuleIndexInfo",1000,600);
335 hgpmod->Draw();
336 hgpmod->GetXaxis()->SetTitle("SDD Module Id");
337 hgpmod->GetYaxis()->SetTitle("Number of tracks");
338 hmpmod->SetLineColor(2);
339 hmpmod->SetMarkerColor(2);
340 hmpmod->SetMarkerStyle(22);
341 hmpmod->SetMarkerSize(0.5);
342 hmpmod->Draw("psame");
343 hbrmod->SetLineColor(kGreen+1);
344 hbrmod->SetMarkerColor(kGreen+1);
345 hbrmod->SetMarkerStyle(20);
346 hbrmod->SetMarkerSize(0.5);
347 hbrmod->Draw("same");
348 hskmod->SetLineColor(kYellow);
349 hskmod->Draw("same");
350 hoamod->SetLineColor(4);
351 hoamod->Draw("same");
352 hnrmod->SetLineColor(6);
353 hnrmod->Draw("same");
354 TLatex* t1=new TLatex(0.7,0.85,"Good Point");
355 t1->SetNDC();
356 t1->SetTextColor(1);
357 t1->Draw();
358 TLatex* t2=new TLatex(0.7,0.8,"Missing Point");
359 t2->SetNDC();
360 t2->SetTextColor(2);
361 t2->Draw();
362 TLatex* t3=new TLatex(0.7,0.75,"Bad Region");
363 t3->SetNDC();
364 t3->SetTextColor(kGreen+1);
365 t3->Draw();
366 ceff0->Update();
367
368 TH1F* heff=new TH1F("heff","",260,239.5,499.5);
369 for(Int_t imod=0; imod<260;imod++){
370 Float_t numer=hgpmod->GetBinContent(imod+1)+hbrmod->GetBinContent(imod+1)+hoamod->GetBinContent(imod+1)+hnrmod->GetBinContent(imod+1);
371 Float_t denom=hapmod->GetBinContent(imod+1);
372 Float_t eff=0.;
373 Float_t erreff=0.;
374 if(denom>0){
375 eff=numer/denom;
376 erreff=TMath::Sqrt(eff*(1-eff)/denom);
377 }
378 heff->SetBinContent(imod+1,eff);
379 heff->SetBinError(imod+1,erreff);
380 }
381
382 printf("---- Modules with efficiency < 90%% ----\n");
383 TCanvas* ceff1=new TCanvas("ceff1","Efficiency",1000,600);
384 heff->Draw();
385 heff->GetXaxis()->SetTitle("SDD Module Id");
386 heff->GetYaxis()->SetTitle("Fraction of tracks with point in good region");
387 for(Int_t ibin=1; ibin<=heff->GetNbinsX(); ibin++){
388 Float_t e=heff->GetBinContent(ibin);
389 if(e<0.9){
390 Int_t iMod=(Int_t)heff->GetBinCenter(ibin);
391 Int_t lay,lad,det;
392 AliITSgeomTGeo::GetModuleId(iMod,lay,lad,det);
393 printf("Module %d - Layer %d Ladder %2d Det %d - Eff. %.3f\n",iMod,lay,lad,det,heff->GetBinContent(ibin));
394 }
395 }
396 ceff1->Update();
397
398 if(hgpxl){
399 hgpxl->SetTitle("");
400 hgpzl->SetTitle("");
401 TCanvas* ceff2=new TCanvas("ceff2","LocalCoord",1000,600);
402 ceff2->Divide(2,1);
403 ceff2->cd(1);
404 hgpxl->Draw();
405 hgpxl->GetXaxis()->SetTitle("Xlocal (cm)");
406 hgpxl->GetYaxis()->SetTitle("Number of tracks");
407 hmpxl->SetLineColor(2);
408 hmpxl->SetMarkerColor(2);
409 hmpxl->SetMarkerStyle(22);
410 hmpxl->SetMarkerSize(0.5);
411 hmpxl->Draw("psame");
412 hbrxl->SetLineColor(kGreen+1);
413 hbrxl->SetMarkerColor(kGreen+1);
414 hbrxl->SetMarkerStyle(20);
415 hbrxl->SetMarkerSize(0.5);
416 hbrxl->Draw("same");
417 t1->Draw();
418 t2->Draw();
419 t3->Draw();
420 ceff2->cd(2);
421 hgpzl->Draw();
422 hgpzl->GetXaxis()->SetTitle("Zlocal (cm)");
423 hgpzl->GetYaxis()->SetTitle("Number of tracks");
424 hmpzl->SetLineColor(2);
425 hmpzl->SetMarkerColor(2);
426 hmpzl->SetMarkerStyle(22);
427 hmpzl->SetMarkerSize(0.5);
428 hmpzl->Draw("psame");
429 hbrzl->SetLineColor(kGreen+1);
430 hbrzl->SetMarkerColor(kGreen+1);
431 hbrzl->SetMarkerStyle(20);
432 hbrzl->SetMarkerSize(0.5);
433 hbrzl->Draw("same");
434 t1->Draw();
435 t2->Draw();
436 t3->Draw();
437 ceff2->Update();
438 }
439
440
441 if(hClSizAn && hClSizTb){
442 TCanvas* ccs=new TCanvas("ccs","Cluster Size",1200,600);
443 ccs->Divide(2,1);
444 ccs->cd(1);
445 gPad->SetLogz();
446 gPad->SetRightMargin(0.12);
447 hClSizAn->Draw("colz");
448 hClSizAn->GetXaxis()->SetTitle("Drift Time (ns)");
449 hClSizAn->GetYaxis()->SetTitle("Cluster Size - Anodes");
450 ccs->cd(2);
451 gPad->SetLogz();
452 gPad->SetRightMargin(0.12);
453 hClSizTb->Draw("colz");
454 hClSizTb->GetXaxis()->SetTitle("Drift Time (ns)");
455 hClSizTb->GetYaxis()->SetTitle("Cluster Size - Time Bins");
456 ccs->Update();
457 }
458
459 TH1F* htimR=(TH1F*)l->FindObject("hDrTimRP");
460 TH1F* htimT=(TH1F*)l->FindObject("hDrTimTPAll");
461 TH1F* htimTe=(TH1F*)l->FindObject("hDrTimTPExtra");
462 TH1F* htimTne=(TH1F*)l->FindObject("hDrTimTPNoExtra");
463 htimR->Rebin(4);
464 htimT->Rebin(4);
465 htimTe->Rebin(4);
466 htimTne->Rebin(4);
467 htimR->SetLineWidth(2);
468 htimT->SetLineWidth(2);
469 htimTe->SetLineWidth(2);
470 htimTne->SetLineWidth(2);
471
472 TCanvas* ctim=new TCanvas("ctim","DriftTime",1400,600);
473 ctim->Divide(2,1);
474 ctim->cd(1);
475 htimR->Draw();
476 htimR->GetYaxis()->SetTitleOffset(1.2);
477 htimR->GetXaxis()->SetTitle("Drift Time (ns)");
478 htimR->GetYaxis()->SetTitle("RecPoints");
479 ctim->cd(2);
480 htimT->Draw();
481 htimTe->SetLineColor(2);
482 htimTe->Draw("same");
483 htimTne->SetLineColor(4);
484 htimTne->Draw("same");
485 htimT->GetXaxis()->SetTitle("Drift Time (ns)");
486 htimT->GetYaxis()->SetTitle("TrackPoints");
487 htimT->GetYaxis()->SetTitleOffset(1.2);
488 TLatex* ta=new TLatex(0.5,0.85,"All Clusters");
489 ta->SetNDC();
490 ta->SetTextColor(1);
491 ta->Draw();
492 TLatex* te=new TLatex(0.5,0.8,"Extra Clusters");
493 te->SetNDC();
494 te->SetTextColor(2);
495 te->Draw();
496 TLatex* tn=new TLatex(0.5,0.75,"Non-Extra Clusters");
497 tn->SetNDC();
498 tn->SetTextColor(4);
499 tn->Draw();
500 ctim->Update();
501
502 TCanvas* cdedx=new TCanvas("cdedx","dedx",1400,600);
503 cdedx->Divide(3,1);
504 cdedx->cd(1);
505 gPad->SetLogz();
506 hdedx3->Draw("col");
507 hdedx3->GetXaxis()->SetTitle("P (GeV/c)");
508 hdedx3->GetYaxis()->SetTitle("dE/dx (keV/300 #mum) Layer 3");
509 hdedx3->GetYaxis()->SetTitleOffset(1.25);
510 cdedx->cd(2);
511 gPad->SetLogz();
512 hdedx4->Draw("col");
513 hdedx4->GetXaxis()->SetTitle("P (GeV/c)");
514 hdedx4->GetYaxis()->SetTitle("dE/dx (keV/300 #mum) Layer 4");
515 hdedx4->GetYaxis()->SetTitleOffset(1.25);
516 cdedx->cd(3);
517 gPad->SetLogz();
518 hdedxmod->Draw("col");
519 hdedxmod->GetXaxis()->SetTitle("SDD Module Id");
520 hdedxmod->GetYaxis()->SetTitle("dE/dx (keV/300 #mum)");
521 hdedxmod->GetYaxis()->SetTitleOffset(1.25);
522 cdedx->Update();
523
524 printf("---- dE/dx vs.DriftTime ----\n");
525 TCanvas* csig=new TCanvas("csig","dedx vs. DriftTime",1000,700);
526 csig->Divide(4,2);
527 TH1F* hSigTim[8];
528 TGraphErrors* gmpv=new TGraphErrors(0);
529 TGraphErrors* gsigg=new TGraphErrors(0);
530 TGraphErrors* gsigl=new TGraphErrors(0);
531 gmpv->SetTitle("");
532 gsigg->SetTitle("");
533 gsigl->SetTitle("");
534 Int_t iPoint=0;
535 TF1 *lfun = new TF1("LangausFun",LangausFun,50.,300.,4);
536 for(Int_t it=0; it<8; it++){
537 hSigTim[it]=(TH1F*)l->FindObject(Form("hSigTimeInt%d",it));
538 csig->cd(it+1);
539 hSigTim[it]->Draw();
540 if(hSigTim[it]->GetEntries()>200){
541 lfun->SetLineWidth(2);
542 lfun->SetParameter(0,5.);
543 lfun->SetParameter(1,80.);
544 lfun->SetParameter(2,hSigTim[it]->GetEntries()/10.);
545 lfun->SetParameter(3,10.);
546 lfun->SetParLimits(3,0.,20);
547
548 hSigTim[it]->Fit("LangausFun","QLR");
549 hSigTim[it]->GetXaxis()->SetTitle(Form("dE/dx, time interval %d",it+1));
550 hSigTim[it]->GetYaxis()->SetTitle("Events");
551 Float_t mpv=lfun->GetParameter(1);
552 Float_t empv=lfun->GetParError(1);
553 Float_t sig=lfun->GetParameter(3);
554 Float_t esig=lfun->GetParError(3);
555 Float_t sigl=lfun->GetParameter(0);
556 Float_t esigl=lfun->GetParError(0);
557 gmpv->SetPoint(iPoint,(Float_t)it,mpv);
558 gmpv->SetPointError(iPoint,0.,empv);
559 gsigg->SetPoint(iPoint,(Float_t)it,sig);
560 gsigg->SetPointError(iPoint,0.,esig);
561 gsigl->SetPoint(iPoint,(Float_t)it,sigl);
562 gsigl->SetPointError(iPoint,0.,esigl);
563 ++iPoint;
564 gPad->Update();
565 printf("Bin %d - MPV=%.3f \t SigmaLandau=%.3f \t SigmaGaus=%.3f\n",it,mpv,sigl,sig);
566 }
567 }
568
569 TCanvas* cpars=new TCanvas("cpars","Params",800,900);
570 cpars->Divide(1,3,0.01,0.);
571 cpars->cd(1);
572 gPad->SetLeftMargin(0.14);
573 gPad->SetFrameLineWidth(2);
574 gPad->SetTickx();
575 gPad->SetTicky();
576 gmpv->SetMarkerStyle(20);
577 // gmpv->SetMinimum(0);
578 // gmpv->SetMaximum(120);
579 gmpv->GetXaxis()->SetLimits(-0.2,6.8);
580 gmpv->Draw("AP");
581 // gmpv->GetXaxis()->SetTitle("Drift Time interval number");
582 gmpv->GetYaxis()->SetTitle("Landau MPV (keV)");
583 gmpv->GetXaxis()->SetTitleSize(0.05);
584 gmpv->GetYaxis()->SetTitleSize(0.05);
585 gmpv->GetYaxis()->SetTitleOffset(1.2);
586 cpars->cd(2);
587 gPad->SetLeftMargin(0.14);
588 gPad->SetFrameLineWidth(2);
589 gPad->SetTickx();
590 gPad->SetTicky();
591 gsigl->SetMarkerStyle(20);
592 gsigl->GetXaxis()->SetLimits(-0.2,6.8);
593 gsigl->Draw("AP");
594 // gsigl->GetXaxis()->SetTitle("Drift Time interval number");
595 gsigl->GetYaxis()->SetTitle("#sigma_{Landau} (keV)");
596 gsigl->GetXaxis()->SetTitleSize(0.05);
597 gsigl->GetYaxis()->SetTitleSize(0.05);
598 gsigl->GetYaxis()->SetTitleOffset(1.2);
599 cpars->cd(3);
600 gPad->SetLeftMargin(0.14);
601 gPad->SetFrameLineWidth(2);
602 gPad->SetTickx();
603 gPad->SetTicky();
604 gsigg->SetMarkerStyle(20);
605 gsigg->GetXaxis()->SetLimits(-0.2,6.8);
606 gsigg->Draw("AP");
607 gsigg->GetXaxis()->SetTitle("Drift Time interval number");
608 gsigg->GetYaxis()->SetTitle("#sigma_{Gauss} (keV)");
609 gsigg->GetXaxis()->SetTitleSize(0.05);
610 gsigg->GetYaxis()->SetTitleSize(0.05);
611 gsigg->GetYaxis()->SetTitleOffset(1.2);
612 cpars->Update();
613
614}
615
616
617Double_t LangausFun(Double_t *x, Double_t *par) {
618
619 //Fit parameters:
620 //par[0]=Width (scale) parameter of Landau density
621 //par[1]=Most Probable (MP, location) parameter of Landau density
622 //par[2]=Total area (integral -inf to inf, normalization constant)
623 //par[3]=Width (sigma) of convoluted Gaussian function
624 //
625 //In the Landau distribution (represented by the CERNLIB approximation),
626 //the maximum is located at x=-0.22278298 with the location parameter=0.
627 //This shift is corrected within this function, so that the actual
628 //maximum is identical to the MP parameter.
629
630 // Numeric constants
631 Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
632 Double_t mpshift = -0.22278298; // Landau maximum location
633
634 // Control constants
635 Double_t np = 100.0; // number of convolution steps
636 Double_t sc = 5.0; // convolution extends to +-sc Gaussian sigmas
637
638 // Variables
639 Double_t xx;
640 Double_t mpc;
641 Double_t fland;
642 Double_t sum = 0.0;
643 Double_t xlow,xupp;
644 Double_t step;
645 Double_t i;
646
647
648 // MP shift correction
649 mpc = par[1] - mpshift * par[0];
650
651 // Range of convolution integral
652 xlow = x[0] - sc * par[3];
653 xupp = x[0] + sc * par[3];
654
655 step = (xupp-xlow) / np;
656
657 // Convolution integral of Landau and Gaussian by sum
658 for(i=1.0; i<=np/2; i++) {
659 xx = xlow + (i-.5) * step;
660 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
661 sum += fland * TMath::Gaus(x[0],xx,par[3]);
662
663 xx = xupp - (i-.5) * step;
664 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
665 sum += fland * TMath::Gaus(x[0],xx,par[3]);
666 }
667
668 return (par[2] * step * sum * invsq2pi / par[3]);
669}
670