Added two missing includes to allow macro compilation (thanks to Laurent for remarkin...
[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 fgkNumOfLDCs = 8;      //number of SSD LDCs
26 static const Int_t fgkNumOfDDLs = 16;      //number of SSD DDLs
27 static const Int_t fgkSSDMODULES = 1698;      //total number of SSD modules
28 static const Int_t fgkSSDLADDERSLAYER5 = 34; //ladders on layer 5
29 static const Int_t fgkSSDLADDERSLAYER6 = 38; //ladders on layer 6
30 static const Int_t fgkSSDMODULESPERLADDERLAYER5 = 22; //modules per ladder - layer 5                           
31 static const Int_t fgkSSDMODULESPERLADDERLAYER6 = 25; //modules per ladder - layer 6                          
32 static const Int_t fgkSSDMODULESLAYER5 = 748; //total number of SSD modules - layer5                           
33 static const Int_t fgkSSDMODULESLAYER6 = 950; //total number of SSD modules - layer6                         
34 static const Int_t fgkNumberOfPSideStrips = 768; //number of P-side strips
35 //============================================================//
36
37 TList *initCM();
38 void makeCM(const char* filename, Int_t nEvents, TList *list);
39
40 //__________________________________________________________//
41 void readSSDCommonMode(const char* filename = "raw.root",
42                        Int_t nEvents = -1) {
43   //Reads the CM pseudo-channels and produces the CM map for both layers 
44   //and for p and n-side.
45   gStyle->SetPalette(1,0);
46
47   TList *list = initCM();
48   //list->ls();
49   Printf("CM histograms: %d",list->GetEntries());
50   makeCM(filename,nEvents,list);
51
52  
53 //__________________________________________________________//
54 TList *initCM() {
55   //Initializes the histograms and returns the TList object
56   TList *list = new TList();
57
58   Int_t gLayer = 0,gLadder = 0, gModule = 0;
59   Int_t gHistCounter = 0;
60   TString gTitle;
61   TH1F *gHistSSDCMModule[2*fgkSSDMODULES]; 
62   for(Int_t iModule = 500; iModule < fgkSSDMODULES + 500; iModule++) {
63     AliITSgeomTGeo::GetModuleId(iModule,gLayer,gLadder,gModule);
64     gTitle = "SSD_CM_PSide_Layer"; gTitle += gLayer;
65     gTitle += "_Ladder"; gTitle += gLadder;
66     gTitle += "_Module"; gTitle += gModule;
67     gHistSSDCMModule[gHistCounter] = new TH1F(gTitle.Data(),gTitle.Data(),
68                                               100,-50.,50.);
69     gHistSSDCMModule[gHistCounter]->GetXaxis()->SetTitle("CM");
70     list->Add(gHistSSDCMModule[gHistCounter]);
71     gHistCounter += 1;
72   }
73   for(Int_t iModule = 500; iModule < fgkSSDMODULES + 500; iModule++) {
74     AliITSgeomTGeo::GetModuleId(iModule,gLayer,gLadder,gModule);
75     gTitle = "SSD_CM_NSide_Layer"; gTitle += gLayer;
76     gTitle += "_Ladder"; gTitle += gLadder;
77     gTitle += "_Module"; gTitle += gModule;
78     gHistSSDCMModule[gHistCounter] = new TH1F(gTitle.Data(),gTitle.Data(),
79                                               100,-50.,50.);
80     gHistSSDCMModule[gHistCounter]->GetXaxis()->SetTitle("CM");
81     list->Add(gHistSSDCMModule[gHistCounter]);
82     gHistCounter += 1;
83   }
84
85   return list;
86 }
87
88 //__________________________________________________________//
89 void makeCM(const char* filename, Int_t nEvents, TList *list) {
90   //Function to read the CM values
91   Int_t gStripNumber = 0;
92   Int_t gLayer = 0,gLadder = 0, gModule = 0;
93   
94   //==================================================//
95   AliCDBManager *fCDBManager = AliCDBManager::Instance();
96   fCDBManager->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
97   Int_t runNumber = atoi(gSystem->Getenv("DATE_RUN_NUMBER"));
98   if(!runNumber) 
99     Printf("DATE_RUN_NUMBER not defined!!!\n");
100   
101   fCDBManager->SetRun(runNumber);
102   AliCDBEntry *geomGRP = fCDBManager->Get("GRP/Geometry/Data");
103   if(!geomGRP) cout<<"GRP geometry not found!!!"<<endl;;
104   //==================================================//
105     
106   //==================================================//
107   TH2F *fHistPSideCMMapLayer5 = new TH2F("fHistPSideCMMapLayer5",
108                                          "Layer 5 - P side;N_{module};N_{ladder}",
109                                          22,1,23,
110                                          34,500,534);
111   fHistPSideCMMapLayer5->GetXaxis()->SetTitleColor(1);
112   fHistPSideCMMapLayer5->GetZaxis()->SetRangeUser(0.,20.);
113   fHistPSideCMMapLayer5->SetStats(kFALSE);
114   fHistPSideCMMapLayer5->GetYaxis()->SetTitleOffset(1.8);
115   fHistPSideCMMapLayer5->GetXaxis()->SetNdivisions(22);
116   fHistPSideCMMapLayer5->GetYaxis()->SetNdivisions(34);
117   fHistPSideCMMapLayer5->GetXaxis()->SetLabelSize(0.03);
118   fHistPSideCMMapLayer5->GetYaxis()->SetLabelSize(0.03);
119   fHistPSideCMMapLayer5->GetZaxis()->SetTitleOffset(1.6);
120   fHistPSideCMMapLayer5->GetZaxis()->SetTitle("RMS(CM) (p-side)");
121
122   TH2F *fHistNSideCMMapLayer5 = new TH2F("fHistNSideCMMapLayer5",
123                                          "Layer 5 - N side;N_{module};N_{ladder}",
124                                          22,1,23,
125                                          34,500,534);
126   fHistNSideCMMapLayer5->GetXaxis()->SetTitleColor(1);
127   fHistNSideCMMapLayer5->GetZaxis()->SetRangeUser(0.,20.);
128   fHistNSideCMMapLayer5->SetStats(kFALSE);
129   fHistNSideCMMapLayer5->GetYaxis()->SetTitleOffset(1.8);
130   fHistNSideCMMapLayer5->GetXaxis()->SetNdivisions(22);
131   fHistNSideCMMapLayer5->GetYaxis()->SetNdivisions(34);
132   fHistNSideCMMapLayer5->GetXaxis()->SetLabelSize(0.03);
133   fHistNSideCMMapLayer5->GetYaxis()->SetLabelSize(0.03);
134   fHistNSideCMMapLayer5->GetZaxis()->SetTitleOffset(1.6);
135   fHistNSideCMMapLayer5->GetZaxis()->SetTitle("RMS(CM) (n-side)");
136   
137   TH2F *fHistPSideCMMapLayer6 = new TH2F("fHistPSideCMMapLayer6",
138                                          "Layer 6 - P side;N_{module};N_{ladder}",
139                                          25,1,26,
140                                          38,600,638);
141   fHistPSideCMMapLayer6->GetXaxis()->SetTitleColor(1);
142   fHistPSideCMMapLayer6->GetZaxis()->SetRangeUser(0.,20.);
143   fHistPSideCMMapLayer6->SetStats(kFALSE);
144   fHistPSideCMMapLayer6->GetYaxis()->SetTitleOffset(1.8);
145   fHistPSideCMMapLayer6->GetXaxis()->SetNdivisions(25);
146   fHistPSideCMMapLayer6->GetYaxis()->SetNdivisions(38);
147   fHistPSideCMMapLayer6->GetXaxis()->SetLabelSize(0.03);
148   fHistPSideCMMapLayer6->GetYaxis()->SetLabelSize(0.03);
149   fHistPSideCMMapLayer6->GetZaxis()->SetTitleOffset(1.6);
150   fHistPSideCMMapLayer6->GetZaxis()->SetTitle("RMS(CM) (p-side)");
151
152   TH2F *fHistNSideCMMapLayer6 = new TH2F("fHistNSideCMMapLayer6",
153                                          "Layer 6 - N side;N_{module};N_{ladder}",
154                                          25,1,26,
155                                          38,600,638);
156   fHistNSideCMMapLayer6->GetXaxis()->SetTitleColor(1);
157   fHistNSideCMMapLayer6->GetZaxis()->SetRangeUser(0.,20.);
158   fHistNSideCMMapLayer6->SetStats(kFALSE);
159   fHistNSideCMMapLayer6->GetYaxis()->SetTitleOffset(1.8);
160   fHistNSideCMMapLayer6->GetXaxis()->SetNdivisions(25);
161   fHistNSideCMMapLayer6->GetYaxis()->SetNdivisions(38);
162   fHistNSideCMMapLayer6->GetXaxis()->SetLabelSize(0.03);
163   fHistNSideCMMapLayer6->GetYaxis()->SetLabelSize(0.03);
164   fHistNSideCMMapLayer6->GetZaxis()->SetTitleOffset(1.6);
165   fHistNSideCMMapLayer6->GetZaxis()->SetTitle("RMS(CM) (n-side)");
166   //==================================================//
167
168   TChain *chain = new TChain("RAW");
169   chain->Add(filename);
170   Int_t nTotalEvents = chain->GetEntries();
171   if(nEvents == -1) nEvents = nTotalEvents;
172   
173   AliRawReaderRoot *rawReader = new AliRawReaderRoot(filename);
174   Int_t iEvent = 0;
175   Int_t fSSDEvent = 0;
176   for(iEvent = 0; iEvent < nEvents; iEvent++) {
177     cout<<"Event: "<<iEvent+1<<"/"<<nEvents<<endl;
178     rawReader->Select("ITSSSD",-1,-1);  
179     rawReader->Reset(); 
180     rawReader->NextEvent();   
181     
182     if(rawReader->GetType() != 7) continue;
183     
184     fSSDEvent += 1;
185     AliITSRawStreamSSD gSSDStream(rawReader);    
186     while (gSSDStream.Next()) {
187       if(gSSDStream.GetModuleID() < 0) continue;
188       AliITSgeomTGeo::GetModuleId(gSSDStream.GetModuleID(),gLayer,gLadder,gModule);
189       //if(gSSDStream.GetModuleID() != 500) continue;
190       //Printf("Module id: %d - Layer: %d - Ladder: %d - Module: %d",gSSDStream.GetModuleID(),gLayer,gLadder,gModule);
191
192       if(gSSDStream.GetStrip() >= 0) continue;
193       gStripNumber = (gSSDStream.GetSideFlag() == 0) ? gSSDStream.GetStrip() : -gSSDStream.GetStrip() + 2*fgkNumberOfPSideStrips;
194       //Printf("Module id: %d - Strip: %d - strip number: %d - Signal: %lf",gSSDStream.GetModuleID(),gSSDStream.GetStrip(),gStripNumber,signal);
195       if(TMath::Abs(gSSDStream.GetStrip()) < 7)
196         ((TH1*)list->At(gSSDStream.GetModuleID()-500))->Fill(gSSDStream.GetSignal());
197       if(TMath::Abs(gSSDStream.GetStrip()) > 6)
198         ((TH1*)list->At(1698+gSSDStream.GetModuleID()-500))->Fill(gSSDStream.GetSignal());
199       
200     }//streamer loop
201   }//event loop
202   
203   //compute the rms of the CM values
204   Double_t rmsPsideCM = 0.0, rmsNsideCM = 0.0;
205   for(Int_t i = 0; i < 1698; i++) {
206     rmsPsideCM = 0.0; rmsNsideCM = 0.0;
207     AliITSgeomTGeo::GetModuleId(i+500,gLayer,gLadder,gModule);
208     //Printf("%s - %s",((TH1*)list->At(i))->GetName(),
209     //((TH1*)list->At(1698+i))->GetName());
210     rmsPsideCM = ((TH1*)list->At(i))->GetRMS();
211     rmsNsideCM = ((TH1*)list->At(1698+i))->GetRMS();
212     if(rmsPsideCM == 0) rmsPsideCM = 0.001;
213     if(rmsNsideCM == 0) rmsNsideCM = 0.001;
214     //Printf("rmsPside: %lf - rmsNside: %lf",rmsPsideCM,rmsNsideCM);
215     if(gLayer == 5) {
216       fHistPSideCMMapLayer5->SetBinContent(gModule,gLadder,rmsPsideCM);
217       fHistNSideCMMapLayer5->SetBinContent(gModule,gLadder,rmsNsideCM);
218     }
219      if(gLayer == 6) {
220       fHistPSideCMMapLayer6->SetBinContent(gModule,gLadder,rmsPsideCM);
221       fHistNSideCMMapLayer6->SetBinContent(gModule,gLadder,rmsNsideCM);
222     }
223     }
224
225   TFile *foutput = TFile::Open("SSD.CM.root","recreate");
226   list->Write();
227   fHistPSideCMMapLayer5->Write(); fHistNSideCMMapLayer5->Write();
228   fHistPSideCMMapLayer6->Write(); fHistNSideCMMapLayer6->Write();
229   foutput->Close();
230 }
231   
232 //__________________________________________________________//
233 void drawSSDCM(const char* filename = "SSD.CM.root") {
234   gStyle->SetPalette(1,0);
235
236   TFile *f = TFile::Open(filename);  
237
238   TCanvas *c1 = new TCanvas("c1","CM values",0,0,700,650);
239   c1->SetFillColor(10); c1->SetHighLightColor(10); c1->Divide(2,2);
240   c1->cd(1);
241   TH1F *fHistPSideCMMapLayer5 = dynamic_cast<TH1F *>(f->Get("fHistPSideCMMapLayer5"));
242   fHistPSideCMMapLayer5->Draw("colz");
243   c1->cd(2);
244   TH1F *fHistNSideCMMapLayer5 = dynamic_cast<TH1F *>(f->Get("fHistNSideCMMapLayer5"));
245   fHistNSideCMMapLayer5->Draw("colz");
246   c1->cd(3);
247   TH1F *fHistPSideCMMapLayer6 = dynamic_cast<TH1F *>(f->Get("fHistPSideCMMapLayer6"));
248   fHistPSideCMMapLayer6->Draw("colz");
249   c1->cd(4);
250   TH1F *fHistNSideCMMapLayer6 = dynamic_cast<TH1F *>(f->Get("fHistNSideCMMapLayer6"));
251   fHistNSideCMMapLayer6->Draw("colz");
252 }