Update of calibration task and replace AliTRDrunCalib.C by AddTaskTRDCalib.C (Raphaelle)
[u/mrichter/AliRoot.git] / TRD / Macros / AliTRDanalyzeTestBeam.C
1
2 void 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
64   AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
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 //}