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