Change to TRDdigitizer
[u/mrichter/AliRoot.git] / TRD / AliTRDanalyzeDigits.C
index d40abcbc7963024c7e0012f36b07c8ddff40bec7..596778c3cb4aa5b388387ac00c7e5282d72bc3ab 100644 (file)
@@ -6,95 +6,185 @@ Int_t AliTRDanalyzeDigits()
 
   Int_t rc = 0;
 
 
   Int_t rc = 0;
 
+  const Int_t kNpla = AliTRDgeometry::Nplan();
+
   if (!gAlice) {
     cout << "<AliTRDanalyzeDigits> No AliRun object found" << endl;
     rc = 1;
     return rc;
   }
   if (!gAlice) {
     cout << "<AliTRDanalyzeDigits> No AliRun object found" << endl;
     rc = 1;
     return rc;
   }
-  gAlice->GetEvent(0);
+  Int_t nPart = gAlice->GetEvent(0);
 
   // Get the pointer to the TRD detector 
 
   // Get the pointer to the TRD detector 
-  AliTRD *TRD = (AliTRD *) gAlice->GetDetector("TRD");
-  if (!TRD) {
+  AliTRD *trd = (AliTRD *) gAlice->GetDetector("TRD");
+  if (!trd) {
     cout << "<AliTRDanalyzeDigits> No TRD detector found" << endl;
     rc = 2;
     return rc;
   }
 
     cout << "<AliTRDanalyzeDigits> No TRD detector found" << endl;
     rc = 2;
     return rc;
   }
 
+  // Get the digitizer object
+  TFile *file = (TFile *) gROOT->GetListOfFiles()->FindObject("TRD_test.root");
+  AliTRDdigitizer *digitizer = (AliTRDdigitizer *) file->Get("TRDdigitizer");
+  if (!digitizer) {
+    cout << "<AliTRDanalyzeDigits> No digitizer object found" << endl;
+    rc = 3;
+    return rc;
+  }
+
   // Define the histograms
   // Define the histograms
-  TH1F *hAmp = new TH1F("hAmp","Amplitude of the digits",256,-0.5,255.5);
+  Int_t adcRange = ((Int_t) digitizer->GetADCoutRange());
+  TH1F *hAmpAll   = new TH1F("hAmpAll"  ,"Amplitude of the digits (all)"
+                            ,adcRange+1,-0.5,((Float_t) adcRange)+0.5);
+  TH1F *hAmpEl    = new TH1F("hAmpEl"   ,"Amplitude of the digits (electrons)"
+                            ,adcRange+1,-0.5,((Float_t) adcRange)+0.5);
+  TH1F *hAmpPi    = new TH1F("hAmpPi"   ,"Amplitude of the digits (pions)"
+                            ,adcRange+1,-0.5,((Float_t) adcRange)+0.5);
+  TH1F *hAmpNoise = new TH1F("hAmpNoise","Amplitude of the digits (noise)"
+                            ,5,-0.5,4.5);
 
   // Get the pointer to the geometry object
 
   // Get the pointer to the geometry object
-  AliTRDgeometry *TRDgeometry;
-  if (TRD) {
-    TRDgeometry = TRD->GetGeometry();
+  AliTRDgeometry *geo;
+  if (trd) {
+    geo = trd->GetGeometry();
   }
   else {
     cout << "<AliTRDanalyzeDigits> No TRD geometry found" << endl;
   }
   else {
     cout << "<AliTRDanalyzeDigits> No TRD geometry found" << endl;
-    rc = 3;
+    rc = 4;
     return rc;
   }
 
   // Create the digits manager
     return rc;
   }
 
   // Create the digits manager
-  AliTRDdigitsManager *DigitsManager = new AliTRDdigitsManager();
+  AliTRDdigitsManager *digitsManager = new AliTRDdigitsManager();
+
+  digitsManager->Open("TRD_test.root");
 
   // Read the digits from the file
 
   // Read the digits from the file
-  if (!(DigitsManager->ReadDigits())) {
+  if (!(digitsManager->ReadDigits())) {
     cout << "<AliTRDanalyzeDigits> Cannot read the digits" << endl;
     cout << "<AliTRDanalyzeDigits> Cannot read the digits" << endl;
-    rc = 4;
+    rc = 5;
     return rc;
   }
     return rc;
   }
-  
-  // Define the detector matrix for one chamber
-  Int_t iSec = TRD->GetSensSector();
-  Int_t iCha = TRD->GetSensChamber();
-  Int_t iPla = 0;
-  Int_t  rowMax = TRDgeometry->GetRowMax(iPla,iCha,iSec);
-  Int_t  colMax = TRDgeometry->GetColMax(iPla);
-  Int_t timeMax = TRDgeometry->GetTimeMax();
-  cout << "<AliTRDanalyzeDigits> Geometry: rowMax = "  <<  rowMax
-                                      << " colMax = "  <<  colMax
-                                      << " timeMax = " << timeMax << endl;
-  AliTRDmatrix *TRDmatrix = new AliTRDmatrix(rowMax,colMax,timeMax,iSec,iCha,iPla);
-  // Get the detector number
-  Int_t iDet = TRDgeometry->GetDetector(iPla,iCha,iSec);
-  cout << "<AliTRDanalyzeDigits> iSec = " << iSec
-                            << " iCha = " << iCha 
-                            << " iPla = " << iPla 
-                            << " iDet = " << iDet << endl;
-
-  // Loop through the detector pixel
+
+  AliTRDmatrix *matrix;
+
   Int_t countDigits = 0;
   Int_t countDigits = 0;
-  for (Int_t time = 0; time < timeMax; time++) {
-    for (Int_t  col = 0;  col <  colMax;  col++) {
-      for (Int_t  row = 0;  row <  rowMax;  row++) {
+  Int_t iSec        = trd->GetSensSector();
+  Int_t iCha        = trd->GetSensChamber();
+  Int_t timeMax     = geo->GetTimeTotal();
 
 
-        AliTRDdigit *Digit = DigitsManager->GetDigit(row,col,time,iDet);
-        Int_t amp = Digit->GetAmp();
+  TProfile *hAmpTimeEl = new TProfile("hAmpTimeEl","Amplitude of the digits (electrons)"
+                                     ,timeMax,-0.5,((Double_t) timeMax)-0.5);
+  TProfile *hAmpTimePi = new TProfile("hAmpTimePi","Amplitude of the digits (pions)"
+                                     ,timeMax,-0.5,((Double_t) timeMax)-0.5);
 
 
-        if (amp > 0) {
-          countDigits++;
-          hAmp->Fill(amp);
-          TRDmatrix->SetSignal(row,col,time,amp);
-       }
+  // Loop over all planes
+  for (Int_t iPla = 0; iPla < kNpla; iPla++) {
 
 
-        delete Digit;
+    Int_t iDet   = geo->GetDetector(iPla,iCha,iSec);
+    Int_t rowMax = geo->GetRowMax(iPla,iCha,iSec);
+    Int_t colMax = geo->GetColMax(iPla);
+  
+    if (iPla == 0) {
+      matrix = new AliTRDmatrix(rowMax,colMax,timeMax,iSec,iCha,iPla);
+    }
 
 
+    // Loop through the detector pixel
+    for (Int_t time = 0; time < timeMax; time++) {
+      for (Int_t  col = 0;  col <  colMax;  col++) {
+        for (Int_t  row = 0;  row <  rowMax;  row++) {
+
+          AliTRDdigit *digit    = digitsManager->GetDigit(row,col,time,iDet);
+          Int_t        amp      = digit->GetAmp();
+          Int_t        track0   = digitsManager->GetTrack(0,row,col,time,iDet);
+          Int_t        track1   = digitsManager->GetTrack(1,row,col,time,iDet);
+          TParticle   *particle = 0;
+          if (track0 > -1) {
+            particle = gAlice->Particle(track0);
+         }
+
+          if (amp > 0) {
+            countDigits++;
+            if (iPla == 0) {
+              matrix->SetSignal(row,col,time,amp);
+           }
+         }
+
+         // Total spectrum
+          hAmpAll->Fill(amp);
+
+         // Noise spectrum
+          if (track0 < 0) {
+            hAmpNoise->Fill(amp);
+         }          
+
+         // Electron digit
+          if ((particle) && (particle->GetPdgCode() ==   11) && (track1 < 0)) {
+            hAmpEl->Fill(amp);
+            hAmpTimeEl->Fill(time,amp);
+         }
+
+          // Pion digit
+          if ((particle) && (particle->GetPdgCode() == -211) && (track1 < 0)) {
+            hAmpPi->Fill(amp);
+            hAmpTimePi->Fill(time,amp);
+         }
+
+          delete digit;
+
+        }
       }
     }
       }
     }
+
   }
 
   cout << "<AliTRDanalyzeDigits> Found " << countDigits << " digits in total" << endl;
 
   // Display the detector matrix
   }
 
   cout << "<AliTRDanalyzeDigits> Found " << countDigits << " digits in total" << endl;
 
   // Display the detector matrix
-  TRDmatrix->Draw();
-  TRDmatrix->ProjRow();
-  TRDmatrix->ProjCol();
-  TRDmatrix->ProjTime();
-
-  TCanvas *cDigits = new TCanvas("cDigits","AliTRDanalyzeDigits",50,50,600,600);
+  matrix->Draw();
+  matrix->ProjRow();
+  matrix->ProjCol();
+  matrix->ProjTime();
+
+  TCanvas *cDigits = new TCanvas("cDigits","AliTRDanalyzeDigits",50,50,600,800);
+  cDigits->Divide(2,3);
+  cDigits->cd(1);
+  gPad->SetLogy();
+  hAmpAll->SetXTitle("Amplitude (ADC-channels)");
+  hAmpAll->SetYTitle("Entries");
+  hAmpAll->Draw();
+  cDigits->cd(2);
+  gPad->SetLogy();
+  hAmpNoise->SetXTitle("Amplitude (ADC-channels)");
+  hAmpNoise->SetYTitle("Entries");
+  hAmpNoise->Draw();
+  cDigits->cd(3);
+  gPad->SetLogy();
+  hAmpEl->SetXTitle("Amplitude (ADC-channels)");
+  hAmpEl->SetYTitle("Entries");
+  hAmpEl->Draw();
+  cDigits->cd(4);
   gPad->SetLogy();
   gPad->SetLogy();
-  hAmp->Draw();
+  hAmpPi->SetXTitle("Amplitude (ADC-channels)");
+  hAmpPi->SetYTitle("Entries");
+  hAmpPi->Draw();
+  cDigits->cd(5);
+  hAmpTimeEl->SetXTitle("Timebin number");
+  hAmpTimeEl->SetYTitle("Mean amplitude");
+  hAmpTimeEl->Draw("HIST");
+  cDigits->cd(6);
+  hAmpTimePi->SetXTitle("Timebin number");
+  hAmpTimePi->SetYTitle("Mean amplitude");
+  hAmpTimePi->Draw("HIST");
+
+  TFile *fileOut = new TFile("digits_test.root","RECREATE");
+  hAmpAll->Write();
+  hAmpNoise->Write();
+  hAmpEl->Write();
+  hAmpPi->Write();
+  hAmpTimeEl->Write();
+  hAmpTimePi->Write();
+  fileOut->Close();
 
   return rc;
 
 
   return rc;