declaration for StepManager()
[u/mrichter/AliRoot.git] / TRD / AliTRDanalyzeDigits.C
CommitLineData
793ff80c 1Int_t AliTRDanalyzeDigits()
2{
3 //
4 // Analyzes the digits
5 //
6
7 Int_t rc = 0;
8
f9428ca8 9 const Int_t kNpla = AliTRDgeometry::Nplan();
10
e2c86a4a 11 if (!gAlice) {
12 cout << "<AliTRDanalyzeDigits> No AliRun object found" << endl;
13 rc = 1;
14 return rc;
15 }
4a0fe73c 16 Int_t nPart = gAlice->GetEvent(0);
e2c86a4a 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
4a0fe73c 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
e2c86a4a 35 // Define the histograms
4a0fe73c 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)"
f9428ca8 44 ,5,-0.5,4.5);
e2c86a4a 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;
4a0fe73c 53 rc = 4;
e2c86a4a 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;
4a0fe73c 63 rc = 5;
e2c86a4a 64 return rc;
65 }
4a0fe73c 66
f9428ca8 67 AliTRDmatrix *TRDmatrix;
4a0fe73c 68
f9428ca8 69 Int_t countDigits = 0;
70 Int_t iSec = TRD->GetSensSector();
71 Int_t iCha = TRD->GetSensChamber();
72 Int_t timeMax = TRDgeometry->GetTimeTotal();
4a0fe73c 73
f9428ca8 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);
4a0fe73c 78
f9428ca8 79 // Loop over all planes
80 for (Int_t iPla = 0; iPla < kNpla; iPla++) {
e2c86a4a 81
f9428ca8 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 }
e2c86a4a 89
f9428ca8 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 }
e2c86a4a 135 }
136 }
f9428ca8 137
e2c86a4a 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
f9428ca8 148 TCanvas *cDigits = new TCanvas("cDigits","AliTRDanalyzeDigits",50,50,600,800);
149 cDigits->Divide(2,3);
4a0fe73c 150 cDigits->cd(1);
151 gPad->SetLogy();
f9428ca8 152 hAmpAll->SetXTitle("Amplitude (ADC-channels)");
153 hAmpAll->SetYTitle("Entries");
4a0fe73c 154 hAmpAll->Draw();
155 cDigits->cd(2);
156 gPad->SetLogy();
f9428ca8 157 hAmpNoise->SetXTitle("Amplitude (ADC-channels)");
158 hAmpNoise->SetYTitle("Entries");
4a0fe73c 159 hAmpNoise->Draw();
160 cDigits->cd(3);
161 gPad->SetLogy();
f9428ca8 162 hAmpEl->SetXTitle("Amplitude (ADC-channels)");
163 hAmpEl->SetYTitle("Entries");
4a0fe73c 164 hAmpEl->Draw();
165 cDigits->cd(4);
e2c86a4a 166 gPad->SetLogy();
f9428ca8 167 hAmpPi->SetXTitle("Amplitude (ADC-channels)");
168 hAmpPi->SetYTitle("Entries");
4a0fe73c 169 hAmpPi->Draw();
f9428ca8 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();
e2c86a4a 187
793ff80c 188 return rc;
189
190}