TRD/AliTRDanalyzeDigits.C

   1 Int_t AliTRDanalyzeDigits()
   2 {
   3   //
   4   // Analyzes the digits
   5   //
   6
   7   Int_t rc = 0;
   8
   9   const Int_t kNpla = AliTRDgeometry::Nplan();
  10
  11   if (!gAlice) {
  12     cout << "<AliTRDanalyzeDigits> No AliRun object found" << endl;
  13     rc = 1;
  14     return rc;
  15   }
  16   Int_t nPart = gAlice->GetEvent(0);
  17
  18   // Get the pointer to the TRD detector
  19   AliTRD *TRD = (AliTRD *) gAlice->GetDetector("TRD");
  20   if (!TRD) {
  21     cout << "<AliTRDanalyzeDigits> No TRD detector found" << endl;
  22     rc = 2;
  23     return rc;
  24   }
  25
  26   // Get the digitizer object
  27   TFile *file = (TFile *) gROOT->GetListOfFiles()->FindObject("TRD_test.root");
  28   AliTRDdigitizer *Digitizer = (AliTRDdigitizer *) file->Get("digitizer");
  29   if (!Digitizer) {
  30     cout << "<AliTRDanalyzeDigits> No digitizer object found" << endl;
  31     rc = 3;
  32     return rc;
  33   }
  34
  35   // Define the histograms
  36   Int_t adcRange = ((Int_t) Digitizer->GetADCoutRange());
  37   TH1F *hAmpAll   = new TH1F("hAmpAll"  ,"Amplitude of the digits (all)"
  38                             ,adcRange+1,-0.5,((Float_t) adcRange)+0.5);
  39   TH1F *hAmpEl    = new TH1F("hAmpEl"   ,"Amplitude of the digits (electrons)"
  40                             ,adcRange+1,-0.5,((Float_t) adcRange)+0.5);
  41   TH1F *hAmpPi    = new TH1F("hAmpPi"   ,"Amplitude of the digits (pions)"
  42                             ,adcRange+1,-0.5,((Float_t) adcRange)+0.5);
  43   TH1F *hAmpNoise = new TH1F("hAmpNoise","Amplitude of the digits (noise)"
  44                             ,5,-0.5,4.5);
  45
  46   // Get the pointer to the geometry object
  47   AliTRDgeometry *TRDgeometry;
  48   if (TRD) {
  49     TRDgeometry = TRD->GetGeometry();
  50   }
  51   else {
  52     cout << "<AliTRDanalyzeDigits> No TRD geometry found" << endl;
  53     rc = 4;
  54     return rc;
  55   }
  56
  57   // Create the digits manager
  58   AliTRDdigitsManager *DigitsManager = new AliTRDdigitsManager();
  59
  60   // Read the digits from the file
  61   if (!(DigitsManager->ReadDigits())) {
  62     cout << "<AliTRDanalyzeDigits> Cannot read the digits" << endl;
  63     rc = 5;
  64     return rc;
  65   }
  66
  67   AliTRDmatrix *TRDmatrix;
  68
  69   Int_t countDigits = 0;
  70   Int_t iSec        = TRD->GetSensSector();
  71   Int_t iCha        = TRD->GetSensChamber();
  72   Int_t timeMax     = TRDgeometry->GetTimeTotal();
  73
  74   TProfile *hAmpTimeEl = new TProfile("hAmpTimeEl","Amplitude of the digits (electrons)"
  75                                       ,timeMax,-0.5,((Double_t) timeMax)-0.5);
  76   TProfile *hAmpTimePi = new TProfile("hAmpTimePi","Amplitude of the digits (pions)"
  77                                       ,timeMax,-0.5,((Double_t) timeMax)-0.5);
  78
  79   // Loop over all planes
  80   for (Int_t iPla = 0; iPla < kNpla; iPla++) {
  81
  82     Int_t iDet   = TRDgeometry->GetDetector(iPla,iCha,iSec);
  83     Int_t rowMax = TRDgeometry->GetRowMax(iPla,iCha,iSec);
  84     Int_t colMax = TRDgeometry->GetColMax(iPla);
  85
  86     if (iPla == 0) {
  87       TRDmatrix = new AliTRDmatrix(rowMax,colMax,timeMax,iSec,iCha,iPla);
  88     }
  89
  90     // Loop through the detector pixel
  91     for (Int_t time = 0; time < timeMax; time++) {
  92       for (Int_t  col = 0;  col <  colMax;  col++) {
  93         for (Int_t  row = 0;  row <  rowMax;  row++) {
  94
  95           AliTRDdigit *Digit = DigitsManager->GetDigit(row,col,time,iDet);
  96           Int_t amp = Digit->GetAmp();
  97
  98           Int_t track0 = DigitsManager->GetTrack(0,row,col,time,iDet);
  99           Int_t track1 = DigitsManager->GetTrack(1,row,col,time,iDet);
 100           TParticle *Part = 0;
 101           if (track0 > -1) {
 102             Part = gAlice->Particle(track0);
 103           }
 104
 105           if (amp > 0) {
 106             countDigits++;
 107             if (iPla == 0) {
 108               TRDmatrix->SetSignal(row,col,time,amp);
 109             }
 110           }
 111
 112           // Total spectrum
 113           hAmpAll->Fill(amp);
 114
 115           // Noise spectrum
 116           if (track0 < 0) {
 117             hAmpNoise->Fill(amp);
 118           }
 119
 120           // Electron digit
 121           if ((Part) && (Part->GetPdgCode() ==   11) && (track1 < 0)) {
 122             hAmpEl->Fill(amp);
 123             hAmpTimeEl->Fill(time,amp);
 124           }
 125
 126           // Pion digit
 127           if ((Part) && (Part->GetPdgCode() == -211) && (track1 < 0)) {
 128             hAmpPi->Fill(amp);
 129             hAmpTimePi->Fill(time,amp);
 130           }
 131
 132           delete Digit;
 133
 134         }
 135       }
 136     }
 137
 138   }
 139
 140   cout << "<AliTRDanalyzeDigits> Found " << countDigits << " digits in total" << endl;
 141
 142   // Display the detector matrix
 143   TRDmatrix->Draw();
 144   TRDmatrix->ProjRow();
 145   TRDmatrix->ProjCol();
 146   TRDmatrix->ProjTime();
 147
 148   TCanvas *cDigits = new TCanvas("cDigits","AliTRDanalyzeDigits",50,50,600,800);
 149   cDigits->Divide(2,3);
 150   cDigits->cd(1);
 151   gPad->SetLogy();
 152   hAmpAll->SetXTitle("Amplitude (ADC-channels)");
 153   hAmpAll->SetYTitle("Entries");
 154   hAmpAll->Draw();
 155   cDigits->cd(2);
 156   gPad->SetLogy();
 157   hAmpNoise->SetXTitle("Amplitude (ADC-channels)");
 158   hAmpNoise->SetYTitle("Entries");
 159   hAmpNoise->Draw();
 160   cDigits->cd(3);
 161   gPad->SetLogy();
 162   hAmpEl->SetXTitle("Amplitude (ADC-channels)");
 163   hAmpEl->SetYTitle("Entries");
 164   hAmpEl->Draw();
 165   cDigits->cd(4);
 166   gPad->SetLogy();
 167   hAmpPi->SetXTitle("Amplitude (ADC-channels)");
 168   hAmpPi->SetYTitle("Entries");
 169   hAmpPi->Draw();
 170   cDigits->cd(5);
 171   hAmpTimeEl->SetXTitle("Timebin number");
 172   hAmpTimeEl->SetYTitle("Mean amplitude");
 173   hAmpTimeEl->Draw("HIST");
 174   cDigits->cd(6);
 175   hAmpTimePi->SetXTitle("Timebin number");
 176   hAmpTimePi->SetYTitle("Mean amplitude");
 177   hAmpTimePi->Draw("HIST");
 178
 179   TFile *fileOut = new TFile("digits_test.root","RECREATE");
 180   hAmpAll->Write();
 181   hAmpNoise->Write();
 182   hAmpEl->Write();
 183   hAmpPi->Write();
 184   hAmpTimeEl->Write();
 185   hAmpTimePi->Write();
 186   fileOut->Close();
 187
 188   return rc;
 189
 190 }