New plots in macro for trending SDD QA-train output (F. Prino)
[u/mrichter/AliRoot.git] / ITS / macrosSDD / TrendQAtrainSDD.C
CommitLineData
341adf04 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 <TSystem.h>
14#include <TNtuple.h>
15#include <TCanvas.h>
16#include <TStyle.h>
17#include <TLatex.h>
18#include <TLegend.h>
19#include <TLegendEntry.h>
20#include "AliITSgeomTGeo.h"
21#endif
22
23Double_t LangausFun(Double_t *x, Double_t *par);
24void MakePlots(TString ntupleFileName);
25
26
27void TrendQAtrainSDD(TString period,
28 TString recoPass,
e1a5dc0a 29 TString qaTrain1,
30 TString qaTrain2,
341adf04 31 Int_t firstRun=0,
32 Int_t lastRun=999999999,
33 TString fileName="QAresults.root"){
34
35 gStyle->SetOptStat(0);
36
37
38 TGrid::Connect("alien:");
39 Int_t year=0;
40 if(period.Contains("LHC09")) year=2009;
41 else if(period.Contains("LHC10")) year=2010;
42
e1a5dc0a 43 TString outFilNam=Form("TrendingSDD_%s_%s_%s.root",period.Data(),recoPass.Data(),qaTrain1.Data());
341adf04 44
45
e1a5dc0a 46 const Int_t nVariables=27;
47 TNtuple* ntsdd=new TNtuple("ntsdd","SDD trending","nrun:fracTrackWithClu1:errfracTrackWithClu1:fracTrackWithClu2:errfracTrackWithClu2:fracTrackWithClu3:errfracTrackWithClu3:fracTrackWithClu4:errfracTrackWithClu4:fracTrackWithClu5:errfracTrackWithClu5:fracTrackWithClu6:errfracTrackWithClu6:meanTrPts3:errmeanTrPts3:meanTrPts4:errmeanTrPts4:minDrTime:errminDrTime:meanDrTime:errmeanDrTime:fracExtra:errfracExtra:meandEdxTB0:errmeandEdxTB0:meandEdxTB5:errmeandEdxTB5");
341adf04 48 Float_t xnt[nVariables];
49
50 TBits* readRun=new TBits(999999);
51 readRun->ResetAllBits();
52 if(!gSystem->Exec(Form("ls -l %s > /dev/null 2>&1",outFilNam.Data()))){
53 TFile* oldfil=new TFile(outFilNam.Data());
54 TNtuple* ntmp=(TNtuple*)oldfil->Get("ntsdd");
55 Bool_t isOK=kFALSE;
56 if(ntmp){
57 if(ntmp->GetNvar()==ntsdd->GetNvar()){
58 isOK=kTRUE;
59 TObjArray* arr1=(TObjArray*)ntsdd->GetListOfBranches();
60 TObjArray* arr2=(TObjArray*)ntmp->GetListOfBranches();
61 for(Int_t iV=0; iV<ntmp->GetNvar(); iV++){
62 TString vnam1=arr1->At(iV)->GetName();
63 TString vnam2=arr2->At(iV)->GetName();
64 if(vnam1!=vnam2) isOK=kFALSE;
65 ntmp->SetBranchAddress(vnam2.Data(),&xnt[iV]);
66 }
67 if(isOK){
68 for(Int_t nE=0; nE<ntmp->GetEntries(); nE++){
69 ntmp->GetEvent(nE);
70 Int_t theRun=(Int_t)(xnt[0]+0.0001);
71 readRun->SetBitNumber(theRun);
72 ntsdd->Fill(xnt);
73 }
74 }
75 }
76 }
77 if(!isOK){
78 printf("Ntuple in local file not OK -> will be recreated\n");
79 }
80 oldfil->Close();
81 delete oldfil;
82 }
83
84 if(!gGrid||!gGrid->IsConnected()) {
85 printf("gGrid not found! exit macro\n");
86 return;
87 }
88
89 TString path=Form("/alice/data/%d/%s/",year,period.Data(),recoPass.Data());
90 TGridResult *gr = gGrid->Query(path,fileName);
91 Int_t nFiles = gr->GetEntries();
92 printf("================>%d files found\n", nFiles);
93 if (nFiles < 1) return;
94
95
96 if (nFiles > 1){
97 for (Int_t iFil = 0; iFil <nFiles ; iFil++) {
98 TString fileNameLong=Form("%s",gr->GetKey(iFil,"turl"));
99 if(!fileNameLong.Contains(recoPass.Data())) continue;
e1a5dc0a 100 if(!fileNameLong.Contains(qaTrain1.Data()) &&
101 !fileNameLong.Contains(qaTrain2.Data())) continue;
102 if(!fileNameLong.Contains(Form("%s/%s",qaTrain1.Data(),fileName.Data())) &&
103 !fileNameLong.Contains(Form("%s/%s",qaTrain2.Data(),fileName.Data()))) continue;
341adf04 104 TString runNumber=fileNameLong;
105 runNumber.ReplaceAll(Form("alien:///alice/data/%d/%s/",year,period.Data()),"");
106 runNumber.Remove(9,runNumber.Sizeof());
107
108 Int_t iRun=atoi(runNumber.Data());
109 if(readRun->TestBitNumber(iRun)){
110 printf("Run %d aleady in local ntuple -> skipping it\n",iRun);
111 continue;
112 }
113 if(iRun<firstRun) continue;
114 if(iRun>lastRun) continue;
115
116 TFile* f=TFile::Open(fileNameLong.Data());
117
118 TDirectoryFile* df=(TDirectoryFile*)f->Get("SDD_Performance");
119 if(!df){
120 printf("Run %d SDD_Performance MISSING -> Exit\n",iRun);
121 continue;
122 }
123 TList* l=(TList*)df->Get("coutputRP");
124 if(!df){
125 printf("Run %d coutputRP TList MISSING -> Exit\n",iRun);
126 continue;
e1a5dc0a 127 }
128 TH1F* hcllay=(TH1F*)l->FindObject("hCluInLay");
129 Float_t fracT[6]={0.,0.,0.,0.,0.,0.};
130 Float_t efracT[6]={0.,0.,0.,0.,0.,0.};
131 if(hcllay->GetBinContent(1)>0){
132 for(Int_t iLay=0; iLay<6; iLay++){
133 fracT[iLay]=hcllay->GetBinContent(iLay+2)/hcllay->GetBinContent(1);
134 efracT[iLay]=TMath::Sqrt(fracT[iLay]*(1-fracT[iLay])/hcllay->GetBinContent(1));
135 }
341adf04 136 }
137 TH1F* hmodT=(TH1F*)l->FindObject("hTPMod");
138 TH1F* hgamod=(TH1F*)l->FindObject("hGAMod");
139 Int_t bestMod=0;
140 for(Int_t iMod=0; iMod<260;iMod++){
141 Int_t gda=(Int_t)hgamod->GetBinContent(iMod+1);
142 if(gda>bestMod) bestMod=gda;
143 }
144 Int_t nChunks=1;
145 if(bestMod>512){
146 nChunks=(Int_t)(bestMod/512.+0.5);
147 }
148 hgamod->Scale(1./nChunks);
149
150 TH1F* hev=(TH1F*)l->FindObject("hNEvents");
151 Int_t nTotEvents=hev->GetBinContent(2);
152 Int_t nTrigEvents=hev->GetBinContent(3);
153 Int_t nEvents=nTotEvents;
154 printf("Run %d Number of Events = %d Triggered=%d\n",iRun,nTotEvents,nTrigEvents);
155 if(nTrigEvents>0){
156 nEvents=nTrigEvents;
157 }
158
159 Int_t nModGood3=0;
160 Int_t nModGood4=0;
161 Int_t nModBadAn=0;
162 Float_t sumtp3=0;
163 Float_t sumtp4=0;
164 Float_t sumEtp3=0;
165 Float_t sumEtp4=0;
166 for(Int_t iMod=0; iMod<260; iMod++){
167 Float_t tps=hmodT->GetBinContent(iMod+1);
168 Float_t ga=hgamod->GetBinContent(iMod+1);
169 if(ga<500) nModBadAn++;
170 Float_t tpsN=0.;
171 Float_t etpsN=0.;
172 if(ga>0){
173 tpsN=tps/ga/(Float_t)nEvents;
174 etpsN=TMath::Sqrt(tps)/ga/(Float_t)nEvents;
175 if(iMod<84){
176 sumtp3+=tpsN;
177 sumEtp3+=(etpsN*etpsN);
178 nModGood3++;
179 }
180 else{
181 sumtp4+=tpsN;
182 sumEtp4+=(etpsN*etpsN);
183 nModGood4++;
184 }
185 }
186 }
187
188 TH1F* htimT=(TH1F*)l->FindObject("hDrTimTPAll");
189 TH1F* htimTe=(TH1F*)l->FindObject("hDrTimTPExtra");
190
191 Double_t fracExtra=0.;
192 Double_t errFracExtra=0.;
193 if(htimT->GetEntries()>0){
194 fracExtra=htimTe->GetEntries()/htimT->GetEntries();
195 errFracExtra=TMath::Sqrt(htimTe->GetEntries())/htimT->GetEntries();
196 }
197 Double_t averPoints=0.;
198 Double_t cntBins=0.;
199 for(Int_t iBin=1; iBin<=htimT->GetNbinsX(); iBin++){
200 Float_t tim=htimT->GetBinCenter(iBin);
201 if(tim>2000. && tim<4000.){
202 averPoints+=htimT->GetBinContent(iBin);
203 cntBins+=1;
204 }
205 }
206 Double_t minTime=-999.;
207 Double_t errMinTime=0.;
208 if(cntBins>0){
209 averPoints/=cntBins;
210 for(Int_t iBin=1; iBin<=htimT->GetNbinsX(); iBin++){
211 if(htimT->GetBinContent(iBin)>0.5*averPoints){
212 minTime=htimT->GetBinCenter(iBin);
213 errMinTime=0.5*htimT->GetBinWidth(iBin);
214 break;
215 }
216 }
217 }
218 TH1F* hSigTim0=(TH1F*)l->FindObject("hSigTimeInt0");
219 TH1F* hSigTim5=(TH1F*)l->FindObject("hSigTimeInt5");
220
e1a5dc0a 221 Int_t index=0;
222 xnt[index++]=(Float_t)iRun;
223 xnt[index++]=fracT[0];
224 xnt[index++]=efracT[0];
225 xnt[index++]=fracT[1];
226 xnt[index++]=efracT[1];
227 xnt[index++]=fracT[2];
228 xnt[index++]=efracT[2];
229 xnt[index++]=fracT[3];
230 xnt[index++]=efracT[3];
231 xnt[index++]=fracT[4];
232 xnt[index++]=efracT[4];
233 xnt[index++]=fracT[5];
234 xnt[index++]=efracT[5];
235 xnt[index++]=sumtp3/nModGood3;
236 xnt[index++]=TMath::Sqrt(sumEtp3)/nModGood3;
237 xnt[index++]=sumtp4/nModGood4;
238 xnt[index++]=TMath::Sqrt(sumEtp4)/nModGood4;
239 xnt[index++]=minTime;
240 xnt[index++]=errMinTime;
241 xnt[index++]=htimT->GetMean();
242 xnt[index++]=htimT->GetMeanError();
243 xnt[index++]=fracExtra;
244 xnt[index++]=errFracExtra;
245 xnt[index++]=hSigTim0->GetMean();
246 xnt[index++]=hSigTim0->GetMeanError();
247 xnt[index++]=hSigTim5->GetMean();
248 xnt[index++]=hSigTim5->GetMeanError();
341adf04 249
250 ntsdd->Fill(xnt);
251 }
252 }
253
254 TFile* outfil=new TFile(outFilNam.Data(),"recreate");
255 outfil->cd();
256 ntsdd->Write();
257 outfil->Close();
258
259 MakePlots(outFilNam);
260
261}
262
263void MakePlots(TString ntupleFileName){
264 TFile* fil=new TFile(ntupleFileName.Data(),"read");
265 if(!fil){
266 printf("File with ntuple does not exist\n");
267 return;
268 }
269 TNtuple* ntsdd=(TNtuple*)fil->Get("ntsdd");
270
271 Float_t nrun;
272 Float_t meanTrPts3,errmeanTrPts3,meanTrPts4,errmeanTrPts4;
273 Float_t minDrTime,errminDrTime,meanDrTime,errmeanDrTime;
e1a5dc0a 274 Float_t fracTrackWithClu1,fracTrackWithClu2,errfracTrackWithClu1,errfracTrackWithClu2;
275 Float_t fracTrackWithClu3,fracTrackWithClu4,errfracTrackWithClu3,errfracTrackWithClu4;
276 Float_t fracTrackWithClu5,fracTrackWithClu6,errfracTrackWithClu5,errfracTrackWithClu6;
341adf04 277 Float_t fracExtra,errfracExtra;
278 Float_t meandEdxTB0,errmeandEdxTB0,meandEdxTB5,errmeandEdxTB5;
279
280 ntsdd->SetBranchAddress("nrun",&nrun);
e1a5dc0a 281 ntsdd->SetBranchAddress("fracTrackWithClu1",&fracTrackWithClu1);
282 ntsdd->SetBranchAddress("errfracTrackWithClu1",&errfracTrackWithClu1);
283 ntsdd->SetBranchAddress("fracTrackWithClu2",&fracTrackWithClu2);
284 ntsdd->SetBranchAddress("errfracTrackWithClu2",&errfracTrackWithClu2);
285 ntsdd->SetBranchAddress("fracTrackWithClu3",&fracTrackWithClu3);
286 ntsdd->SetBranchAddress("errfracTrackWithClu3",&errfracTrackWithClu3);
287 ntsdd->SetBranchAddress("fracTrackWithClu4",&fracTrackWithClu4);
288 ntsdd->SetBranchAddress("errfracTrackWithClu4",&errfracTrackWithClu4);
289 ntsdd->SetBranchAddress("fracTrackWithClu5",&fracTrackWithClu5);
290 ntsdd->SetBranchAddress("errfracTrackWithClu5",&errfracTrackWithClu5);
291 ntsdd->SetBranchAddress("fracTrackWithClu6",&fracTrackWithClu6);
292 ntsdd->SetBranchAddress("errfracTrackWithClu6",&errfracTrackWithClu6);
293
341adf04 294 ntsdd->SetBranchAddress("meanTrPts3",&meanTrPts3);
295 ntsdd->SetBranchAddress("errmeanTrPts3",&errmeanTrPts3);
296 ntsdd->SetBranchAddress("meanTrPts4",&meanTrPts4);
297 ntsdd->SetBranchAddress("errmeanTrPts4",&errmeanTrPts4);
298 ntsdd->SetBranchAddress("minDrTime",&minDrTime);
299 ntsdd->SetBranchAddress("errminDrTime",&errminDrTime);
300 ntsdd->SetBranchAddress("meanDrTime",&meanDrTime);
301 ntsdd->SetBranchAddress("errmeanDrTime",&errmeanDrTime);
302 ntsdd->SetBranchAddress("fracExtra",&fracExtra);
303 ntsdd->SetBranchAddress("errfracExtra",&errfracExtra);
304 ntsdd->SetBranchAddress("meandEdxTB0",&meandEdxTB0);
305 ntsdd->SetBranchAddress("errmeandEdxTB0",&errmeandEdxTB0);
306 ntsdd->SetBranchAddress("meandEdxTB5",&meandEdxTB5);
307 ntsdd->SetBranchAddress("errmeandEdxTB5",&errmeandEdxTB5);
308
309 TH1F* histotrp3=new TH1F("histotrp3","",(Int_t)ntsdd->GetEntries(),0.,ntsdd->GetEntries());
310 TH1F* histotrp4=new TH1F("histotrp4","",(Int_t)ntsdd->GetEntries(),0.,ntsdd->GetEntries());
311 TH1F* histominTime=new TH1F("histominTime","",(Int_t)ntsdd->GetEntries(),0.,ntsdd->GetEntries());
312 TH1F* histomeanTime=new TH1F("histomeanTime","",(Int_t)ntsdd->GetEntries(),0.,ntsdd->GetEntries());
313 TH1F* histofracExtra=new TH1F("histofracExtra","",(Int_t)ntsdd->GetEntries(),0.,ntsdd->GetEntries());
314 TH1F* histodEdxTB0=new TH1F("histodEdxTB0","",(Int_t)ntsdd->GetEntries(),0.,ntsdd->GetEntries());
315 TH1F* histodEdxTB5=new TH1F("histodEdxTB5","",(Int_t)ntsdd->GetEntries(),0.,ntsdd->GetEntries());
e1a5dc0a 316 TH1F* histoTrackClu1=new TH1F("histoTrackClu1","",(Int_t)ntsdd->GetEntries(),0.,ntsdd->GetEntries());
317 TH1F* histoTrackClu2=new TH1F("histoTrackClu2","",(Int_t)ntsdd->GetEntries(),0.,ntsdd->GetEntries());
318 TH1F* histoTrackClu3=new TH1F("histoTrackClu3","",(Int_t)ntsdd->GetEntries(),0.,ntsdd->GetEntries());
319 TH1F* histoTrackClu4=new TH1F("histoTrackClu4","",(Int_t)ntsdd->GetEntries(),0.,ntsdd->GetEntries());
320 TH1F* histoTrackClu5=new TH1F("histoTrackClu5","",(Int_t)ntsdd->GetEntries(),0.,ntsdd->GetEntries());
321 TH1F* histoTrackClu6=new TH1F("histoTrackClu6","",(Int_t)ntsdd->GetEntries(),0.,ntsdd->GetEntries());
322
341adf04 323
324 for(Int_t i=0; i<ntsdd->GetEntries();i++){
325 ntsdd->GetEvent(i);
e1a5dc0a 326 histoTrackClu1->SetBinContent(i+1,fracTrackWithClu1);
327 histoTrackClu1->SetBinError(i+1,errfracTrackWithClu1);
328 histoTrackClu1->GetXaxis()->SetBinLabel(i+1,Form("%d",(Int_t)nrun));
329 histoTrackClu2->SetBinContent(i+1,fracTrackWithClu2);
330 histoTrackClu2->SetBinError(i+1,errfracTrackWithClu2);
331 histoTrackClu2->GetXaxis()->SetBinLabel(i+1,Form("%d",(Int_t)nrun));
332 histoTrackClu3->SetBinContent(i+1,fracTrackWithClu3);
333 histoTrackClu3->SetBinError(i+1,errfracTrackWithClu3);
334 histoTrackClu3->GetXaxis()->SetBinLabel(i+1,Form("%d",(Int_t)nrun));
335 histoTrackClu4->SetBinContent(i+1,fracTrackWithClu4);
336 histoTrackClu4->SetBinError(i+1,errfracTrackWithClu4);
337 histoTrackClu4->GetXaxis()->SetBinLabel(i+1,Form("%d",(Int_t)nrun));
338 histoTrackClu5->SetBinContent(i+1,fracTrackWithClu5);
339 histoTrackClu5->SetBinError(i+1,errfracTrackWithClu5);
340 histoTrackClu5->GetXaxis()->SetBinLabel(i+1,Form("%d",(Int_t)nrun));
341 histoTrackClu6->SetBinContent(i+1,fracTrackWithClu6);
342 histoTrackClu6->SetBinError(i+1,errfracTrackWithClu6);
343 histoTrackClu6->GetXaxis()->SetBinLabel(i+1,Form("%d",(Int_t)nrun));
341adf04 344 histominTime->SetBinContent(i+1,minDrTime);
345 histominTime->SetBinError(i+1,errminDrTime);
346 histominTime->GetXaxis()->SetBinLabel(i+1,Form("%d",(Int_t)nrun));
347 histomeanTime->SetBinContent(i+1,meanDrTime);
348 histomeanTime->SetBinError(i+1,errmeanDrTime);
349 histomeanTime->GetXaxis()->SetBinLabel(i+1,Form("%d",(Int_t)nrun));
350 histotrp3->SetBinContent(i+1,meanTrPts3);
351 histotrp3->SetBinError(i+1,errmeanTrPts3);
352 histotrp3->GetXaxis()->SetBinLabel(i+1,Form("%d",(Int_t)nrun));
353 histotrp4->SetBinContent(i+1,meanTrPts4);
354 histotrp4->SetBinError(i+1,errmeanTrPts3);
355 histotrp4->GetXaxis()->SetBinLabel(i+1,Form("%d",(Int_t)nrun));
356 histofracExtra->SetBinContent(i+1,fracExtra);
357 histofracExtra->SetBinError(i+1,errfracExtra);
358 histofracExtra->GetXaxis()->SetBinLabel(i+1,Form("%d",(Int_t)nrun));
359 histodEdxTB0->SetBinContent(i+1,meandEdxTB0);
360 histodEdxTB0->SetBinError(i+1,errmeandEdxTB0);
361 histodEdxTB0->GetXaxis()->SetBinLabel(i+1,Form("%d",(Int_t)nrun));
362 histodEdxTB5->SetBinContent(i+1,meandEdxTB5);
363 histodEdxTB5->SetBinError(i+1,errmeandEdxTB5);
364 histodEdxTB5->GetXaxis()->SetBinLabel(i+1,Form("%d",(Int_t)nrun));
365 }
366
367 gStyle->SetOptStat(0);
368
369
370 TCanvas* c1=new TCanvas("c1","Occupancy");
371 histotrp3->SetLineColor(1);
372 histotrp3->SetMarkerStyle(20);
373 histotrp3->Draw();
374 histotrp3->SetMinimum(0.);
375 histotrp4->SetLineColor(2);
376 histotrp4->SetMarkerColor(2);
377 histotrp4->SetMarkerStyle(22);
378 histotrp4->Draw("same");
379 histotrp3->GetYaxis()->SetTitle("Track Point Occupancy");
380 TLegend* leg=new TLegend(0.7,0.15,0.88,0.35);
381 TLegendEntry* ent=leg->AddEntry(histotrp3,"Layer3","PL");
382 ent=leg->AddEntry(histotrp4,"Layer4","PL");
383 ent->SetTextColor(histotrp4->GetMarkerColor());
384 leg->SetFillStyle(0);
385 leg->Draw();
386 c1->Update();
387
388 TCanvas* c2=new TCanvas("c2","DriftTime",600,900);
389 c2->Divide(1,2);
390 c2->cd(1);
391 histominTime->Draw();
392 histominTime->SetMinimum(250);
393 histominTime->SetMaximum(1000);
394 histominTime->GetYaxis()->SetTitle("Minimum Drift Time (ns)");
395 c2->cd(2);
396 histomeanTime->Draw();
397 histomeanTime->GetYaxis()->SetTitle("Average Drift Time (ns)");
398 c2->Update();
399
400 TCanvas* c3=new TCanvas("c3","ExtraClusters");
401 histofracExtra->Draw();
402 histofracExtra->GetYaxis()->SetTitle("Fraction of Extra Clusters");
403 c3->Update();
404
405
406 TCanvas* c4=new TCanvas("c4","Charge");
407 histodEdxTB0->SetLineColor(1);
408 histodEdxTB0->SetMarkerStyle(20);
409 histodEdxTB0->Draw();
410 histodEdxTB0->SetMinimum(0.);
411 histodEdxTB5->SetLineColor(4);
412 histodEdxTB5->SetMarkerColor(4);
413 histodEdxTB5->SetMarkerStyle(23);
414 histodEdxTB5->Draw("same");
e1a5dc0a 415 histodEdxTB0->GetYaxis()->SetTitle("<dE/dx> (keV/300 #mum)");
341adf04 416 TLegend* leg2=new TLegend(0.6,0.15,0.88,0.35);
417 ent=leg2->AddEntry(histodEdxTB0,"Small drift time","PL");
418 ent=leg2->AddEntry(histodEdxTB5,"Large drift time","PL");
419 ent->SetTextColor(histodEdxTB5->GetMarkerColor());
420 leg2->SetFillStyle(0);
421 leg2->Draw();
422 c4->Update();
423
e1a5dc0a 424 TCanvas* c5=new TCanvas("c5","TrackWithSDD");
425 histoTrackClu3->Draw();
426 histoTrackClu3->SetLineColor(1);
427 histoTrackClu3->SetMarkerStyle(20);
428 histoTrackClu3->Draw();
429 histoTrackClu3->SetMinimum(0.);
430 histoTrackClu4->SetLineColor(2);
431 histoTrackClu4->SetMarkerColor(2);
432 histoTrackClu4->SetMarkerStyle(22);
433 histoTrackClu4->Draw("same");
434 histoTrackClu1->SetLineColor(kGray+1);
435 histoTrackClu1->SetMarkerColor(kGray+1);
436 histoTrackClu1->SetMarkerStyle(24);
437 histoTrackClu1->Draw("same");
438 histoTrackClu2->SetLineColor(kGray+2);
439 histoTrackClu2->SetMarkerColor(kGray+2);
440 histoTrackClu2->SetMarkerStyle(26);
441 histoTrackClu2->Draw("same");
442 histoTrackClu5->SetLineColor(4);
443 histoTrackClu5->SetMarkerColor(4);
444 histoTrackClu5->SetMarkerStyle(29);
445 histoTrackClu5->Draw("same");
446 histoTrackClu6->SetLineColor(kBlue+1);
447 histoTrackClu6->SetMarkerColor(kBlue+1);
448 histoTrackClu6->SetMarkerStyle(30);
449 histoTrackClu6->Draw("same");
450 histoTrackClu3->GetYaxis()->SetTitle("Fraction of Tracks with Point in SDD Layers");
451 TLegend* leg3=new TLegend(0.7,0.15,0.88,0.35);
452 ent=leg3->AddEntry(histoTrackClu1,"Layer1","PL");
453 ent->SetTextColor(histoTrackClu1->GetMarkerColor());
454 ent=leg3->AddEntry(histoTrackClu2,"Layer2","PL");
455 ent->SetTextColor(histoTrackClu2->GetMarkerColor());
456 ent=leg3->AddEntry(histoTrackClu3,"Layer3","PL");
457 ent->SetTextColor(histoTrackClu3->GetMarkerColor());
458 ent=leg3->AddEntry(histoTrackClu4,"Layer4","PL");
459 ent->SetTextColor(histoTrackClu4->GetMarkerColor());
460 ent=leg3->AddEntry(histoTrackClu5,"Layer5","PL");
461 ent->SetTextColor(histoTrackClu5->GetMarkerColor());
462 ent=leg3->AddEntry(histoTrackClu6,"Layer6","PL");
463 ent->SetTextColor(histoTrackClu6->GetMarkerColor());
464
465 leg3->SetFillStyle(0);
466 leg3->Draw();
467 c5->Update();
341adf04 468
469}
470
471Double_t LangausFun(Double_t *x, Double_t *par) {
472
473 //Fit parameters:
474 //par[0]=Width (scale) parameter of Landau density
475 //par[1]=Most Probable (MP, location) parameter of Landau density
476 //par[2]=Total area (integral -inf to inf, normalization constant)
477 //par[3]=Width (sigma) of convoluted Gaussian function
478 //
479 //In the Landau distribution (represented by the CERNLIB approximation),
480 //the maximum is located at x=-0.22278298 with the location parameter=0.
481 //This shift is corrected within this function, so that the actual
482 //maximum is identical to the MP parameter.
483
484 // Numeric constants
485 Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
486 Double_t mpshift = -0.22278298; // Landau maximum location
487
488 // Control constants
489 Double_t np = 100.0; // number of convolution steps
490 Double_t sc = 5.0; // convolution extends to +-sc Gaussian sigmas
491
492 // Variables
493 Double_t xx;
494 Double_t mpc;
495 Double_t fland;
496 Double_t sum = 0.0;
497 Double_t xlow,xupp;
498 Double_t step;
499 Double_t i;
500
501
502 // MP shift correction
503 mpc = par[1] - mpshift * par[0];
504
505 // Range of convolution integral
506 xlow = x[0] - sc * par[3];
507 xupp = x[0] + sc * par[3];
508
509 step = (xupp-xlow) / np;
510
511 // Convolution integral of Landau and Gaussian by sum
512 for(i=1.0; i<=np/2; i++) {
513 xx = xlow + (i-.5) * step;
514 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
515 sum += fland * TMath::Gaus(x[0],xx,par[3]);
516
517 xx = xupp - (i-.5) * step;
518 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
519 sum += fland * TMath::Gaus(x[0],xx,par[3]);
520 }
521
522 return (par[2] * step * sum * invsq2pi / par[3]);
523}