Update of 2 SSD macros (Panos)
[u/mrichter/AliRoot.git] / ITS / readSSDCommonMode.C
1 #if !defined(__CINT__) || defined(__MAKECINT__)
2 #include "TString.h"
3 #include "TH1.h"
4 #include "TH2.h"
5 #include "TFile.h"
6 #include "TList.h"
7 #include "TChain.h"
8 #include "TObject.h"
9 #include "TCanvas.h"
10 #include "TStyle.h"
11 #include "TSystem.h"
12 #include "Riostream.h"
13
14 #include "AliITSgeomTGeo.h"
15 #include "AliLog.h"
16 #include "AliRawReaderRoot.h"
17 #include "AliITSRawStreamSSD.h"
18 #include "AliCDBManager.h"
19 #include "AliCDBEntry.h"
20 #endif
21
22 /*  $Id:  $    */
23
24 //============================================================//
25 static const Int_t fgkNumOfChips = 6;      //number of SSD chips per module per side
26 static const Int_t fgkNumOfLDCs = 8;      //number of SSD LDCs
27 static const Int_t fgkNumOfDDLs = 16;      //number of SSD DDLs
28 static const Int_t fgkSSDMODULES = 1698;      //total number of SSD modules
29 static const Int_t fgkSSDLADDERSLAYER5 = 34; //ladders on layer 5
30 static const Int_t fgkSSDLADDERSLAYER6 = 38; //ladders on layer 6
31 static const Int_t fgkSSDMODULESPERLADDERLAYER5 = 22; //modules per ladder - layer 5                           
32 static const Int_t fgkSSDMODULESPERLADDERLAYER6 = 25; //modules per ladder - layer 6                          
33 static const Int_t fgkSSDMODULESLAYER5 = 748; //total number of SSD modules - layer5                           
34 static const Int_t fgkSSDMODULESLAYER6 = 950; //total number of SSD modules - layer6                         
35 static const Int_t fgkNumberOfPSideStrips = 768; //number of P-side strips
36 //============================================================//
37
38 TList *initCM();
39 void makeCM(const char* filename, Int_t nEvents, TList *list);
40
41 //__________________________________________________________//
42 void readSSDCommonMode(const char* filename = "raw.root",
43                        Int_t nEvents = -1) {
44   //Reads the CM pseudo-channels and produces the CM map for both layers 
45   //and for p and n-side.
46   gStyle->SetPalette(1,0);
47
48   TList *list = initCM();
49   //list->ls();
50   Printf("CM histograms: %d",list->GetEntries());
51   makeCM(filename,nEvents,list);
52
53  
54 //__________________________________________________________//
55 TList *initCM() {
56   //Initializes the histograms and returns the TList object
57   TList *list = new TList();
58
59   Int_t gLayer = 0,gLadder = 0, gModule = 0;
60   Int_t gHistCounter = 0;
61   TString gTitle;
62   TH1F *gHistSSDCMModule[2*fgkSSDMODULES][fgkNumOfChips]; 
63   for(Int_t iModule = 500; iModule < fgkSSDMODULES + 500; iModule++) {
64     AliITSgeomTGeo::GetModuleId(iModule,gLayer,gLadder,gModule);
65     for(Int_t iChip = 0; iChip < fgkNumOfChips; iChip++) {
66       gTitle = "SSD_CM_PSide_Layer"; gTitle += gLayer;
67       gTitle += "_Ladder"; gTitle += gLadder;
68       gTitle += "_Module"; gTitle += gModule;
69       gTitle += "_Chip"; gTitle += iChip+1;
70       gHistSSDCMModule[gHistCounter][iChip] = new TH1F(gTitle.Data(),
71                                                        gTitle.Data(),
72                                                        100,-50.,50.);
73       gHistSSDCMModule[gHistCounter][iChip]->GetXaxis()->SetTitle("CM");
74       list->Add(gHistSSDCMModule[gHistCounter][iChip]);
75     }
76     gHistCounter += 1;
77   }
78   for(Int_t iModule = 500; iModule < fgkSSDMODULES + 500; iModule++) {
79     AliITSgeomTGeo::GetModuleId(iModule,gLayer,gLadder,gModule);
80     for(Int_t iChip = 0; iChip < fgkNumOfChips; iChip++) {
81       gTitle = "SSD_CM_NSide_Layer"; gTitle += gLayer;
82       gTitle += "_Ladder"; gTitle += gLadder;
83       gTitle += "_Module"; gTitle += gModule;
84       gTitle += "_Chip"; gTitle += iChip+1;
85       gHistSSDCMModule[gHistCounter][iChip] = new TH1F(gTitle.Data(),
86                                                        gTitle.Data(),
87                                                        100,-50.,50.);
88     gHistSSDCMModule[gHistCounter][iChip]->GetXaxis()->SetTitle("CM");
89     list->Add(gHistSSDCMModule[gHistCounter][iChip]);
90     }
91     gHistCounter += 1;
92   }
93
94   return list;
95 }
96
97 //__________________________________________________________//
98 void makeCM(const char* filename, Int_t nEvents, TList *list) {
99   //Function to read the CM values
100   Int_t gStripNumber = 0;
101   Int_t gLayer = 0,gLadder = 0, gModule = 0;
102   
103   //==================================================//
104   AliCDBManager *fCDBManager = AliCDBManager::Instance();
105   fCDBManager->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
106   Int_t runNumber = atoi(gSystem->Getenv("DATE_RUN_NUMBER"));
107   if(!runNumber) 
108     Printf("DATE_RUN_NUMBER not defined!!!\n");
109   
110   fCDBManager->SetRun(runNumber);
111   AliCDBEntry *geomGRP = fCDBManager->Get("GRP/Geometry/Data");
112   if(!geomGRP) cout<<"GRP geometry not found!!!"<<endl;;
113   //==================================================//
114     
115   //==================================================//
116   TH2F *fHistPSideMeanCMMapLayer5 = new TH2F("fHistPSideMeanCMMapLayer5",
117                                          "Layer 5 - P side;N_{module};N_{ladder}",
118                                          22,1,23,
119                                          34,500,534);
120   fHistPSideMeanCMMapLayer5->GetXaxis()->SetTitleColor(1);
121   fHistPSideMeanCMMapLayer5->GetZaxis()->SetRangeUser(0.,20.);
122   fHistPSideMeanCMMapLayer5->SetStats(kFALSE);
123   fHistPSideMeanCMMapLayer5->GetYaxis()->SetTitleOffset(1.8);
124   fHistPSideMeanCMMapLayer5->GetXaxis()->SetNdivisions(22);
125   fHistPSideMeanCMMapLayer5->GetYaxis()->SetNdivisions(34);
126   fHistPSideMeanCMMapLayer5->GetXaxis()->SetLabelSize(0.03);
127   fHistPSideMeanCMMapLayer5->GetYaxis()->SetLabelSize(0.03);
128   fHistPSideMeanCMMapLayer5->GetZaxis()->SetTitleOffset(1.6);
129   fHistPSideMeanCMMapLayer5->GetZaxis()->SetTitle("RMS(CM) (p-side)");
130
131   TH2F *fHistNSideMeanCMMapLayer5 = new TH2F("fHistNSideMeanCMMapLayer5",
132                                          "Layer 5 - N side;N_{module};N_{ladder}",
133                                          22,1,23,
134                                          34,500,534);
135   fHistNSideMeanCMMapLayer5->GetXaxis()->SetTitleColor(1);
136   fHistNSideMeanCMMapLayer5->GetZaxis()->SetRangeUser(0.,20.);
137   fHistNSideMeanCMMapLayer5->SetStats(kFALSE);
138   fHistNSideMeanCMMapLayer5->GetYaxis()->SetTitleOffset(1.8);
139   fHistNSideMeanCMMapLayer5->GetXaxis()->SetNdivisions(22);
140   fHistNSideMeanCMMapLayer5->GetYaxis()->SetNdivisions(34);
141   fHistNSideMeanCMMapLayer5->GetXaxis()->SetLabelSize(0.03);
142   fHistNSideMeanCMMapLayer5->GetYaxis()->SetLabelSize(0.03);
143   fHistNSideMeanCMMapLayer5->GetZaxis()->SetTitleOffset(1.6);
144   fHistNSideMeanCMMapLayer5->GetZaxis()->SetTitle("RMS(CM) (n-side)");
145   
146   TH2F *fHistPSideMeanCMMapLayer6 = new TH2F("fHistPSideMeanCMMapLayer6",
147                                          "Layer 6 - P side;N_{module};N_{ladder}",
148                                          25,1,26,
149                                          38,600,638);
150   fHistPSideMeanCMMapLayer6->GetXaxis()->SetTitleColor(1);
151   fHistPSideMeanCMMapLayer6->GetZaxis()->SetRangeUser(0.,20.);
152   fHistPSideMeanCMMapLayer6->SetStats(kFALSE);
153   fHistPSideMeanCMMapLayer6->GetYaxis()->SetTitleOffset(1.8);
154   fHistPSideMeanCMMapLayer6->GetXaxis()->SetNdivisions(25);
155   fHistPSideMeanCMMapLayer6->GetYaxis()->SetNdivisions(38);
156   fHistPSideMeanCMMapLayer6->GetXaxis()->SetLabelSize(0.03);
157   fHistPSideMeanCMMapLayer6->GetYaxis()->SetLabelSize(0.03);
158   fHistPSideMeanCMMapLayer6->GetZaxis()->SetTitleOffset(1.6);
159   fHistPSideMeanCMMapLayer6->GetZaxis()->SetTitle("RMS(CM) (p-side)");
160
161   TH2F *fHistNSideMeanCMMapLayer6 = new TH2F("fHistNSideMeanCMMapLayer6",
162                                          "Layer 6 - N side;N_{module};N_{ladder}",
163                                          25,1,26,
164                                          38,600,638);
165   fHistNSideMeanCMMapLayer6->GetXaxis()->SetTitleColor(1);
166   fHistNSideMeanCMMapLayer6->GetZaxis()->SetRangeUser(0.,20.);
167   fHistNSideMeanCMMapLayer6->SetStats(kFALSE);
168   fHistNSideMeanCMMapLayer6->GetYaxis()->SetTitleOffset(1.8);
169   fHistNSideMeanCMMapLayer6->GetXaxis()->SetNdivisions(25);
170   fHistNSideMeanCMMapLayer6->GetYaxis()->SetNdivisions(38);
171   fHistNSideMeanCMMapLayer6->GetXaxis()->SetLabelSize(0.03);
172   fHistNSideMeanCMMapLayer6->GetYaxis()->SetLabelSize(0.03);
173   fHistNSideMeanCMMapLayer6->GetZaxis()->SetTitleOffset(1.6);
174   fHistNSideMeanCMMapLayer6->GetZaxis()->SetTitle("RMS(CM) (n-side)");
175
176   //____________________________________________________________//
177   TH2F *fHistPSideRMSCMMapLayer5 = new TH2F("fHistPSideRMSCMMapLayer5",
178                                          "Layer 5 - P side;N_{module};N_{ladder}",
179                                          22,1,23,
180                                          34,500,534);
181   fHistPSideRMSCMMapLayer5->GetXaxis()->SetTitleColor(1);
182   fHistPSideRMSCMMapLayer5->GetZaxis()->SetRangeUser(0.,20.);
183   fHistPSideRMSCMMapLayer5->SetStats(kFALSE);
184   fHistPSideRMSCMMapLayer5->GetYaxis()->SetTitleOffset(1.8);
185   fHistPSideRMSCMMapLayer5->GetXaxis()->SetNdivisions(22);
186   fHistPSideRMSCMMapLayer5->GetYaxis()->SetNdivisions(34);
187   fHistPSideRMSCMMapLayer5->GetXaxis()->SetLabelSize(0.03);
188   fHistPSideRMSCMMapLayer5->GetYaxis()->SetLabelSize(0.03);
189   fHistPSideRMSCMMapLayer5->GetZaxis()->SetTitleOffset(1.6);
190   fHistPSideRMSCMMapLayer5->GetZaxis()->SetTitle("RMS(CM) (p-side)");
191
192   TH2F *fHistNSideRMSCMMapLayer5 = new TH2F("fHistNSideRMSCMMapLayer5",
193                                          "Layer 5 - N side;N_{module};N_{ladder}",
194                                          22,1,23,
195                                          34,500,534);
196   fHistNSideRMSCMMapLayer5->GetXaxis()->SetTitleColor(1);
197   fHistNSideRMSCMMapLayer5->GetZaxis()->SetRangeUser(0.,20.);
198   fHistNSideRMSCMMapLayer5->SetStats(kFALSE);
199   fHistNSideRMSCMMapLayer5->GetYaxis()->SetTitleOffset(1.8);
200   fHistNSideRMSCMMapLayer5->GetXaxis()->SetNdivisions(22);
201   fHistNSideRMSCMMapLayer5->GetYaxis()->SetNdivisions(34);
202   fHistNSideRMSCMMapLayer5->GetXaxis()->SetLabelSize(0.03);
203   fHistNSideRMSCMMapLayer5->GetYaxis()->SetLabelSize(0.03);
204   fHistNSideRMSCMMapLayer5->GetZaxis()->SetTitleOffset(1.6);
205   fHistNSideRMSCMMapLayer5->GetZaxis()->SetTitle("RMS(CM) (n-side)");
206   
207   TH2F *fHistPSideRMSCMMapLayer6 = new TH2F("fHistPSideRMSCMMapLayer6",
208                                          "Layer 6 - P side;N_{module};N_{ladder}",
209                                          25,1,26,
210                                          38,600,638);
211   fHistPSideRMSCMMapLayer6->GetXaxis()->SetTitleColor(1);
212   fHistPSideRMSCMMapLayer6->GetZaxis()->SetRangeUser(0.,20.);
213   fHistPSideRMSCMMapLayer6->SetStats(kFALSE);
214   fHistPSideRMSCMMapLayer6->GetYaxis()->SetTitleOffset(1.8);
215   fHistPSideRMSCMMapLayer6->GetXaxis()->SetNdivisions(25);
216   fHistPSideRMSCMMapLayer6->GetYaxis()->SetNdivisions(38);
217   fHistPSideRMSCMMapLayer6->GetXaxis()->SetLabelSize(0.03);
218   fHistPSideRMSCMMapLayer6->GetYaxis()->SetLabelSize(0.03);
219   fHistPSideRMSCMMapLayer6->GetZaxis()->SetTitleOffset(1.6);
220   fHistPSideRMSCMMapLayer6->GetZaxis()->SetTitle("RMS(CM) (p-side)");
221
222   TH2F *fHistNSideRMSCMMapLayer6 = new TH2F("fHistNSideRMSCMMapLayer6",
223                                          "Layer 6 - N side;N_{module};N_{ladder}",
224                                          25,1,26,
225                                          38,600,638);
226   fHistNSideRMSCMMapLayer6->GetXaxis()->SetTitleColor(1);
227   fHistNSideRMSCMMapLayer6->GetZaxis()->SetRangeUser(0.,20.);
228   fHistNSideRMSCMMapLayer6->SetStats(kFALSE);
229   fHistNSideRMSCMMapLayer6->GetYaxis()->SetTitleOffset(1.8);
230   fHistNSideRMSCMMapLayer6->GetXaxis()->SetNdivisions(25);
231   fHistNSideRMSCMMapLayer6->GetYaxis()->SetNdivisions(38);
232   fHistNSideRMSCMMapLayer6->GetXaxis()->SetLabelSize(0.03);
233   fHistNSideRMSCMMapLayer6->GetYaxis()->SetLabelSize(0.03);
234   fHistNSideRMSCMMapLayer6->GetZaxis()->SetTitleOffset(1.6);
235   fHistNSideRMSCMMapLayer6->GetZaxis()->SetTitle("RMS(CM) (n-side)");
236   //==================================================//
237
238   TChain *chain = new TChain("RAW");
239   chain->Add(filename);
240   Int_t nTotalEvents = chain->GetEntries();
241   if(nEvents == -1) nEvents = nTotalEvents;
242   
243   AliRawReaderRoot *rawReader = new AliRawReaderRoot(filename);
244   Int_t iEvent = 0;
245   Int_t fSSDEvent = 0;
246   for(iEvent = 0; iEvent < nEvents; iEvent++) {
247     cout<<"Event: "<<iEvent+1<<"/"<<nEvents<<endl;
248     rawReader->Select("ITSSSD",-1,-1);  
249     rawReader->Reset(); 
250     rawReader->NextEvent();   
251     
252     if(rawReader->GetType() != 7) continue;
253     
254     fSSDEvent += 1;
255     AliITSRawStreamSSD gSSDStream(rawReader);    
256     while (gSSDStream.Next()) {
257       if(gSSDStream.GetModuleID() < 0) continue;
258       AliITSgeomTGeo::GetModuleId(gSSDStream.GetModuleID(),gLayer,gLadder,gModule);
259       //if(gSSDStream.GetModuleID() != 500) continue;
260       //Printf("Module id: %d - Layer: %d - Ladder: %d - Module: %d",gSSDStream.GetModuleID(),gLayer,gLadder,gModule);
261
262       if(gSSDStream.GetStrip() >= 0) continue;
263       gStripNumber = (gSSDStream.GetSideFlag() == 0) ? gSSDStream.GetStrip() : -gSSDStream.GetStrip() + 2*fgkNumberOfPSideStrips;
264       //Printf("Module id: %d - Strip: %d - strip number: %d - Signal: %lf",gSSDStream.GetModuleID(),gSSDStream.GetStrip(),gStripNumber,signal);
265       if(TMath::Abs(gSSDStream.GetStrip()) < 7)
266         ((TH1*)list->At((gSSDStream.GetModuleID()-500)*fgkNumOfChips+gSSDStream.GetStrip()))->Fill(gSSDStream.GetSignal());
267       if(TMath::Abs(gSSDStream.GetStrip()) > 6)
268         ((TH1*)list->At(fgkSSDMODULES*fgkNumOfChips+(gSSDStream.GetModuleID()-500)*fgkNumOfChips+gSSDStream.GetStrip()-fgkNumOfChips))->Fill(gSSDStream.GetSignal());
269     }//streamer loop
270   }//event loop
271   
272   //compute the rms of the CM values
273   TH1F *gHistCMDummy = new TH1F("gHistCMDummy","",100,-50.,50.);
274   Double_t meanPsideCM = 0.0, meanNsideCM = 0.0;
275   Double_t rmsPsideCM = 0.0, rmsNsideCM = 0.0;
276   for(Int_t iModule = 0; iModule < fgkSSDMODULES; iModule++) {
277       meanPsideCM = 0.0; meanNsideCM = 0.0;
278     rmsPsideCM = 0.0; rmsNsideCM = 0.0;
279     AliITSgeomTGeo::GetModuleId(iModule+500,gLayer,gLadder,gModule);
280   
281     gHistCMDummy->Clear();
282     for(Int_t iChip = 0; iChip < fgkNumOfChips; iChip++) {
283       //cout<<"Name: "<<dynamic_cast<TH1*>(list->At(iModule*fgkNumOfChips+iChip))->GetName()<<endl;
284       gHistCMDummy->Add((TH1*)list->At(iModule*fgkNumOfChips+iChip));
285     } 
286     meanPsideCM = TMath::Abs(gHistCMDummy->GetMean());
287     rmsPsideCM = gHistCMDummy->GetRMS();
288
289     gHistCMDummy->Clear();
290     for(Int_t iChip = 0; iChip < fgkNumOfChips; iChip++) {
291       //cout<<"Name: "<<dynamic_cast<TH1*>(list->At(fgkSSDMODULES*fgkNumOfChips+iModule*fgkNumOfChips+iChip))->GetName()<<endl;
292       gHistCMDummy->Add((TH1*)list->At(fgkSSDMODULES*fgkNumOfChips+iModule*fgkNumOfChips+iChip));
293     }
294     meanNsideCM = TMath::Abs(gHistCMDummy->GetMean());
295     rmsNsideCM = gHistCMDummy->GetRMS();
296
297     if(meanPsideCM == 0) meanPsideCM = 0.001;
298     if(meanNsideCM == 0) meanNsideCM = 0.001;
299     if(rmsPsideCM == 0) rmsPsideCM = 0.001;
300     if(rmsNsideCM == 0) rmsNsideCM = 0.001;
301  
302     if(gLayer == 5) {
303       fHistPSideMeanCMMapLayer5->SetBinContent(gModule,gLadder,meanPsideCM);
304       fHistNSideMeanCMMapLayer5->SetBinContent(gModule,gLadder,meanNsideCM);
305       fHistPSideRMSCMMapLayer5->SetBinContent(gModule,gLadder,rmsPsideCM);
306       fHistNSideRMSCMMapLayer5->SetBinContent(gModule,gLadder,rmsNsideCM);
307     }
308      if(gLayer == 6) {
309       fHistPSideMeanCMMapLayer6->SetBinContent(gModule,gLadder,meanPsideCM);
310       fHistNSideMeanCMMapLayer6->SetBinContent(gModule,gLadder,meanNsideCM);
311       fHistPSideRMSCMMapLayer6->SetBinContent(gModule,gLadder,rmsPsideCM);
312       fHistNSideRMSCMMapLayer6->SetBinContent(gModule,gLadder,rmsNsideCM);
313     }
314     }
315
316   TFile *foutput = TFile::Open("SSD.CM.root","recreate");
317   list->Write();
318   fHistPSideMeanCMMapLayer5->Write(); fHistNSideMeanCMMapLayer5->Write();
319   fHistPSideMeanCMMapLayer6->Write(); fHistNSideMeanCMMapLayer6->Write();
320   fHistPSideRMSCMMapLayer5->Write(); fHistNSideRMSCMMapLayer5->Write();
321   fHistPSideRMSCMMapLayer6->Write(); fHistNSideRMSCMMapLayer6->Write();
322   foutput->Close();
323 }
324   
325 //__________________________________________________________//
326 void drawSSDCM(const char* filename = "SSD.CM.root") {
327   gStyle->SetPalette(1,0);
328
329   TFile *f = TFile::Open(filename);  
330
331   TCanvas *c1 = new TCanvas("c1","Mean of CM values",0,0,700,650);
332   c1->SetFillColor(10); c1->SetHighLightColor(10); c1->Divide(2,2);
333   c1->cd(1);
334   TH2F *fHistPSideMeanCMMapLayer5 = dynamic_cast<TH2F *>(f->Get("fHistPSideMeanCMMapLayer5"));
335   fHistPSideMeanCMMapLayer5->Draw("colz");
336   c1->cd(2);
337   TH2F *fHistNSideMeanCMMapLayer5 = dynamic_cast<TH2F *>(f->Get("fHistNSideMeanCMMapLayer5"));
338   fHistNSideMeanCMMapLayer5->Draw("colz");
339   c1->cd(3);
340   TH2F *fHistPSideMeanCMMapLayer6 = dynamic_cast<TH2F *>(f->Get("fHistPSideMeanCMMapLayer6"));
341   fHistPSideMeanCMMapLayer6->Draw("colz");
342   c1->cd(4);
343   TH2F *fHistNSideMeanCMMapLayer6 = dynamic_cast<TH2F *>(f->Get("fHistNSideMeanCMMapLayer6"));
344   fHistNSideMeanCMMapLayer6->Draw("colz");
345
346   TCanvas *c2 = new TCanvas("c2","RMS of CM values",100,100,700,650);
347   c2->SetFillColor(10); c2->SetHighLightColor(10); c2->Divide(2,2);
348   c2->cd(1);
349   TH2F *fHistPSideRMSCMMapLayer5 = dynamic_cast<TH2F *>(f->Get("fHistPSideRMSCMMapLayer5"));
350   fHistPSideRMSCMMapLayer5->Draw("colz");
351   c2->cd(2);
352   TH2F *fHistNSideRMSCMMapLayer5 = dynamic_cast<TH2F *>(f->Get("fHistNSideRMSCMMapLayer5"));
353   fHistNSideRMSCMMapLayer5->Draw("colz");
354   c2->cd(3);
355   TH2F *fHistPSideRMSCMMapLayer6 = dynamic_cast<TH2F *>(f->Get("fHistPSideRMSCMMapLayer6"));
356   fHistPSideRMSCMMapLayer6->Draw("colz");
357   c2->cd(4);
358   TH2F *fHistNSideRMSCMMapLayer6 = dynamic_cast<TH2F *>(f->Get("fHistNSideRMSCMMapLayer6"));
359   fHistNSideRMSCMMapLayer6->Draw("colz");
360 }