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("TRDdigitizer");
  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 *geo;
  48   if (trd) {
  49     geo = 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   digitsManager->Open("TRD_test.root");
  61
  62   // Read the digits from the file
  63   if (!(digitsManager->ReadDigits())) {
  64     cout << "<AliTRDanalyzeDigits> Cannot read the digits" << endl;
  65     rc = 5;
  66     return rc;
  67   }
  68
  69   AliTRDmatrix *matrix;
  70
  71   Int_t countDigits = 0;
  72   Int_t iSec        = trd->GetSensSector();
  73   Int_t iCha        = trd->GetSensChamber();
  74   Int_t timeMax     = geo->GetTimeTotal();
  75
  76   TProfile *hAmpTimeEl = new TProfile("hAmpTimeEl","Amplitude of the digits (electrons)"
  77                                       ,timeMax,-0.5,((Double_t) timeMax)-0.5);
  78   TProfile *hAmpTimePi = new TProfile("hAmpTimePi","Amplitude of the digits (pions)"
  79                                       ,timeMax,-0.5,((Double_t) timeMax)-0.5);
  80
  81   // Loop over all planes
  82   for (Int_t iPla = 0; iPla < kNpla; iPla++) {
  83
  84     Int_t iDet   = geo->GetDetector(iPla,iCha,iSec);
  85     Int_t rowMax = geo->GetRowMax(iPla,iCha,iSec);
  86     Int_t colMax = geo->GetColMax(iPla);
  87
  88     if (iPla == 0) {
  89       matrix = new AliTRDmatrix(rowMax,colMax,timeMax,iSec,iCha,iPla);
  90     }
  91
  92     // Loop through the detector pixel
  93     for (Int_t time = 0; time < timeMax; time++) {
  94       for (Int_t  col = 0;  col <  colMax;  col++) {
  95         for (Int_t  row = 0;  row <  rowMax;  row++) {
  96
  97           AliTRDdigit *digit    = digitsManager->GetDigit(row,col,time,iDet);
  98           Int_t        amp      = digit->GetAmp();
  99           Int_t        track0   = digitsManager->GetTrack(0,row,col,time,iDet);
 100           Int_t        track1   = digitsManager->GetTrack(1,row,col,time,iDet);
 101           TParticle   *particle = 0;
 102           if (track0 > -1) {
 103             particle = gAlice->Particle(track0);
 104           }
 105
 106           if (amp > 0) {
 107             countDigits++;
 108             if (iPla == 0) {
 109               matrix->SetSignal(row,col,time,amp);
 110             }
 111           }
 112
 113           // Total spectrum
 114           hAmpAll->Fill(amp);
 115
 116           // Noise spectrum
 117           if (track0 < 0) {
 118             hAmpNoise->Fill(amp);
 119           }
 120
 121           // Electron digit
 122           if ((particle) && (particle->GetPdgCode() ==   11) && (track1 < 0)) {
 123             hAmpEl->Fill(amp);
 124             hAmpTimeEl->Fill(time,amp);
 125           }
 126
 127           // Pion digit
 128           if ((particle) && (particle->GetPdgCode() == -211) && (track1 < 0)) {
 129             hAmpPi->Fill(amp);
 130             hAmpTimePi->Fill(time,amp);
 131           }
 132
 133           delete digit;
 134
 135         }
 136       }
 137     }
 138
 139   }
 140
 141   cout << "<AliTRDanalyzeDigits> Found " << countDigits << " digits in total" << endl;
 142
 143   // Display the detector matrix
 144   matrix->Draw();
 145   matrix->ProjRow();
 146   matrix->ProjCol();
 147   matrix->ProjTime();
 148
 149   TCanvas *cDigits = new TCanvas("cDigits","AliTRDanalyzeDigits",50,50,600,800);
 150   cDigits->Divide(2,3);
 151   cDigits->cd(1);
 152   gPad->SetLogy();
 153   hAmpAll->SetXTitle("Amplitude (ADC-channels)");
 154   hAmpAll->SetYTitle("Entries");
 155   hAmpAll->Draw();
 156   cDigits->cd(2);
 157   gPad->SetLogy();
 158   hAmpNoise->SetXTitle("Amplitude (ADC-channels)");
 159   hAmpNoise->SetYTitle("Entries");
 160   hAmpNoise->Draw();
 161   cDigits->cd(3);
 162   gPad->SetLogy();
 163   hAmpEl->SetXTitle("Amplitude (ADC-channels)");
 164   hAmpEl->SetYTitle("Entries");
 165   hAmpEl->Draw();
 166   cDigits->cd(4);
 167   gPad->SetLogy();
 168   hAmpPi->SetXTitle("Amplitude (ADC-channels)");
 169   hAmpPi->SetYTitle("Entries");
 170   hAmpPi->Draw();
 171   cDigits->cd(5);
 172   hAmpTimeEl->SetXTitle("Timebin number");
 173   hAmpTimeEl->SetYTitle("Mean amplitude");
 174   hAmpTimeEl->Draw("HIST");
 175   cDigits->cd(6);
 176   hAmpTimePi->SetXTitle("Timebin number");
 177   hAmpTimePi->SetYTitle("Mean amplitude");
 178   hAmpTimePi->Draw("HIST");
 179
 180   TFile *fileOut = new TFile("digits_test.root","RECREATE");
 181   hAmpAll->Write();
 182   hAmpNoise->Write();
 183   hAmpEl->Write();
 184   hAmpPi->Write();
 185   hAmpTimeEl->Write();
 186   hAmpTimePi->Write();
 187   fileOut->Close();
 188
 189   return rc;
 190
 191 }