New plots in macro for trending SDD QA-train output (F. Prino)
[u/mrichter/AliRoot.git] / ITS / macrosSDD / TrendQAtrainSDD.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 <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
23 Double_t LangausFun(Double_t *x, Double_t *par);
24 void MakePlots(TString ntupleFileName);
25
26
27 void TrendQAtrainSDD(TString period,
28                      TString recoPass,
29                      TString qaTrain1,
30                      TString qaTrain2,
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
43   TString outFilNam=Form("TrendingSDD_%s_%s_%s.root",period.Data(),recoPass.Data(),qaTrain1.Data());
44
45
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");
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;
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;
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;
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         }
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
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();
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
263 void 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;
274   Float_t fracTrackWithClu1,fracTrackWithClu2,errfracTrackWithClu1,errfracTrackWithClu2;
275   Float_t fracTrackWithClu3,fracTrackWithClu4,errfracTrackWithClu3,errfracTrackWithClu4;
276   Float_t fracTrackWithClu5,fracTrackWithClu6,errfracTrackWithClu5,errfracTrackWithClu6;
277   Float_t fracExtra,errfracExtra;
278   Float_t meandEdxTB0,errmeandEdxTB0,meandEdxTB5,errmeandEdxTB5;
279
280   ntsdd->SetBranchAddress("nrun",&nrun);
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
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());
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
323
324   for(Int_t i=0; i<ntsdd->GetEntries();i++){
325     ntsdd->GetEvent(i);
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));
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");
415   histodEdxTB0->GetYaxis()->SetTitle("<dE/dx> (keV/300 #mum)");
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
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();
468
469 }
470
471 Double_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 }