Change to TRDdigitizer
[u/mrichter/AliRoot.git] / 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 }