TRD/Macros/AliTRDanalyseRecPoints.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 #include <TRD/AliTRDcluster.h>
   3 #include <AliRunLoader.h>
   4 #include <TCanvas.h>
   5 #include <TF1.h>
   6 #include <TH1.h>
   7 #include <TH2.h>
   8 #include <TLine.h>
   9 #include <TMath.h>
  10 #include <TObjArray.h>
  11 #include <TPad.h>
  12 #include <TROOT.h>
  13 #include <TStyle.h>
  14 #include <TSystem.h>
  15 #include <TTree.h>
  16 #endif
  17
  18 //#define EVENT_BY_EVENT      // Show amplitude event-by-event
  19 //#define EVENT_BY_EVENT2     // Show tracks event-by-event
  20 //#define PRINT_OTHER         // Display various additional histograms
  21
  22 // Number of noisy detectors
  23 const Int_t numNoisyDets = 10;
  24
  25 // Holds the detector number of detectors that produce too much noise and
  26 // are therefore excluded from the calculation
  27 const Int_t exDets[numNoisyDets] = {25, 511, 515, 538, 523, 249, 294, 27, 24, 247};
  28
  29 Bool_t isBadCluster(Int_t det, Int_t Q, Int_t PadCol, Int_t PadRow)
  30 {
  31   if (Q < 20) return 1;
  32
  33   // Exclude bad cols
  34   if (PadCol == 142 || PadCol == 136 || PadCol == 66 || PadCol == 65 || PadCol <= 2) return 1;
  35   //if (PadCol == 142 || PadCol == 66 || PadCol == 65 || PadCol <= 2) return 1;
  36   //if (PadCol == 142 || PadCol <= 2) return 1;
  37
  38   // Exclude noisy detectors
  39   for (Int_t indDet = 0; indDet < numNoisyDets; indDet++)
  40   {
  41     if (det == exDets[indDet])
  42     {
  43       // Exclude several regions of bad detectors
  44       if ((PadCol >= 134 && PadRow <= 5) || (PadCol == 141 && PadRow >= 10)) return 1;
  45       //if (PadRow == 15 || (PadCol >= 134 && PadRow <= 11) || (PadCol >= 119 && PadCol <= 123 && PadRow == 11)) return 1;
  46
  47       break;
  48     }
  49   }
  50
  51   // Cluster seems to be clean
  52   return 0;
  53 }
  54
  55
  56 //=====================================
  57
  58 void AliTRDanalyseRecPoints (const char *filename="galice.root")
  59 {
  60   gROOT->SetStyle("Plain");
  61   gStyle->SetPalette(1);
  62
  63   // open run loader and load gAlice, kinematics and header
  64   AliRunLoader* runLoader = AliRunLoader::Open(filename);
  65   if (!runLoader)
  66   {
  67     printf("Error: readKine:\nGetting run loader from file \"%s\" failed", filename);
  68     return;
  69   }
  70
  71   runLoader->LoadHeader();
  72   runLoader->LoadRecPoints("TRD");
  73   TObjArray *module = new TObjArray();
  74
  75   runLoader->CdGAFile();
  76
  77
  78   // Flag indicating whether the signal comes from a noisy detector
  79   Bool_t fromNoisyDet = 0;
  80
  81   // Counts the number of clusters per event
  82   Int_t clusters = 0;
  83
  84   // Threshold for the amplitude
  85   const Int_t kQThreshold = 30;
  86
  87   // Selection criterions for clusters
  88   const Int_t kThresholdLower = 30;
  89   const Int_t kThresholdUpper = 500;
  90
  91   // Number of events passing the criterion
  92   Int_t numPassed = 0;
  93
  94   // hist
  95 #ifdef PRINT_OTHER
  96   TH1D *mNCls = new TH1D("NCls", ";number of clusters", 400, 0, 400);
  97   TH1D *mResY = new TH1D("resY", ";sigma Y (cm)", 100, 0, 1);
  98   TH1D *mResZ = new TH1D("resZ", ";sigma Z (cm)", 100, 1, 3);
  99
 100   TH1D *mNPads = new TH1D("NPads", ";number of pads", 10, -0.5, 9.5);
 101   TH1D *mCenter = new TH1D("center", ";center", 100, -1, 1);
 102   TH1D *mTime = new TH1D("timeBin", ";time bin", 25, -0.5, 24.5);
 103   TH1D *mPlane = new TH1D("plane", ";plane", 6, -0.5, 5.5);
 104
 105   TH1D *mX = new TH1D("x", "; X (cm)", 100, 0, 3.5);
 106   TH1D *mY = new TH1D("y", "; y / chamber width (%)", 100, -1, 1);
 107   TH1D *mZ = new TH1D("z", "; Z (cm)", 200, -400, 400);
 108
 109   TH2D *mPlaneY = new TH2D("planeY", ";plane;Y (cm);#Clusters", 6, -0.5, 5.5, 100, -60, 60);
 110   TH2D *mYX = new TH2D("yx", ";Y (cm);X (cm);#Clusters", 100, -60, 60, 120, 280, 400);
 111 #endif
 112
 113   TH1D *mChrg = new TH1D("chrg", ";Charge(no pad-filter)", 200, 0, 200);
 114   TH1D *mChrg2 = new TH1D("chrg2", ";Charge(pads<9)", 200, 0, 200);
 115   TH1D *mChrg3 = new TH1D("chrg3", ";Charge(2<pads<9)", 200, 0, 200);
 116
 117   TH2D *mQvsPads = new TH2D("QvsPads", ";Number of pads;Q;#Clusters", 15, -0.5, 14.5, 100, 0, 1000);
 118
 119   // There are 540 detectors
 120   TH2D *mQvsDetector = new TH2D("QvsDet", ";Number of detector;Q;#Clusters", 540, -0.5, 539.5, 100, 0, 1000);
 121
 122   // Supermoduleid=int(detectorNum / 30)
 123   TH2D *mQvsSModId = new TH2D("QvsSModId", ";Supermodule id;Q;#Clusters", 18, -0.5, 17.5, 100, 0, 1000);
 124   // Histograms for the already installed supermodules
 125   TH1D *mSMod0 = new TH1D("SMod0", ";Charge", 200, 0, 200);
 126   TH1D *mSMod8 = new TH1D("SMod8", ";Charge", 200, 0, 200);
 127   TH1D *mSMod9 = new TH1D("SMod9", ";Charge", 200, 0, 200);
 128   TH1D *mSMod17 = new TH1D("SMod17", ";Charge", 200, 0, 200);
 129
 130   // Index + 2 corresponds to the number of activated pads
 131   TH1D *mPads[14];
 132   char name[8];
 133   for (Int_t ind = 0; ind < 14; ind++)
 134   {
 135     sprintf(name,"%dPads",ind + 2);
 136     mPads[ind] = new TH1D(name, ";Charge", 200, 0, 200);
 137   }
 138
 139   // Amplitude distribution will be displayed event by event for the first 100 events
 140   TH1D *mChrgEBE = new TH1D("chrgEBE", ";Charge", 200, 0, 200);
 141
 142   // Keep track of the number of clusters per event
 143   TH1D *mClusters = new TH1D("Clusters", ";Number of clusters per event", 300, -0.5, 299.5);
 144
 145   // Keep track of the number of clusters per chamber with Q < kQThreshold
 146   char clusTitle[55];
 147   sprintf(clusTitle,";Number of clusters with Q < %d vs detector", kQThreshold);
 148   TH1D *mClusCham = new TH1D("ClusCham", clusTitle, 540, -0.5, 539.5);
 149   // Same but this time noisy detectors excluded
 150   sprintf(clusTitle,";Number of clusters with Q < %d vs detector", kQThreshold);
 151   TH1D *mClusChamEx = new TH1D("ClusChamEx", clusTitle, 540, -0.5, 539.5);
 152
 153   TH2D *mColRow = new TH2D("ColRow", ";Column;Row;#Clusters", 144, -0.5, 143.5, 18, -0.5, 17.5);
 154   TH2D *mColRowFiltered = new TH2D("ColRowFiltered", ";Column;Row;#Clusters", 144, -0.5, 143.5, 18, -0.5, 17.5);
 155
 156   TH2D *mColTime = new TH2D("ColTime", ";TimeBin;Column", 18, -0.5, 17.5, 144, -0.5, 143.5);
 157   TH2D *mColTimeFiltered = new TH2D("ColTimeFiltered", ";TimeBin;Column;#Clusters", 18, -0.5, 17.5, 144, -0.5, 143.5);
 158
 159   // YX histograms for the already installed supermodules
 160   TH2D *mTrSMod0 = new TH2D("TrSMod0", ";Col;time+layer*20;#Clusters", 144, -0.5, 143.5, 120, -0.5, 119.5);
 161   TH2D *mTrSMod8 = new TH2D("TrSMod8", ";Col;time+layer*20;#Clusters", 144, -0.5, 143.5, 120, -0.5, 119.5);
 162   TH2D *mTrSMod9 = new TH2D("TrSMod9", ";Col;time+layer*20;#Clusters", 144, -0.5, 143.5, 120, -0.5, 119.5);
 163   TH2D *mTrSMod17 = new TH2D("TrSMod17", ";Col;time+layer*20;#Clusters", 144, -0.5, 143.5, 120, -0.5, 119.5);
 164
 165
 166 #ifdef EVENT_BY_EVENT
 167   TCanvas* eventc = new TCanvas("eventc","Amplitude distribution for event 0",700,500);
 168   char title[40];
 169 #endif
 170
 171 #ifdef EVENT_BY_EVENT2
 172   TCanvas* sTrc = new TCanvas("sTrc","Track -  Supermodules",800,800);
 173   sTrc->SetLogz();
 174   sTrc->Divide(2,2);
 175   char title2[40];
 176 #endif
 177
 178
 179   int nEvents = runLoader->GetNumberOfEvents();
 180
 181   AliTRDcluster *cls = 0;
 182
 183   for(Int_t ev = 0; ev < nEvents - 1; ev++)
 184   {
 185     clusters = 0;
 186
 187     TTree *tree = runLoader->GetTreeR("TRD", 0);
 188     tree->SetBranchAddress("TRDcluster", &module);
 189
 190     int N = tree->GetEntries();
 191
 192     // Check number of clusters
 193     for(Int_t ind = 0; ind < N; ind++)
 194     {
 195       tree->GetEntry(ind);
 196       Int_t m = module->GetEntries();
 197
 198       for (Int_t j = 0; j < m; j++)
 199       {
 200         if (cls != 0) delete cls;
 201         cls = (AliTRDcluster*)module->At(j);
 202         if (!isBadCluster(cls->GetDetector(), (Int_t)cls->GetQ(), cls->GetPadCol(), cls->GetPadRow())) clusters++;
 203         else
 204               {
 205           // Uncomment the following lines AND comment the next if-clause to display the histograms for the
 206           // EXCLUDED clusters
 207
 208                 //mQvsPads->Fill(cls->GetNPads(),cls->GetQ());
 209                 //mQvsDetector->Fill(cls->GetDetector(),cls->GetQ());
 210                 //mQvsSModId->Fill((Int_t)(cls->GetDetector() / 30.),cls->GetQ());
 211                 //mColTime->Fill(cls->GetLocalTimeBin(), cls->GetPadCol());
 212           //mColRow->Fill(cls->GetPadCol(), cls->GetPadRow());
 213                 //mClusCham->Fill(cls->GetDetector());
 214           //mYX->Fill(cls->GetY(), cls->GetX());
 215         }
 216       }
 217     }
 218
 219     // Process event, if number of clusters is in the correct range
 220     if (clusters >= kThresholdLower && clusters <= kThresholdUpper)
 221     {
 222       clusters = 0;
 223       numPassed++;
 224
 225       for (Int_t i = 0; i < N; i++)
 226       {
 227               tree->GetEntry(i);
 228               int m = module->GetEntries();
 229
 230 #ifdef PRINT_OTHER
 231               mNCls->Fill(m);
 232 #endif
 233
 234               for(Int_t j = 0; j < m; j++)
 235               {
 236           fromNoisyDet = 0;
 237                 clusters++;
 238
 239           if (cls != 0) delete cls;
 240                 cls = (AliTRDcluster*)module->At(j);
 241
 242                 // Filter the clusters of the events, too
 243           if (cls->GetQ() < 20) continue;
 244
 245                 mColTime->Fill(cls->GetLocalTimeBin(), cls->GetPadCol());
 246
 247           if (cls->GetQ() <= kQThreshold)
 248           {
 249             mClusCham->Fill(cls->GetDetector());
 250                 }
 251
 252           // Exclude noisy detectors
 253           for (Int_t indDet = 0; indDet < numNoisyDets; indDet++)
 254                 {
 255                   if (cls->GetDetector() == exDets[indDet])
 256                   {
 257                   fromNoisyDet = 1;
 258                     break;
 259                   }
 260                 }
 261
 262           // Exclude I
 263           if (cls->GetPadCol() == 142 || cls->GetPadCol() == 136 || cls->GetPadCol() == 66 || cls->GetPadCol() == 65 || cls->GetPadCol() <= 2) continue;
 264           //if (cls->GetPadCol() == 142 || cls->GetPadCol() == 66 || cls->GetPadCol() == 65 || cls->GetPadCol() <= 2) continue;
 265           //if (cls->GetPadCol() == 142 || cls->GetPadCol() <= 2) continue;
 266
 267                 if (fromNoisyDet)
 268                 {
 269                   mColRow->Fill(cls->GetPadCol(), cls->GetPadRow());
 270
 271                   // Exclude II
 272                   //if (cls->GetPadRow() == 15 || (cls->GetPadCol() >= 134 && cls->GetPadRow() <= 11) || (cls->GetPadCol() >= 119 && cls->GetPadCol() <= 123 && cls->GetPadRow() == 11)) continue;
 273
 274             if ((cls->GetPadCol() >= 134 && cls->GetPadRow() <= 5) || (cls->GetPadCol() == 141 && cls->GetPadRow() >= 10)) continue;
 275
 276             mColRowFiltered->Fill(cls->GetPadCol(), cls->GetPadRow());
 277                 }
 278
 279           mColTimeFiltered->Fill(cls->GetLocalTimeBin(), cls->GetPadCol());
 280
 281           if (cls->GetQ() <= kQThreshold)
 282           {
 283             mClusChamEx->Fill(cls->GetDetector());
 284                 }
 285
 286 #ifdef PRINT_OTHER
 287                 mX->Fill(cls->GetX());
 288                 mZ->Fill(cls->GetZ());
 289           mYX->Fill(cls->GetY(), cls->GetX());
 290                 mResY->Fill(TMath::Sqrt(cls->GetSigmaY2()));
 291                 mResZ->Fill(TMath::Sqrt(cls->GetSigmaZ2()));
 292           mNPads->Fill(cls->GetNPads());
 293           mCenter->Fill(cls->GetCenter());
 294                 mTime->Fill(cls->GetLocalTimeBin());
 295 #endif
 296
 297                 if (cls->GetNPads() > 2 && cls->GetNPads() < 9) mChrg3->Fill(cls->GetQ());
 298           if (cls->GetNPads() < 9) mChrg2->Fill(cls->GetQ());
 299           mChrg->Fill(cls->GetQ());
 300
 301                 mQvsPads->Fill(cls->GetNPads(),cls->GetQ());
 302                 mQvsDetector->Fill(cls->GetDetector(),cls->GetQ());
 303                 mQvsSModId->Fill((Int_t)(cls->GetDetector() / 30.),cls->GetQ());
 304
 305           Int_t superModule = (Int_t)(cls->GetDetector() / 30.);
 306
 307           // Position in supermodule = cls->GetDetector() - superModule * 30
 308           // There are 6 layers per supermodule:
 309           Int_t layer = (cls->GetDetector() - superModule * 30) % 6;
 310                 switch (superModule)
 311                 {
 312                   case 0:
 313                     mSMod0->Fill(cls->GetQ());
 314                     mTrSMod0->Fill(cls->GetPadCol(), cls->GetLocalTimeBin() + layer * 20);
 315                     break;
 316                   case 8:
 317                     mSMod8->Fill(cls->GetQ());
 318                     mTrSMod8->Fill(cls->GetPadCol(), cls->GetLocalTimeBin() + layer * 20);
 319                     break;
 320                   case 9:
 321                     mSMod9->Fill(cls->GetQ());
 322                     mTrSMod9->Fill(cls->GetPadCol(), cls->GetLocalTimeBin() + layer * 20);
 323                     break;
 324                   case 17:
 325                     mSMod17->Fill(cls->GetQ());
 326                     mTrSMod17->Fill(cls->GetPadCol(), cls->GetLocalTimeBin() + layer * 20);
 327                     break;
 328                   default:
 329                     break;
 330                 }
 331
 332                 // Only consider events with maximal 15 pads here
 333                 if (cls->GetNPads() - 2 < 14 && cls->GetNPads() - 2 >= 0)
 334                 {
 335                   mPads[cls->GetNPads() - 2]->Fill(cls->GetQ());
 336                 }
 337
 338                 mChrgEBE->Fill(cls->GetQ());
 339               }
 340       }
 341
 342
 343 #ifdef EVENT_BY_EVENT
 344       // Event by event display for the first 100 events
 345       if (numPassed < 100)
 346       {
 347         sprintf(title,"Charge distribution for event %d", ev);
 348         eventc->SetTitle(title);
 349         mChrgEBE->Draw();
 350         eventc->Update();
 351         gSystem->Sleep(400);
 352         mChrgEBE->Reset();
 353       }
 354 #endif
 355
 356 #ifdef EVENT_BY_EVENT2
 357       // Event by event display
 358       sprintf(title2,"Supermodule tracks for event %d", ev);
 359       sTrc->SetTitle(title2);
 360       sTrc->cd(1);
 361       mTrSMod0->Draw("colz");
 362       sTrc->cd(2);
 363       mTrSMod8->Draw("colz");
 364       sTrc->cd(3);
 365       mTrSMod9->Draw("colz");
 366       sTrc->cd(4);
 367       mTrSMod17->Draw("colz");
 368       sTrc->Update();
 369       gSystem->Sleep(400);
 370       mTrSMod0->Reset();
 371       mTrSMod8->Reset();
 372       mTrSMod9->Reset();
 373       mTrSMod17->Reset();
 374 #endif
 375
 376
 377       mClusters->Fill(clusters);
 378     }
 379
 380     runLoader->GetNextEvent();
 381   }
 382
 383 #ifdef PRINT_OTHER
 384   new TCanvas();
 385   gPad->SetLogy();
 386   mNCls->Draw();
 387
 388   new TCanvas();
 389   gPad->SetLogy();
 390   mResY->Draw();
 391
 392   new TCanvas();
 393   gPad->SetLogy();
 394   mResZ->Draw();
 395 #endif
 396
 397   new TCanvas("chrgc", "Charge(no pad-filter)", 700, 500);
 398   mChrg->Draw();
 399   // Try landau-fit
 400   TF1* lFit = new TF1("lFit", "landau", 20, 200);
 401   cout << "\n\nFit (no pad-filter):\n";
 402   mChrg->Fit(lFit, "", "", 20, 200);
 403
 404   new TCanvas("chrg2c", "Charge(pads<9)", 700, 500);
 405   mChrg2->Draw();
 406   // Try landau-fit
 407   TF1* lFit2 = new TF1("lFit2", "landau", 20, 200);
 408   cout << "\n\nFit (pads<9):\n";
 409   mChrg2->Fit(lFit2, "", "", 20, 200);
 410
 411   new TCanvas("chrg3c", "Charge(2<pads<9)", 700, 500);
 412   mChrg3->Draw();
 413   // Try landau-fit
 414   TF1* lFit3 = new TF1("lFit3", "landau", 20, 200);
 415   cout << "\n\nFit (2<pads<9):\n";
 416   mChrg3->Fit(lFit3, "", "", 20, 200);
 417
 418 #ifdef PRINT_OTHER
 419   new TCanvas();
 420   mNPads->Draw();
 421
 422   new TCanvas();
 423   mTime->Draw();
 424
 425   new TCanvas();
 426   mX->Draw();
 427
 428   new TCanvas();
 429   mY->Draw();
 430
 431   new TCanvas();
 432   mZ->Draw();
 433
 434   new TCanvas();
 435   mCenter->Draw();
 436
 437   new TCanvas();
 438   mYX->Draw("colz");
 439
 440   new TCanvas();
 441   mPlane->Draw();
 442
 443   new TCanvas();
 444   mPlaneY->Draw("colz");
 445 #endif
 446
 447   new TCanvas();
 448   gPad->SetLogz();
 449   mQvsPads->Draw("colz");
 450
 451   new TCanvas();
 452   gPad->SetLogz();
 453   mQvsDetector->Draw("colz");
 454   TLine* tl = 0;
 455   for (Int_t num = 30; num < 540; num += 30)
 456   {
 457     // Draw vertical lines every 30 bins to sketch the layer modules
 458     tl = new TLine(num - 0.5, 0, num - 0.5, 1000);
 459     tl->Draw();
 460   }
 461
 462   new TCanvas();
 463   gPad->SetLogz();
 464   mQvsSModId->Draw("colz");
 465
 466   TCanvas* sc = new TCanvas("sc","Supermodules", 800, 800);
 467   sc->SetLogz();
 468   sc->Divide(2,2);
 469   sc->cd(1);
 470   mSMod0->Draw();
 471   sc->cd(2);
 472   mSMod8->Draw();
 473   sc->cd(3);
 474   mSMod9->Draw();
 475   sc->cd(4);
 476   mSMod17->Draw();
 477
 478
 479   TCanvas* padc = new TCanvas("padc","Pads", 1000, 1000);
 480   padc->Divide(4,4);
 481   for (Int_t ind = 0; ind < 14; ind++)
 482   {
 483     padc->cd(ind + 1);
 484     mPads[ind]->Draw();
 485   }
 486
 487   new TCanvas();
 488   mClusters->Draw();
 489
 490   TCanvas* clusChamc = new TCanvas("clusChamc", "Clusters against detector with threshold", 800, 800);
 491   clusChamc->Divide(1,2);
 492   clusChamc->cd(1);
 493   mClusCham->Draw();
 494   clusChamc->cd(2);
 495   mClusChamEx->Draw();
 496
 497   TCanvas* noisyc = new TCanvas("noisyc", "Col/Row of noisy detectors");
 498   noisyc->Divide(1,2);
 499   noisyc->cd(1);
 500   gPad->SetLogz();
 501   mColRow->Draw("colz");
 502   noisyc->cd(2);
 503   gPad->SetLogz();
 504   mColRowFiltered->Draw("colz");
 505
 506   TCanvas* ColTimec = new TCanvas("ColTimec", "Col/Time");
 507   ColTimec->Divide(1,2);
 508   ColTimec->cd(1);
 509   gPad->SetLogz();
 510   mColTime->Draw("colz");
 511   ColTimec->cd(2);
 512   gPad->SetLogz();
 513   mColTimeFiltered->Draw("colz");
 514
 515   cout << numPassed << " out of " << nEvents << " events passed the criterion (number of 'good' clusters >= ";
 516   cout << kThresholdLower << " && clusters <= " << kThresholdUpper << ").\n";
 517 }