Change of module numbers for display purposes
[u/mrichter/AliRoot.git] / ITS / macrosSDD / AnalyzeSDDInjectorsAllMod.C
1 #if !defined(__CINT__) || defined(__MAKECINT__)
2 #include <TH2F.h>
3 #include <TCanvas.h>
4 #include <TGraphErrors.h>
5 #include <TStopwatch.h>
6 #include <TStyle.h>
7 #include <TLatex.h>
8 #include <TFile.h>
9 #include <TMath.h>
10 #include <TNtuple.h>
11 #include <TGrid.h>
12 #include <TF1.h>
13 #include <TLine.h>
14 #include "AliRawReader.h"
15 #include "AliRawReaderDate.h"
16 #include "AliRawReaderRoot.h"
17 #include "AliITSOnlineSDDInjectors.h"
18 #include "AliITSRawStreamSDD.h"
19 #include "AliITSRawStreamSDDCompressed.h"
20 #include "AliITSDDLModuleMapSDD.h"
21 #endif
22
23 // Macro for the analysis of PULSER runs (equivalent to ITSSDDINJda.cxx)
24 // Two functions named AnalyzeSDDInjectorsAllModules: 
25 // The first is for analyzing a local raw data file and takes as agrument the file name.
26 // The second is for running on ALIEN
27 // All DDLs are analyzed, the argument nDDL selects the DDL to be plotted
28 // Origin: F. Prino (prino@to.infn.it)
29
30
31 void AnalyzeSDDInjectorsAllMod(Char_t *datafil, 
32                                Int_t adcfreq=20, 
33                                Int_t nDDL=0, 
34                                Int_t firstEv=18, 
35                                Int_t lastEv=30,
36                                Int_t jpad=20, 
37                                Int_t statuscut=7){
38
39
40   const Int_t kTotDDL=24;
41   const Int_t kModPerDDL=12;
42   const Int_t kSides=2;
43   Bool_t writtenoutput=kFALSE;
44
45   AliITSDDLModuleMapSDD* dmap=new AliITSDDLModuleMapSDD();
46   dmap->SetJun09Map();
47
48   TNtuple* ntsp=new TNtuple("ntsp","","mod:sid:an:stat:vall:errvall:v23:v13:v12:c1:c2:c3");
49   Float_t xnt[12];
50   TH2F** histo = new TH2F*[kTotDDL*kModPerDDL*kSides];
51   Int_t nWrittenEv[kTotDDL*kModPerDDL*kSides];
52   TGraphErrors** gvel = new TGraphErrors*[kTotDDL*kModPerDDL*kSides];
53   AliITSOnlineSDDInjectors **anal=new AliITSOnlineSDDInjectors*[kTotDDL*kModPerDDL*kSides];
54   AliITSOnlineSDDInjectors **anal23=new AliITSOnlineSDDInjectors*[kTotDDL*kModPerDDL*kSides];
55   AliITSOnlineSDDInjectors **anal13=new AliITSOnlineSDDInjectors*[kTotDDL*kModPerDDL*kSides];
56   AliITSOnlineSDDInjectors **anal12=new AliITSOnlineSDDInjectors*[kTotDDL*kModPerDDL*kSides];
57   TH1F** hvdriftl=new TH1F*[260];  
58   TH1F** hvdriftr=new TH1F*[260];  
59   Char_t hisnam[20];
60   for(Int_t idet=0; idet<260;idet++){
61     sprintf(hisnam,"vdriftl%03d",idet);
62     hvdriftl[idet]=new TH1F(hisnam,"",500,5.5,8.0);
63     sprintf(hisnam,"vdriftr%03d",idet);
64     hvdriftr[idet]=new TH1F(hisnam,"",500,5.5,8.0);
65   }
66   for(Int_t iddl=0; iddl<kTotDDL;iddl++){
67     for(Int_t imod=0; imod<kModPerDDL;imod++){
68       for(Int_t isid=0;isid<kSides;isid++){
69         Int_t index=kSides*(kModPerDDL*iddl+imod)+isid;
70         sprintf(hisnam,"h%02dc%02ds%d",iddl,imod,isid);
71
72         histo[index]=new TH2F(hisnam,"",256,-0.5,255.5,256,-0.5,255.5);
73         anal[index]=new AliITSOnlineSDDInjectors(iddl,imod,isid);
74         if(adcfreq==40) anal[index]->Set40MHzConfig();
75         else anal[index]->Set20MHzConfig();
76
77         anal23[index]=new AliITSOnlineSDDInjectors(iddl,imod,isid);
78         if(adcfreq==40) anal23[index]->Set40MHzConfig();
79         else anal23[index]->Set20MHzConfig();
80         anal23[index]->SetUseLine(0,kFALSE);
81
82         anal13[index]=new AliITSOnlineSDDInjectors(iddl,imod,isid);
83         if(adcfreq==40) anal13[index]->Set40MHzConfig();
84         else anal13[index]->Set20MHzConfig();
85         anal13[index]->SetUseLine(1,kFALSE);
86
87         anal12[index]=new AliITSOnlineSDDInjectors(iddl,imod,isid);
88         if(adcfreq==40) anal12[index]->Set40MHzConfig();
89         else anal12[index]->Set20MHzConfig();
90         anal12[index]->SetUseLine(2,kFALSE);
91
92         nWrittenEv[index]=0;
93       }
94     }
95   }
96   TGraphErrors *gvvsmod0=new TGraphErrors(0);
97   TGraphErrors *gvvsmod1=new TGraphErrors(0);
98   TGraphErrors *gtvsmod0=new TGraphErrors(0);
99   TGraphErrors *gtvsmod1=new TGraphErrors(0);
100   Float_t gvmin=6.0, gvmax=7.5;
101   Float_t gtmin=288., gtmax=308.;
102   TH1F* hanst=new TH1F("hanst","",8,-0.5,7.5);
103   TH1F* hpad7l=new TH1F("hpad7l","",33,-0.5,32.5);
104   TH1F* hpad7r=new TH1F("hpad7r","",33,-0.5,32.5);
105
106   TCanvas* c0 = new TCanvas("c0","Event display",900,900);
107   gStyle->SetPalette(1);
108   TCanvas* c1 = new TCanvas("c1","Drift Speed vs. anode",900,900);
109   Char_t text[50];
110   UInt_t timeSt=0;
111
112   Int_t iev=firstEv;
113   AliRawReader *rd; 
114   if(strstr(datafil,".root")!=0){
115     rd=new AliRawReaderRoot(datafil,iev);
116   }else{
117     rd=new AliRawReaderDate(datafil,iev);
118   }
119
120   Char_t gname[15];
121   TF1 *funz=new TF1("funz","[0]+[1]*x+[2]*x*x+[3]*x*x*x",0.,255.);
122   TLatex *t0=new TLatex();
123   t0->SetNDC();
124   t0->SetTextSize(0.06);
125   t0->SetTextColor(4);
126   Int_t readEv=0;
127   do{
128     c0->Clear();
129     c0->Divide(4,6,0.001,0.001);
130     c1->Clear();
131     c1->Divide(4,6,0.001,0.001);
132     printf("Event # %d\n",iev);
133     timeSt=rd->GetTimestamp();
134     rd->Reset();
135     for(Int_t iddl=0; iddl<kTotDDL;iddl++){
136       for(Int_t imod=0; imod<kModPerDDL;imod++){
137         for(Int_t isid=0;isid<kSides;isid++){
138           Int_t index=kSides*(kModPerDDL*iddl+imod)+isid;
139           histo[index]->Reset();
140         }
141       }
142     }
143
144
145     UChar_t cdhAttr=AliITSRawStreamSDD::ReadBlockAttributes(rd);
146     UInt_t amSamplFreq=AliITSRawStreamSDD::ReadAMSamplFreqFromCDH(cdhAttr);
147     AliITSRawStream* s=AliITSRawStreamSDD::CreateRawStreamSDD(rd,cdhAttr);
148     if(!writtenoutput){
149       printf("Use %s raw stream, sampling frequency %d MHz\n",s->ClassName(),amSamplFreq);
150       writtenoutput=kTRUE;
151     }
152     while(s->Next()){
153       Int_t iDDL=rd->GetDDLID();
154       Int_t iCarlos=s->GetCarlosId();
155       if(s->IsCompletedModule()) continue;
156       if(s->IsCompletedDDL()) continue;
157       if(iDDL>=0 && iDDL<kTotDDL){ 
158         Int_t index=kSides*(kModPerDDL*iDDL+iCarlos)+s->GetChannel(); 
159         histo[index]->Fill(s->GetCoord2(),s->GetCoord1(),s->GetSignal());
160       }
161     }
162     
163     for(Int_t iddl=0; iddl<kTotDDL;iddl++){
164       for(Int_t imod=0; imod<kModPerDDL;imod++){
165         for(Int_t isid=0;isid<kSides;isid++){
166           Int_t index=kSides*(kModPerDDL*iddl+imod)+isid;
167           anal[index]->AnalyzeEvent(histo[index]); 
168           anal23[index]->AnalyzeEvent(histo[index]); 
169           anal13[index]->AnalyzeEvent(histo[index]); 
170           anal12[index]->AnalyzeEvent(histo[index]); 
171           nWrittenEv[index]++;
172           Int_t iMod=dmap->GetModuleNumber(iddl,imod);
173           if(iMod!=-1){     
174             for(Int_t ipad=0;ipad<33;ipad++){
175               Int_t st=anal[index]->GetInjPadStatus(ipad);
176               hanst->Fill(st);
177               if(anal[index]->GetInjPadStatus(ipad)>=statuscut){
178                 if(isid==0) hpad7l->Fill(ipad);
179                 if(isid==1) hpad7r->Fill(ipad);
180               }
181               xnt[0]=(Float_t)iMod;
182               xnt[1]=(Float_t)isid;
183               xnt[2]=(Float_t)anal[index]->GetAnodeNumber(ipad);
184               xnt[3]=(Float_t)st;
185               xnt[4]=anal[index]->GetDriftSpeed(ipad);
186               xnt[5]=anal[index]->GetDriftSpeedErr(ipad);
187               xnt[6]=anal23[index]->GetDriftSpeed(ipad);
188               xnt[7]=anal13[index]->GetDriftSpeed(ipad);
189               xnt[8]=anal12[index]->GetDriftSpeed(ipad);
190               xnt[9]=anal[index]->GetCentroid(ipad,0);
191               xnt[10]=anal[index]->GetCentroid(ipad,1);
192               xnt[11]=anal[index]->GetCentroid(ipad,2);
193               ntsp->Fill(xnt);
194             }
195             if(anal[index]->GetInjPadStatus(jpad)>=statuscut){
196               Float_t vel=anal[index]->GetDriftSpeed(jpad);
197               if(isid==0) hvdriftl[iMod-240]->Fill(vel);
198               if(isid==1) hvdriftr[iMod-240]->Fill(vel);
199             }
200           }
201           if(iddl==nDDL){
202             Int_t index2=kSides*imod+isid;
203             c0->cd(index2+1);
204             histo[index]->SetMaximum(100.);
205             histo[index]->DrawCopy("colz");
206             sprintf(text,"DDL %d channel %d Side %d",nDDL,imod,isid);
207             t0->DrawLatex(0.15,0.92,text);
208             c0->Update();
209             c1->cd(index2+1);
210             gvel[index]=anal[index]->GetDriftSpeedGraph();
211             gvel[index]->SetMarkerStyle(20);
212             gvel[index]->SetTitle("");
213             sprintf(gname,"gvel%dev%d",index,iev);
214             gvel[index]->SetName(gname);
215             //      gvel[index]->SetMinimum(0);
216             //gvel[index]->SetMaximum(200);
217
218             gvel[index]->GetXaxis()->SetLimits(0,256);
219             gvel[index]->GetXaxis()->SetTitle("Anode");
220             gvel[index]->GetYaxis()->SetTitle("Drift vel.");
221             gvel[index]->GetXaxis()->SetTitleSize(0.07);
222             gvel[index]->GetYaxis()->SetTitleSize(0.07);
223             gvel[index]->GetXaxis()->SetTitleOffset(0.6);
224             gvel[index]->GetYaxis()->SetTitleOffset(0.6);
225             if(gvel[index]->GetN()>0) gvel[index]->Draw("AP");
226             Double_t *param=anal[index]->GetDriftSpeedFitParam();
227             funz->SetParameters(param[0],param[1],param[2],param[3]);
228             funz->SetLineColor(2);
229             funz->DrawCopy("LSAME");
230             t0->DrawLatex(0.15,0.92,text);
231             c1->Update();
232           }
233         }
234       }
235     }
236     delete s;
237     iev++;
238     readEv++;
239     printf(" --- OK\n");
240   }while(rd->NextEvent()&&iev<=lastEv);
241   printf("Total number of events = %d\n",readEv);
242   Float_t nfac=1./(Float_t)readEv/33./520.;
243   hanst->Scale(nfac);
244   nfac=1./(Float_t)readEv;
245   hpad7l->Scale(nfac);
246   hpad7r->Scale(nfac);
247
248   TFile *outfil1=new TFile("DriftSpeedVsAnode.root","recreate");  
249   for(Int_t iddl=0; iddl<kTotDDL;iddl++){
250     for(Int_t imod=0; imod<kModPerDDL;imod++){
251       for(Int_t isid=0;isid<kSides;isid++){
252         Int_t index=kSides*(kModPerDDL*iddl+imod)+isid;
253         anal[index]->FitMeanDriftSpeedVsAnode();
254         anal[index]->WriteToASCII(0,timeSt,0);
255         anal[index]->WriteInjectorStatusToASCII();
256         anal[index]->WriteToROOT(outfil1);
257       }
258     }
259   }
260   outfil1->Close();
261
262   Int_t ipt0=0, ipt1=0;
263   Float_t Edrift=(1800-45)/291/0.012;  
264   TFile *outfil=new TFile("DriftSpeedHistos.root","recreate");
265   ntsp->Write();
266   for(Int_t iMod=0; iMod<260; iMod++){
267     outfil->cd();
268     //    hvdriftl[iMod]->Write();    
269     //    hvdriftr[iMod]->Write();
270     Float_t modid=iMod+240;
271     if(hvdriftl[iMod]->GetEntries()>0){
272       Float_t avevell=hvdriftl[iMod]->GetMean();
273       Float_t rmsvell=hvdriftl[iMod]->GetRMS();
274       if(avevell > 5.5 && avevell < 8.5){
275         gvvsmod0->SetPoint(ipt0,modid,avevell);
276         gvvsmod0->SetPointError(ipt0,0,rmsvell);
277         Float_t mob=avevell*1.E5/Edrift;  
278         Float_t temper=293.15*TMath::Power((mob/1350.),-1/2.4); 
279         gtvsmod0->SetPoint(ipt0,modid,temper);
280         ++ipt0;
281       }
282     }
283     if(hvdriftr[iMod]->GetEntries()>0){
284       Float_t avevelr=hvdriftr[iMod]->GetMean();
285       Float_t rmsvelr=hvdriftr[iMod]->GetRMS();
286       if(avevelr > 5.5 && avevelr < 8.5){
287         gvvsmod1->SetPoint(ipt1,modid,avevelr);
288         gvvsmod1->SetPointError(ipt1,0,rmsvelr);
289         Float_t mob=avevelr*1.E5/Edrift;
290         Float_t temper=293.15*TMath::Power((mob/1350.),-1./2.4); 
291         gtvsmod1->SetPoint(ipt1,modid,temper);
292         ++ipt1;
293       }
294     }
295   }
296   gvvsmod0->SetName("gvvsmod0");
297   gvvsmod1->SetName("gvvsmod1");
298   gtvsmod0->SetName("gtvsmod0");
299   gtvsmod1->SetName("gtvsmod1");
300   outfil->cd();
301   gvvsmod0->Write();
302   gvvsmod1->Write();
303   gtvsmod0->Write();
304   gtvsmod1->Write();
305   outfil->Close();
306
307   TCanvas* c8=new TCanvas("c8","Drift Speed vs. mod");
308   gvvsmod0->SetTitle("");
309   gvvsmod1->SetTitle("");
310   gvvsmod0->SetMarkerStyle(20);
311   gvvsmod1->SetMarkerStyle(21);
312   gvvsmod1->SetMarkerColor(2);
313   gvvsmod0->Draw("AP");
314   gvvsmod0->SetMinimum(gvmin);
315   gvvsmod0->SetMaximum(gvmax);
316   gvvsmod0->GetXaxis()->SetTitle("Module Number");
317   Char_t title[50];
318   sprintf(title,"Vdrift at injector pad %d",jpad);
319   gvvsmod0->GetYaxis()->SetTitle(title);  
320   gvvsmod1->Draw("PSAME");
321   TLatex* tleft=new TLatex(0.7,0.82,"Side 0");
322   tleft->SetNDC();
323   tleft->SetTextColor(1);
324   tleft->Draw();
325   TLatex* tright=new TLatex(0.7,0.75,"Side 1");
326   tright->SetNDC();
327   tright->SetTextColor(2);
328   tright->Draw();
329
330   TLine *lin=new TLine(323,gvmin,323,gvmax);
331   lin->SetLineColor(4);
332   lin->Draw();
333   c8->Update();
334   c8->SaveAs("VdriftVsMod.gif");
335
336   TCanvas* c8t=new TCanvas("c8t","Temeprature vs. mod");
337   gtvsmod0->SetTitle("");
338   gtvsmod1->SetTitle("");
339   gtvsmod0->SetMarkerStyle(20);
340   gtvsmod1->SetMarkerStyle(21);
341   gtvsmod1->SetMarkerColor(2);
342   gtvsmod0->Draw("AP");
343   gtvsmod0->SetMinimum(gtmin);
344   gtvsmod0->SetMaximum(gtmax);
345   gtvsmod0->GetXaxis()->SetTitle("Module Number");
346   sprintf(title,"Estimated Temperature (K)");
347   gtvsmod0->GetYaxis()->SetTitle(title);  
348   gtvsmod1->Draw("PSAME");
349   tleft->Draw();
350   tright->Draw();
351   TLine *lint=new TLine(323,gtmin,323,gtmax);
352   lint->SetLineColor(4);
353   lint->Draw();
354   c8t->Update();
355   c8t->SaveAs("TempVsMod.gif");
356
357   TCanvas* c9=new TCanvas("c9","Injector status");
358   hanst->SetStats(0);
359   hanst->Draw();
360   hanst->GetXaxis()->SetTitle("Injector pad status");
361   hanst->GetXaxis()->CenterTitle();
362   c9->SaveAs("InjStatus.gif");
363
364 //   TCanvas* c10=new TCanvas("c10","Pad status 7",1200,600);
365 //   hpad7l->SetStats(0);
366 //   hpad7r->SetStats(0);
367 //   c10->Divide(2,1);
368 //   c10->cd(1);
369 //   hpad7l->Draw(); 
370 //   hpad7l->GetXaxis()->SetTitle("Side Left -- Pad number");
371 //   hpad7l->GetXaxis()->CenterTitle();
372 //   hpad7l->GetYaxis()->SetTitle("Number of status 7");
373 //   c10->cd(2);
374 //   hpad7r->Draw();
375 //   hpad7r->GetXaxis()->SetTitle("Side Right -- Pad number");
376 //   hpad7r->GetXaxis()->CenterTitle();
377 //   hpad7r->GetYaxis()->SetTitle("Number of status 7");
378 //   printf("Side 0, maximum pad=%d\n",hpad7l->GetMaximumBin());
379 //   printf("Side 1, maximum pad=%d\n",hpad7r->GetMaximumBin());
380   
381
382 }
383
384 void AnalyzeSDDInjectorsAllMod(Int_t nrun, Int_t n2, Int_t year=2010, Char_t* dir="LHC10b_SDD",
385                                Int_t adcfreq=20, 
386                                Int_t nDDL=0, 
387                                Int_t firstEv=18, 
388                                Int_t lastEv=25,
389                                Int_t jpad=20, 
390                                Int_t statuscut=7){
391
392   TGrid::Connect("alien:",0,0,"t");
393   Char_t filnam[200];
394   sprintf(filnam,"alien:///alice/data/%d/%s/%09d/raw/%02d%09d%03d.10.root",year,dir,nrun,year-2000,nrun,n2);
395   printf("Open file %s\n",filnam);
396   AnalyzeSDDInjectorsAllMod(filnam,adcfreq,nDDL,firstEv,lastEv,jpad,statuscut);
397 }
398
399
400