]> 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
Overlaps corrected, new shape of sectors
[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
193da5be 22void PlotOutputQAtrainSDD(TString option="local",
7541e187 23 Int_t nRun=0,
193da5be 24 TString period="LHC11a",
9a6276c2 25 TString recoPass="pass2_with_SDD",
193da5be 26 TString qaTrain="",
27 TString fileName="QAresults.root"){
7541e187 28
29 gStyle->SetOptStat(0);
30 // gStyle->SetOptTitle(0);
31 TFile *f;
32 TString path;
193da5be 33 Int_t year=2011;
34 if(period.Contains("LHC10")) year=2010;
35 else if(period.Contains("LHC09")) year=2009;
7541e187 36
37 if(option.Contains("local")){
38 f=new TFile(fileName.Data());
5fcb23cb 39 printf("Opened file %s\n",f->GetName());
7541e187 40 }else{
41 TGrid::Connect("alien:");
9a6276c2 42 path=Form("/alice/data/%d/%s/%09d/ESDs/%s",year,period.Data(),nRun,recoPass.Data());
193da5be 43 printf("search in path %s\n",path.Data());
7541e187 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);
193da5be 55 if(qaTrain.Length()>0){
56 Int_t found=kFALSE;
7541e187 57 for (Int_t iFil = 0; iFil <nFiles ; iFil++) {
193da5be 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;
7541e187 66 }
67 }
193da5be 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"));
7541e187 93 }
7541e187 94 f=TFile::Open(fileName.Data());
95 }
5fcb23cb 96
7541e187 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{
193da5be 159 printf("No request on the trigger done when running the task\n");
7541e187 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.;
193da5be 196 if(ga>0 && nEvents>0){
7541e187 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 }
193da5be 221 if(nEvents<1) return;
222
7541e187 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 }
c74ed8cf 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 }
7541e187 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
c74ed8cf 332
7541e187 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 }
c74ed8cf 350
7541e187 351 if(hcllay){
352 Double_t norm=hcllay->GetBinContent(1);
c74ed8cf 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 }
7541e187 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);
f9430df2 405 TH1F* hfracskip=new TH1F("hfracskip","",260,239.5,499.5);
406
7541e187 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);
f9430df2 409 Float_t denom=hapmod->GetBinContent(imod+1)-hskmod->GetBinContent(imod+1);
7541e187 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);
f9430df2 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
7541e187 429 }
430
431 printf("---- Modules with efficiency < 90%% ----\n");
f9430df2 432 TCanvas* ceff1=new TCanvas("ceff1","Efficiency",1000,1000);
433 ceff1->Divide(1,2);
434 ceff1->cd(1);
7541e187 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 }
f9430df2 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
7541e187 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);
d489963e 607 Float_t maxf=lfun->GetMaximumX(70.,90.);
608 printf("mpv=%f maxfunc=%f\n",mpv,maxf);
7541e187 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);
d489963e 634 gmpv->SetMinimum(70);
635 gmpv->SetMaximum(90);
7541e187 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);
d489963e 650 gsigl->SetMinimum(7.);
651 gsigl->SetMaximum(11.);
7541e187 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);
d489963e 665 gsigg->SetMinimum(3.);
666 gsigg->SetMaximum(7.);
7541e187 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
678Double_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