AliPhysics  b76e98e (b76e98e)
anaTreeV3.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 class TCalRun : public TObject {
25 public:
26  TCalRun() : fLedM(-1), fLedR(-1), fMonM(-1), fMonR(-1), fSMT(-1), fLength(0), fRunNo(-1), fBad(0) {;}
27  virtual ~TCalRun() {;}
28  Double32_t fLedM; //[0,0,16] led mean
29  Double32_t fLedR; //[0,0,16] led rms
30  Double32_t fMonM; //[0,0,16] mon mean
31  Double32_t fMonR; //[0,0,16] mon rms
32  Double32_t fSMT; //[0,0,16] sm T
33  UInt_t fLength; // run length in seconds
34  Int_t fRunNo; // run number
35  Short_t fBad; // bad cell
36  ClassDef(TCalRun, 1); // CalRun class
37 };
38 
39 class TCalCellInfo : public TObject {
40 public:
41  TCalCellInfo() : fCellId(-1), fRuns("TCalRun") {}
42  virtual ~TCalCellInfo() {;}
43  Int_t fCellId; // cell IDs
44  TClonesArray fRuns;
45  ClassDef(TCalCellInfo, 1); // CalCellInfo class
46 };
47 
48 Bool_t writeDetailed = kFALSE; // switch on detailed information storing
49 Int_t debugInfo = 1; // enable different levels of debuggin information
51 
52 void anaTreeV3(
53  const char *ifile ="treefile.root",
54  const char *ofile ="outhist.root",
55  Bool_t applyDurLimit = kFALSE, // bolean to switch on duration limit
56  Float_t durMin = 0, // minimum length of a run in minutes
57  Float_t durMax = 10000, // maximum length of a run in minutes
58  Bool_t appBC = kFALSE, // boolean to switch on bad channel
59  Int_t referenceRun = 286313, // define reference run to which all are being calibrated
60  Int_t nBinsT = 1000,
61  TString iterFile = "",
62  Int_t singleSM = -1,
63  Int_t nBins2D = 500
64  )
65 {
66  // Load EMCAL geometry for reference run
67  AliEMCALGeometry *g = AliEMCALGeometry::GetInstanceFromRunNumber(referenceRun);
68  const Int_t kSM = g->GetNumberOfSuperModules();
69  const Int_t kNcells = g->GetNCells();
70 
71  // Return info about optional settings
72  if (applyDurLimit){
73  cout << "INFO: only runs with a length of " << durMin << " to " << durMax << " minutes will be considered in the analysis" << endl;
74  }
75  if (appBC){
76  cout << "INFO: will be using bad channel map" << endl;
77  }
78  Bool_t enable2DWriting = kFALSE;
79  TH1D* referenceLEDDiffLEDMon = NULL;
80  if (iterFile.CompareTo("") != 0){
81  enable2DWriting = kTRUE;
82  TFile* fileIter = new TFile(iterFile);
83  referenceLEDDiffLEDMon = (TH1D*)fileIter->Get("LedDiffLedMonVsCellID");
84  }
85  Bool_t runSingleSM = kFALSE;
86  if (singleSM != -1 && (singleSM >= 0 && singleSM < 21))
87  runSingleSM = kTRUE;
88 
89  // initialize info from tree created by $ALICE_PHYSICS/PWGPP/EMCAL/TeCMacros/createTree.C
90  TCalInfo *info = 0;
91  TFile *in = TFile::Open(ifile,"read");
92  TTree *treeRuns = (TTree*)in->Get("tcal");
93  treeRuns->SetBranchAddress("event",&info);
94  treeRuns->Branch("event", &info, 32000, 99);
95  Int_t Nev = treeRuns->GetEntries();
96 
97  TCalCellInfo *cellinfo = NULL;
98  TTree *treeCells = (TTree*)in->Get("tcalcell");
99  treeCells->SetBranchAddress("cells",&cellinfo);
100  treeCells->Branch("cells", &cellinfo, 32000, 99);
101  Int_t nCellsTree = treeCells->GetEntries();
102 
103  Int_t minRunNo = 1e6;
104  Int_t maxRunNo = -1;
105  for (Int_t i=0;i<Nev;++i) {
106  treeRuns->GetEvent(i);
107  if (info->fRunNo < minRunNo) minRunNo = info->fRunNo;
108  if (info->fRunNo > maxRunNo) maxRunNo = info->fRunNo;
109  }
110  cout << minRunNo << "\t" << maxRunNo << endl;
111 
112  TProfile *gLedVsT[kSM];
113  TProfile *gLedMonVsT[kSM];
114  TProfile *gRatVsT[kSM];
115  TProfile* gCellIdVsRat = new TProfile("", "Led/LedMon run; Cell ID", kNcells+1, -0.5, kNcells+1-0.5);
116  gCellIdVsRat->SetName("LedDiffLedMonVsCellID");
117  TH2D* gCellIdVsLed = new TH2D("", "Led run; Cell ID", kNcells+1, -0.5, kNcells+1-0.5, 1000, 0, 1000);
118  gCellIdVsLed->SetName("LedVsCellID");
119  TH2D* gCellIdVsLedMon = new TH2D("", "LedMon run; Cell ID", kNcells+1, -0.5, kNcells+1-0.5, 1000, 0, 1000);
120  gCellIdVsLedMon->SetName("LedMonVsCellID");
121  TH1F* hAverageT[kSM];
122  TH1F* hAverageTSorted[kSM];
123  TH2F* hAverageTPerSM = new TH2F("","T per SM; SM; T", kSM, -0.5, kSM-0.5, nBinsT, 15, 40);
124  hAverageTPerSM->SetName("MeanSMTemperature");
125  hAverageTPerSM->Sumw2();
126  const char* opt = "S"; //"S" for spread
127  for (Int_t i=0;i<kSM;++i) {
128  if (runSingleSM && i != singleSM )
129  continue;
130  gLedVsT[i] = new TProfile("","Led info;T;",nBinsT,17, 27);
131  gLedVsT[i]->SetName(Form("ledsm%d",i));
132  gLedMonVsT[i] = new TProfile("","Led info;T;",nBinsT,17, 27);
133  gLedMonVsT[i]->SetName(Form("ledmonsm%d",i));
134  gRatVsT[i] = new TProfile("","Led/LedMon;T",nBinsT,17, 27);
135  gRatVsT[i]->SetName(Form("ledovermonsm%d",i));
136  hAverageT[i] = new TH1F ("",Form("T SM %d ; run ID; T",i),Nev,0.5,Nev+0.5);
137  hAverageT[i]->SetName(Form("TAverageSM%dvsRunId",i));
138  hAverageTSorted[i] = new TH1F ("",Form("T SM %d ; run number; T",i),maxRunNo-minRunNo+1,minRunNo-0.5,maxRunNo+0.5);
139  hAverageTSorted[i]->SetName(Form("TAverageSM%dvsRunNumber",i));
140  }
141 
142  TH1F* hSensorsT[160];
143  TH1F* hSensorsTSorted[160];
144  for (Int_t sens = 0; sens< 160; sens++){
145  hSensorsT[sens] = new TH1F ("",Form("T sensor %d ; run ID; T",sens),Nev,0.5,Nev+0.5);
146  hSensorsT[sens]->SetName(Form("Tsensor%dvsRunId",sens));
147  hSensorsTSorted[sens] = new TH1F ("",Form("T sensor %d ; run number; T",sens),maxRunNo-minRunNo+1,minRunNo-0.5,maxRunNo+0.5);
148  hSensorsTSorted[sens]->SetName(Form("Tsensor%dvsRunNumber",sens));
149  }
150 
151  Bool_t isRefRun = kFALSE;
152  Bool_t hadRefRun = kFALSE;
153 
154  cout << "-> there are " << Nev << " contained in this tree, starting to analyse them" << endl;
155 
156  for (Int_t i=0;i<Nev;++i) {
157  treeRuns->GetEvent(i);
158  UInt_t deltaTimeS = ((UInt_t)info->fLastTime-(UInt_t)info->fFirstTime)/10; // run duration in seconds
159  Float_t deltaTime = ((Float_t)info->fLastTime-(Float_t)info->fFirstTime)/60.; // run duration in minutes
160  if (i%50 == 0)
161  cout << "starting with run " << i << "/" << Nev << endl;
162  cout << info->fRunNo << "\t" << info->fFirstTime << "\t"<<info->fLastTime << "\t"<< deltaTime<< endl;
163 
164  if (applyDurLimit){
165  if (deltaTime < durMin || deltaTime > durMax){
166  cout << "INFO: skipped run due to mismatch in run length" << endl;
167  continue;
168  }
169  }
170 
171  TClonesArray &sms = info->fSMs;
172  for (Int_t sm=0;sm<sms.GetEntries();++sm) {
173  TCalSM *smInfot = static_cast<TCalSM*>(sms.At(sm));
174  if (runSingleSM && sm != singleSM )
175  continue;
176  hAverageT[sm]->SetBinContent(i+1,smInfot->fAvgT);
177  hAverageTSorted[sm]->SetBinContent(hAverageTSorted[sm]->FindBin(info->fRunNo),smInfot->fAvgT);
178  }
179  for (Int_t sm=0;sm<sms.GetEntries();++sm) {
180  TCalSM *smInfot = static_cast<TCalSM*>(sms.At(sm));
181  hAverageTPerSM->Fill(sm,smInfot->fAvgT,deltaTime/60);
182  hSensorsT[sm*8]->SetBinContent(i+1,smInfot->fT1);
183  hSensorsT[sm*8+1]->SetBinContent(i+1,smInfot->fT2);
184  hSensorsT[sm*8+2]->SetBinContent(i+1,smInfot->fT3);
185  hSensorsT[sm*8+3]->SetBinContent(i+1,smInfot->fT4);
186  hSensorsT[sm*8+4]->SetBinContent(i+1,smInfot->fT5);
187  hSensorsT[sm*8+5]->SetBinContent(i+1,smInfot->fT6);
188  hSensorsT[sm*8+6]->SetBinContent(i+1,smInfot->fT7);
189  hSensorsT[sm*8+7]->SetBinContent(i+1,smInfot->fT8);
190  hSensorsTSorted[sm*8]->SetBinContent(hSensorsTSorted[sm*8]->FindBin(info->fRunNo),smInfot->fT1);
191  hSensorsTSorted[sm*8+1]->SetBinContent(hSensorsTSorted[sm*8]->FindBin(info->fRunNo),smInfot->fT2);
192  hSensorsTSorted[sm*8+2]->SetBinContent(hSensorsTSorted[sm*8]->FindBin(info->fRunNo),smInfot->fT3);
193  hSensorsTSorted[sm*8+3]->SetBinContent(hSensorsTSorted[sm*8]->FindBin(info->fRunNo),smInfot->fT4);
194  hSensorsTSorted[sm*8+4]->SetBinContent(hSensorsTSorted[sm*8]->FindBin(info->fRunNo),smInfot->fT5);
195  hSensorsTSorted[sm*8+5]->SetBinContent(hSensorsTSorted[sm*8]->FindBin(info->fRunNo),smInfot->fT6);
196  hSensorsTSorted[sm*8+6]->SetBinContent(hSensorsTSorted[sm*8]->FindBin(info->fRunNo),smInfot->fT7);
197  hSensorsTSorted[sm*8+7]->SetBinContent(hSensorsTSorted[sm*8]->FindBin(info->fRunNo),smInfot->fT8);
198  }
199  }
200 
201  TFile *out = NULL;
202  if (runSingleSM)
203  out = TFile::Open(ofile,"update");
204  else
205  out = TFile::Open(ofile,"recreate");
206  for (Int_t i=0;i<kSM;++i) {
207  if (runSingleSM && i != singleSM )
208  continue;
209  hAverageT[i]->Write();
210  hAverageTSorted[i]->Write();
211  }
212  hAverageTPerSM->Write(hAverageTPerSM->GetName(),TObject::kOverwrite);
213  for (Int_t sens = 0; sens < 160; sens++){
214  hSensorsT[sens]->Write(hSensorsT[sens]->GetName(),TObject::kOverwrite);
215  hSensorsTSorted[sens]->Write(hSensorsTSorted[sens]->GetName(),TObject::kOverwrite);
216  }
217 
218 
219  Int_t smcurr = -1;
220  Int_t cellIDFirstInSM = 0;
221 
222  for (Int_t j = 0; j < nCellsTree; j++){
223  treeCells->GetEvent(j);
224 
225  if (cellinfo->fCellId > kNcells)
226  continue;
227  if (cellinfo->fCellId != j)
228  cout << "WARNING: cellID and counter are out of sync" << endl;
229 
230  if (cellinfo->fCellId%50 == 0)
231  cout << "next 50 cells: " << cellinfo->fCellId << "/" << kNcells<< endl;
232 
233  Int_t iTower = -1, iIphi = -1, iIeta = -1, sm=-1;
234  g->GetCellIndex(cellinfo->fCellId,sm,iTower,iIphi,iIeta);
235 
236  TRandom3 ledRandom(1289790);
237  TRandom3 ledMonRandom(909088);
238 
239  TClonesArray &runs = cellinfo->fRuns;
240 
241  if (runSingleSM && sm != singleSM )
242  continue;
243 
244  if (smcurr != sm){
245  smcurr = sm;
246  out->mkdir(Form("cellsInSM%d",sm));
247  out->cd(Form("cellsInSM%d",sm));
248  if (debugInfo) cout << "SM: \t" << sm-1 << "\t"<< cellIDFirstInSM << "\t" << j << "\t" << j -cellIDFirstInSM<< endl;
249  cellIDFirstInSM = j;
250  }
251 
252  vector<Double_t> vecCellInfo[4];// = new vector<Double_t>[4]; // T x led x ledM x RunNo
253  Int_t nEntriesCell = 0;
254 
255  TProfile *gRatCellVsT = new TProfile("",Form("Led/LedMon cell ID%i ;T;",j),nBinsT,17, 27,opt);
256  gRatCellVsT->SetName(Form("ledovermonCell%d",j));
257 
258  for (Int_t r=0;r<runs.GetEntries();++r) {
259  TCalRun *run = static_cast<TCalRun*>(runs.At(r));
260 
261 // if (r%50 == 0)
262 // cout << "starting with run " << r << "/" << Nev << endl;
263 
264  Int_t badcell = run->fBad;
265  Double_t ledM = run->fLedM;
266  Double_t ledR = run->fLedR;
267  Double_t monM = run->fMonM;
268  Double_t monR = run->fMonR;
269  Double_t smT = run->fSMT;
270  ULong_t deltaTimeS= (ULong_t)run->fLength/10;
271  Int_t runNo = run->fRunNo;
272 
273  if (runSingleSM && sm != singleSM )
274  continue;
275 
276  if (appBC && badcell > 0){
277  if (debugInfo > 1) cout << "found bad cell for " << cellinfo->fCellId << "\t " << runNo << endl;
278  continue;
279  }
280  if ((smT<5)||(smT>45))
281  continue;
282 
283  for ( UInt_t s = 0; s < deltaTimeS; s++){
284  Double_t ledcurrent = 0;
285  Double_t ledMoncurrent = 0;
286  if (!((ledM<=0)||(ledR<=0)) ){
287  ledcurrent = ledRandom.Gaus(ledM,ledR);
288  gLedVsT[sm]->Fill(smT,ledcurrent);
289  gCellIdVsLed->Fill(cellinfo->fCellId,ledcurrent);
290  }
291  if (!((monM<=0)||(monR<=0)) ){
292  ledMoncurrent = ledMonRandom.Gaus(monM,monR);
293  gLedMonVsT[sm]->Fill(smT,ledMoncurrent);
294  gCellIdVsLedMon->Fill(cellinfo->fCellId,ledMoncurrent);
295  }
296  if ( (ledcurrent<=0) || (ledMoncurrent<=0) )
297  continue;
298  gCellIdVsRat->Fill(cellinfo->fCellId,ledcurrent/ledMoncurrent);
299  gRatVsT[sm]->Fill(smT,ledcurrent/ledMoncurrent);
300  gRatCellVsT->Fill(smT,ledcurrent/ledMoncurrent);
301 
302  vecCellInfo[0].push_back(smT);
303  vecCellInfo[1].push_back(ledcurrent);
304  vecCellInfo[2].push_back(ledMoncurrent);
305  vecCellInfo[3].push_back(runNo);
306  nEntriesCell++;
307  }
308  }
309  if (nEntriesCell > 0){
310  TGraph* graphRatVsT = new TGraph(nEntriesCell);
311  graphRatVsT->SetName(Form("graphRatvsT_Cell%d", cellinfo->fCellId));
312  TGraph* graphLedVsT = new TGraph(nEntriesCell);
313  graphLedVsT->SetName(Form("graphLedvsT_Cell%d", cellinfo->fCellId));
314  TGraph* graphLedMonVsT = new TGraph(nEntriesCell);
315  graphLedMonVsT->SetName(Form("graphLedMonvsT_Cell%d", cellinfo->fCellId));
316  TGraph* graphRatVsRun = new TGraph(nEntriesCell);
317  graphRatVsRun->SetName(Form("graphRatvsRun_Cell%d", cellinfo->fCellId));
318  for (Int_t k = 0; k< nEntriesCell; k++){
319  graphRatVsT->SetPoint(k, vecCellInfo[0].at(k), vecCellInfo[1].at(k)/vecCellInfo[2].at(k));
320  graphLedVsT->SetPoint(k, vecCellInfo[0].at(k), vecCellInfo[1].at(k));
321  graphLedMonVsT->SetPoint(k, vecCellInfo[0].at(k), vecCellInfo[2].at(k));
322  graphRatVsRun->SetPoint(k, vecCellInfo[3].at(k), vecCellInfo[1].at(k)/vecCellInfo[2].at(k));
323  }
324  graphRatVsT->Sort();
325  graphLedVsT->Sort();
326  graphLedMonVsT->Sort();
327  graphRatVsRun->Sort();
328 
329  if (gRatCellVsT->GetEntries() > 0)
330  gRatCellVsT->Write();
331  if (graphRatVsT){
332  graphRatVsT->Write();
333  }
334  if (graphRatVsRun){
335  graphRatVsRun->Write();
336  }
337  if (graphLedVsT){
338  graphLedVsT->Write();
339  }
340  if (graphLedMonVsT){
341  graphLedMonVsT->Write();
342  }
343  delete graphLedVsT;
344  delete graphRatVsT;
345  delete graphRatVsRun;
346  delete graphLedMonVsT;
347  }
348  delete gRatCellVsT;
349  }
350 
351  for (Int_t i=0;i<kSM;++i) {
352  if (runSingleSM && i != singleSM )
353  continue;
354  gLedVsT[i]->Write();
355  gLedMonVsT[i]->Write();
356  gRatVsT[i]->Write();
357  hAverageT[i]->Write();
358  hAverageTSorted[i]->Write();
359  }
360  if (runSingleSM){
361  if (out->Get("LedDiffLedMonVsCellID")){
362  TProfile* tempgCellIdVsRat = (TProfile*)out->Get("LedDiffLedMonVsCellID");
363  TH2D* tempgCellIdVsLed = (TH2D*)out->Get("LedVsCellID");
364  TH2D* tempgCellIdVsLedMon = (TH2D*)out->Get("LedMonVsCellID");
365  gCellIdVsRat->Add(gCellIdVsRat);
366  gCellIdVsLed->Add(tempgCellIdVsLed);
367  gCellIdVsLedMon->Add(tempgCellIdVsLedMon);
368  gCellIdVsRat->Write(gCellIdVsRat->GetName(),TObject::kOverwrite);
369  gCellIdVsLed->Write(gCellIdVsLed->GetName(),TObject::kOverwrite);
370  gCellIdVsLedMon->Write(gCellIdVsLedMon->GetName(),TObject::kOverwrite);
371  } else {
372  gCellIdVsRat->Write();
373  gCellIdVsLed->Write();
374  gCellIdVsLedMon->Write();
375  }
376  } else {
377  hAverageTPerSM->Write();
378  gCellIdVsRat->Write();
379  gCellIdVsLed->Write();
380  gCellIdVsLedMon->Write();
381  }
382 }
Int_t fRunNo
Definition: anaTreeV3.C:34
Double32_t fT2
Definition: createTree.C:46
Short_t fBad
Definition: anaTreeV3.C:35
double Double_t
Definition: External.C:58
Double32_t fMonR
Definition: anaTreeV3.C:31
Definition: External.C:236
virtual ~TCalRun()
Definition: anaTreeV3.C:27
UInt_t fFirstTime
Definition: createTree.C:64
Bool_t writeDetailed
Definition: anaTreeV3.C:48
Double32_t fT6
Definition: createTree.C:50
Double32_t fT3
Definition: createTree.C:47
TCalRun()
Definition: anaTreeV3.C:26
Double32_t fT1
Definition: createTree.C:45
Double32_t fLedR
Definition: anaTreeV3.C:29
TClonesArray fRuns
Definition: anaTreeV3.C:44
int Int_t
Definition: External.C:63
unsigned int UInt_t
Definition: External.C:33
float Float_t
Definition: External.C:68
Definition: External.C:228
Definition: External.C:212
Double32_t fSMT
Definition: anaTreeV3.C:32
unsigned long ULong_t
Definition: External.C:38
Double32_t fT5
Definition: createTree.C:49
short Short_t
Definition: External.C:23
UInt_t fLastTime
Definition: createTree.C:65
Double32_t fMonM
Definition: anaTreeV3.C:30
Double32_t fT4
Definition: createTree.C:48
Int_t fRunNo
Definition: createTree.C:62
Int_t debugInfo
Definition: anaTreeV3.C:49
TClonesArray fSMs
Definition: createTree.C:69
Int_t testCells
Definition: anaTreeV3.C:50
virtual ~TCalCellInfo()
Definition: anaTreeV3.C:42
Double32_t fAvgT
Definition: createTree.C:44
Double32_t fT8
Definition: createTree.C:52
bool Bool_t
Definition: External.C:53
UInt_t fLength
Definition: anaTreeV3.C:33
Double32_t fLedM
Definition: anaTreeV3.C:28
Double32_t fT7
Definition: createTree.C:51
void anaTreeV3(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, TString iterFile="", Int_t singleSM=-1, Int_t nBins2D=500)
Definition: anaTreeV3.C:52
Int_t fCellId
Definition: anaTreeV3.C:43