ITS/macrosSDD/PlotOutputQAtrainSDD.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 <TCanvas.h>
  14 #include <TStyle.h>
  15 #include <TLatex.h>
  16 #include "AliITSgeomTGeo.h"
  17 #endif
  18
  19 Double_t LangausFun(Double_t *x, Double_t *par);
  20
  21
  22 void PlotOutputQAtrainSDD(TString option="local",
  23                           Int_t nRun=0,
  24                           TString period="LHC11a",
  25                           TString recoPass="pass2_with_SDD",
  26                           TString qaTrain="",
  27                           TString fileName="QAresults.root"){
  28
  29   gStyle->SetOptStat(0);
  30   //  gStyle->SetOptTitle(0);
  31   TFile *f;
  32   TString path;
  33   Int_t year=2011;
  34   if(period.Contains("LHC10")) year=2010;
  35   else if(period.Contains("LHC09")) year=2009;
  36
  37   if(option.Contains("local")){
  38     f=new TFile(fileName.Data());
  39     printf("Opened file %s\n",f->GetName());
  40   }else{
  41     TGrid::Connect("alien:");
  42     path=Form("/alice/data/%d/%s/%09d/ESDs/%s",year,period.Data(),nRun,recoPass.Data());
  43     printf("search in path %s\n",path.Data());
  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);
  55     if(qaTrain.Length()>0){
  56       Int_t found=kFALSE;
  57       for (Int_t iFil = 0; iFil <nFiles ; iFil++) {
  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;
  66         }
  67       }
  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"));
  93     }
  94     f=TFile::Open(fileName.Data());
  95   }
  96
  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{
 159     printf("No request on the trigger done when running the task\n");
 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.;
 196     if(ga>0 && nEvents>0){
 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   }
 221   if(nEvents<1) return;
 222
 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   }
 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   }
 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
 332
 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   }
 350
 351   if(hcllay){
 352     Double_t norm=hcllay->GetBinContent(1);
 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     }
 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);
 405   TH1F* hfracskip=new TH1F("hfracskip","",260,239.5,499.5);
 406
 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);
 409     Float_t denom=hapmod->GetBinContent(imod+1)-hskmod->GetBinContent(imod+1);
 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);
 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
 429   }
 430
 431   printf("---- Modules with efficiency < 90%% ----\n");
 432   TCanvas* ceff1=new TCanvas("ceff1","Efficiency",1000,1000);
 433   ceff1->Divide(1,2);
 434   ceff1->cd(1);
 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   }
 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
 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);
 607       Float_t empv=lfun->GetParError(1);
 608       Float_t sig=lfun->GetParameter(3);
 609       Float_t esig=lfun->GetParError(3);
 610       Float_t sigl=lfun->GetParameter(0);
 611       Float_t esigl=lfun->GetParError(0);
 612       gmpv->SetPoint(iPoint,(Float_t)it,mpv);
 613       gmpv->SetPointError(iPoint,0.,empv);
 614       gsigg->SetPoint(iPoint,(Float_t)it,sig);
 615       gsigg->SetPointError(iPoint,0.,esig);
 616       gsigl->SetPoint(iPoint,(Float_t)it,sigl);
 617       gsigl->SetPointError(iPoint,0.,esigl);
 618       ++iPoint;
 619       gPad->Update();
 620       printf("Bin %d - MPV=%.3f  \t SigmaLandau=%.3f  \t SigmaGaus=%.3f\n",it,mpv,sigl,sig);
 621     }
 622   }
 623
 624   TCanvas* cpars=new TCanvas("cpars","Params",800,900);
 625   cpars->Divide(1,3,0.01,0.);
 626   cpars->cd(1);
 627   gPad->SetLeftMargin(0.14);
 628   gPad->SetFrameLineWidth(2);
 629   gPad->SetTickx();
 630   gPad->SetTicky();
 631   gmpv->SetMarkerStyle(20);
 632   //  gmpv->SetMinimum(0);
 633   //  gmpv->SetMaximum(120);
 634   gmpv->GetXaxis()->SetLimits(-0.2,6.8);
 635   gmpv->Draw("AP");
 636   //  gmpv->GetXaxis()->SetTitle("Drift Time interval number");
 637   gmpv->GetYaxis()->SetTitle("Landau MPV (keV)");
 638   gmpv->GetXaxis()->SetTitleSize(0.05);
 639   gmpv->GetYaxis()->SetTitleSize(0.05);
 640   gmpv->GetYaxis()->SetTitleOffset(1.2);
 641   cpars->cd(2);
 642   gPad->SetLeftMargin(0.14);
 643   gPad->SetFrameLineWidth(2);
 644   gPad->SetTickx();
 645   gPad->SetTicky();
 646   gsigl->SetMarkerStyle(20);
 647   gsigl->GetXaxis()->SetLimits(-0.2,6.8);
 648   gsigl->Draw("AP");
 649   //  gsigl->GetXaxis()->SetTitle("Drift Time interval number");
 650   gsigl->GetYaxis()->SetTitle("#sigma_{Landau} (keV)");
 651   gsigl->GetXaxis()->SetTitleSize(0.05);
 652   gsigl->GetYaxis()->SetTitleSize(0.05);
 653   gsigl->GetYaxis()->SetTitleOffset(1.2);
 654   cpars->cd(3);
 655   gPad->SetLeftMargin(0.14);
 656   gPad->SetFrameLineWidth(2);
 657   gPad->SetTickx();
 658   gPad->SetTicky();
 659   gsigg->SetMarkerStyle(20);
 660   gsigg->GetXaxis()->SetLimits(-0.2,6.8);
 661   gsigg->Draw("AP");
 662   gsigg->GetXaxis()->SetTitle("Drift Time interval number");
 663   gsigg->GetYaxis()->SetTitle("#sigma_{Gauss} (keV)");
 664   gsigg->GetXaxis()->SetTitleSize(0.05);
 665   gsigg->GetYaxis()->SetTitleSize(0.05);
 666   gsigg->GetYaxis()->SetTitleOffset(1.2);
 667   cpars->Update();
 668
 669 }
 670
 671
 672 Double_t LangausFun(Double_t *x, Double_t *par) {
 673
 674   //Fit parameters:
 675   //par[0]=Width (scale) parameter of Landau density
 676   //par[1]=Most Probable (MP, location) parameter of Landau density
 677   //par[2]=Total area (integral -inf to inf, normalization constant)
 678   //par[3]=Width (sigma) of convoluted Gaussian function
 679   //
 680   //In the Landau distribution (represented by the CERNLIB approximation),
 681   //the maximum is located at x=-0.22278298 with the location parameter=0.
 682   //This shift is corrected within this function, so that the actual
 683   //maximum is identical to the MP parameter.
 684
 685   // Numeric constants
 686   Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
 687   Double_t mpshift  = -0.22278298;       // Landau maximum location
 688
 689   // Control constants
 690   Double_t np = 100.0;      // number of convolution steps
 691   Double_t sc =   5.0;      // convolution extends to +-sc Gaussian sigmas
 692
 693   // Variables
 694   Double_t xx;
 695   Double_t mpc;
 696   Double_t fland;
 697   Double_t sum = 0.0;
 698   Double_t xlow,xupp;
 699   Double_t step;
 700   Double_t i;
 701
 702
 703   // MP shift correction
 704   mpc = par[1] - mpshift * par[0];
 705
 706   // Range of convolution integral
 707   xlow = x[0] - sc * par[3];
 708   xupp = x[0] + sc * par[3];
 709
 710   step = (xupp-xlow) / np;
 711
 712   // Convolution integral of Landau and Gaussian by sum
 713   for(i=1.0; i<=np/2; i++) {
 714     xx = xlow + (i-.5) * step;
 715     fland = TMath::Landau(xx,mpc,par[0]) / par[0];
 716     sum += fland * TMath::Gaus(x[0],xx,par[3]);
 717
 718     xx = xupp - (i-.5) * step;
 719     fland = TMath::Landau(xx,mpc,par[0]) / par[0];
 720     sum += fland * TMath::Gaus(x[0],xx,par[3]);
 721   }
 722
 723   return (par[2] * step * sum * invsq2pi / par[3]);
 724 }
 725