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                           TString fileName="QAresults.root",
  24                           Int_t nRun=0,
  25                           Int_t year=2010,
  26                           TString period="LHC10e"){
  27
  28   gStyle->SetOptStat(0);
  29   //  gStyle->SetOptTitle(0);
  30   TFile *f;
  31   TString path;
  32
  33   if(option.Contains("local")){
  34     f=new TFile(fileName.Data());
  35   }else{
  36     TGrid::Connect("alien:");
  37     path=Form("/alice/data/%d/%s/%09d/ESDs/",year,period.Data(),nRun);
  38     if(!gGrid||!gGrid->IsConnected()) {
  39       printf("gGrid not found! exit macro\n");
  40       return;
  41     }
  42     TGridResult *gr = gGrid->Query(path,fileName);
  43     Int_t nFiles = gr->GetEntries();
  44     if (nFiles < 1){
  45       printf("QA file for run %d not found\n",nRun);
  46       return;
  47     }
  48     printf("================>%d files found\n", nFiles);
  49     Int_t theFile=0;
  50     Int_t maxVer=0;
  51     if (nFiles > 1){
  52       for (Int_t iFil = 0; iFil <nFiles ; iFil++) {
  53         fileName=Form("%s",gr->GetKey(iFil,"turl"));
  54         TString cutFilename=fileName.ReplaceAll("/QAresults.root","");
  55         Int_t last=cutFilename.Sizeof()-3;
  56         cutFilename.Remove(0,last);
  57         Int_t qaver=atoi(cutFilename.Data());
  58         if(qaver>maxVer){
  59           maxVer=qaver;
  60           theFile=iFil;
  61         }
  62       }
  63     }
  64     fileName=Form("%s",gr->GetKey(theFile,"turl"));
  65     f=TFile::Open(fileName.Data());
  66   }
  67
  68   TDirectoryFile* df=(TDirectoryFile*)f->Get("SDD_Performance");
  69   if(!df){
  70     printf("SDD_Performance MISSING -> Exit\n");
  71     return;
  72   }
  73   TList* l=(TList*)df->Get("coutputRP");
  74   if(!df){
  75     printf("coutputRP TList MISSING -> Exit\n");
  76     return;
  77   }
  78
  79   TH1F* hcllay=(TH1F*)l->FindObject("hCluInLay");
  80   TH1F* hapmod=(TH1F*)l->FindObject("hAllPmod");
  81   TH1F* hgpmod=(TH1F*)l->FindObject("hGoodPmod");
  82   TH1F* hmpmod=(TH1F*)l->FindObject("hMissPmod");
  83   TH1F* hbrmod=(TH1F*)l->FindObject("hBadRegmod");
  84   TH1F* hskmod=(TH1F*)l->FindObject("hSkippedmod");
  85   TH1F* hoamod=(TH1F*)l->FindObject("hOutAccmod");
  86   TH1F* hnrmod=(TH1F*)l->FindObject("hNoRefitmod");
  87
  88   //  TH1F* hapxl=(TH1F*)l->FindObject("hAllPxloc");
  89   TH1F* hgpxl=(TH1F*)l->FindObject("hGoodPxloc");
  90   TH1F* hmpxl=(TH1F*)l->FindObject("hMissPxloc");
  91   TH1F* hbrxl=(TH1F*)l->FindObject("hBadRegxloc");
  92   //  TH1F* hapzl=(TH1F*)l->FindObject("hAllPzloc");
  93   TH1F* hgpzl=(TH1F*)l->FindObject("hGoodPzloc");
  94   TH1F* hmpzl=(TH1F*)l->FindObject("hMissPzloc");
  95   TH1F* hbrzl=(TH1F*)l->FindObject("hBadRegzloc");
  96
  97   TH2F* hClSizAn=(TH2F*)l->FindObject("hCluSizAn");
  98   TH2F* hClSizTb=(TH2F*)l->FindObject("hCluSizTb");
  99
 100   TH2F* hdedx3=(TH2F*)l->FindObject("hdEdxL3VsP");
 101   TH2F* hdedx4=(TH2F*)l->FindObject("hdEdxL4VsP");
 102   TH2F* hdedxmod=(TH2F*)l->FindObject("hdEdxVsMod");
 103
 104
 105   TH1F* hmodR=(TH1F*)l->FindObject("hRPMod");
 106   TH1F* hmodT=(TH1F*)l->FindObject("hTPMod");
 107   TH1F* hgamod=(TH1F*)l->FindObject("hGAMod");
 108
 109   TH2F* h2dmodR3=new TH2F("h2dmodR3","Rec Points, Layer 3",6,0.5,6.5,14,0.5,14.5);
 110   TH2F* h2dmodR4=new TH2F("h2dmodR4","Rec Points, Layer 4",8,0.5,8.5,22,0.5,22.5);
 111   TH2F* h2dmodT3=new TH2F("h2dmodT3","Track Points, Layer 3",6,0.5,6.5,14,0.5,14.5);
 112   TH2F* h2dmodT4=new TH2F("h2dmodT4","Track Points, Layer 4",8,0.5,8.5,22,0.5,22.5);
 113   TH2F* h2dmodR3N=new TH2F("h2dmodR3N","Rec Points/GoodAnode/Event, Layer 3",6,0.5,6.5,14,0.5,14.5);
 114   TH2F* h2dmodR4N=new TH2F("h2dmodR4N","Rec Points/GoodAnode/Event, Layer 4",8,0.5,8.5,22,0.5,22.5);
 115   TH2F* h2dmodT3N=new TH2F("h2dmodT3N","Track Points/GoodAnode/Event, Layer 3",6,0.5,6.5,14,0.5,14.5);
 116   TH2F* h2dmodT4N=new TH2F("h2dmodT4N","Track Points/GoodAnode/Event, Layer 4",8,0.5,8.5,22,0.5,22.5);
 117   TH1F* hmodRN=new TH1F("hmodRN","Normalized Rec Points per Module",260,239.5,499.5);
 118   TH1F* hmodTN=new TH1F("hmodTN","Normalized Track Points per Module",260,239.5,499.5);
 119
 120   TH1F* hev=(TH1F*)l->FindObject("hNEvents");
 121   Int_t nTotEvents=hev->GetBinContent(2);
 122   Int_t nTrigEvents=hev->GetBinContent(3);
 123   Int_t nEvents=nTotEvents;
 124   printf("---- Statistics ----\n");
 125   printf("Number of Events = %d\n",nTotEvents);
 126   if(nTrigEvents>0){
 127     printf("Number of Triggered Events = %d\n",nTrigEvents);
 128     nEvents=nTrigEvents;
 129   }else{
 130     printf("No request on the trigger done whenrunning the task\n");
 131   }
 132   Int_t bestMod=0;
 133   for(Int_t iMod=0; iMod<260;iMod++){
 134     Int_t gda=(Int_t)hgamod->GetBinContent(iMod+1);
 135     if(gda>bestMod) bestMod=gda;
 136   }
 137   Int_t nChunks=1;
 138   if(bestMod>512){
 139     nChunks=(Int_t)(bestMod/512.+0.5);
 140   }
 141   printf("Chunks merged = %d\n",nChunks);
 142   hgamod->Scale(1./nChunks);
 143   TCanvas* cgan=new TCanvas("cgan","Good Anodes");
 144   cgan->SetTickx();
 145   cgan->SetTicky();
 146   hgamod->SetMarkerStyle(20);
 147   hgamod->SetMarkerSize(0.6);
 148   hgamod->SetMinimum(-10.);
 149   hgamod->Draw("P");
 150   hgamod->GetXaxis()->SetTitle("SDD Module Id");
 151   hgamod->GetYaxis()->SetTitle("Number of good anodes");
 152   cgan->Update();
 153
 154   printf("---- Modules with > 2%% of bad anodes ----\n");
 155   for(Int_t iMod=0; iMod<260; iMod++){
 156     Int_t idMod=iMod+240;
 157     Float_t rps=hmodR->GetBinContent(iMod+1);
 158     Float_t tps=hmodT->GetBinContent(iMod+1);
 159     Float_t ga=hgamod->GetBinContent(iMod+1);
 160     if(ga<500){
 161       printf("Module %d - Good Anodes = %d\n",idMod,(Int_t)ga);
 162     }
 163     Float_t rpsN=0.;
 164     Float_t tpsN=0.;
 165     Float_t erpsN=0.;
 166     Float_t etpsN=0.;
 167     if(ga>0){
 168       rpsN=rps/ga/(Float_t)nEvents;
 169       tpsN=tps/ga/(Float_t)nEvents;
 170       erpsN=TMath::Sqrt(rps)/ga/(Float_t)nEvents;
 171       etpsN=TMath::Sqrt(tps)/ga/(Float_t)nEvents;
 172     }
 173     hmodRN->SetBinContent(iMod+1,rpsN);
 174     hmodTN->SetBinContent(iMod+1,tpsN);
 175     hmodRN->SetBinError(iMod+1,erpsN);
 176     hmodTN->SetBinError(iMod+1,etpsN);
 177     Int_t iLay,iLad,iDet;
 178     AliITSgeomTGeo::GetModuleId(idMod,iLay,iLad,iDet);
 179     if(iLay==3){
 180       h2dmodR3->SetBinContent(iDet,iLad,rps);
 181       h2dmodT3->SetBinContent(iDet,iLad,tps);
 182       h2dmodR3N->SetBinContent(iDet,iLad,rpsN);
 183       h2dmodT3N->SetBinContent(iDet,iLad,tpsN);
 184     }
 185     else if(iLay==4){
 186       h2dmodR4->SetBinContent(iDet,iLad,rps);
 187       h2dmodT4->SetBinContent(iDet,iLad,tps);
 188       h2dmodR4N->SetBinContent(iDet,iLad,rpsN);
 189       h2dmodT4N->SetBinContent(iDet,iLad,tpsN);
 190     }
 191   }
 192
 193   gStyle->SetPalette(1);
 194
 195   if(hmodR->GetEntries()>0){
 196     TCanvas* cmodR=new TCanvas("cmodR","RecPoint Occup",1200,1200);
 197     cmodR->Divide(2,3);
 198     cmodR->cd(1);
 199     gPad->SetLeftMargin(0.14);
 200     hmodR->Draw();
 201     hmodR->GetXaxis()->SetTitle("SDD Module Id");
 202     hmodR->GetYaxis()->SetTitle("RecPoints");
 203     hmodR->GetYaxis()->SetTitleOffset(1.55);
 204     cmodR->cd(2);
 205     gPad->SetLeftMargin(0.14);
 206     hmodRN->Draw("E");
 207     hmodRN->GetXaxis()->SetTitle("SDD Module Id");
 208     hmodRN->GetYaxis()->SetTitle("RecPoints/GoodAnode/Event");
 209     hmodRN->GetYaxis()->SetTitleOffset(1.55);
 210     cmodR->cd(3);
 211     gPad->SetLeftMargin(0.14);
 212     h2dmodR3->Draw("colz");
 213     h2dmodR3->GetXaxis()->SetTitle("Detector");
 214     h2dmodR3->GetYaxis()->SetTitle("Ladder");
 215     cmodR->cd(4);
 216     gPad->SetLeftMargin(0.14);
 217     h2dmodR3N->Draw("colz");
 218     h2dmodR3N->GetXaxis()->SetTitle("Detector");
 219     h2dmodR3N->GetYaxis()->SetTitle("Ladder");
 220     cmodR->cd(5);
 221     gPad->SetLeftMargin(0.14);
 222     h2dmodR4->Draw("colz");
 223     h2dmodR4->GetXaxis()->SetTitle("Detector");
 224     h2dmodR4->GetYaxis()->SetTitle("Ladder");
 225     cmodR->cd(6);
 226     gPad->SetLeftMargin(0.14);
 227     gPad->SetLeftMargin(0.14);
 228     h2dmodR4N->Draw("colz");
 229     h2dmodR4N->GetXaxis()->SetTitle("Detector");
 230     h2dmodR4N->GetYaxis()->SetTitle("Ladder");
 231     cmodR->Update();
 232   }
 233
 234   TCanvas* cmodT=new TCanvas("cmodT","TrackPoint Occup",1200,1200);
 235   cmodT->Divide(2,3);
 236   cmodT->cd(1);
 237   hmodT->Draw();
 238   hmodT->GetXaxis()->SetTitle("SDD Module Id");
 239   hmodT->GetYaxis()->SetTitle("TrackPoints");
 240   hmodT->GetYaxis()->SetTitleOffset(1.4);
 241   cmodT->cd(2);
 242   gPad->SetLeftMargin(0.14);
 243   hmodTN->Draw("E");
 244   hmodTN->GetXaxis()->SetTitle("SDD Module Id");
 245   hmodTN->GetYaxis()->SetTitle("TrackPoints");
 246   hmodTN->GetYaxis()->SetTitleOffset(1.4);
 247   cmodT->cd(3);
 248   gPad->SetLeftMargin(0.14);
 249   h2dmodT3->Draw("colz");
 250   h2dmodT3->GetXaxis()->SetTitle("Detector");
 251   h2dmodT3->GetYaxis()->SetTitle("Ladder");
 252   cmodT->cd(4);
 253   gPad->SetLeftMargin(0.14);
 254   h2dmodT3N->Draw("colz");
 255   h2dmodT3N->GetXaxis()->SetTitle("Detector");
 256   h2dmodT3N->GetYaxis()->SetTitle("Ladder");
 257   cmodT->cd(5);
 258   gPad->SetLeftMargin(0.14);
 259   h2dmodT4->Draw("colz");
 260   h2dmodT4->GetXaxis()->SetTitle("Detector");
 261   h2dmodT4->GetYaxis()->SetTitle("Ladder");
 262   cmodT->cd(6);
 263   gPad->SetLeftMargin(0.14);
 264   h2dmodT4N->Draw("colz");
 265   h2dmodT4N->GetXaxis()->SetTitle("Detector");
 266   h2dmodT4N->GetYaxis()->SetTitle("Ladder");
 267   cmodT->Update();
 268
 269   TH1F* htplad3=(TH1F*)l->FindObject("hTPLad3");
 270   TH1F* htplad4=(TH1F*)l->FindObject("hTPLad4");
 271   TH1F* hgalad3=(TH1F*)l->FindObject("hGALad3");
 272   TH1F* hgalad4=(TH1F*)l->FindObject("hGALad4");
 273   TH1F* hnormOcc3=new TH1F("hnormOcc3","",14,-0.5,13.5);
 274   TH1F* hnormOcc4=new TH1F("hnormOcc4","",22,-0.5,21.5);
 275   Bool_t tpok=kFALSE;
 276   for(Int_t ilad=0;ilad<14;ilad++){
 277     Float_t occ=0.;
 278     Float_t eocc=0.;
 279     Int_t gd3=hgalad3->GetBinContent(ilad+1);
 280     if(gd3>0){
 281       occ=(Float_t)htplad3->GetBinContent(ilad+1)/(Float_t)gd3/(Float_t)nEvents;
 282       eocc=TMath::Sqrt((Float_t)htplad3->GetBinContent(ilad+1))/(Float_t)gd3/(Float_t)nEvents;
 283     }
 284     hnormOcc3->SetBinContent(ilad+1,occ);
 285     hnormOcc3->SetBinError(ilad+1,eocc);
 286   }
 287   for(Int_t ilad=0;ilad<22;ilad++){
 288     Float_t occ=0.;
 289     Float_t eocc=0.;
 290     Int_t gd4=hgalad4->GetBinContent(ilad+1);
 291     if(gd4>0){
 292       occ=(Float_t)htplad4->GetBinContent(ilad+1)/(Float_t)gd4/(Float_t)nEvents;
 293       eocc=TMath::Sqrt((Float_t)htplad4->GetBinContent(ilad+1))/(Float_t)gd4/(Float_t)nEvents;
 294     }
 295     hnormOcc4->SetBinContent(ilad+1,occ);
 296     hnormOcc4->SetBinError(ilad+1,eocc);
 297   }
 298
 299
 300   if(tpok){
 301     TCanvas* cn0=new TCanvas("cn0","Normalized Ladder Occupancy",1400,600);
 302     cn0->Divide(2,1);
 303     cn0->cd(1);
 304     gPad->SetLeftMargin(0.14);
 305     hnormOcc3->Draw();
 306     hnormOcc3->GetXaxis()->SetTitle("Ladder number (layer 3)");
 307     hnormOcc3->GetYaxis()->SetTitle("TrackPoints/GoodAnodes/Events");
 308     hnormOcc3->GetYaxis()->SetTitleOffset(1.35);
 309     cn0->cd(2);
 310     gPad->SetLeftMargin(0.14);
 311     hnormOcc4->Draw();
 312     hnormOcc4->GetXaxis()->SetTitle("Ladder number (layer 4)");
 313     hnormOcc4->GetYaxis()->SetTitle("TrackPoints/GoodAnode/Events");
 314     hnormOcc4->GetYaxis()->SetTitleOffset(1.35);
 315     cn0->Update();
 316   }
 317   if(hcllay){
 318     Double_t norm=hcllay->GetBinContent(1);
 319     hcllay->Scale(1./norm);
 320     hcllay->SetTitle("");
 321     hcllay->GetXaxis()->SetRange(2,7);
 322     hcllay->SetMinimum(0.);
 323     hcllay->SetMaximum(1.1);
 324     hcllay->SetMarkerStyle(23);
 325     TCanvas* ceffL=new TCanvas("ceffL","PointPerLayer",1000,600);
 326     ceffL->SetGridy();
 327     hcllay->Draw();
 328     hcllay->GetXaxis()->SetTitle("Layer");
 329     hcllay->GetYaxis()->SetTitle("Fraction of tracks with point in layer");
 330     ceffL->Update();
 331   }
 332
 333   hgpmod->SetTitle("");
 334   TCanvas* ceff0=new TCanvas("ceff0","ModuleIndexInfo",1000,600);
 335   hgpmod->Draw();
 336   hgpmod->GetXaxis()->SetTitle("SDD Module Id");
 337   hgpmod->GetYaxis()->SetTitle("Number of tracks");
 338   hmpmod->SetLineColor(2);
 339   hmpmod->SetMarkerColor(2);
 340   hmpmod->SetMarkerStyle(22);
 341   hmpmod->SetMarkerSize(0.5);
 342   hmpmod->Draw("psame");
 343   hbrmod->SetLineColor(kGreen+1);
 344   hbrmod->SetMarkerColor(kGreen+1);
 345   hbrmod->SetMarkerStyle(20);
 346   hbrmod->SetMarkerSize(0.5);
 347   hbrmod->Draw("same");
 348   hskmod->SetLineColor(kYellow);
 349   hskmod->Draw("same");
 350   hoamod->SetLineColor(4);
 351   hoamod->Draw("same");
 352   hnrmod->SetLineColor(6);
 353   hnrmod->Draw("same");
 354   TLatex* t1=new TLatex(0.7,0.85,"Good Point");
 355   t1->SetNDC();
 356   t1->SetTextColor(1);
 357   t1->Draw();
 358   TLatex* t2=new TLatex(0.7,0.8,"Missing Point");
 359   t2->SetNDC();
 360   t2->SetTextColor(2);
 361   t2->Draw();
 362   TLatex* t3=new TLatex(0.7,0.75,"Bad Region");
 363   t3->SetNDC();
 364   t3->SetTextColor(kGreen+1);
 365   t3->Draw();
 366   ceff0->Update();
 367
 368   TH1F* heff=new TH1F("heff","",260,239.5,499.5);
 369   for(Int_t imod=0; imod<260;imod++){
 370     Float_t numer=hgpmod->GetBinContent(imod+1)+hbrmod->GetBinContent(imod+1)+hoamod->GetBinContent(imod+1)+hnrmod->GetBinContent(imod+1);
 371     Float_t denom=hapmod->GetBinContent(imod+1);
 372     Float_t eff=0.;
 373     Float_t erreff=0.;
 374     if(denom>0){
 375       eff=numer/denom;
 376       erreff=TMath::Sqrt(eff*(1-eff)/denom);
 377     }
 378     heff->SetBinContent(imod+1,eff);
 379     heff->SetBinError(imod+1,erreff);
 380   }
 381
 382   printf("---- Modules with efficiency < 90%% ----\n");
 383   TCanvas* ceff1=new TCanvas("ceff1","Efficiency",1000,600);
 384   heff->Draw();
 385   heff->GetXaxis()->SetTitle("SDD Module Id");
 386   heff->GetYaxis()->SetTitle("Fraction of tracks with point in good region");
 387   for(Int_t ibin=1; ibin<=heff->GetNbinsX(); ibin++){
 388     Float_t e=heff->GetBinContent(ibin);
 389     if(e<0.9){
 390       Int_t iMod=(Int_t)heff->GetBinCenter(ibin);
 391       Int_t lay,lad,det;
 392       AliITSgeomTGeo::GetModuleId(iMod,lay,lad,det);
 393       printf("Module %d - Layer %d Ladder %2d Det %d  -   Eff. %.3f\n",iMod,lay,lad,det,heff->GetBinContent(ibin));
 394     }
 395   }
 396   ceff1->Update();
 397
 398   if(hgpxl){
 399     hgpxl->SetTitle("");
 400     hgpzl->SetTitle("");
 401     TCanvas* ceff2=new TCanvas("ceff2","LocalCoord",1000,600);
 402     ceff2->Divide(2,1);
 403     ceff2->cd(1);
 404     hgpxl->Draw();
 405     hgpxl->GetXaxis()->SetTitle("Xlocal (cm)");
 406     hgpxl->GetYaxis()->SetTitle("Number of tracks");
 407     hmpxl->SetLineColor(2);
 408     hmpxl->SetMarkerColor(2);
 409     hmpxl->SetMarkerStyle(22);
 410     hmpxl->SetMarkerSize(0.5);
 411     hmpxl->Draw("psame");
 412     hbrxl->SetLineColor(kGreen+1);
 413     hbrxl->SetMarkerColor(kGreen+1);
 414     hbrxl->SetMarkerStyle(20);
 415     hbrxl->SetMarkerSize(0.5);
 416     hbrxl->Draw("same");
 417     t1->Draw();
 418     t2->Draw();
 419     t3->Draw();
 420     ceff2->cd(2);
 421     hgpzl->Draw();
 422     hgpzl->GetXaxis()->SetTitle("Zlocal (cm)");
 423     hgpzl->GetYaxis()->SetTitle("Number of tracks");
 424     hmpzl->SetLineColor(2);
 425     hmpzl->SetMarkerColor(2);
 426     hmpzl->SetMarkerStyle(22);
 427     hmpzl->SetMarkerSize(0.5);
 428     hmpzl->Draw("psame");
 429     hbrzl->SetLineColor(kGreen+1);
 430     hbrzl->SetMarkerColor(kGreen+1);
 431     hbrzl->SetMarkerStyle(20);
 432     hbrzl->SetMarkerSize(0.5);
 433     hbrzl->Draw("same");
 434     t1->Draw();
 435     t2->Draw();
 436     t3->Draw();
 437     ceff2->Update();
 438   }
 439
 440
 441   if(hClSizAn && hClSizTb){
 442     TCanvas* ccs=new TCanvas("ccs","Cluster Size",1200,600);
 443     ccs->Divide(2,1);
 444     ccs->cd(1);
 445     gPad->SetLogz();
 446     gPad->SetRightMargin(0.12);
 447     hClSizAn->Draw("colz");
 448     hClSizAn->GetXaxis()->SetTitle("Drift Time (ns)");
 449     hClSizAn->GetYaxis()->SetTitle("Cluster Size - Anodes");
 450     ccs->cd(2);
 451     gPad->SetLogz();
 452     gPad->SetRightMargin(0.12);
 453     hClSizTb->Draw("colz");
 454     hClSizTb->GetXaxis()->SetTitle("Drift Time (ns)");
 455     hClSizTb->GetYaxis()->SetTitle("Cluster Size - Time Bins");
 456     ccs->Update();
 457   }
 458
 459   TH1F* htimR=(TH1F*)l->FindObject("hDrTimRP");
 460   TH1F* htimT=(TH1F*)l->FindObject("hDrTimTPAll");
 461   TH1F* htimTe=(TH1F*)l->FindObject("hDrTimTPExtra");
 462   TH1F* htimTne=(TH1F*)l->FindObject("hDrTimTPNoExtra");
 463   htimR->Rebin(4);
 464   htimT->Rebin(4);
 465   htimTe->Rebin(4);
 466   htimTne->Rebin(4);
 467   htimR->SetLineWidth(2);
 468   htimT->SetLineWidth(2);
 469   htimTe->SetLineWidth(2);
 470   htimTne->SetLineWidth(2);
 471
 472   TCanvas* ctim=new TCanvas("ctim","DriftTime",1400,600);
 473   ctim->Divide(2,1);
 474   ctim->cd(1);
 475   htimR->Draw();
 476   htimR->GetYaxis()->SetTitleOffset(1.2);
 477   htimR->GetXaxis()->SetTitle("Drift Time (ns)");
 478   htimR->GetYaxis()->SetTitle("RecPoints");
 479   ctim->cd(2);
 480   htimT->Draw();
 481   htimTe->SetLineColor(2);
 482   htimTe->Draw("same");
 483   htimTne->SetLineColor(4);
 484   htimTne->Draw("same");
 485   htimT->GetXaxis()->SetTitle("Drift Time (ns)");
 486   htimT->GetYaxis()->SetTitle("TrackPoints");
 487   htimT->GetYaxis()->SetTitleOffset(1.2);
 488   TLatex* ta=new TLatex(0.5,0.85,"All Clusters");
 489   ta->SetNDC();
 490   ta->SetTextColor(1);
 491   ta->Draw();
 492   TLatex* te=new TLatex(0.5,0.8,"Extra Clusters");
 493   te->SetNDC();
 494   te->SetTextColor(2);
 495   te->Draw();
 496   TLatex* tn=new TLatex(0.5,0.75,"Non-Extra Clusters");
 497   tn->SetNDC();
 498   tn->SetTextColor(4);
 499   tn->Draw();
 500   ctim->Update();
 501
 502   TCanvas* cdedx=new TCanvas("cdedx","dedx",1400,600);
 503   cdedx->Divide(3,1);
 504   cdedx->cd(1);
 505   gPad->SetLogz();
 506   hdedx3->Draw("col");
 507   hdedx3->GetXaxis()->SetTitle("P (GeV/c)");
 508   hdedx3->GetYaxis()->SetTitle("dE/dx (keV/300 #mum) Layer 3");
 509   hdedx3->GetYaxis()->SetTitleOffset(1.25);
 510   cdedx->cd(2);
 511   gPad->SetLogz();
 512   hdedx4->Draw("col");
 513   hdedx4->GetXaxis()->SetTitle("P (GeV/c)");
 514   hdedx4->GetYaxis()->SetTitle("dE/dx (keV/300 #mum) Layer 4");
 515   hdedx4->GetYaxis()->SetTitleOffset(1.25);
 516   cdedx->cd(3);
 517   gPad->SetLogz();
 518   hdedxmod->Draw("col");
 519   hdedxmod->GetXaxis()->SetTitle("SDD Module Id");
 520   hdedxmod->GetYaxis()->SetTitle("dE/dx (keV/300 #mum)");
 521   hdedxmod->GetYaxis()->SetTitleOffset(1.25);
 522   cdedx->Update();
 523
 524   printf("---- dE/dx vs.DriftTime ----\n");
 525   TCanvas* csig=new TCanvas("csig","dedx vs. DriftTime",1000,700);
 526   csig->Divide(4,2);
 527   TH1F* hSigTim[8];
 528   TGraphErrors* gmpv=new TGraphErrors(0);
 529   TGraphErrors* gsigg=new TGraphErrors(0);
 530   TGraphErrors* gsigl=new TGraphErrors(0);
 531   gmpv->SetTitle("");
 532   gsigg->SetTitle("");
 533   gsigl->SetTitle("");
 534   Int_t iPoint=0;
 535   TF1 *lfun = new TF1("LangausFun",LangausFun,50.,300.,4);
 536   for(Int_t it=0; it<8; it++){
 537     hSigTim[it]=(TH1F*)l->FindObject(Form("hSigTimeInt%d",it));
 538     csig->cd(it+1);
 539     hSigTim[it]->Draw();
 540     if(hSigTim[it]->GetEntries()>200){
 541       lfun->SetLineWidth(2);
 542       lfun->SetParameter(0,5.);
 543       lfun->SetParameter(1,80.);
 544       lfun->SetParameter(2,hSigTim[it]->GetEntries()/10.);
 545       lfun->SetParameter(3,10.);
 546       lfun->SetParLimits(3,0.,20);
 547
 548       hSigTim[it]->Fit("LangausFun","QLR");
 549       hSigTim[it]->GetXaxis()->SetTitle(Form("dE/dx, time interval %d",it+1));
 550       hSigTim[it]->GetYaxis()->SetTitle("Events");
 551       Float_t mpv=lfun->GetParameter(1);
 552       Float_t empv=lfun->GetParError(1);
 553       Float_t sig=lfun->GetParameter(3);
 554       Float_t esig=lfun->GetParError(3);
 555       Float_t sigl=lfun->GetParameter(0);
 556       Float_t esigl=lfun->GetParError(0);
 557       gmpv->SetPoint(iPoint,(Float_t)it,mpv);
 558       gmpv->SetPointError(iPoint,0.,empv);
 559       gsigg->SetPoint(iPoint,(Float_t)it,sig);
 560       gsigg->SetPointError(iPoint,0.,esig);
 561       gsigl->SetPoint(iPoint,(Float_t)it,sigl);
 562       gsigl->SetPointError(iPoint,0.,esigl);
 563       ++iPoint;
 564       gPad->Update();
 565       printf("Bin %d - MPV=%.3f  \t SigmaLandau=%.3f  \t SigmaGaus=%.3f\n",it,mpv,sigl,sig);
 566     }
 567   }
 568
 569   TCanvas* cpars=new TCanvas("cpars","Params",800,900);
 570   cpars->Divide(1,3,0.01,0.);
 571   cpars->cd(1);
 572   gPad->SetLeftMargin(0.14);
 573   gPad->SetFrameLineWidth(2);
 574   gPad->SetTickx();
 575   gPad->SetTicky();
 576   gmpv->SetMarkerStyle(20);
 577   //  gmpv->SetMinimum(0);
 578   //  gmpv->SetMaximum(120);
 579   gmpv->GetXaxis()->SetLimits(-0.2,6.8);
 580   gmpv->Draw("AP");
 581   //  gmpv->GetXaxis()->SetTitle("Drift Time interval number");
 582   gmpv->GetYaxis()->SetTitle("Landau MPV (keV)");
 583   gmpv->GetXaxis()->SetTitleSize(0.05);
 584   gmpv->GetYaxis()->SetTitleSize(0.05);
 585   gmpv->GetYaxis()->SetTitleOffset(1.2);
 586   cpars->cd(2);
 587   gPad->SetLeftMargin(0.14);
 588   gPad->SetFrameLineWidth(2);
 589   gPad->SetTickx();
 590   gPad->SetTicky();
 591   gsigl->SetMarkerStyle(20);
 592   gsigl->GetXaxis()->SetLimits(-0.2,6.8);
 593   gsigl->Draw("AP");
 594   //  gsigl->GetXaxis()->SetTitle("Drift Time interval number");
 595   gsigl->GetYaxis()->SetTitle("#sigma_{Landau} (keV)");
 596   gsigl->GetXaxis()->SetTitleSize(0.05);
 597   gsigl->GetYaxis()->SetTitleSize(0.05);
 598   gsigl->GetYaxis()->SetTitleOffset(1.2);
 599   cpars->cd(3);
 600   gPad->SetLeftMargin(0.14);
 601   gPad->SetFrameLineWidth(2);
 602   gPad->SetTickx();
 603   gPad->SetTicky();
 604   gsigg->SetMarkerStyle(20);
 605   gsigg->GetXaxis()->SetLimits(-0.2,6.8);
 606   gsigg->Draw("AP");
 607   gsigg->GetXaxis()->SetTitle("Drift Time interval number");
 608   gsigg->GetYaxis()->SetTitle("#sigma_{Gauss} (keV)");
 609   gsigg->GetXaxis()->SetTitleSize(0.05);
 610   gsigg->GetYaxis()->SetTitleSize(0.05);
 611   gsigg->GetYaxis()->SetTitleOffset(1.2);
 612   cpars->Update();
 613
 614 }
 615
 616
 617 Double_t LangausFun(Double_t *x, Double_t *par) {
 618
 619   //Fit parameters:
 620   //par[0]=Width (scale) parameter of Landau density
 621   //par[1]=Most Probable (MP, location) parameter of Landau density
 622   //par[2]=Total area (integral -inf to inf, normalization constant)
 623   //par[3]=Width (sigma) of convoluted Gaussian function
 624   //
 625   //In the Landau distribution (represented by the CERNLIB approximation),
 626   //the maximum is located at x=-0.22278298 with the location parameter=0.
 627   //This shift is corrected within this function, so that the actual
 628   //maximum is identical to the MP parameter.
 629
 630   // Numeric constants
 631   Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
 632   Double_t mpshift  = -0.22278298;       // Landau maximum location
 633
 634   // Control constants
 635   Double_t np = 100.0;      // number of convolution steps
 636   Double_t sc =   5.0;      // convolution extends to +-sc Gaussian sigmas
 637
 638   // Variables
 639   Double_t xx;
 640   Double_t mpc;
 641   Double_t fland;
 642   Double_t sum = 0.0;
 643   Double_t xlow,xupp;
 644   Double_t step;
 645   Double_t i;
 646
 647
 648   // MP shift correction
 649   mpc = par[1] - mpshift * par[0];
 650
 651   // Range of convolution integral
 652   xlow = x[0] - sc * par[3];
 653   xupp = x[0] + sc * par[3];
 654
 655   step = (xupp-xlow) / np;
 656
 657   // Convolution integral of Landau and Gaussian by sum
 658   for(i=1.0; i<=np/2; i++) {
 659     xx = xlow + (i-.5) * step;
 660     fland = TMath::Landau(xx,mpc,par[0]) / par[0];
 661     sum += fland * TMath::Gaus(x[0],xx,par[3]);
 662
 663     xx = xupp - (i-.5) * step;
 664     fland = TMath::Landau(xx,mpc,par[0]) / par[0];
 665     sum += fland * TMath::Gaus(x[0],xx,par[3]);
 666   }
 667
 668   return (par[2] * step * sum * invsq2pi / par[3]);
 669 }
 670