TRD/qaRec/AliTRDcheckDET.cxx

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