Correct formula for phi, in case of current local X <0, in GetLocalXat and GetPhiZat
[u/mrichter/AliRoot.git] / ITS / AnalyzeSDDNoiseAllMod.C
1 #if !defined(__CINT__) || defined(__MAKECINT__)
2 #include <TH2F.h>
3 #include <TCanvas.h>
4 #include <TStopwatch.h>
5 #include <TStyle.h>
6 #include <TGrid.h>
7 #include <TSystem.h>
8 #include <TLatex.h>
9 #include <TFile.h>
10 #include "AliRawReader.h"
11 #include "AliRawReaderDate.h"
12 #include "AliRawReaderRoot.h"
13 #include "AliITSOnlineSDDBase.h"
14 #include "AliITSOnlineSDDCMN.h"
15 #include "AliITSRawStreamSDD.h"
16 #include "TPaveStats.h"
17 #endif
18
19 // Macro for the analysis of PEDESTAL runs (equivalent to ITSSDDBASda.cxx)
20 // Two functions named AnalyzeSDDNoiseAllModules: 
21 // The first is for analyzing a local raw data file and takes as agrument the file name.
22 // The second is for running on ALIEN
23 // All DDLs are analyzed, the argument nDDL selects the DDL to be plotted
24 // Origin: F. Prino (prino@to.infn.it)
25
26 void AnalyzeSDDNoiseAllMod(Char_t *datafil, Int_t nDDL, Int_t firstEv=10, Int_t lastEv=12){
27
28   const Int_t kTotDDL=24;
29   const Int_t kModPerDDL=12;
30   const Int_t kSides=2;
31
32   AliITSOnlineSDDBase **base=new AliITSOnlineSDDBase*[kTotDDL*kModPerDDL*kSides];
33   TH2F **histo=new TH2F*[kTotDDL*kModPerDDL*kSides];
34
35   Char_t hisnam[20];
36   for(Int_t iddl=0; iddl<kTotDDL;iddl++){
37     for(Int_t imod=0; imod<kModPerDDL;imod++){
38       for(Int_t isid=0;isid<kSides;isid++){
39         Int_t index=kSides*(kModPerDDL*iddl+imod)+isid;
40         base[index]=new AliITSOnlineSDDBase(iddl,imod,isid);
41         sprintf(hisnam,"h%02dc%02ds%d",iddl,imod,isid);
42         histo[index]=new TH2F(hisnam,"",128,-0.5,127.5,256,-0.5,255.5);
43       }
44     }
45   }
46
47
48   TCanvas* c0 = new TCanvas("c0","Ev Display",900,900);
49   gStyle->SetPalette(1);
50   Char_t text[50];
51
52   Int_t iev=firstEv;
53   AliRawReader *rd; 
54   if(strstr(datafil,".root")!=0){
55     rd=new AliRawReaderRoot(datafil,iev);
56   }else{
57     rd=new AliRawReaderDate(datafil,iev);
58   }
59   TLatex *t0=new TLatex();
60   t0->SetNDC();
61   t0->SetTextSize(0.06);
62   t0->SetTextColor(4);
63
64   do{
65     c0->Clear();
66     c0->Divide(4,6,0.001,0.001);
67     printf("Event # %d ",iev);
68     rd->Reset();
69     for(Int_t iddl=0; iddl<kTotDDL;iddl++){
70       for(Int_t imod=0; imod<kModPerDDL;imod++){
71         for(Int_t isid=0;isid<kSides;isid++){
72           Int_t index=kSides*(kModPerDDL*iddl+imod)+isid;
73           histo[index]->Reset();
74         }
75       }
76     }
77     AliITSRawStreamSDD s(rd);
78     while(s.Next()){
79       Int_t iDDL=rd->GetDDLID();
80       Int_t iCarlos=s.GetCarlosId();
81       if(s.IsCompletedModule()) continue;
82       if(s.IsCompletedDDL()) continue;
83       if(iDDL>=0 && iDDL<kTotDDL){ 
84         Int_t index=kSides*(kModPerDDL*iDDL+iCarlos)+s.GetChannel(); 
85         histo[index]->Fill(s.GetCoord2(),s.GetCoord1(),s.GetSignal());
86       }
87     }
88     iev++;
89     for(Int_t iddl=0; iddl<kTotDDL;iddl++){
90       for(Int_t imod=0; imod<kModPerDDL;imod++){
91         for(Int_t isid=0;isid<kSides;isid++){
92           Int_t index=kSides*(kModPerDDL*iddl+imod)+isid;
93           base[index]->AddEvent(histo[index]);
94           if(iddl==nDDL){
95             Int_t index2=kSides*imod+isid;
96             c0->cd(index2+1);
97             histo[index]->DrawCopy("colz");
98             sprintf(text,"DDL %d channel %d Side %d",nDDL,imod,isid);
99             t0->DrawLatex(0.15,0.92,text);
100             c0->Update();
101           }
102         }
103       }
104     }
105     printf(" --- OK\n");
106   }while(rd->NextEvent()&&iev<=lastEv);
107
108   for(Int_t iddl=0; iddl<kTotDDL;iddl++){
109     for(Int_t imod=0; imod<kModPerDDL;imod++){
110       for(Int_t isid=0;isid<kSides;isid++){
111         Int_t index=kSides*(kModPerDDL*iddl+imod)+isid;
112         base[index]->ValidateAnodes();
113         base[index]->WriteToASCII(); // fondamentale!!!!!!!!!
114         delete base[index];
115       }
116     }
117   }
118   delete rd;
119   delete [] base;
120   
121   printf("Start second analysis for Common Mode correction\n");
122   AliITSOnlineSDDCMN **corr=new AliITSOnlineSDDCMN*[kTotDDL*kModPerDDL*kSides];
123   Bool_t isFilled[kTotDDL*kModPerDDL*kSides];
124
125   for(Int_t iddl=0; iddl<kTotDDL;iddl++){
126     for(Int_t imod=0; imod<kModPerDDL;imod++){
127       for(Int_t isid=0;isid<kSides;isid++){
128         Int_t index=kSides*(kModPerDDL*iddl+imod)+isid;
129         corr[index]=new AliITSOnlineSDDCMN(iddl,imod,isid);
130         isFilled[index]=0;
131       }
132     }
133   }
134
135   iev=firstEv;
136   AliRawReader *rd2; 
137   if(strstr(datafil,".root")!=0){
138     rd2=new AliRawReaderRoot(datafil,iev);
139   }else{
140     rd2=new AliRawReaderDate(datafil,iev);
141   }
142   do{
143     c0->Clear();
144     c0->Divide(4,6,0.001,0.001);
145     printf("Event # %d ",iev);
146     rd2->Reset();
147     for(Int_t iddl=0; iddl<kTotDDL;iddl++){
148       for(Int_t imod=0; imod<kModPerDDL;imod++){
149         for(Int_t isid=0;isid<kSides;isid++){
150           Int_t index=kSides*(kModPerDDL*iddl+imod)+isid;
151           histo[index]->Reset();
152         }
153       }
154     }
155     
156     AliITSRawStreamSDD s(rd2);
157     while(s.Next()){
158       Int_t iDDL=rd2->GetDDLID();
159       Int_t iCarlos=s.GetCarlosId();
160       if(s.IsCompletedModule()) continue;
161       if(s.IsCompletedDDL()) continue;
162       if(iDDL>=0 && iDDL<kTotDDL){ 
163         Int_t index=kSides*(kModPerDDL*iDDL+iCarlos)+s.GetChannel(); 
164         histo[index]->Fill(s.GetCoord2(),s.GetCoord1(),s.GetSignal());
165         isFilled[index]=1;
166       }
167     }
168     iev++;
169     for(Int_t iddl=0; iddl<kTotDDL;iddl++){
170       for(Int_t imod=0; imod<kModPerDDL;imod++){
171         for(Int_t isid=0;isid<kSides;isid++){
172           Int_t index=kSides*(kModPerDDL*iddl+imod)+isid;
173           if(isFilled[index]) corr[index]->AddEvent(histo[index]);
174           if(iddl==nDDL){
175             Int_t index2=kSides*imod+isid;
176             c0->cd(index2+1);
177             histo[index]->DrawCopy("colz");
178             sprintf(text,"DDL %d channel %d Side %d",nDDL,imod,isid);
179             t0->DrawLatex(0.15,0.92,text);
180             c0->Update();
181           }
182         }
183       }
184     }
185     printf(" --- OK\n");
186   }while(rd2->NextEvent()&&iev<=lastEv);
187
188   TH1F *htotbas=new TH1F("htotbas","",100,0.,150.);
189   TH1F *htotbaseq=new TH1F("htotbaseq","",100,0.,150.);
190   TH1F *htotnoise=new TH1F("htotnoise","",100,0.,10.);
191   TH1F *htotnoisecorr=new TH1F("htotnoisecorr","",100,0.,10.);
192   TH1F *hstatus=new TH1F("hstatus","",2,-0.5,1.5);
193
194   TFile *outfil=new TFile("SDDbase-results.root","recreate");
195   for(Int_t iddl=0; iddl<kTotDDL;iddl++){
196     for(Int_t imod=0; imod<kModPerDDL;imod++){
197       for(Int_t isid=0;isid<kSides;isid++){
198         Int_t index=kSides*(kModPerDDL*iddl+imod)+isid;
199         if(isFilled[index]){
200           corr[index]->ValidateAnodes();
201           corr[index]->WriteToASCII();
202           corr[index]->WriteToROOT(outfil);
203           for(Int_t ian=0; ian<256;ian++){
204             Float_t basl=corr[index]->GetAnodeBaseline(ian);
205             Float_t basleq=corr[index]->GetAnodeEqualizedBaseline(ian);
206             Float_t noi=corr[index]->GetAnodeRawNoise(ian);
207             Float_t cornoi=corr[index]->GetAnodeCorrNoise(ian);
208             Int_t anstatus=corr[index]->IsAnodeGood(ian);
209             hstatus->Fill(anstatus);
210             htotbas->Fill(basl);
211             htotbaseq->Fill(basleq);
212             htotnoise->Fill(noi);
213             htotnoisecorr->Fill(cornoi);
214           }
215         }
216       }
217     }
218   }
219   outfil->Close();
220   // Draw Statistics of baselines and noise
221   TCanvas *call=new TCanvas("call","General stats",700,700);
222   call->Divide(2,2);
223   call->cd(1);
224   htotbas->Draw();
225   htotbas->GetXaxis()->SetTitle("Baselines");
226   htotbas->GetXaxis()->SetTitleSize(0.07);
227   htotbas->GetXaxis()->SetTitleOffset(0.6);
228   call->cd(2);
229   htotbaseq->Draw();
230   htotbaseq->GetXaxis()->SetTitle("Baselines after equalization");
231   htotbaseq->GetXaxis()->SetTitleSize(0.07);
232   htotbaseq->GetXaxis()->SetTitleOffset(0.6);
233   call->cd(3);
234   htotnoisecorr->SetLineColor(2);
235   htotnoisecorr->Draw();
236   call->Update();
237   TPaveStats *st1=(TPaveStats*)htotnoisecorr->GetListOfFunctions()->FindObject("stats");
238   st1->SetY1NDC(0.51);
239   st1->SetY2NDC(0.7);
240   htotnoisecorr->GetXaxis()->SetTitle("Noise");
241   htotnoisecorr->GetXaxis()->SetTitleSize(0.07);
242   htotnoisecorr->GetXaxis()->SetTitleOffset(0.6);
243   htotnoise->Draw("SAMES");
244   call->Update();
245   TPaveStats *st2=(TPaveStats*)htotnoise->GetListOfFunctions()->FindObject("stats");
246   st2->SetY1NDC(0.71);
247   st2->SetY2NDC(0.9);
248   
249   call->cd(4);
250   hstatus->Draw();
251   hstatus->GetXaxis()->SetTitle("Anode Status (0=bad 1=good)");
252   hstatus->GetXaxis()->SetTitleSize(0.07);
253   hstatus->GetXaxis()->SetTitleOffset(0.6);
254   call->Update();
255   call->SaveAs("GenStatsPedestal.gif");
256
257   // Draw baselines and noise for all modules of the selected DDL
258
259   TH1F** hbas = new TH1F*[kSides*kModPerDDL];
260   TH1F** hrawn = new TH1F*[kSides*kModPerDDL];
261   TH1F** hcorrn = new TH1F*[kSides*kModPerDDL];
262   TH1F** hdbas = new TH1F*[kSides*kModPerDDL];
263   TH1F** hdrawn = new TH1F*[kSides*kModPerDDL];
264   TH1F** hdcorrn = new TH1F*[kSides*kModPerDDL];
265   TCanvas *c1=new TCanvas("c1","DDL: Baselines vs anode",900,900);
266   c1->SetBottomMargin(0.14);
267   c1->Divide(4,6,0.001,0.001);
268   TCanvas *c2=new TCanvas("c2","DDL: Noise vs anode",900,900);
269   c2->SetBottomMargin(0.14);
270   c2->Divide(4,6,0.001,0.001);
271   TCanvas *c3=new TCanvas("c3","DDL: Baselines distr",900,900);
272   c3->SetBottomMargin(0.14);
273   c3->Divide(4,6,0.001,0.001);
274   TCanvas *c4=new TCanvas("c4","DDL: Noise Distr",900,900);
275   c4->SetBottomMargin(0.14);
276   c4->Divide(4,6,0.001,0.001);
277   TLatex *t1=new TLatex(0.15,0.2,"Raw Noise");
278   t1->SetNDC();
279   t1->SetTextSize(0.05);
280   TLatex *t2=new TLatex(0.4,0.2,"Corrected Noise");
281   t2->SetNDC();
282   t2->SetTextSize(0.05);
283   t2->SetTextColor(2);
284   TLatex *t3=new TLatex();
285   t3->SetNDC();
286   t3->SetTextSize(0.06);
287   t3->SetTextColor(4);
288
289   for(Int_t imod=0; imod<kModPerDDL;imod++){
290     for(Int_t isid=0;isid<kSides;isid++){
291       Int_t index1=kSides*(kModPerDDL*nDDL+imod)+isid;
292       Int_t index2=kSides*imod+isid;
293       sprintf(text,"DDL %d channel %d Side %d",nDDL,imod,isid);
294       hbas[index2]=corr[index1]->GetBaselineAnodeHisto();
295       hrawn[index2]=corr[index1]->GetRawNoiseAnodeHisto();
296       hcorrn[index2]=corr[index1]->GetCorrNoiseAnodeHisto();
297       hdbas[index2]=corr[index1]->GetBaselineHisto();
298       hdrawn[index2]=corr[index1]->GetRawNoiseHisto();
299       hdcorrn[index2]=corr[index1]->GetCorrNoiseHisto();
300       c1->cd(index2+1);
301       hbas[index2]->Draw();
302       hbas[index2]->SetMinimum(0);
303       hbas[index2]->SetMaximum(75);
304       hbas[index2]->GetXaxis()->SetTitle("Anode");
305       hbas[index2]->GetYaxis()->SetTitle("Baseline");
306       hbas[index2]->GetXaxis()->SetTitleSize(0.07);
307       hbas[index2]->GetYaxis()->SetTitleSize(0.07);
308       hbas[index2]->GetXaxis()->SetTitleOffset(0.6);
309       hbas[index2]->GetYaxis()->SetTitleOffset(0.7);
310       t3->DrawLatex(0.15,0.92,text);
311       c1->Update();
312
313
314       c2->cd(index2+1); 
315       hrawn[index2]->SetMinimum(1.);
316       hrawn[index2]->SetMaximum(6.);
317       hrawn[index2]->Draw();
318       hrawn[index2]->GetXaxis()->SetTitle("Anode");
319       hrawn[index2]->GetYaxis()->SetTitle("Noise");
320       hrawn[index2]->GetXaxis()->SetTitleSize(0.07);
321       hrawn[index2]->GetYaxis()->SetTitleSize(0.07);
322       hrawn[index2]->GetXaxis()->SetTitleOffset(0.6);
323       hrawn[index2]->GetYaxis()->SetTitleOffset(0.7);
324       gStyle->SetOptStat(0);
325       hrawn[index2]->SetStats(0);
326       hcorrn[index2]->SetLineColor(2);
327       hcorrn[index2]->Draw("SAME");
328       t1->Draw();
329       t2->Draw();
330       t3->DrawLatex(0.15,0.92,text);
331       c2->Update();
332
333       c3->cd(index2+1);
334       hdbas[index2]->Draw();
335       hdbas[index2]->GetXaxis()->SetTitle("Baseline");
336       hdbas[index2]->GetXaxis()->SetTitleSize(0.07);
337       hdbas[index2]->GetXaxis()->SetTitleOffset(0.6);
338       t3->DrawLatex(0.15,0.92,text);
339       c3->Update();
340
341       c4->cd(index2+1); 
342       hdrawn[index2]->Draw();
343       hdrawn[index2]->GetXaxis()->SetTitle("Noise");
344       hdrawn[index2]->GetXaxis()->SetTitleSize(0.07);
345       hdrawn[index2]->GetXaxis()->SetTitleOffset(0.6);
346       hdcorrn[index2]->SetLineColor(2);
347       hdcorrn[index2]->Draw("SAME");
348       t1->Draw();
349       t2->Draw();
350       t3->DrawLatex(0.15,0.92,text);
351       c4->Update();
352     }
353   }
354
355   c1->SaveAs("Baselines.gif");
356   c2->SaveAs("Noise.gif");
357   c3->SaveAs("BaselinesDist.gif");
358   c4->SaveAs("NoiseDist.gif");
359   
360   Char_t delfil[100];
361   for(Int_t iddl=0; iddl<kTotDDL;iddl++){
362     for(Int_t imod=0; imod<kModPerDDL;imod++){
363       for(Int_t isid=0;isid<kSides;isid++){
364         sprintf(delfil,"rm SDDbase_step1_ddl%02dc%02d_sid%d.data",iddl,imod,isid);
365         gSystem->Exec(delfil);
366       }
367     }
368   }
369
370 }
371
372 void AnalyzeSDDNoiseAllMod(Int_t nrun, Int_t n2, Char_t* dir="LHC08d_SDD",Int_t nDDL=0, Int_t firstEv=15, Int_t lastEv=18){
373   TGrid::Connect("alien:",0,0,"t");
374   Char_t filnam[200];
375   sprintf(filnam,"alien:///alice/data/2008/%s/%09d/raw/08%09d%03d.10.root",dir,nrun,nrun,n2);
376   printf("Open file %s\n",filnam);
377   AnalyzeSDDNoiseAllMod(filnam,nDDL,firstEv,lastEv);
378 }