Update of calibration task and replace AliTRDrunCalib.C by AddTaskTRDCalib.C (Raphaelle)
[u/mrichter/AliRoot.git] / TRD / Macros / AliTRDanalyzeTestBeam.C
CommitLineData
1a6c0f30 1
2void AliTRDanalyzeTestBeam(Int_t run, Int_t begin, Int_t end) {
3
4 gROOT->SetStyle("Plain");
5 //gStyle->SetPadTopMargin(0.02);
6 //gStyle->SetPadRightMargin(0.02);
7 //gStyle->SetPadLeftMargin(0.1);
8
9 gStyle->SetPalette(1);
10 gStyle->SetOptStat(0);
11 gStyle->SetOptDate();
12
13 TGaxis::SetMaxDigits(3);
14
15 TStopwatch st;
16 st.Start();
17
18 const Int_t N = 640;
19
20 // declare histograms
21 const int nbins = 100;
22 const int nstart = 0;
23
24 TH1D *mSi1L = new TH1D("mSi1L", ";number of Si1 fired pads", 50, nstart-0.5, nstart+nbins-0.5);
25 TH1D *mSi2L = new TH1D("mS21L", ";number of Si2 fired pads", 50, nstart-0.5, nstart+nbins-0.5);
26
27 TProfile *mSi1ChP = new TProfile("mSi1ChP",";Si1 pad number;signal", 1280, -0.5, 1279.5, 0, 200, "s");
28 TProfile *mSi2ChP = new TProfile("mS21ChP",";Si2 pad number;signal", 1280, -0.5, 1279.5, 0, 200, "s");
29
30 TH1D *mSi1N = new TH1D("mSi1N", "Noise Dist Si1;ADC", 100, 0, 50);
31 TH1D *mSi2N = new TH1D("mSi2N", "Noise Dist Si2;ADC", 100, 0, 50);
32
33 TH1D *mSiCh[4];
34 mSiCh[0] = new TH1D("mSi1ChX", ";Si1X pad amplitude (ADC)", 250, -0.5, 249.5);
35 mSiCh[1] = new TH1D("mSi1ChY", ";Si1Y pad amplitude (ADC)", 250, -0.5, 249.5);
36 mSiCh[2] = new TH1D("mSi2ChX", ";Si2X pad amplitude (ADC)", 250, -0.5, 249.5);
37 mSiCh[3] = new TH1D("mSi2ChY", ";Si2Y pad amplitude (ADC)", 250, -0.5, 249.5);
38
39 TH1D *mSiFullCh[4];
40 mSiFullCh[0] = new TH1D("mSi1fChX", "Si1X;max amplitude (ADC)", 300, -0.5, 299.5);
41 mSiFullCh[1] = new TH1D("mSi1fChY", "Si1Y;max amplitude (ADC)", 300, -0.5, 299.5);
42 mSiFullCh[2] = new TH1D("mSi2fChX", "Si2X;max amplitude (ADC)", 300, -0.5, 299.5);
43 mSiFullCh[3] = new TH1D("mSi2fChY", "Si2Y;max amplitude (ADC)", 300, -0.5, 299.5);
44
45 TH2D *mPos[2];
46 mPos[0] = new TH2D("posSi1", ";x Si1 (mm);y Si1 (mm)", 128, 0, 32, 128, 0, 32);
47 mPos[1] = new TH2D("posSi2", ";x Si2 (mm);y Si2 (mm)", 128, 0, 32, 128, 0, 32);
48
49 TH2D *mCor[2];
50 mCor[0] = new TH2D("corX", ";Si1 X (mm);Si2 X (mm)", 128, 0, 32, 128, 0, 32);
51 mCor[1] = new TH2D("corY", ";Si1 Y (mm);Si2 Y (mm)", 128, 0, 32, 128, 0, 32);
52
53 TH2D *mChCor[2];
54 mChCor[0] = new TH2D("ChCorSi1", ";Si1 amp X;Si1 amp Y", 100, 0, 200, 100, 0, 200);
55 mChCor[1] = new TH2D("ChCorSi2", ";Si2 amp X;Si2 amp Y", 100, 0, 200, 100, 0, 200);
56
57 gStyle->SetOptStat(11);
58 TH1D *mPb = new TH1D("mPb", ";Amp Pb (ADC)", 150, -0.5, 4499.5);
59 TH1D *mCher = new TH1D("mCher", ";Amp Cherenkov (ADC)", 150, -0.5, 4499.5);
60 TH2D *mPbCher = new TH2D("mPbCher", ";amp Cherenkov;amp Pb", 150, -0.5, 4499.5, 150, -0.5, 4599.5);
61 // gStyle->SetOptStat(0);
62
63 // needed by the AliTRDRawStreamTB
162637e4 64 AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
1a6c0f30 65 AliTRDcalibDB *calib = AliTRDcalibDB::Instance();
66 calib->SetRun(0);
67
68 // TRD data monitoring
69 TH1D *mDet = new TH1D("mDet", ";chamber", 20, -0.5, 19.5);
70 TH1D *mROB = new TH1D("mROB", ";rob", 20, -0.5, 19.5);
71 TH1D *mTRDsig = new TH1D("mTRDsig", ";trdSignal", 100, -0.5, 299.5);
72
73 //AliLog::SetClassDebugLevel("AliTRDRawStreamTB", 10);
74
75 int counter = 0;
76 // for(Int_t run = 365; run < 386; run++) {
77 //for(Int_t run = 369; run < 382; run++) {
78 //for(int run = 387; run < 389; run++) {
79
80 // if (run == 389) continue;
81 cout << run << endl;
82
83 for(Int_t fn=begin; fn<end; fn++) {
84
85 // connect to data
86 //const char *base="/Users/radomski/data/1GeV/";
87 const char *base = "./";
88 const char *filename = Form("%s/run%d_gdc_daq09.%03d.raw",base,run,fn);
89 if (gSystem->AccessPathName(filename)) continue;
90 cout << filename << endl;
91
92 AliTRDtestBeam *data = new AliTRDtestBeam(filename);
93
94 // process data
95 while (data->NextEvent()) {
96
97 if (!(counter%1000)) cout << "Event = " << counter << endl;
98 counter++;
99
100 /*
101 AliTRDRawStreamTB *tb = data->GetTRDrawStream();
102 while(tb->Next()) {
103 mROB->Fill(tb->GetROB());
104 mDet->Fill(tb->GetDet());
105 int *sig = tb->GetSignals();
106 mTRDsig->Fill(sig[0]);
107 mTRDsig->Fill(sig[1]);
108 mTRDsig->Fill(sig[2]);
109 }
110 delete tb;
111 */
112
113 mCher->Fill(data->GetCher());
114 mPb->Fill(data->GetPb());
115 mPbCher->Fill(data->GetCher(), data->GetPb());
116
117 mSi1L->Fill(data->GetNSi1());
118 mSi2L->Fill(data->GetNSi2());
119
120 for(int i=0; i<data->GetNSi1(); i++) {
121 Int_t q = data->GetSi1Charge(i);
122 Int_t a = data->GetSi1Address(i);
123 if (a == 0) continue; // noisy channels
124 mSi1ChP->Fill(a, q);
125 if (a < N) mSiCh[0]->Fill(q);
126 else mSiCh[1]->Fill(q);
127 }
128
129 for(int i=0; i<data->GetNSi2(); i++) {
130 Int_t q = data->GetSi2Charge(i);
131 Int_t a = data->GetSi2Address(i);
132 if (a == 0 || a == 1279) continue; // noisy channels
133 mSi2ChP->Fill(a, q);
134 if (a < N) mSiCh[2]->Fill(q);
135 else mSiCh[3]->Fill(q);
136 }
137
138 mSiFullCh[0]->Fill(data->GetQx(0));
139 mSiFullCh[1]->Fill(data->GetQy(0));
140 mSiFullCh[2]->Fill(data->GetQx(1));
141 mSiFullCh[3]->Fill(data->GetQy(1));
142
143 for(int k=0; k<2; k++)
144 mChCor[k]->Fill(data->GetQx(k), data->GetQy(k));
145
146 /*
147 if (data->GetQx(0) < 20) continue;
148 if (data->GetQx(1) < 20) continue;
149 if (data->GetQy(0) < 20) continue;
150 if (data->GetQy(1) < 20) continue;
151 */
152
153 for(int k=0; k<2; k++)
154 mPos[k]->Fill(data->GetX(k), data->GetY(k));
155
156 mCor[0]->Fill(data->GetX(0), data->GetX(1));
157 mCor[1]->Fill(data->GetY(0), data->GetY(1));
158 }
159
160 delete data;
161 }
162 //}
163
164 // process histograms
165 for(int i=1; i<1281; i++) mSi1N->Fill(mSi1ChP->GetBinError(i));
166 for(int i=1; i<1281; i++) mSi2N->Fill(mSi2ChP->GetBinError(i));
167
168 // display
169 cout << "Number of Events = " << counter << endl;
170
171 /**/
172 TCanvas *c = new TCanvas("siliconSignal", "silicon signal");
173 c->Divide(2,2, 0.01, 0.01);
174 c->cd(1);
175 mSi1L->Draw();
176 c->cd(2);
177 mSi2L->Draw();
178 c->cd(3);
179 mSi1ChP->Draw();
180 c->cd(4);
181 mSi2ChP->Draw();
182 /* */
183
184 /**/
185 c = new TCanvas();
186 c->Divide(2,2,0.01,0.01);
187 c->cd(1);
188 mSi1N->Draw();
189 c->cd(2);
190 mSi2N->Draw();
191 /**/
192
193 /**/
194 // pads
195 c = new TCanvas("siPads", "Silicon Pads");
196 c->Divide(2,2, 0.01, 0.01);
197 for(int i=0; i<4; i++) {
198 c->cd(i+1);
199 gPad->SetLogy();
200 mSiCh[i]->Draw();
201 }
202
203 // clusters
204 c = new TCanvas("siCluster", "silicon clusters");
205 c->Divide(2,2, 0.01, 0.01);
206 for(int i=0; i<4; i++) {
207 c->cd(i+1);
208 gPad->SetLogy();
209 mSiFullCh[i]->Draw();
210 }
211
212 // position and correlation
213 c = new TCanvas("siPosition", "reconstructed position");
214 c->Divide(2,2, 0.01, 0.01);
215 for(int i=0; i<2; i++) {
216 c->cd(1+i);
217 mPos[i]->Draw("col");
218 c->cd(3+i);
219 mCor[i]->Draw("col");
220 }
221
222 c = new TCanvas("siCharge", "si charge correlation");
223 c->Divide(2,2, 0.01, 0.01);
224 c->cd(1);
225 gPad->SetLogz();
226 mChCor[0]->Draw("col");
227 c->cd(2);
228 gPad->SetLogz();
229 mChCor[1]->Draw("col");
230 /**/
231
232 new TCanvas();
233 gPad->SetLogy();
234 mCher->Draw();
235
236 // electron sample
237 int bin = mCher->FindBin(500.);
238 double ef = (mCher->Integral(bin, 151)/ mCher->GetSum());
239 TLatex *l = new TLatex(2e3, 0.02*mCher->GetSum(), Form("Electron fraction = %.2f ", ef));
240 l->Draw();
241
242 new TCanvas();
243 gPad->SetLogy();
244 mPb->Draw();
245
246 new TCanvas();
247 gPad->SetLogz();
248 mPbCher->Draw("colz");
249
250 /*
251 c = new TCanvas();
252 c->Divide(2,2,0.01, 0.01);
253 c->cd(1);
254 gPad->SetLogy();
255 mTRDsig->Draw();
256
257 c->cd(2);
258 mROB->Draw();
259
260 c->cd(3);
261 mDet->Draw();
262 */
263
264 st.Stop();
265 st.Print();
266}
267
268//Int_t SelectEvent() {
269//}