AliPhysics  59e0e03 (59e0e03)
anaTree.C
Go to the documentation of this file.
1 #if !defined(__CINT__) || defined(__MAKECINT__)
2 #include <AliEMCALGeometry.h>
3 #include <AliOADBContainer.h>
4 
5 #include <TCanvas.h>
6 #include <TClonesArray.h>
7 #include <TDatime.h>
8 #include <TFile.h>
9 #include <TGraphErrors.h>
10 #include <TGrid.h>
11 #include <TH2F.h>
12 #include <TLegend.h>
13 #include <TMap.h>
14 #include <TNtuple.h>
15 #include <TObject.h>
16 #include <TProfile.h>
17 #include <TProfile2D.h>
18 #include <TSystem.h>
19 #include <TString.h>
20 #include "createTree.C"
21 #endif
22 
23 Bool_t writeDetailed = kFALSE; // switch on detailed information storing
24 Int_t debugInfo = 1; // enable different levels of debuggin information
25 Bool_t useWeight = kTRUE;
26 Int_t nBinsT = 1000;
27 
28 void anaTree(
29  const char *ifile ="treefile.root",
30  const char *ofile ="outhist.root",
31  Bool_t applyDurLimit = kFALSE, // bolean to switch on duration limit
32  Float_t durMin = 0, // minimum length of a run in minutes
33  Float_t durMax = 10000, // maximum length of a run in minutes
34  Bool_t appBC = kFALSE, // boolean to switch on bad channel
35  Int_t referenceRun = -1 // define reference run to which all are being calibrated
36  )
37 {
38  // Load EMCAL geometry for run 2
39  AliEMCALGeometry *g=AliEMCALGeometry::GetInstance("EMCAL_COMPLETE12SMV1_DCAL_8SM");
40  const Int_t kSM=g->GetNumberOfSuperModules();
41  const Int_t kNcells=g->GetNCells();
42  const Int_t gain = 1;
43 
44  // Return info about optional settings
45  if (applyDurLimit){
46  cout << "INFO: only runs with a length of " << durMin << " to " << durMax << " minutes will be considered in the analysis" << endl;
47  }
48  if (appBC){
49  cout << "INFO: will be using bad channel map" << endl;
50  }
51  if (referenceRun != -1){
52  cout << "INFO: will be storing ratios from run: " << referenceRun << " as reference numbers" << endl;
53  }
54 
55  // initialize info from tree created by $ALICE_PHYSICS/PWGPP/EMCAL/TeCMacros/createTree.C
56  TCalInfo *info = 0;
57  TFile *in = TFile::Open(ifile,"read");
58  TTree *tt = (TTree*)in->Get("tcal");
59  tt->SetBranchAddress("event",&info);
60  tt->Branch("event", &info, 32000, 99);
61  Int_t Nev=tt->GetEntries();
62 
63  Int_t minRunNo = 1e6;
64  Int_t maxRunNo = -1;
65  for (Int_t i=0;i<Nev;++i) {
66  tt->GetEvent(i);
67  if (info->fRunNo < minRunNo) minRunNo = info->fRunNo;
68  if (info->fRunNo > maxRunNo) maxRunNo = info->fRunNo;
69  }
70 
71  TProfile *gLedVsT[20];
72  TProfile *gLedMonVsT[20];
73  TH2F* hLedVsLength[20];
74  TH2F* hLedMonVsLength[20];
75  TProfile *gRatVsT[20];
76  TH1D* hRefRunCellIdVsRat = new TH1D("", "Led/LedMon Ref run; Cell ID", kNcells+1, -0.5, kNcells+1-0.5);
77  hRefRunCellIdVsRat->SetName("ReferenceRunRatios");
78  TH1D* hRefSMVsT = new TH1D("", "T Ref run; SM", 20, -0.5, 20-0.5);
79  hRefSMVsT->SetName("ReferenceRunTemperatures");
80  TH1F* hAverageT[20];
81  TH1F* hAverageTSorted[20];
82 
83  const char* opt = ""; //"S" for spread
84  for (Int_t i=0;i<20;++i) {
85  gLedVsT[i] = new TProfile("","Led info;T;",nBinsT,15,40,opt);
86  gLedVsT[i]->SetName(Form("ledsm%d",i));
87  gLedMonVsT[i] = new TProfile("","Led info;T;",nBinsT,15,40,opt);
88  gLedMonVsT[i]->SetName(Form("ledmonsm%d",i));
89  gRatVsT[i] = new TProfile("","Led/LedMon;T",nBinsT,15,40,opt);
90  gRatVsT[i]->SetName(Form("ledovermonsm%d",i));
91  hLedVsLength[i] = new TH2F ("","Led info; t [h];",500,0,20, 20000, 0, 20000);
92  hLedVsLength[i]->SetName(Form("ledVsLength%d",i));
93  hLedMonVsLength[i] = new TH2F ("","Led Mon info; t [h];",500,0,20, 20000, 0, 40000);
94  hLedMonVsLength[i]->SetName(Form("ledMonVsLength%d",i));
95  hAverageT[i] = new TH1F ("",Form("T SM %d ; run ID; T",i),Nev,0.5,Nev+0.5);
96  hAverageT[i]->SetName(Form("TAverageSM%dvsRunId",i));
97  hAverageTSorted[i] = new TH1F ("",Form("T SM %d ; run number; T",i),maxRunNo-minRunNo+1,minRunNo-0.5,maxRunNo+0.5);
98  hAverageTSorted[i]->SetName(Form("TAverageSM%dvsRunNumber",i));
99  }
100 
101  TH1F* hSensorsT[160];
102  TH1F* hSensorsTSorted[160];
103  for (Int_t sens = 0; sens< 160; sens++){
104  hSensorsT[sens] = new TH1F ("",Form("T sensor %d ; run ID; T",sens),Nev,0.5,Nev+0.5);
105  hSensorsT[sens]->SetName(Form("Tsensor%dvsRunId",sens));
106  hSensorsTSorted[sens] = new TH1F ("",Form("T sensor %d ; run number; T",sens),maxRunNo-minRunNo+1,minRunNo-0.5,maxRunNo+0.5);
107  hSensorsTSorted[sens]->SetName(Form("Tsensor%dvsRunNumber",sens));
108  }
109 
110  TProfile *gLedCellVsT[kNcells+1];
111  TProfile *gLedCellRMSDiffMeanVsT[kNcells+1];
112  TProfile *gLedMonCellVsT[kNcells+1];
113  TProfile *gLedMonCellRMSDiffMeanVsT[kNcells+1];
114  TProfile *gRatCellVsT[kNcells+1];
115  TProfile *gRatECellVsT[kNcells+1];
116 
117  cout << "Initializing cell histos" << endl;
118  for (Int_t j=0;j<kNcells+1;++j) {
119  if (j%500 == 0) cout << "-->next 500: " << j << endl;
120  if (writeDetailed){
121  gLedCellVsT[j] = new TProfile("",Form("Led info cell ID%i ;T;",j),nBinsT,15,40,opt);
122  gLedCellVsT[j]->SetName(Form("ledCell%d",j));
123  gLedMonCellVsT[j] = new TProfile("",Form("LedMon info cell ID%i ;T;",j),nBinsT,15,40,opt);
124  gLedMonCellVsT[j]->SetName(Form("ledMonCell%d",j));
125  gLedCellRMSDiffMeanVsT[j] = new TProfile("",Form("Led rms/mean info cell ID%i ;T;",j),nBinsT,15,40,opt);
126  gLedCellRMSDiffMeanVsT[j]->SetName(Form("ledCellRMSDiffMean%d",j));
127  gLedMonCellRMSDiffMeanVsT[j] = new TProfile("",Form("LedMon rms/mean info cell ID%i ;T;",j),nBinsT,15,40,opt);
128  gLedMonCellRMSDiffMeanVsT[j]->SetName(Form("ledMonRMSDiffMeanCell%d",j));
129  }
130  gRatCellVsT[j] = new TProfile("",Form("Led/LedMon cell ID%i ;T;",j),nBinsT,15,40,opt);
131  gRatCellVsT[j]->SetName(Form("ledovermonCell%d",j));
132  gRatECellVsT[j] = new TProfile("",Form("Led/LedMon Error cell ID%i ;T;",j),nBinsT,15,40,opt);
133  gRatECellVsT[j]->SetName(Form("ledovermonECell%d",j));
134  }
135  cout << "-> done initializing histos" << endl;
136 
137  Bool_t isRefRun = kFALSE;
138  Bool_t hadRefRun = kFALSE;
139 
140  for (Int_t i=0;i<Nev;++i) {
141  tt->GetEvent(i);
142  Float_t deltaTime = ((Float_t)info->fLastTime-(Float_t)info->fFirstTime)/60.; // run duration in minutes
143  cout << info->fRunNo << "\t" << info->fFirstTime << "\t"<<info->fLastTime << "\t"<< deltaTime<< endl;
144 
145  if (referenceRun != -1 && referenceRun == info->fRunNo){
146  isRefRun = kTRUE;
147  hadRefRun = kTRUE;
148  } else {
149  isRefRun = kFALSE;
150  }
151 
152  if (applyDurLimit){
153  if (deltaTime < durMin || deltaTime > durMax){
154  cout << "INFO: skipped run due to mismatch in run length" << endl;
155  continue;
156  }
157  }
158 
159  TClonesArray &sms = info->fSMs;
160  for (Int_t sm=0;sm<sms.GetEntries();++sm) {
161  TCalSM *smInfot = static_cast<TCalSM*>(sms.At(sm));
162  hAverageT[sm]->SetBinContent(i+1,smInfot->fAvgT);
163  hAverageTSorted[sm]->SetBinContent(hAverageTSorted[sm]->FindBin(info->fRunNo),smInfot->fAvgT);
164  hSensorsT[sm*8]->SetBinContent(i+1,smInfot->fT1);
165  hSensorsT[sm*8+1]->SetBinContent(i+1,smInfot->fT2);
166  hSensorsT[sm*8+2]->SetBinContent(i+1,smInfot->fT3);
167  hSensorsT[sm*8+3]->SetBinContent(i+1,smInfot->fT4);
168  hSensorsT[sm*8+4]->SetBinContent(i+1,smInfot->fT5);
169  hSensorsT[sm*8+5]->SetBinContent(i+1,smInfot->fT6);
170  hSensorsT[sm*8+6]->SetBinContent(i+1,smInfot->fT7);
171  hSensorsT[sm*8+7]->SetBinContent(i+1,smInfot->fT8);
172  hSensorsTSorted[sm*8]->SetBinContent(hSensorsTSorted[sm*8]->FindBin(info->fRunNo),smInfot->fT1);
173  hSensorsTSorted[sm*8+1]->SetBinContent(hSensorsTSorted[sm*8]->FindBin(info->fRunNo),smInfot->fT2);
174  hSensorsTSorted[sm*8+2]->SetBinContent(hSensorsTSorted[sm*8]->FindBin(info->fRunNo),smInfot->fT3);
175  hSensorsTSorted[sm*8+3]->SetBinContent(hSensorsTSorted[sm*8]->FindBin(info->fRunNo),smInfot->fT4);
176  hSensorsTSorted[sm*8+4]->SetBinContent(hSensorsTSorted[sm*8]->FindBin(info->fRunNo),smInfot->fT5);
177  hSensorsTSorted[sm*8+5]->SetBinContent(hSensorsTSorted[sm*8]->FindBin(info->fRunNo),smInfot->fT6);
178  hSensorsTSorted[sm*8+6]->SetBinContent(hSensorsTSorted[sm*8]->FindBin(info->fRunNo),smInfot->fT7);
179  hSensorsTSorted[sm*8+7]->SetBinContent(hSensorsTSorted[sm*8]->FindBin(info->fRunNo),smInfot->fT8);
180  }
181 
182 
183  TClonesArray &cells = info->fCells;
184  for (Int_t j=0;j<cells.GetEntries();++j) {
185  TCalCell *cell = static_cast<TCalCell*>(cells.At(j));
186  Int_t cellID = cell->fId;
187  Int_t sm = cell->fSM;
188  Int_t badcell = cell->fBad;
189  Double_t ledM = cell->fLedM;
190  Double_t ledR = cell->fLedR;
191  Double_t monM = cell->fMonM;
192  Double_t monR = cell->fMonR;
193  Double_t locT = cell->fLocT;
194  Double_t smT = cell->fSMT;
195 
196  if (appBC && badcell > 0){
197  if (debugInfo > 1) cout << "found bad cell for " << cellID << "\t " << info->fRunNo << endl;
198  continue;
199  }
200 
201  Double_t T = smT; // use local or SM T, change here
202  if ((T<5)||(T>45))
203  continue;
204 
205  if (isRefRun){
206  if (hRefSMVsT->GetBinContent(sm+1) < 5) hRefSMVsT->SetBinContent(sm+1,smT);
207  }
208 
209  if ((ledM<=0)||(ledR<=0))
210  continue;
211  if ((monM<=0)||(monR<=0))
212  continue;
213  Double_t w = 1./(ledR*ledR);
214  if (!useWeight) w = 1;
215  gLedVsT[sm]->Fill(T,ledM,w);
216  if (writeDetailed){
217  gLedCellVsT[cellID]->Fill(smT,ledM,w);
218  gLedCellRMSDiffMeanVsT[cellID]->Fill(smT,ledR/ledM,1);
219  }
220  Double_t w3 = 1./(monR*monR);
221  if (!useWeight) w3 = 1;
222  gLedMonVsT[sm]->Fill(T,monM,w3);
223  if (writeDetailed){
224  gLedMonCellVsT[cellID]->Fill(smT,monM,w3);
225  gLedMonCellRMSDiffMeanVsT[cellID]->Fill(smT,monR/monM,1);
226  }
227  Double_t ratErr = ledM/monM * TMath::Sqrt((ledR*ledR)/(ledM*ledM)+(monR*monR)/(monM*monM));
228  Double_t w2=1./(ratErr*ratErr);
229  if (!useWeight) w2 = 1;
230  gRatVsT[sm]->Fill(T,ledM/monM,w2);
231  gRatCellVsT[cellID]->Fill(smT,ledM/monM,w2);
232  gRatECellVsT[cellID]->Fill(smT,ratErr,w2);
233  hLedVsLength[sm]->Fill(deltaTime/60,ledM,w);
234  hLedMonVsLength[sm]->Fill(deltaTime/60,monM,w);
235  if (isRefRun){
236  hRefRunCellIdVsRat->SetBinContent(cellID+1,ledM/monM);
237  hRefRunCellIdVsRat->SetBinError(cellID+1,ratErr);
238  }
239  }
240  }
241 
242  TFile *out = TFile::Open(ofile,"recreate");
243  if (!out)
244  out = TFile::Open("dummyfile.root","update");
245  for (Int_t i=0;i<20;++i) {
246  gLedVsT[i]->Write();
247  gLedMonVsT[i]->Write();
248  gRatVsT[i]->Write();
249  hLedVsLength[i]->Write();
250  hLedMonVsLength[i]->Write();
251  hAverageT[i]->Write();
252  hAverageTSorted[i]->Write();
253  }
254  if (hadRefRun){
255  hRefRunCellIdVsRat->Write();
256  hRefSMVsT->Write();
257  }
258  for (Int_t sens = 0; sens < 160; sens++){
259  hSensorsT[sens]->Write();
260  hSensorsTSorted[sens]->Write();
261  }
262 
263  Int_t smcurr = -1;
264  Int_t cellIDFirstInSM = 0;
265  for (Int_t j=0;j<kNcells+1;++j) {
266  Int_t iTower = -1, iIphi = -1, iIeta = -1, sm=-1;
267  g->GetCellIndex(j,sm,iTower,iIphi,iIeta);
268  if (smcurr != sm){
269  smcurr = sm;
270  out->mkdir(Form("cellsInSM%d",sm));
271  out->cd(Form("cellsInSM%d",sm));
272  if (debugInfo) cout << "SM: \t" << sm-1 << "\t"<< cellIDFirstInSM << "\t" << j << "\t" << j -cellIDFirstInSM<< endl;
273  cellIDFirstInSM = j;
274  }
275 
276  if (writeDetailed) {
277  if (gLedCellVsT[j]->GetEntries() > 0)
278  gLedCellVsT[j]->Write();
279  if (gLedMonCellVsT[j]->GetEntries() > 0)
280  gLedMonCellVsT[j]->Write();
281  if (gLedCellRMSDiffMeanVsT[j]->GetEntries() > 0)
282  gLedCellRMSDiffMeanVsT[j]->Write();
283  if (gLedMonCellRMSDiffMeanVsT[j]->GetEntries() > 0)
284  gLedMonCellRMSDiffMeanVsT[j]->Write();
285  }
286  if (gRatCellVsT[j]->GetEntries() > 0)
287  gRatCellVsT[j]->Write();
288  if (gRatECellVsT[j]->GetEntries() > 0)
289  gRatECellVsT[j]->Write();
290  }
291 }
Double32_t fT2
Definition: createTree.C:46
double Double_t
Definition: External.C:58
Double32_t fSMT
Definition: createTree.C:35
Short_t fBad
Definition: createTree.C:29
void anaTree(const char *ifile="treefile.root", const char *ofile="outhist.root", Bool_t applyDurLimit=kFALSE, Float_t durMin=0, Float_t durMax=10000, Bool_t appBC=kFALSE, Int_t referenceRun=-1)
Definition: anaTree.C:28
Definition: External.C:236
Bool_t writeDetailed
Definition: anaTree.C:23
UInt_t fFirstTime
Definition: createTree.C:64
Double32_t fT6
Definition: createTree.C:50
Short_t fId
Definition: createTree.C:25
Double32_t fT3
Definition: createTree.C:47
Double32_t fT1
Definition: createTree.C:45
TClonesArray fCells
Definition: createTree.C:70
Int_t nBinsT
Definition: anaTree.C:26
Double32_t fMonM
Definition: createTree.C:32
UShort_t T(UShort_t m, UShort_t t)
Definition: RingBits.C:60
int Int_t
Definition: External.C:63
float Float_t
Definition: External.C:68
Short_t fSM
Definition: createTree.C:26
Double32_t fLedR
Definition: createTree.C:31
Definition: External.C:212
Double32_t fT5
Definition: createTree.C:49
Bool_t useWeight
Definition: anaTree.C:25
Int_t debugInfo
Definition: anaTree.C:24
UInt_t fLastTime
Definition: createTree.C:65
Double32_t fT4
Definition: createTree.C:48
Double32_t fLocT
Definition: createTree.C:34
Int_t fRunNo
Definition: createTree.C:62
TClonesArray fSMs
Definition: createTree.C:69
Double32_t fMonR
Definition: createTree.C:33
Double32_t fAvgT
Definition: createTree.C:44
Double32_t fT8
Definition: createTree.C:52
bool Bool_t
Definition: External.C:53
Double32_t fLedM
Definition: createTree.C:30
Double32_t fT7
Definition: createTree.C:51