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