TRD/qaRec/AliTRDcheckDetector.cxx

   1 #include <TAxis.h>
   2 #include <TCanvas.h>
   3 #include <TFile.h>
   4 #include <TH1F.h>
   5 #include <TGaxis.h>
   6 #include <TGraph.h>
   7 #include <TLegend.h>
   8 #include <TMath.h>
   9 #include <TMap.h>
  10 #include <TObjArray.h>
  11 #include <TObject.h>
  12 #include <TObjString.h>
  13
  14 #include <TPad.h>
  15 #include <TProfile.h>
  16 #include <TProfile2D.h>
  17 #include <TROOT.h>
  18
  19 #include "AliLog.h"
  20 #include "AliTRDcluster.h"
  21 #include "AliESDHeader.h"
  22 #include "AliESDRun.h"
  23 #include "AliESDtrack.h"
  24 #include "AliTRDgeometry.h"
  25 #include "AliTRDpadPlane.h"
  26 #include "AliTRDSimParam.h"
  27 #include "AliTRDseedV1.h"
  28 #include "AliTRDtrackV1.h"
  29 #include "AliTRDtrackerV1.h"
  30 #include "AliTRDReconstructor.h"
  31 #include "AliTrackReference.h"
  32 #include "AliTrackPointArray.h"
  33 #include "AliTracker.h"
  34 #include "TTreeStream.h"
  35
  36 #include "AliTRDtrackInfo/AliTRDtrackInfo.h"
  37 #include "AliTRDtrackInfo/AliTRDeventInfo.h"
  38 #include "AliTRDcheckDetector.h"
  39
  40 #include <cstdio>
  41 #include <iostream>
  42
  43 ////////////////////////////////////////////////////////////////////////////
  44 //                                                                        //
  45 //  Reconstruction QA                                                     //
  46 //                                                                        //
  47 //  Task doing basic checks for tracking and detector performance         //
  48 //                                                                        //
  49 //  Authors:                                                              //
  50 //    Anton Andronic <A.Andronic@gsi.de>                                  //
  51 //    Alexandru Bercuci <A.Bercuci@gsi.de>                                //
  52 //    Markus Fasel <M.Fasel@gsi.de>                                       //
  53 //                                                                        //
  54 ////////////////////////////////////////////////////////////////////////////
  55
  56 //_______________________________________________________
  57 AliTRDcheckDetector::AliTRDcheckDetector():
  58   AliTRDrecoTask("DetChecker", "Basic Detector Checker")
  59   ,fEventInfo(0x0)
  60   ,fTriggerNames(0x0)
  61   ,fReconstructor(0x0)
  62   ,fGeo(0x0)
  63 {
  64   //
  65   // Default constructor
  66   //
  67   DefineInput(1,AliTRDeventInfo::Class());
  68   fReconstructor = new AliTRDReconstructor;
  69   fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
  70   fGeo = new AliTRDgeometry;
  71   InitFunctorList();
  72 }
  73
  74 //_______________________________________________________
  75 AliTRDcheckDetector::~AliTRDcheckDetector(){
  76   //
  77   // Destructor
  78   //
  79   if(fTriggerNames) delete fTriggerNames;
  80   delete fReconstructor;
  81   delete fGeo;
  82 }
  83
  84 //_______________________________________________________
  85 void AliTRDcheckDetector::ConnectInputData(Option_t *opt){
  86   //
  87   // Connect the Input data with the task
  88   //
  89   AliTRDrecoTask::ConnectInputData(opt);
  90   fEventInfo = dynamic_cast<AliTRDeventInfo *>(GetInputData(1));
  91 }
  92
  93 //_______________________________________________________
  94 void AliTRDcheckDetector::CreateOutputObjects(){
  95   //
  96   // Create Output Objects
  97   //
  98   OpenFile(0,"RECREATE");
  99   fContainer = Histos();
 100   if(!fTriggerNames) fTriggerNames = new TMap();
 101 }
 102
 103 //_______________________________________________________
 104 void AliTRDcheckDetector::Exec(Option_t *opt){
 105   //
 106   // Execution function
 107   // Filling TRD quality histos
 108   //
 109   if(!HasMCdata() && fEventInfo->GetEventHeader()->GetEventType() != 7) return; // For real data we select only physical events
 110   AliTRDrecoTask::Exec(opt);
 111   Int_t nTracks = 0;            // Count the number of tracks per event
 112   Int_t triggermask = fEventInfo->GetEventHeader()->GetTriggerMask();
 113   TString triggername =  fEventInfo->GetRunInfo()->GetFiredTriggerClasses(triggermask);
 114   if(fDebugLevel > 6)printf("Trigger cluster: %d, Trigger class: %s\n", triggermask, triggername.Data());
 115   dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTrigger))->Fill(triggermask);
 116   for(Int_t iti = 0; iti < fTracks->GetEntriesFast(); iti++){
 117     if(!fTracks->UncheckedAt(iti)) continue;
 118     AliTRDtrackInfo *fTrackInfo = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(iti));
 119     if(!fTrackInfo->GetTrack()) continue;
 120     nTracks++;
 121   }
 122   if(nTracks){
 123     dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTriggerTracks))->Fill(triggermask);
 124     dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNtracksEvent))->Fill(nTracks);
 125   }
 126   if(triggermask <= 20 && !fTriggerNames->FindObject(Form("%d", triggermask))){
 127     fTriggerNames->Add(new TObjString(Form("%d", triggermask)), new TObjString(triggername));
 128     // also set the label for both histograms
 129     TH1 *histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTriggerTracks));
 130     histo->GetXaxis()->SetBinLabel(histo->FindBin(triggermask), triggername);
 131     histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTrigger));
 132     histo->GetXaxis()->SetBinLabel(histo->FindBin(triggermask), triggername);
 133   }
 134   PostData(0, fContainer);
 135 }
 136
 137 //_______________________________________________________
 138 void AliTRDcheckDetector::Terminate(Option_t *){
 139   //
 140   // Terminate function
 141   //
 142 }
 143
 144 //_______________________________________________________
 145 Bool_t AliTRDcheckDetector::PostProcess(){
 146   //
 147   // Do Postprocessing (for the moment set the number of Reference histograms)
 148   //
 149
 150   TH1 * h = 0x0;
 151
 152   // Calculate of the trigger clusters purity
 153   h = dynamic_cast<TH1F *>(fContainer->FindObject("hEventsTrigger"));
 154   TH1F *h1 = dynamic_cast<TH1F *>(fContainer->FindObject("hEventsTriggerTracks"));
 155   h1->Divide(h);
 156   Float_t purities[20], val = 0;
 157   TString triggernames[20];
 158   Int_t nTriggerClasses = 0;
 159   for(Int_t ibin = 1; ibin <= h->GetNbinsX(); ibin++){
 160     if((val = h1->GetBinContent(ibin))){
 161       purities[nTriggerClasses] = val;
 162       triggernames[nTriggerClasses] = h1->GetXaxis()->GetBinLabel(ibin);
 163       nTriggerClasses++;
 164     }
 165   }
 166   h = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kTriggerPurity));
 167   TAxis *ax = h->GetXaxis();
 168   for(Int_t itrg = 0; itrg < nTriggerClasses; itrg++){
 169     h->Fill(itrg, purities[itrg]);
 170     ax->SetBinLabel(itrg+1, triggernames[itrg].Data());
 171   }
 172   ax->SetRangeUser(-0.5, nTriggerClasses+.5);
 173   h->GetYaxis()->SetRangeUser(0,1);
 174
 175   fNRefFigures = 14;
 176
 177   return kTRUE;
 178 }
 179
 180 //_______________________________________________________
 181 Bool_t AliTRDcheckDetector::GetRefFigure(Int_t ifig){
 182   //
 183   // Setting Reference Figures
 184   //
 185   switch(ifig){
 186   case kNclustersTrack:
 187     ((TH1F*)fContainer->FindObject("hNcls"))->Draw("pl");
 188     return kTRUE;
 189   case kNclustersTracklet:
 190     ((TH1F*)fContainer->FindObject("hNclTls"))->Draw("pc");
 191     return kTRUE;
 192   case kNtrackletsTrack:
 193     MakePlotNTracklets();
 194     return kTRUE;
 195   case kNtrackletsCross:
 196     if(!MakeBarPlot((TH1F*)fContainer->FindObject("hNtlsCross"), kRed)) break;
 197     return kTRUE;
 198   case kNtrackletsFindable:
 199     if(!MakeBarPlot((TH1F*)fContainer->FindObject("hNtlsFindable"), kGreen)) break;
 200     return kTRUE;
 201   case kNtracksEvent:
 202     ((TH1F*)fContainer->FindObject("hNtrks"))->Draw("pl");
 203     return kTRUE;
 204   case kNtracksSector:
 205     if(!MakeBarPlot((TH1F*)fContainer->FindObject("hNtrksSector"), kGreen)) break;
 206     return kTRUE;
 207   case kChi2:
 208     ((TH1F*)((TObjArray*)fContainer->FindObject("Chi2"))->At(0))->Draw("");
 209     return kTRUE;
 210   case kPH:
 211     MakePlotPulseHeight();
 212     return kTRUE;
 213   case kChargeCluster:
 214     ((TH1F*)fContainer->FindObject("hQcl"))->Draw("c");
 215     return kTRUE;
 216   case kChargeTracklet:
 217     ((TH1F*)fContainer->FindObject("hQtrklt"))->Draw("c");
 218     return kTRUE;
 219   case kNeventsTrigger:
 220     ((TH1F*)fContainer->FindObject("hEventsTrigger"))->Draw("");
 221     return kTRUE;
 222   case kNeventsTriggerTracks:
 223     ((TH1F*)fContainer->FindObject("hEventsTriggerTracks"))->Draw("");
 224     return kTRUE;
 225   case kTriggerPurity:
 226     if(!MakeBarPlot((TH1F*)fContainer->FindObject("hTriggerPurity"), kGreen)) break;
 227     break;
 228   default:
 229     break;
 230   }
 231   AliInfo(Form("Reference plot [%d] missing result", ifig));
 232   return kFALSE;
 233 }
 234
 235 //_______________________________________________________
 236 TObjArray *AliTRDcheckDetector::Histos(){
 237   //
 238   // Create QA histograms
 239   //
 240   if(fContainer) return fContainer;
 241
 242   fContainer = new TObjArray(20);
 243   //fContainer->SetOwner(kTRUE);
 244
 245   // Register Histograms
 246   TH1 * h = 0x0;
 247   if(!(h = (TH1F *)gROOT->FindObject("hNcls"))){
 248     h = new TH1F("hNcls", "N_{clusters} / track", 181, -0.5, 180.5);
 249     h->GetXaxis()->SetTitle("N_{clusters}");
 250     h->GetYaxis()->SetTitle("Entries");
 251   } else h->Reset();
 252   fContainer->AddAt(h, kNclustersTrack);
 253
 254   if(!(h = (TH1F *)gROOT->FindObject("hNclTls"))){
 255     h = new TH1F("hNclTls","N_{clusters} / tracklet", 51, -0.5, 50.5);
 256     h->GetXaxis()->SetTitle("N_{clusters}");
 257     h->GetYaxis()->SetTitle("Entries");
 258   } else h->Reset();
 259   fContainer->AddAt(h, kNclustersTracklet);
 260
 261   if(!(h = (TH1F *)gROOT->FindObject("hNtls"))){
 262     h = new TH1F("hNtls", "N_{tracklets} / track", AliTRDgeometry::kNlayer, 0.5, 6.5);
 263     h->GetXaxis()->SetTitle("N^{tracklet}");
 264     h->GetYaxis()->SetTitle("freq. [%]");
 265   } else h->Reset();
 266   fContainer->AddAt(h, kNtrackletsTrack);
 267
 268   if(!(h = (TH1F *)gROOT->FindObject("htlsSTA"))){
 269     h = new TH1F("hNtlsSTA", "N_{tracklets} / track (Stand Alone)", AliTRDgeometry::kNlayer, 0.5, 6.5);
 270     h->GetXaxis()->SetTitle("N^{tracklet}");
 271     h->GetYaxis()->SetTitle("freq. [%]");
 272   }
 273   fContainer->AddAt(h, kNtrackletsSTA);
 274
 275   if(!(h = (TH1F *)gROOT->FindObject("htlsBAR"))){
 276     h = new TH1F("hNtlsBAR", "N_{tracklets} / track (Barrel)", AliTRDgeometry::kNlayer, 0.5, 6.5);
 277     h->GetXaxis()->SetTitle("N^{tracklet}");
 278     h->GetYaxis()->SetTitle("freq. [%]");
 279   }
 280   fContainer->AddAt(h, kNtrackletsBAR);
 281
 282   //
 283   if(!(h = (TH1F *)gROOT->FindObject("hNtlsCross"))){
 284     h = new TH1F("hNtlsCross", "N_{tracklets}^{cross} / track", 7, -0.5, 6.5);
 285     h->GetXaxis()->SetTitle("n_{row cross}");
 286     h->GetYaxis()->SetTitle("freq. [%]");
 287   } else h->Reset();
 288   fContainer->AddAt(h, kNtrackletsCross);
 289
 290   if(!(h = (TH1F *)gROOT->FindObject("hNtlsFindable"))){
 291     h = new TH1F("hNtlsFindable", "Found/Findable Tracklets" , 101, -0.005, 1.005);
 292     h->GetXaxis()->SetTitle("r [a.u]");
 293     h->GetYaxis()->SetTitle("Entries");
 294   } else h->Reset();
 295   fContainer->AddAt(h, kNtrackletsFindable);
 296
 297   if(!(h = (TH1F *)gROOT->FindObject("hNtrks"))){
 298     h = new TH1F("hNtrks", "N_{tracks} / event", 100, 0, 100);
 299     h->GetXaxis()->SetTitle("N_{tracks}");
 300     h->GetYaxis()->SetTitle("Entries");
 301   } else h->Reset();
 302   fContainer->AddAt(h, kNtracksEvent);
 303
 304   if(!(h = (TH1F *)gROOT->FindObject("hNtrksSector"))){
 305     h = new TH1F("hNtrksSector", "N_{tracks} / sector", AliTRDgeometry::kNsector, -0.5, 17.5);
 306     h->GetXaxis()->SetTitle("sector");
 307     h->GetYaxis()->SetTitle("freq. [%]");
 308   } else h->Reset();
 309   fContainer->AddAt(h, kNtracksSector);
 310
 311   // <PH> histos
 312   TObjArray *arr = new TObjArray(2);
 313   arr->SetOwner(kTRUE);  arr->SetName("<PH>");
 314   fContainer->AddAt(arr, kPH);
 315   if(!(h = (TH1F *)gROOT->FindObject("hPHt"))){
 316     h = new TProfile("hPHt", "<PH>", 31, -0.5, 30.5);
 317     h->GetXaxis()->SetTitle("Time / 100ns");
 318     h->GetYaxis()->SetTitle("<PH> [a.u]");
 319   } else h->Reset();
 320   arr->AddAt(h, 0);
 321   if(!(h = (TH1F *)gROOT->FindObject("hPHx")))
 322     h = new TProfile("hPHx", "<PH>", 31, -0.08, 4.88);
 323   else h->Reset();
 324   arr->AddAt(h, 1);
 325
 326   // Chi2 histos
 327   arr = new TObjArray(2);
 328   arr->SetOwner(kTRUE); arr->SetName("Chi2");
 329   fContainer->AddAt(arr, kChi2);
 330   if(!(h = (TH1F *)gROOT->FindObject("hChi2")))
 331     h = new TH1F("hChi2", "#Chi2", 200, 0, 20);
 332   else h->Reset();
 333   arr->AddAt(h, 0);
 334   if(!(h = (TH1F *)gROOT->FindObject("hChi2n")))
 335     h = new TH1F("hChi2n", "Norm. Chi2 (tracklets)", 50, 0, 5);
 336   else h->Reset();
 337   arr->AddAt(h, 1);
 338
 339
 340   if(!(h = (TH1F *)gROOT->FindObject("hQcl"))){
 341     h = new TH1F("hQcl", "Q_{cluster}", 200, 0, 1200);
 342     h->GetXaxis()->SetTitle("Q_{cluster} [a.u.]");
 343     h->GetYaxis()->SetTitle("Entries");
 344   }else h->Reset();
 345   fContainer->AddAt(h, kChargeCluster);
 346
 347   if(!(h = (TH1F *)gROOT->FindObject("hQtrklt"))){
 348     h = new TH1F("hQtrklt", "Q_{tracklet}", 6000, 0, 6000);
 349     h->GetXaxis()->SetTitle("Q_{tracklet} [a.u.]");
 350     h->GetYaxis()->SetTitle("Entries");
 351   }else h->Reset();
 352   fContainer->AddAt(h, kChargeTracklet);
 353
 354
 355   if(!(h = (TH1F *)gROOT->FindObject("hEventsTrigger")))
 356     h = new TH1F("hEventsTrigger", "Trigger Class", 100, 0, 100);
 357   else h->Reset();
 358   fContainer->AddAt(h, kNeventsTrigger);
 359
 360   if(!(h = (TH1F *)gROOT->FindObject("hEventsTriggerTracks")))
 361     h = new TH1F("hEventsTriggerTracks", "Trigger Class (Tracks)", 100, 0, 100);
 362   else h->Reset();
 363   fContainer->AddAt(h, kNeventsTriggerTracks);
 364
 365   if(!(h = (TH1F *)gROOT->FindObject("hTriggerPurity"))){
 366     h = new TH1F("hTriggerPurity", "Trigger Purity", 10, -0.5, 9.5);
 367     h->GetXaxis()->SetTitle("Trigger Cluster");
 368     h->GetYaxis()->SetTitle("freq.");
 369   } else h->Reset();
 370   fContainer->AddAt(h, kTriggerPurity);
 371
 372   return fContainer;
 373 }
 374
 375 /*
 376 * Plotting Functions
 377 */
 378
 379 //_______________________________________________________
 380 TH1 *AliTRDcheckDetector::PlotNClustersTracklet(const AliTRDtrackV1 *track){
 381   //
 382   // Plot the mean number of clusters per tracklet
 383   //
 384   if(track) fTrack = track;
 385   if(!fTrack){
 386     AliWarning("No Track defined.");
 387     return 0x0;
 388   }
 389   TH1 *h = 0x0;
 390   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNclustersTracklet)))){
 391     AliWarning("No Histogram defined.");
 392     return 0x0;
 393   }
 394   AliTRDseedV1 *tracklet = 0x0;
 395   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
 396     if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
 397     h->Fill(tracklet->GetN2());
 398   }
 399   return h;
 400 }
 401
 402 //_______________________________________________________
 403 TH1 *AliTRDcheckDetector::PlotNClustersTrack(const AliTRDtrackV1 *track){
 404   //
 405   // Plot the number of clusters in one track
 406   //
 407   if(track) fTrack = track;
 408   if(!fTrack){
 409     AliWarning("No Track defined.");
 410     return 0x0;
 411   }
 412   TH1 *h = 0x0;
 413   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNclustersTrack)))){
 414     AliWarning("No Histogram defined.");
 415     return 0x0;
 416   }
 417
 418   Int_t nclusters = 0;
 419   AliTRDseedV1 *tracklet = 0x0;
 420   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
 421     if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
 422     nclusters += tracklet->GetN();
 423     if(fDebugLevel > 2){
 424       Int_t crossing = Int_t(tracklet->IsRowCross());
 425       Int_t detector = tracklet->GetDetector();
 426       Float_t theta = TMath::ATan(tracklet->GetZref(1));
 427       Float_t phi = TMath::ATan(tracklet->GetYref(1));
 428       Float_t momentum = 0.;
 429       Int_t pdg = 0;
 430       Int_t kinkIndex = fESD ? fESD->GetKinkIndex() : 0;
 431       UShort_t TPCncls = fESD ? fESD->GetTPCncls() : 0;
 432       if(fMC){
 433         if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
 434         pdg = fMC->GetPDG();
 435       }
 436       (*fDebugStream) << "NClustersTrack"
 437         << "Detector="  << detector
 438         << "crossing="  << crossing
 439         << "momentum="  << momentum
 440         << "pdg="                               << pdg
 441         << "theta="                     << theta
 442         << "phi="                               << phi
 443         << "kinkIndex=" << kinkIndex
 444         << "TPCncls="           << TPCncls
 445         << "nclusters=" << nclusters
 446         << "\n";
 447     }
 448   }
 449   h->Fill(nclusters);
 450   return h;
 451 }
 452
 453
 454 //_______________________________________________________
 455 TH1 *AliTRDcheckDetector::PlotNTrackletsTrack(const AliTRDtrackV1 *track){
 456   //
 457   // Plot the number of tracklets
 458   //
 459   if(track) fTrack = track;
 460   if(!fTrack){
 461     AliWarning("No Track defined.");
 462     return 0x0;
 463   }
 464   TH1 *h = 0x0, *hMethod = 0x0;
 465   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsTrack)))){
 466     AliWarning("No Histogram defined.");
 467     return 0x0;
 468   }
 469   Int_t status = fESD->GetStatus();
 470 /*  printf("in/out/refit/pid: TRD[%d|%d|%d|%d]\n", status &AliESDtrack::kTRDin ? 1 : 0, status &AliESDtrack::kTRDout ? 1 : 0, status &AliESDtrack::kTRDrefit ? 1 : 0, status &AliESDtrack::kTRDpid ? 1 : 0);*/
 471   if((status & AliESDtrack::kTRDin) != 0){
 472     // Full BarrelTrack
 473     if(!(hMethod = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsBAR))))
 474       AliWarning("Method: Barrel.  Histogram not processed!");
 475   } else {
 476     // Stand alone Track
 477     if(!(hMethod = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsSTA))))
 478       AliWarning("Method: StandAlone.  Histogram not processed!");
 479   }
 480   Int_t nTracklets = fTrack->GetNumberOfTracklets();
 481   h->Fill(nTracklets);
 482   hMethod->Fill(nTracklets);
 483   if(fDebugLevel > 3){
 484     if(nTracklets == 1){
 485       // If we have one Tracklet, check in which layer this happens
 486       Int_t layer = -1;
 487       AliTRDseedV1 *tracklet = 0x0;
 488       for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
 489         if((tracklet = fTrack->GetTracklet(il)) && tracklet->IsOK()){layer =  il; break;}
 490       }
 491       (*fDebugStream) << "NTrackletsTrack"
 492         << "Layer=" << layer
 493         << "\n";
 494     }
 495   }
 496   return h;
 497 }
 498
 499
 500 //_______________________________________________________
 501 TH1 *AliTRDcheckDetector::PlotNTrackletsRowCross(const AliTRDtrackV1 *track){
 502   //
 503   // Plot the number of tracklets
 504   //
 505   if(track) fTrack = track;
 506   if(!fTrack){
 507     AliWarning("No Track defined.");
 508     return 0x0;
 509   }
 510   TH1 *h = 0x0;
 511   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsCross)))){
 512     AliWarning("No Histogram defined.");
 513     return 0x0;
 514   }
 515
 516   Int_t ncross = 0;
 517   AliTRDseedV1 *tracklet = 0x0;
 518   for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
 519     if(!(tracklet = fTrack->GetTracklet(il)) || !tracklet->IsOK()) continue;
 520
 521     if(tracklet->IsRowCross()) ncross++;
 522   }
 523   h->Fill(ncross);
 524   return h;
 525 }
 526
 527 //_______________________________________________________
 528 TH1 *AliTRDcheckDetector::PlotFindableTracklets(const AliTRDtrackV1 *track){
 529   //
 530   // Plots the ratio of number of tracklets vs.
 531   // number of findable tracklets
 532   //
 533   // Findable tracklets are defined as track prolongation
 534   // to layer i does not hit the dead area +- epsilon
 535   //
 536   // In order to check whether tracklet hist active area in Layer i,
 537   // the track is refitted and the fitted position + an uncertainty
 538   // range is compared to the chamber border (also with a different
 539   // uncertainty)
 540   //
 541   // For the track fit two cases are distinguished:
 542   // If the track is a stand alone track (defined by the status bit
 543   // encoding, then the track is fitted with the tilted Rieman model
 544   // Otherwise the track is fitted with the Kalman fitter in two steps:
 545   // Since the track parameters are give at the outer point, we first
 546   // fit in direction inwards. Afterwards we fit again in direction outwards
 547   // to extrapolate the track to layers which are not reached by the track
 548   // For the Kalman model, the radial track points have to be shifted by
 549   // a distance epsilon in the direction that we want to fit
 550   //
 551   const Float_t epsilon = 0.01;   // dead area tolerance
 552   const Float_t epsilon_R = 1;    // shift in radial direction of the anode wire position (Kalman filter only)
 553   const Float_t delta_y = 0.7;    // Tolerance in the track position in y-direction
 554   const Float_t delta_z = 7.0;    // Tolerance in the track position in z-direction (Padlength)
 555   Double_t x_anode[AliTRDgeometry::kNlayer] = {300.2, 312.8, 325.4, 338.0, 350.6, 363.2}; // Take the default X0
 556
 557   if(track) fTrack = track;
 558   if(!fTrack){
 559     AliWarning("No Track defined.");
 560     return 0x0;
 561   }
 562   TH1 *h = 0x0;
 563   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsFindable)))){
 564     AliWarning("No Histogram defined.");
 565     return 0x0;
 566   }
 567   Int_t nFound = 0, nFindable = 0;
 568   Int_t stack = -1;
 569   Double_t ymin = 0., ymax = 0., zmin = 0., zmax = 0.;
 570   Double_t y = 0., z = 0.;
 571   AliTRDseedV1 *tracklet = 0x0;
 572   AliTRDpadPlane *pp;
 573   for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
 574     if((tracklet = fTrack->GetTracklet(il)) && tracklet->IsOK()){
 575       tracklet->SetReconstructor(fReconstructor);
 576       nFound++;
 577     }
 578   }
 579   // 2 Different cases:
 580   // 1st stand alone: here we cannot propagate, but be can do a Tilted Rieman Fit
 581   // 2nd barrel track: here we propagate the track to the layers
 582   AliTrackPoint points[6];
 583   Float_t xyz[3];
 584   memset(xyz, 0, sizeof(Float_t) * 3);
 585   if(((fESD->GetStatus() & AliESDtrack::kTRDout) > 0) && !((fESD->GetStatus() & AliESDtrack::kTRDin) > 0)){
 586     // stand alone track
 587     for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
 588       xyz[0] = x_anode[il];
 589       points[il].SetXYZ(xyz);
 590     }
 591     AliTRDtrackerV1::FitRiemanTilt(const_cast<AliTRDtrackV1 *>(fTrack), 0x0, kTRUE, 6, points);
 592   } else {
 593     // barrel track
 594     //
 595     // 2 Steps:
 596     // -> Kalman inwards
 597     // -> Kalman outwards
 598     AliTRDtrackV1 copy_track(*fTrack);  // Do Kalman on a (non-constant) copy of the track
 599     AliTrackPoint points_inward[6], points_outward[6];
 600     for(Int_t il = AliTRDgeometry::kNlayer; il--;){
 601       // In order to avoid complications in the Kalman filter if the track points have the same radial
 602       // position like the tracklets, we have to shift the radial postion of the anode wire by epsilon
 603       // in the direction we want to go
 604       // The track points have to be in reverse order for the Kalman Filter inwards
 605       xyz[0] = x_anode[AliTRDgeometry::kNlayer - il - 1] - epsilon_R;
 606       points_inward[il].SetXYZ(xyz);
 607       xyz[0] = x_anode[il] + epsilon_R;
 608       points_outward[il].SetXYZ(xyz);
 609     }
 610     /*for(Int_t ipt = 0; ipt < AliTRDgeometry::kNlayer; ipt++)
 611       printf("%d. X = %f\n", ipt, points[ipt].GetX());*/
 612     // Kalman inwards
 613     AliTRDtrackerV1::FitKalman(&copy_track, 0x0, kFALSE, 6, points_inward);
 614     memcpy(points, points_inward, sizeof(AliTrackPoint) * 6); // Preliminary store the inward results in the Array points
 615     // Kalman outwards
 616     AliTRDtrackerV1::FitKalman(&copy_track, 0x0, kTRUE, 6, points_inward);
 617     memcpy(points, points_outward, sizeof(AliTrackPoint) * AliTRDgeometry::kNlayer);
 618   }
 619   for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
 620     y = points[il].GetY();
 621     z = points[il].GetZ();
 622     if((stack = fGeo->GetStack(z, il)) < 0) continue; // Not findable
 623     pp = fGeo->GetPadPlane(il, stack);
 624     ymin = pp->GetCol0() + epsilon;
 625     ymax = pp->GetColEnd() - epsilon;
 626     zmin = pp->GetRowEnd() + epsilon;
 627     zmax = pp->GetRow0() - epsilon;
 628     // ignore y-crossing (material)
 629     if((z + delta_z > zmin && z - delta_z < zmax) && (y + delta_y > ymin && y - delta_y < ymax)) nFindable++;
 630       if(fDebugLevel > 3){
 631         Double_t pos_tracklet[2] = {tracklet ? tracklet->GetYfit(0) : 0, tracklet ? tracklet->GetZfit(0) : 0};
 632         Int_t hasTracklet = tracklet ? 1 : 0;
 633         (*fDebugStream)   << "FindableTracklets"
 634           << "layer="     << il
 635           << "ytracklet=" << pos_tracklet[0]
 636           << "ytrack="    << y
 637           << "ztracklet=" << pos_tracklet[1]
 638           << "ztrack="    << z
 639           << "tracklet="  << hasTracklet
 640           << "\n";
 641       }
 642   }
 643
 644   h->Fill(nFindable > 0 ? TMath::Min(nFound/static_cast<Double_t>(nFindable), 1.) : 1);
 645   if(fDebugLevel > 2) AliInfo(Form("Findable[Found]: %d[%d|%f]", nFindable, nFound, nFound/static_cast<Float_t>(nFindable > 0 ? nFindable : 1)));
 646   return h;
 647 }
 648
 649
 650 //_______________________________________________________
 651 TH1 *AliTRDcheckDetector::PlotChi2(const AliTRDtrackV1 *track){
 652   //
 653   // Plot the chi2 of the track
 654   //
 655   if(track) fTrack = track;
 656   if(!fTrack){
 657     AliWarning("No Track defined.");
 658     return 0x0;
 659   }
 660   TH1 *h = 0x0;
 661   if(!(h = dynamic_cast<TH1F *>(((TObjArray*)(fContainer->At(kChi2)))->At(0)))){
 662     AliWarning("No Histogram defined.");
 663     return 0x0;
 664   }
 665   h->Fill(fTrack->GetChi2());
 666   return h;
 667 }
 668
 669 //_______________________________________________________
 670 TH1 *AliTRDcheckDetector::PlotChi2Norm(const AliTRDtrackV1 *track){
 671   //
 672   // Plot the chi2 of the track
 673   //
 674   if(track) fTrack = track;
 675   if(!fTrack){
 676     AliWarning("No Track defined.");
 677     return 0x0;
 678   }
 679   TH1 *h = 0x0;
 680   if(!(h = dynamic_cast<TH1F *>(((TObjArray*)(fContainer->At(kChi2)))->At(1)))){
 681     AliWarning("No Histogram defined.");
 682     return 0x0;
 683   }
 684   Int_t nTracklets = 0;
 685   AliTRDseedV1 *tracklet = 0x0;
 686   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
 687     if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
 688     nTracklets++;
 689   }
 690   if(!nTracklets) return 0x0;
 691   h->Fill(fTrack->GetChi2()/nTracklets);
 692   return h;
 693 }
 694
 695
 696 //_______________________________________________________
 697 TH1 *AliTRDcheckDetector::PlotPHt(const AliTRDtrackV1 *track){
 698   //
 699   // Plot the average pulse height
 700   //
 701   if(track) fTrack = track;
 702   if(!fTrack){
 703     AliWarning("No Track defined.");
 704     return 0x0;
 705   }
 706   TProfile *h = 0x0;
 707   if(!(h = dynamic_cast<TProfile *>(((TObjArray*)(fContainer->At(kPH)))->At(0)))){
 708     AliWarning("No Histogram defined.");
 709     return 0x0;
 710   }
 711   AliTRDseedV1 *tracklet = 0x0;
 712   AliTRDcluster *c = 0x0;
 713   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
 714     if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK())continue;
 715     Int_t crossing = Int_t(tracklet->IsRowCross());
 716     Int_t detector = tracklet->GetDetector();
 717     tracklet->ResetClusterIter();
 718     while((c = tracklet->NextCluster())){
 719       if(!c->IsInChamber()) continue;
 720       Int_t localtime        = c->GetLocalTimeBin();
 721       Double_t absolute_charge = TMath::Abs(c->GetQ());
 722       h->Fill(localtime, absolute_charge);
 723       if(fDebugLevel > 3){
 724         Double_t distance[2];
 725         GetDistanceToTracklet(distance, tracklet, c);
 726         Float_t theta = TMath::ATan(tracklet->GetZref(1));
 727         Float_t phi = TMath::ATan(tracklet->GetYref(1));
 728         Float_t momentum = 0.;
 729         Int_t pdg = 0;
 730         Int_t kinkIndex = fESD ? fESD->GetKinkIndex() : 0;
 731         UShort_t TPCncls = fESD ? fESD->GetTPCncls() : 0;
 732         if(fMC){
 733           if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
 734           pdg = fMC->GetPDG();
 735         }
 736         (*fDebugStream) << "PHt"
 737           << "Detector="        << detector
 738           << "crossing="        << crossing
 739           << "Timebin="         << localtime
 740           << "Charge="          << absolute_charge
 741           << "momentum="        << momentum
 742           << "pdg="                             << pdg
 743           << "theta="                   << theta
 744           << "phi="                             << phi
 745           << "kinkIndex="       << kinkIndex
 746           << "TPCncls="         << TPCncls
 747           << "dy="        << distance[0]
 748           << "dz="        << distance[1]
 749           << "c.="        << c
 750           << "\n";
 751       }
 752     }
 753   }
 754   return h;
 755 }
 756
 757 //_______________________________________________________
 758 TH1 *AliTRDcheckDetector::PlotPHx(const AliTRDtrackV1 *track){
 759   //
 760   // Plots the average pulse height vs the distance from the anode wire
 761   // (plus const anode wire offset)
 762   //
 763   if(track) fTrack = track;
 764   if(!fTrack){
 765     AliWarning("No Track defined.");
 766     return 0x0;
 767   }
 768   TProfile *h = 0x0;
 769   if(!(h = dynamic_cast<TProfile *>(((TObjArray*)(fContainer->At(kPH)))->At(1)))){
 770     AliWarning("No Histogram defined.");
 771     return 0x0;
 772   }
 773   Float_t offset = .5*AliTRDgeometry::CamHght();
 774   AliTRDseedV1 *tracklet = 0x0;
 775   AliTRDcluster *c = 0x0;
 776   Double_t distance = 0;
 777   Double_t x, y;
 778   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
 779     if(!(tracklet = fTrack->GetTracklet(itl)) || !(tracklet->IsOK())) continue;
 780     tracklet->ResetClusterIter();
 781     while((c = tracklet->NextCluster())){
 782       if(!c->IsInChamber()) continue;
 783       tracklet->GetClusterXY(c, x, y);
 784       distance = tracklet->GetX0() - (c->GetX() + 0.3) + offset;
 785       h->Fill(distance, TMath::Abs(c->GetQ()));
 786     }
 787   }
 788   return h;
 789 }
 790
 791 //_______________________________________________________
 792 TH1 *AliTRDcheckDetector::PlotChargeCluster(const AliTRDtrackV1 *track){
 793   //
 794   // Plot the cluster charge
 795   //
 796   if(track) fTrack = track;
 797   if(!fTrack){
 798     AliWarning("No Track defined.");
 799     return 0x0;
 800   }
 801   TH1 *h = 0x0;
 802   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kChargeCluster)))){
 803     AliWarning("No Histogram defined.");
 804     return 0x0;
 805   }
 806   AliTRDseedV1 *tracklet = 0x0;
 807   AliTRDcluster *c = 0x0;
 808   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
 809     if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK())continue;
 810     for(Int_t itime = 0; itime < AliTRDtrackerV1::GetNTimeBins(); itime++){
 811       if(!(c = tracklet->GetClusters(itime))) continue;
 812       h->Fill(c->GetQ());
 813     }
 814   }
 815   return h;
 816 }
 817
 818 //_______________________________________________________
 819 TH1 *AliTRDcheckDetector::PlotChargeTracklet(const AliTRDtrackV1 *track){
 820   //
 821   // Plot the charge deposit per chamber
 822   //
 823   if(track) fTrack = track;
 824   if(!fTrack){
 825     AliWarning("No Track defined.");
 826     return 0x0;
 827   }
 828   TH1 *h = 0x0;
 829   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kChargeTracklet)))){
 830     AliWarning("No Histogram defined.");
 831     return 0x0;
 832   }
 833   AliTRDseedV1 *tracklet = 0x0;
 834   AliTRDcluster *c = 0x0;
 835   Double_t Qtot = 0;
 836   Int_t nTracklets =fTrack->GetNumberOfTracklets();
 837   for(Int_t itl = 0x0; itl < AliTRDgeometry::kNlayer; itl++){
 838     if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
 839     Qtot = 0.;
 840     for(Int_t ic = AliTRDseedV1::kNclusters; ic--;){
 841       if(!(c = tracklet->GetClusters(ic))) continue;
 842       Qtot += TMath::Abs(c->GetQ());
 843     }
 844     h->Fill(Qtot);
 845     if(fDebugLevel > 3){
 846       Int_t crossing = (Int_t)tracklet->IsRowCross();
 847       Int_t detector = tracklet->GetDetector();
 848       Float_t theta = TMath::ATan(tracklet->GetZfit(1));
 849       Float_t phi = TMath::ATan(tracklet->GetYfit(1));
 850       Float_t momentum = 0.;
 851       Int_t pdg = 0;
 852       Int_t kinkIndex = fESD ? fESD->GetKinkIndex() : 0;
 853       UShort_t TPCncls = fESD ? fESD->GetTPCncls() : 0;
 854       if(fMC){
 855               if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
 856         pdg = fMC->GetPDG();
 857       }
 858       (*fDebugStream) << "ChargeTracklet"
 859         << "Detector="  << detector
 860         << "crossing="  << crossing
 861         << "momentum="  << momentum
 862         << "nTracklets="<< nTracklets
 863         << "pdg="                               << pdg
 864         << "theta="                     << theta
 865         << "phi="                               << phi
 866         << "kinkIndex=" << kinkIndex
 867         << "TPCncls="           << TPCncls
 868         << "QT="        << Qtot
 869         << "\n";
 870     }
 871   }
 872   return h;
 873 }
 874
 875 //_______________________________________________________
 876 TH1 *AliTRDcheckDetector::PlotNTracksSector(const AliTRDtrackV1 *track){
 877   //
 878   // Plot the number of tracks per Sector
 879   //
 880   if(track) fTrack = track;
 881   if(!fTrack){
 882     AliWarning("No Track defined.");
 883     return 0x0;
 884   }
 885   TH1 *h = 0x0;
 886   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtracksSector)))){
 887     AliWarning("No Histogram defined.");
 888     return 0x0;
 889   }
 890
 891   // TODO we should compare with
 892   // sector = Int_t(track->GetAlpha() / AliTRDgeometry::GetAlpha());
 893
 894   AliTRDseedV1 *tracklet = 0x0;
 895   Int_t sector = -1;
 896   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
 897     if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
 898     sector = static_cast<Int_t>(tracklet->GetDetector()/AliTRDgeometry::kNdets);
 899     break;
 900   }
 901   h->Fill(sector);
 902   return h;
 903 }
 904
 905
 906 //________________________________________________________
 907 void AliTRDcheckDetector::SetRecoParam(AliTRDrecoParam *r)
 908 {
 909
 910   fReconstructor->SetRecoParam(r);
 911 }
 912
 913 //________________________________________________________
 914 void AliTRDcheckDetector::GetDistanceToTracklet(Double_t *dist, AliTRDseedV1 *tracklet, AliTRDcluster *c)
 915 {
 916   Float_t x = c->GetX();
 917   dist[0] = c->GetY() - tracklet->GetYat(x);
 918   dist[1] = c->GetZ() - tracklet->GetZat(x);
 919 }
 920
 921 //________________________________________________________
 922 void AliTRDcheckDetector::MakePlotNTracklets(){
 923   //
 924   // Make nice bar plot of the number of tracklets in each method
 925   //
 926   TH1F *hBAR = (TH1F *)fContainer->FindObject("hNtlsBAR");
 927   TH1F *hSTA = (TH1F *)fContainer->FindObject("hNtlsSTA");
 928   TH1F *hCON = (TH1F *)fContainer->FindObject("hNtls");
 929
 930   hBAR->Scale(100./hCON->Integral());
 931   hBAR->SetFillColor(kRed);
 932   hBAR->SetBarWidth(0.2);
 933   hBAR->SetBarOffset(0.2);
 934   hBAR->SetTitle("");
 935   hBAR->SetStats(kFALSE);
 936   hBAR->GetYaxis()->SetRangeUser(0.,40.);
 937   hBAR->GetYaxis()->SetTitleOffset(1.2);
 938   hBAR->Draw("bar1");
 939
 940   hSTA->Scale(100./hCON->Integral());
 941   hSTA->SetFillColor(kBlue);
 942   hSTA->SetBarWidth(0.2);
 943   hSTA->SetBarOffset(0.4);
 944   hSTA->SetTitle("");
 945   hSTA->SetStats(kFALSE);
 946   hSTA->GetYaxis()->SetRangeUser(0.,40.);
 947   hSTA->GetYaxis()->SetTitleOffset(1.2);
 948   hSTA->Draw("bar1same");
 949
 950   hCON->Scale(100./hCON->Integral());
 951   hCON->SetFillColor(kGreen);
 952   hCON->SetBarWidth(0.2);
 953   hCON->SetBarOffset(0.6);
 954   hCON->SetStats(kFALSE);
 955   hCON->GetYaxis()->SetRangeUser(0.,40.);
 956   hCON->GetYaxis()->SetTitleOffset(1.2);
 957   hCON->Draw("bar1same");
 958
 959   TLegend *leg = new TLegend(0.6, 0.75, 0.89, 0.89);
 960   leg->AddEntry(hBAR, "Barrel", "f");
 961   leg->AddEntry(hSTA, "Stand Alone", "f");
 962   leg->AddEntry(hCON, "Convoluted", "f");
 963
 964   leg->Draw();
 965   gPad->Update();
 966 }
 967
 968 //________________________________________________________
 969 void AliTRDcheckDetector::MakePlotPulseHeight(){
 970   //
 971   // Create Plot of the Pluse Height Spectrum
 972   //
 973   TH1 *h, *h1, *h2;
 974   TObjArray *arr = (TObjArray*)fContainer->FindObject("<PH>");
 975   h = (TH1F*)arr->At(0);
 976   h->SetMarkerStyle(24);
 977   h->SetMarkerColor(kBlack);
 978   h->SetLineColor(kBlack);
 979   h->Draw("e1");
 980 //   copy the second histogram in a new one with the same x-dimension as the phs with respect to time
 981   h1 = (TH1F *)arr->At(1);
 982   h2 = new TH1F("hphs1","Average PH", 31, -0.5, 30.5);
 983   for(Int_t ibin = h1->GetXaxis()->GetFirst(); ibin < h1->GetNbinsX(); ibin++)
 984     h2->SetBinContent(ibin, h1->GetBinContent(ibin));
 985   h2->SetMarkerStyle(22);
 986   h2->SetMarkerColor(kBlue);
 987   h2->SetLineColor(kBlue);
 988   h2->Draw("e1same");
 989   gPad->Update();
 990 //   create axis according to the histogram dimensions of the original second histogram
 991   TGaxis *axis = new TGaxis(gPad->GetUxmin(),
 992                     gPad->GetUymax(),
 993                     gPad->GetUxmax(),
 994                     gPad->GetUymax(),
 995                     -0.08, 4.88, 510,"-L");
 996   axis->SetLineColor(kBlue);
 997   axis->SetLabelColor(kBlue);
 998   axis->SetTextColor(kBlue);
 999   axis->SetTitle("x_{0}-x_{c} [cm]");
1000   axis->Draw();
1001 }
1002
1003 //________________________________________________________
1004 Bool_t AliTRDcheckDetector::MakeBarPlot(TH1 *histo, Int_t color){
1005   //
1006   // Draw nice bar plots
1007   //
1008   if(!histo->GetEntries()) return kFALSE;
1009   histo->Scale(100./histo->Integral());
1010   histo->SetFillColor(color);
1011   histo->SetBarOffset(.2);
1012   histo->SetBarWidth(.6);
1013   histo->Draw("bar1");
1014   return kTRUE;
1015 }