Update of the macros
[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("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 }