Changes in digits IO. Add merging of summable digits
[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
abaf1f1d 19 AliTRD *trd = (AliTRD *) gAlice->GetDetector("TRD");
20 if (!trd) {
e2c86a4a 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");
abaf1f1d 28 AliTRDdigitizer *digitizer = (AliTRDdigitizer *) file->Get("digitizer");
29 if (!digitizer) {
4a0fe73c 30 cout << "<AliTRDanalyzeDigits> No digitizer object found" << endl;
31 rc = 3;
32 return rc;
33 }
34
e2c86a4a 35 // Define the histograms
abaf1f1d 36 Int_t adcRange = ((Int_t) digitizer->GetADCoutRange());
4a0fe73c 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
abaf1f1d 47 AliTRDgeometry *geo;
48 if (trd) {
49 geo = trd->GetGeometry();
e2c86a4a 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
abaf1f1d 58 AliTRDdigitsManager *digitsManager = new AliTRDdigitsManager();
59
60 digitsManager->Open("TRD_test.root");
e2c86a4a 61
62 // Read the digits from the file
abaf1f1d 63 if (!(digitsManager->ReadDigits())) {
e2c86a4a 64 cout << "<AliTRDanalyzeDigits> Cannot read the digits" << endl;
4a0fe73c 65 rc = 5;
e2c86a4a 66 return rc;
67 }
4a0fe73c 68
abaf1f1d 69 AliTRDmatrix *matrix;
4a0fe73c 70
f9428ca8 71 Int_t countDigits = 0;
abaf1f1d 72 Int_t iSec = trd->GetSensSector();
73 Int_t iCha = trd->GetSensChamber();
74 Int_t timeMax = geo->GetTimeTotal();
4a0fe73c 75
f9428ca8 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);
4a0fe73c 80
f9428ca8 81 // Loop over all planes
82 for (Int_t iPla = 0; iPla < kNpla; iPla++) {
e2c86a4a 83
abaf1f1d 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);
f9428ca8 87
88 if (iPla == 0) {
abaf1f1d 89 matrix = new AliTRDmatrix(rowMax,colMax,timeMax,iSec,iCha,iPla);
f9428ca8 90 }
e2c86a4a 91
f9428ca8 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
abaf1f1d 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;
f9428ca8 102 if (track0 > -1) {
abaf1f1d 103 particle = gAlice->Particle(track0);
f9428ca8 104 }
105
106 if (amp > 0) {
107 countDigits++;
108 if (iPla == 0) {
abaf1f1d 109 matrix->SetSignal(row,col,time,amp);
f9428ca8 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
abaf1f1d 122 if ((particle) && (particle->GetPdgCode() == 11) && (track1 < 0)) {
f9428ca8 123 hAmpEl->Fill(amp);
124 hAmpTimeEl->Fill(time,amp);
125 }
126
127 // Pion digit
abaf1f1d 128 if ((particle) && (particle->GetPdgCode() == -211) && (track1 < 0)) {
f9428ca8 129 hAmpPi->Fill(amp);
130 hAmpTimePi->Fill(time,amp);
131 }
132
abaf1f1d 133 delete digit;
f9428ca8 134
135 }
e2c86a4a 136 }
137 }
f9428ca8 138
e2c86a4a 139 }
140
141 cout << "<AliTRDanalyzeDigits> Found " << countDigits << " digits in total" << endl;
142
143 // Display the detector matrix
abaf1f1d 144 matrix->Draw();
145 matrix->ProjRow();
146 matrix->ProjCol();
147 matrix->ProjTime();
e2c86a4a 148
f9428ca8 149 TCanvas *cDigits = new TCanvas("cDigits","AliTRDanalyzeDigits",50,50,600,800);
150 cDigits->Divide(2,3);
4a0fe73c 151 cDigits->cd(1);
152 gPad->SetLogy();
f9428ca8 153 hAmpAll->SetXTitle("Amplitude (ADC-channels)");
154 hAmpAll->SetYTitle("Entries");
4a0fe73c 155 hAmpAll->Draw();
156 cDigits->cd(2);
157 gPad->SetLogy();
f9428ca8 158 hAmpNoise->SetXTitle("Amplitude (ADC-channels)");
159 hAmpNoise->SetYTitle("Entries");
4a0fe73c 160 hAmpNoise->Draw();
161 cDigits->cd(3);
162 gPad->SetLogy();
f9428ca8 163 hAmpEl->SetXTitle("Amplitude (ADC-channels)");
164 hAmpEl->SetYTitle("Entries");
4a0fe73c 165 hAmpEl->Draw();
166 cDigits->cd(4);
e2c86a4a 167 gPad->SetLogy();
f9428ca8 168 hAmpPi->SetXTitle("Amplitude (ADC-channels)");
169 hAmpPi->SetYTitle("Entries");
4a0fe73c 170 hAmpPi->Draw();
f9428ca8 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();
e2c86a4a 188
793ff80c 189 return rc;
190
191}