CMS 3D CMS Logo
HIPplots.cc
Go to the documentation of this file.
1 #include "TROOT.h"
2 #include "TAttFill.h"
3 #include "TColor.h"
4 #include "HIPplots.h"
5 #include <string>
6 #include <sstream>
7 #include <cmath>
8 
9 #include "TProfile.h"
10 #include "TPaveStats.h"
11 #include "TList.h"
12 #include "TNtuple.h"
13 #include "TString.h"
14 #include "TTree.h"
15 
16 #include "TMatrixD.h"
17 #include "TVectorD.h"
18 #include "TLegend.h"
19 
20 HIPplots::HIPplots(int IOV, char* path, char* outFile):
21 _IOV(IOV),
22 _path(path),
23 _outFile(outFile)
24 {
25  //Id,ObjId,Nhit,Jtvj,Jtve,AlignableChi2,AlignableNdof
26  _inFile_uservars = Form("%s/IOUserVariables.root", _path.Data());
27  if (!CheckFileExistence(_inFile_uservars)) _inFile_uservars = Form("%s/IOUserVariables_%d.root", _path.Data(), IOV);
28 
29  //aligned absolute positions
30  _inFile_alipos = Form("%s/IOAlignedPositions.root", _path.Data());
31  if (!CheckFileExistence(_inFile_alipos)) _inFile_alipos = Form("%s/IOAlignedPositions_%d.root", _path.Data(), IOV);
32 
33  _inFile_HIPalign = Form("%s/HIPAlignmentAlignables.root", _path.Data());
34  if (!CheckFileExistence(_inFile_HIPalign)) _inFile_HIPalign = Form("%s/HIPAlignmentAlignables_%d.root", _path.Data(), IOV);
35 
36  _inFile_params = Form("%s/IOAlignmentParameters.root", _path.Data()); //Id (raw id), ObjId, Par (alignable movements), CovariantMatrix
37  _inFile_mispos = Form("%s/IOMisalignedPositions.root", _path.Data());
38  _inFile_surveys = Form("%s/HIPSurveyResiduals.root", _path.Data());
39  _inFile_truepos = Form("%s/IOTruePositions.root", _path.Data()); //redundant
40 
41  SetPeakThreshold(8.0);
42  plotbadchi2=true;
43 }
44 
45 
46 TLegend* HIPplots::MakeLegend(double x1,
47  double y1,
48  double x2,
49  double y2)
50 {
51  TLegend* legend = new TLegend(x1, y1, x2, y2, "", "NBNDC");
52  legend->SetNColumns(6);
53  legend->SetFillColor(0);
54  legend->SetBorderSize(0);
55  int COLOR_CODE[6]={ 28, 2, 3, 4, 6, 7 };
56 
57  // TO BE UPDATED FOR PHASE-1
58  TString detNames[6];
59  detNames[0] = TString("PXB");
60  detNames[1] = TString("PXF");
61  detNames[2] = TString("TIB");
62  detNames[3] = TString("TID");
63  detNames[4] = TString("TOB");
64  detNames[5] = TString("TEC");
65 
66  for (unsigned int isublevel = 0; isublevel < 6; isublevel++)
67  {
68  TGraph* g = new TGraph(0);
69  g->SetLineColor(COLOR_CODE[isublevel]);
70  g->SetMarkerColor(COLOR_CODE[isublevel]);
71  g->SetFillColor(COLOR_CODE[isublevel]);
72  g->SetFillColor(1);
73  g->SetMarkerStyle(8);
74  g->SetMarkerSize(1);
75  legend->AddEntry(g, detNames[isublevel], "lp");
76  }
77  return legend;
78 }
79 
80 
81 void HIPplots::extractAlignParams(int currentPar, int minHits, int subDet, int doubleSided){
82 
83  cout << "--- extracting AlignParams ; Par " << currentPar << " ---" << endl << endl; //0-5: u,v,w,alpha,beta,gamma
84  //void ExtractAlignPars(int pp, string SubDetString){
85  // const int par_now=pp;
86  //load first tree
87 
88  // TFile* fp = new TFile( _inFile_params,"READ");
89  // const TList* keysp = fp->GetListOfKeys();
90  // const unsigned int maxIteration = keysp->GetSize() - 1;
91  cout << "Loaded par file -. OK" << endl;
92  TFile* fv = new TFile(_inFile_uservars, "READ");
93  const TList* keysv = fv->GetListOfKeys();
94  const unsigned int maxIteration = keysv->GetSize() - 1;
95 
96  char fileaction[16];
97  if (CheckFileExistence(_outFile))sprintf(fileaction, "UPDATE");
98  else sprintf(fileaction, "NEW");
99 
100  TFile* fout = new TFile(_outFile, fileaction);
101 
102  TTree* tree0 = (TTree*)fv->Get(keysv->At(0)->GetName());
103  //const char* tree0v_Name = keysv->At(0)->GetName();
104  //tree0->AddFriend( tree0v_Name, fv );
105 
106  int nHit;
107  unsigned int detId;
108  double par[6];
109  unsigned int detId0;
110 
111  tree0->SetBranchAddress("Id", &detId0);
112 
113  const int ndets=tree0->GetEntries();
114  TH1D* hpar1[ndets];
115  char ppdirname[16];
116  sprintf(ppdirname, "ShiftsPar%d", currentPar);
117  fout->mkdir(ppdirname);
118  gDirectory->cd(ppdirname);
119  //delare histos
120 
121 
122  TH1D* hpariter[maxIteration];
123  for (unsigned int a = 0; a < maxIteration; a++){
124  char histoname[32], histotitle[32];
125  sprintf(histoname, "Par_%d_Iter_%d", currentPar, a);
126  sprintf(histotitle, "Par %d for iteration #%d", currentPar, a);
127  hpariter[a]=new TH1D(histoname, histoname, 1000, -1.0, 1.0);
128  }
129 
130  for (int i = 0; i < ndets; i++){
131  char histoname[32], histotitle[32];
132  sprintf(histoname, "Par_%d_SiDet_%d", currentPar, i);
133  sprintf(histotitle, "Par %d for detector #%d", currentPar, i);
134  hpar1[i]=new TH1D(histoname, histoname, maxIteration+1, -0.5, maxIteration+0.5);
135  }
136 
137  //file where to save aligned det id list
138  ofstream flist("./List_aligned_dets.txt", ios::out);
139 
140 
141  //Loop on trees
142  int modules_accepted = 0;
143  for (unsigned int iter = 0; iter <= maxIteration; iter++){//loop on i (HIP iterations -> trees in file)
144 
145  // TTree* tmpTree = (TTree*) fp->Get(keysp->At(iter)->GetName() );
146  // const char* tmpTreev = keysv->At(iter)->GetName();
147  TTree* tmpTree = (TTree*)fv->Get(keysv->At(iter)->GetName());
148  // tmpTree->AddFriend( tmpTreev, fv ); //get access to both IOAlignmentParameters and IOUserVariables
149  tmpTree->SetBranchAddress("Id", &detId);
150  tmpTree->SetBranchAddress("Nhit", &nHit); //taken from IOUserVariables
151  tmpTree->SetBranchAddress("Par", &par); //taken from IOAlignmentParameters
152 
153  std::cout << "iteration: " << iter << "..." << std::endl;
154 
155  modules_accepted=0;
156  //loop on entries
157  for (int j = 0; j < tmpTree->GetEntries(); j++){ //loop on j (modules -> entries in each tree)
158 
159  tmpTree->GetEntry(j);
160  bool passSubdetCut = true;
161  if (subDet > 0){ if (GetSubDet(detId) != subDet) passSubdetCut = false; } //find the subDet from module detId
162  if (doubleSided > 0){
163  if (GetSubDet(detId)%2 == 1 && doubleSided == 1 && GetBarrelLayer(detId) < 3) passSubdetCut = false;
164  if (GetSubDet(detId)%2 == 1 && doubleSided == 2 && GetBarrelLayer(detId) > 2) passSubdetCut = false;
165  }
166 
167  if ((nHit >= minHits)&&(passSubdetCut)){ //cut on min nhits per module
168  hpar1[j]->SetBinContent(iter+1, par[currentPar]); //x-iteration, y-movement in this iteration
169  //std::cout << "iteration: " << iter << ",par" << currentPar << "= " << par[currentPar] << std::endl;
170 
171  int COLOR_CODE[6]={ 28, 2, 3, 4, 6, 7 };
172  hpar1[j]->SetLineColor(COLOR_CODE[GetSubDet(detId)-1]);
173  hpar1[j]->SetMarkerColor(COLOR_CODE[GetSubDet(detId)-1]);
174 
175  //hpariter[iter]->Fill(par[currentPar]);
176  //if(currentPar<3&&par[currentPar]>0.5)cout << "LARGE PARCHANGE from DET " << detId << " in subdet " << subDet << ". Par#" << currentPar << " shifted by " << par[currentPar] << " at iter " << iter << endl;
177 
178  if ((iter==maxIteration-1)&&(currentPar==0))flist << detId << endl;
179  modules_accepted++;
180 
181  }
182  }//end loop on entries (i.e. modules)
183  delete tmpTree;
184  }//end loop on iteration
185  std::cout << "Modules accepted: " << modules_accepted << std::endl;
186  std::cout << "Writing..." << std::endl;
187  //Save
188  fout->Write();
189  flist.close();
190  //delete hpar1;
191  std::cout << "Deleting..." << std::endl;
192  delete fout;
193  // std::cout << "Deleting..." << std::endl;
194  // delete fp;
195  //std::cout << "Deleting..." << std::endl;
196  delete fv;
197 }
198 
199 void HIPplots::extractAlignShifts(int currentPar, int minHits, int subDet){
200 
201  cout << "\n--- extracting AlignShifts ; Par " << currentPar << " ---" << endl << endl;
202  TFile* fa = new TFile(_inFile_alipos, "READ");
203  const TList* keysa = fa->GetListOfKeys();
204  const unsigned int maxIteration = keysa->GetSize();
205 
206  TFile* fv = new TFile(_inFile_uservars, "READ");
207  const TList* keysv = fv->GetListOfKeys();
208 
209  // TFile* fm = new TFile( _inFile_mispos,"READ");
210  //TFile* ft = new TFile( _inFile_truepos,"READ");
211 
212  char fileaction[16];
213  if (CheckFileExistence(_outFile))sprintf(fileaction, "UPDATE");
214  else sprintf(fileaction, "NEW");
215  TFile* fout = new TFile(_outFile, fileaction);
216 
217  TTree* tree0 = (TTree*)fa->Get(keysa->At(0)->GetName());
218  const char* tree0v_Name = keysv->At(0)->GetName();
219  tree0->AddFriend(tree0v_Name, fv); //get access to both IOAlignedPositions and IOUserVariables
220 
221  unsigned int detId0;
222 
223  tree0->SetBranchAddress("Id", &detId0);
224 
225  const int ndets=tree0->GetEntries();
226  TH1D* hshift[ndets]; //hshift[i] absolute position change for det i, at all iterations
227  // TH1D* hiter_mis[maxIteration]; //hiter_mis[i] misalignment position change at iteration i (defualt case : misalignment=0)
228  TH1D* hiter_ali[maxIteration]; //hiter_ali[i] absolute position change at iteration i, for all detectors
229  char ppdirname[16];
230  sprintf(ppdirname, "Shifts%d", currentPar);
231  fout->mkdir(ppdirname);
232  gDirectory->cd(ppdirname);
233  //delare histos
234  for (unsigned int a = 0; a < maxIteration; a++){
235  char histoname[64], histotitle[64];
236  // sprintf(histoname,"MisalignedShift_%d_Iter_%d",currentPar,a);
237  // sprintf(histotitle,"Misaligned Shift %d for iteration #%d",currentPar,a);
238  // hiter_mis[a]=new TH1D(histoname,histoname,1000,-1.0,1.0);
239  sprintf(histoname, "AlignedShift_%d_Iter_%d", currentPar, a);
240  sprintf(histotitle, "Aligned Shift %d for iteration #%d", currentPar, a);
241  hiter_ali[a]=new TH1D(histoname, histoname, ndets, 0, ndets);
242  }
243  for (int i = 0; i < ndets; i++){
244  char histoname[64], histotitle[64];
245  sprintf(histoname, "Shift_%d_SiDet_%d", currentPar, i);
246  sprintf(histotitle, "Shift %d for detector #%d", currentPar, i);
247  hshift[i]=new TH1D(histoname, histoname, maxIteration, -0.5, maxIteration-0.5);
248  }
249 
250  //Loop on trees
251  int modules_accepted = 0;
252  int nHit;
253  unsigned int detId;
254  double tpos[3];
255  double apos[3];
256  // double mpos[3];
257  double trot[9];
258  double arot[9];
259  // double mrot[9];
260  double aerr[6];
261 
262  TString detNames[6];
263  detNames[0] = TString("PXB");
264  detNames[1] = TString("PXF");
265  detNames[2] = TString("TIB");
266  detNames[3] = TString("TID");
267  detNames[4] = TString("TOB");
268  detNames[5] = TString("TEC");
269 
270 
271  // TTree* tmpTree_m = (TTree*)fm->Get("AlignablesAbsPos_1"); //misaligned position
272  //TTree* tmpTree_true = (TTree*)ft->Get("AlignablesOrgPos_1");
273  TTree* tmpTree_true = (TTree*)fa->Get("AlignablesAbsPos_0"); //starting position
274 
275  if (currentPar < 3){
276  // tmpTree_m->SetBranchAddress("Pos", &mpos );
277  tmpTree_true->SetBranchAddress("Pos", &tpos);
278  }
279  if (currentPar >= 3){
280  // tmpTree_m->SetBranchAddress("Rot", &mrot );
281  tmpTree_true->SetBranchAddress("Rot", &trot);
282  }
283 
284  for (unsigned int iter = 1; iter < maxIteration; iter++){//loop on i (HIP iterations -> trees in file)
285 
286  TTree* tmpTree = (TTree*)fa->Get(keysa->At(iter)->GetName());
287  const char* tmpTreev = keysv->At(iter)->GetName();
288  tmpTree->AddFriend(tmpTreev, fv);
289  tmpTree->SetBranchAddress("Id", &detId);
290  tmpTree->SetBranchAddress("Nhit", &nHit);
291  tmpTree->SetBranchAddress("ParError", &aerr);
292 
293  if (currentPar < 3){ tmpTree->SetBranchAddress("Pos", &apos); }
294  tmpTree->SetBranchAddress("Rot", &arot);
295 
296  //std::cout << "iteration: " << iter << "..." << std::endl;
297  modules_accepted=0;
298  //loop on entries
299  for (int j = 0; j < tmpTree->GetEntries(); j++){ //loop on j (modules -> entries in each tree)
300  tmpTree_true->GetEntry(j);
301  tmpTree->GetEntry(j);
302  //cout << "det" << j << "," << detId << "," << nHit << endl;
303  // tmpTree_m->GetEntry(j);
304  // tmpTree_true->GetEntry(j);
305 
306 
307  bool passSubdetCut = true;
308  int mysubDet=GetSubDet(detId);
309  if (subDet > 0){ if (GetSubDet(detId) != subDet) passSubdetCut = false; }
310 
311  int COLOR_CODE[6]={ 28, 2, 3, 4, 6, 7 };
312  hshift[j]->SetLineColor(COLOR_CODE[GetSubDet(detId)-1]);
313  hshift[j]->SetMarkerColor(COLOR_CODE[GetSubDet(detId)-1]);
314 
315  if ((nHit >= minHits)&&(passSubdetCut)){
316  //maths
317  if (currentPar < 3){
318  // TVectorD dr_mis(3, mpos);
319  TVectorD dr_ali(3, apos);
320  TVectorD r_true(3, tpos);
321  TMatrixD R_true(3, 3, trot);
322  //std::cout << "dr_ali 00 : " << dr_ali[currentPar] << std::endl;
323  // dr_mis -= r_true;
324  dr_ali -= r_true;
325  //std::cout << "dr_ali 0 : " << dr_ali[currentPar] << std::endl;
326  //to local
327  //if (dr_mis != 0.) dr_mis = R_true* dr_mis;
328  //if (dr_ali != 0.) dr_ali = R_true* dr_ali;
329  //std::cout << "currentPar: " << currentPar << std::endl;
330  //std::cout << "dr_mis: " << dr_mis[currentPar] << std::endl;
331  //std::cout << "dr_ali: " << dr_ali[currentPar] << std::endl;
332  // if(currentPar<3&&dr_ali[currentPar]>1.0)cout << "LARGE SHIFT for DET " << detId << " in subdet " << mysubDet << ". Par#" << currentPar << " shifted by " << dr_ali[currentPar] << " at iter " << iter << endl;
333  hshift[j]->SetBinContent(iter+1, dr_ali[currentPar]); //aligned position - start position
334  // hiter_mis[iter]->SetBinContent(j+1,dr_mis[currentPar]);
335  //if(j=0&&currentPar==5) std::cout << "iter=" << iter << ",dr_ali: " << dr_ali[currentPar] << std::endl;
336  // cout << "iter=" << iter << "dr_ali: " << dr_ali[currentPar] << endl;
337  hiter_ali[iter]->SetBinContent(j+1, dr_ali[currentPar]);
338  // cout << "bin: " << hiter_ali[iter]->GetBinContent(j+1) << endl;
339  // hiter_ali[iter]->SetBinError(j+1,5);
340  }
341  if (currentPar >= 3){
342  // TMatrixD dR_mis(3, 3, mrot);
343  TMatrixD dR_ali(3, 3, arot);
344  TMatrixD R_true(3, 3, trot);
345  // dR_mis = dR_mis* TMatrixD(TMatrixD::kTransposed, R_true);
346  dR_ali = dR_ali* TMatrixD(TMatrixD::kTransposed, R_true);
347  //-std::atan2(dR(2, 1), dR(2, 2)), std::asin(dR(2, 0)), -std::atan2(dR(1, 0), dR(0, 0)));
348  // double dR_mis_euler = 0;
349  double dR_ali_euler = 0;
350  if (currentPar == 3){
351  // dR_mis_euler = -std::atan2(dR_mis(2, 1), dR_mis(2, 2));
352  dR_ali_euler = -std::atan2(dR_ali(2, 1), dR_ali(2, 2));
353  }
354  if (currentPar == 4){
355  // dR_mis_euler = -std::asin(dR_mis(2, 0));
356  dR_ali_euler = -std::asin(dR_ali(2, 0));
357  }
358  if (currentPar == 5){
359  // dR_mis_euler = -std::atan2(dR_mis(1, 0), dR_mis(0, 0));
360  dR_ali_euler = -std::atan2(dR_ali(1, 0), dR_ali(0, 0));
361  }
362  hshift[j]->SetBinContent(iter+1, dR_ali_euler);
363  // hiter_mis[iter]->SetBinContent(j+1,dR_mis_euler);
364  hiter_ali[iter]->SetBinContent(j+1, dR_ali_euler);
365  }
366  modules_accepted++;
367  hiter_ali[iter]->GetXaxis()->SetBinLabel(j+1, detNames[GetSubDet(detId)-1]);
368  hiter_ali[iter]->SetBinError(j+1, aerr[currentPar]);
369  //hiter_ali[iter]->SetBinError(j+1,0.001);
370  hiter_ali[iter]->LabelsOption("u", "x");
371  }
372  }
373  // delete tmpTree;
374  }
375  delete tmpTree_true;
376 
377  std::cout << "Modules accepted: " << modules_accepted << std::endl;
378  std::cout << "Writing..." << std::endl;
379  //Save
380  fout->Write();
381  //delete hpar1;
382  std::cout << "Deleting..." << std::flush;
383  delete fout;
384  std::cout << "Deleting..." << std::flush;
385  delete fa;
386  // delete ft;
387  // delete fm;
388  std::cout << "Deleting..." << std::endl;
389  delete fv;
390 
391 
392 }
393 
394 void HIPplots::plotAlignParams(string ShiftsOrParams, char* plotName){
395 
396 
397  cout << "_|_| plotting AlignParams |_|_" << endl << "---> " << ShiftsOrParams << endl;
398  bool bParams = false;
399  bool bShifts = false;
400  if (ShiftsOrParams == "PARAMS") bParams = true;
401  if (ShiftsOrParams == "SHIFTS") bShifts = true;
402 
403  int i = 0;
404 
405  TFile* f = new TFile(_outFile, "READ");
406  // f->ls();
407 
408  TCanvas* c_params = new TCanvas("can_params", "CAN_PARAMS", 1200, 900);
409  c_params->Divide(3, 2);
410  // cout << "(1) I am in " << gDirectory->GetPath() << endl;
411  TDirectory* d;
412  int ndets = 0;
413  if (bParams){
414  d = (TDirectory*)f->Get("ShiftsPar0");
415  ndets = GetNIterations(d, "Par_0_SiDet_");
416  // std::cout << "INHERE!" << std::endl;
417  }
418  if (bShifts){
419  d = (TDirectory*)f->Get("Shifts0");
420  ndets = GetNIterations(d, "Shift_0_SiDet_");
421  // std::cout << "INHERE!" << std::endl;
422  }
423 
424  for (int iPar = 0; iPar < 6; iPar++){
425 
426  c_params->cd(iPar+1);
427  gPad->SetTopMargin(0.15);
428  gPad->SetBottomMargin(0.15);
429  char ppdirname[32];
430  if (bParams) sprintf(ppdirname, "ShiftsPar%d", iPar);
431  if (bShifts) sprintf(ppdirname, "Shifts%d", iPar);
432  if (iPar > 0)gDirectory->cd("../");
433  gDirectory->cd(ppdirname);
434  // cout << "(2) I am in " << gDirectory->GetPath() << endl;
435 
436  TH1D* hpar1[ndets];
437  char histoname[16];
438  int sampling_ratio=1;
439  int ndets_plotted=(int)ndets/sampling_ratio;
440  cout << "Plotting " << ndets_plotted << " detectors over a total of " << ndets << endl;
441  i=0;
442  double histomax, histomin;
443  if (iPar>=3){ //alpha, beta, gamma
444  histomax=0.5;//in mrad
445  histomin=-0.5;
446  }
447  else if (iPar>=2) { //w
448  if (bShifts){
449  histomax=200.0;//in microns
450  histomin=-200.0;
451  }
452  else{
453  histomax=100.0;//in microns
454  histomin=-100.0;
455  }
456  }
457  else { //u, v
458  if (bShifts){
459  histomax=200.0;//in microns
460  histomin=-200.0;
461  }
462  else{
463  histomax=100.0;//in microns
464  histomin=-100.0;
465  }
466  }
467  while ((i<ndets_plotted) && (i*sampling_ratio<ndets)){
468  if (bParams) sprintf(histoname, "Par_%d_SiDet_%d", iPar, i*sampling_ratio);
469  if (bShifts) sprintf(histoname, "Shift_%d_SiDet_%d", iPar, i*sampling_ratio);
470  hpar1[i]=(TH1D*)gDirectory->Get(histoname);
471 
472  if (iPar>=3)hpar1[i]->Scale(1000.0); //convert from rad to mrad
473  else hpar1[i]->Scale(10000.0); //convert from cm to um
474 
475  hpar1[i]->SetMarkerStyle(7);
476  hpar1[i]->SetStats(0);
477 
478  double tmpmax=hpar1[i]->GetBinContent(hpar1[i]->GetMaximumBin());
479  double tmpmin=hpar1[i]->GetBinContent(hpar1[i]->GetMinimumBin());
480 
481  //Auto-adjust axis range
482  if (tmpmax>histomax)histomax=tmpmax*1.2;
483  if (tmpmin<histomin)histomin=tmpmin*1.2;
484 
485  //if(i%1300==0)cout << "Actual maximum is " << histomax << " Actual minimum is " << histomin << endl;
486  if (i==0){
487 
488  hpar1[i]->SetXTitle("Iteration");
489 
490  if (bParams){
491  if (iPar==0)hpar1[i]->SetYTitle("#delta u param (#mum)");
492  else if (iPar==1)hpar1[i]->SetYTitle("#delta v param (#mum)");
493  else if (iPar==2)hpar1[i]->SetYTitle("#delta w param (#mum)");
494  else if (iPar==3)hpar1[i]->SetYTitle("#delta #alpha param (mrad)");
495  else if (iPar==4)hpar1[i]->SetYTitle("#delta #beta param (mrad)");
496  else if (iPar==5)hpar1[i]->SetYTitle("#delta #gamma param (mrad)");
497  else hpar1[i]->SetYTitle("dunno");
498  }
499  else if (bShifts){
500  if (iPar==0)hpar1[i]->SetYTitle("#delta u shift (#mum)");
501  else if (iPar==1)hpar1[i]->SetYTitle("#delta v shift (#mum)");
502  else if (iPar==2)hpar1[i]->SetYTitle("#delta w shift (#mum)");
503  else if (iPar==3)hpar1[i]->SetYTitle("#delta #alpha shift (mrad)");
504  else if (iPar==4)hpar1[i]->SetYTitle("#delta #beta shift (mrad)");
505  else if (iPar==5)hpar1[i]->SetYTitle("#delta #gamma shift (mrad)");
506  else hpar1[i]->SetYTitle("dunno");
507  }
508 
509  hpar1[i]->GetYaxis()->SetTitleOffset(1.5);
510  //hpar1[i]->SetTitle("Par1: x shift");
511  hpar1[i]->SetTitle("");
512  hpar1[i]->SetMaximum(histomax);
513  hpar1[i]->SetMinimum(histomin);
514  hpar1[i]->Draw("PL");
515  }
516 
517  else hpar1[i]->Draw("PLsame");
518  i++;
519  }//end loop on NDets
520  hpar1[0]->SetMaximum(histomax);
521  hpar1[0]->SetMinimum(histomin);
522  //hpar1[0]->SetLineColor(1);
523  cout << "Plotted " << i << " aligned detectors" << endl;
524  }//end loop on pars
525  c_params->cd();
526  TLegend* legend = MakeLegend(.1, .93, .9, .98);
527  legend->Draw();
528  c_params->SaveAs(plotName);
529  std::cout << "Deleting..." << std::flush;
530  delete c_params;
531  std::cout << "Deleting..." << std::flush;
532  //delete f;
533  std::cout << "Deleting..." << std::endl;
534 
535 }//end PlotAlignPars
536 
537 
538 
539 void HIPplots::plotAlignParamsAtIter(int iter, string ShiftsOrParams, char* plotName){
540 
541  cout << "Welcome to HIPplots::plotAlignParamsAtIter " << iter << endl;
542  bool bParams = false;
543  bool bShifts = false;
544  if (ShiftsOrParams == "PARAMS") bParams = true;
545  else if (ShiftsOrParams == "SHIFTS") bShifts = true;
546  else { cout << "ERROR in plotAliParamsAtIter!!! Wrong input argument: " << ShiftsOrParams << " . Exiting" << endl; return; }
547 
548  int i = 0;
549  TFile* f = new TFile(_outFile, "READ");
550  //f->ls();
551 
552  TCanvas* c_params = new TCanvas("can_params", "CAN_PARAMS", 1200, 700);
553  c_params->Divide(3, 2);
554  //cout << "(1) I am in " << gDirectory->GetPath() << endl;
555 
556  //TDirectory* d = (TDirectory*)f->Get("ShiftsPar0");
557  //const int ndets = GetNIterations(d,"Par_0_SiDet_");
558 
559  for (int iPar = 0; iPar < 6; iPar++){
560 
561  c_params->cd(iPar+1);
562 
563  gPad->SetTopMargin(0.15);
564  gPad->SetBottomMargin(0.15);
565  char ppdirname[32];
566  if (bParams) sprintf(ppdirname, "ShiftsPar%d", iPar);
567  if (bShifts) sprintf(ppdirname, "Shifts%d", iPar);
568  if (iPar > 0)gDirectory->cd("../");
569  gDirectory->cd(ppdirname);
570  //cout << "(2) I am in " << gDirectory->GetPath() << endl;
571 
572  TH1D* hiter;
573  char histoname[16];
574  if (bParams) sprintf(histoname, "Par_%d_Iter_%d", iPar, iter);
575  if (bShifts) sprintf(histoname, "AlignedShift_%d_Iter_%d", iPar, iter);
576  hiter = (TH1D*)gDirectory->Get(histoname);
577 
578 
579  //hiter->SetMarkerStyle(1);
580  hiter->SetStats(1);
581 
582  hiter->SetXTitle("Sub-Det");
583  if (iPar==0)hiter->SetYTitle("#delta u (#mum)");
584  else if (iPar==1)hiter->SetYTitle("#delta v (#mum)");
585  else if (iPar==2)hiter->SetYTitle("#delta w (#mum)");
586  else if (iPar==3)hiter->SetYTitle("#delta #alpha (#murad)");
587  else if (iPar==4)hiter->SetYTitle("#delta #beta (#murad)");
588  else if (iPar==5)hiter->SetYTitle("#delta #gamma (#murad)");
589  else hiter->SetXTitle("dunno");
590  hiter->GetYaxis()->SetTitleOffset(1.5);
591 
592  double histomax, histomin;
593  if (iPar==2||iPar==5) { histomax=200; histomin=-200; }
594  else { histomax=100; histomin=-100; }
595 
596  if (iPar>=3)hiter->Scale(1000000.0); //convert from rad to micro rad
597  else hiter->Scale(10000.0); //convert from cm to um
598  double tmpmax=hiter->GetBinContent(hiter->GetMaximumBin());
599  double tmpmin=hiter->GetBinContent(hiter->GetMinimumBin());
600  if (tmpmax>histomax)histomax=tmpmax*1.2;
601  if (tmpmin<histomin)histomin=tmpmin*1.2;
602  hiter->SetMaximum(histomax);
603  hiter->SetMinimum(histomin);
604  //hiter->SetAxisRange(0, 6, "X");
605  hiter->SetLineColor(1);
606  TColor* col = gROOT->GetColor(kCyan+1);
607  col->SetAlpha(0.3);
608  hiter->SetFillColor(col->GetNumber());
609  // hiter->SetFillColor(kCyan);
610  hiter->SetMarkerStyle(7);
611  hiter->Draw();
612  hiter->Draw("PE1");
613  /*
614  if (bShifts) {
615  TH1D* hiter2;
616  char histoname2[16];
617  sprintf( histoname2, "MisalignedShift_%d_Iter_%d", iPar, iter );
618  hiter2 = (TH1D*) gDirectory->Get(histoname2);
619  hiter2->SetLineColor(1);
620  hiter2->SetFillColor(kRed-9);
621  hiter2->Draw("SAME");
622  }
623 
624  */
625  // cout << "Plotted " << i << " aligned detectors" << endl;
626  }//end loop on pars
627  c_params->SaveAs(plotName);
628  delete c_params;
629  std::cout << "Deleting..." << std::endl;
630  delete f;
631 
632 }//end PlotAlignParamsAtIter
633 
634 
635 
636 void HIPplots::extractAlignableChiSquare(int minHits, int subDet, int doubleSided){
637 
638  //TFile* fp = new TFile( _inFile_params,"READ");
639  // const TList* keysp = fp->GetListOfKeys();
640  // const unsigned int maxIteration = keysp->GetSize();
641  cout << "\n--- Welcome to extractAlignableChiSquare ---" << endl;
642  cout << "\nInput parameters:\n\tMinimum number of hits per alignbale = " << minHits << "\n\tSubdetetctor selection" << subDet << endl;
643 
644  if (minHits<1){
645  cout << "Warning ! Allowing to select modules with NO hits. Chi2 not defined for them. Setting automatically minNhits=1 !!!" << endl;
646  minHits=1;
647  }
648 
649  TFile* fv = new TFile(_inFile_uservars, "READ");
650  const TList* keysv = fv->GetListOfKeys();
651  const unsigned int maxIteration = keysv->GetSize() - 1;
652  cout << "MaxIteration is " << maxIteration << endl;
653 
654  char fileaction[16];
655  if (CheckFileExistence(_outFile))sprintf(fileaction, "UPDATE");
656  else sprintf(fileaction, "NEW");
657  TFile* fout = new TFile(_outFile, fileaction);
658 
659  TTree* tree0 = (TTree*)fv->Get(keysv->At(0)->GetName());
660  unsigned int detId0;
661  tree0->SetBranchAddress("Id", &detId0);
662  const int ndets=tree0->GetEntries();
663  TH1D* halichi2n[ndets];//norm chi2 for each module as a function of iteration
664  TH1D* htotchi2n[maxIteration];//distrib of norm chi2 for all modules at a given iteration
665  TH1D* hprobdist[maxIteration];
666  char ppdirname[16];
667  sprintf(ppdirname, "AlignablesChi2n");
668  fout->mkdir(ppdirname);
669  gDirectory->cd(ppdirname);
670 
671  for (int iter = 1; iter <=(int)maxIteration; iter++){
672  char histoname[64], histotitle[128];
673  sprintf(histoname, "Chi2n_iter%d", iter);
674  sprintf(histotitle, "Distribution of Normalised #chi^{2} for All alignables at iter %d", iter);
675  htotchi2n[iter-1]=new TH1D(histoname, histotitle, 1000, 0.0, 10.0);
676  sprintf(histoname, "ProbChi2_%d", iter);
677  sprintf(histotitle, "Distribution of Prob(#chi^{2},ndof) at iter %d", iter);
678  hprobdist[iter-1]=new TH1D(histoname, histotitle, 100, 0.0, 1.0);
679  }
680  gDirectory->mkdir("AlignablewiseChi2n");
681  gDirectory->cd("AlignablewiseChi2n");
682  //declare histos
683  for (int i = 0; i < (int)ndets; i++){
684  tree0->GetEntry(i);
685  char histoname[64], histotitle[64];
686  sprintf(histoname, "Chi2n_%d", i);
687  sprintf(histotitle, "Normalised #chi^{2} for detector #%d", i);
688  halichi2n[i]=new TH1D(histoname, histotitle, maxIteration, 0.5, maxIteration+0.5);
689  }
690 
691  //Loop on trees and fill histos
692  int modules_accepted = 0;
693  for (unsigned int iter = 1; iter <= maxIteration; iter++){//loop on i (HIP iterations -> trees in file)
694 
695  TTree* tmpTreeUV = (TTree*)fv->Get(keysv->At(iter)->GetName()); //get the UserVariable tree at each iteration
696  cout << "Taking tree " << keysv->At(iter)->GetName() << endl;
697  //tmpTreeUV->GetListOfLeaves()->ls();
698  tmpTreeUV->SetBranchStatus("*", 0);
699  tmpTreeUV->SetBranchStatus("Id", 1);
700  tmpTreeUV->SetBranchStatus("Nhit", 1);
701  tmpTreeUV->SetBranchStatus("AlignableChi2", 1);
702  tmpTreeUV->SetBranchStatus("AlignableNdof", 1);
703  double alichi2=0.0;
704  unsigned int alindof=0;
705  int nHit=0;
706  unsigned int detId=0;
707  tmpTreeUV->SetBranchAddress("AlignableChi2", &alichi2);
708  tmpTreeUV->SetBranchAddress("AlignableNdof", &alindof);
709  tmpTreeUV->SetBranchAddress("Id", &detId);
710  tmpTreeUV->SetBranchAddress("Nhit", &nHit);
711 
712  modules_accepted=0;
713  //loop on entries
714  for (int j = 0; j < tmpTreeUV->GetEntries(); j++){ //loop on j (modules -> entries in each tree)
715  tmpTreeUV->GetEntry(j);
716 
717  bool passSubdetCut = true;
718  if (subDet > 0){ if (GetSubDet(detId) != subDet) passSubdetCut = false; } //get subDet Id from module detId,and plot only the selected subDet
719  if (doubleSided > 0){
720  if (GetSubDet(detId)%2 == 1 && doubleSided == 1 && GetBarrelLayer(detId) < 3) passSubdetCut = false;
721  if (GetSubDet(detId)%2 == 1 && doubleSided == 2 && GetBarrelLayer(detId) > 2) passSubdetCut = false;
722  }
723 
724  if ((nHit >= minHits)&&(passSubdetCut)){ //cut on nhits per module
725  halichi2n[j]->SetBinContent(iter, double(alichi2 / alindof));
726  // halichi2n[j]->SetBinContent(iter,double(alichi2));
727  double prob=TMath::Prob(double(alichi2), int(alindof));
728 
729  htotchi2n[iter-1]->Fill(double(alichi2 / alindof));
730  hprobdist[iter-1]->Fill(prob);
731  modules_accepted++;
732  }
733 
734  }//end loop on j - alignables
735  cout << "alignables accepted at iteration " << iter << " = " << modules_accepted << endl;
736  delete tmpTreeUV;
737  }//end loop on iterations
738 
739  //Ma che e'???
740  /* cout << "Prob for chi2=0,02 ndof=40 -> " << TMath::Prob(0.02,40) << endl;
741  cout << "Prob for chi2=20, ndof=40 -> " << TMath::Prob(20.0,40) << endl;
742  cout << "Prob for chi2=40, ndof=40 -> " << TMath::Prob(40.0,40) << endl;
743  cout << "Prob for chi2=60, ndof=40 -> " << TMath::Prob(60.0,40) << endl;
744  cout << "Prob for chi2=2000, ndof=40 -> " << TMath::Prob(2000.0,40) << endl; */
745 
746  //Save
747  fout->Write();
748  delete fout;
749  delete fv;
750  cout << "Finished extractAlignableChiSquare" << endl;
751 }//end extractAlignableChiSquare
752 
753 
754 void HIPplots::extractSurveyResiduals(int currentPar, int subDet){
755 
756  cout << "\n--- extractSurveyResiduals has been called ---" << endl;
757 
758  TFile* fv = new TFile(_inFile_surveys, "READ");
759  const TList* keysv = fv->GetListOfKeys();
760  const unsigned int maxIteration = keysv->GetSize();
761  cout << "MaxIteration is " << maxIteration << endl;
762 
763  char fileaction[16];
764  if (CheckFileExistence(_outFile))sprintf(fileaction, "UPDATE");
765  else sprintf(fileaction, "NEW");
766  TFile* fout = new TFile(_outFile, fileaction);
767 
768  TTree* tree0 = (TTree*)fv->Get(keysv->At(0)->GetName());
769  unsigned int detId0;
770  tree0->SetBranchAddress("Id", &detId0);
771  const int ndets=tree0->GetEntries();
772  TH1D* hsurvey[ndets];//norm chi2 for each module as a function of iteration
773  TH1D* htotres[maxIteration];//distrib of norm chi2 for all modules at a given iteration
774 
775  char ppdirname[16];
776  sprintf(ppdirname, "SurveyResiduals");
777  fout->mkdir(ppdirname);
778  gDirectory->cd(ppdirname);
779 
780  for (int iter = 1; iter <=(int)maxIteration; iter++){
781  char histoname[64], histotitle[128];
782  sprintf(histoname, "SurveyRes_Par%d_iter%d", currentPar, iter);
783  sprintf(histotitle, "Distribution of survey Residuals for All alignables ; Par=%d ; iter %d", currentPar, iter);
784  htotres[iter-1]=new TH1D(histoname, histotitle, 1000, 0.0, 10.0);
785  }
786  gDirectory->mkdir("SurveyResiduals_alignables");
787  gDirectory->cd("SurveyResiduals_alignables");
788  //declare histos
789  for (int i = 0; i < (int)ndets; i++){
790  tree0->GetEntry(i);
791  char histoname[64], histotitle[64];
792  sprintf(histoname, "SurveyRes_Par%d_%d", currentPar, i);
793  sprintf(histotitle, "Survey residual for detector #%d - Par=%d", i, currentPar);
794  hsurvey[i]=new TH1D(histoname, histotitle, maxIteration, 0.5, maxIteration+0.5);
795  }
796 
797  //Loop on trees and fill histos
798  int modules_accepted = 0;
799  for (unsigned int iter = 1; iter <= maxIteration; iter++){//loop on i (HIP iterations -> trees in file)
800 
801  TTree* tmpTreeUV = (TTree*)fv->Get(keysv->At(iter)->GetName()); //get the UserVariable tree at each iteration
802  cout << "Taking tree " << keysv->At(iter)->GetName() << endl;
803  //tmpTreeUV->GetListOfLeaves()->ls();
804  tmpTreeUV->SetBranchStatus("*", 0);
805  tmpTreeUV->SetBranchStatus("Id", 1);
806  tmpTreeUV->SetBranchStatus("Par", 1);
807 
808  double par[6];
809  unsigned int detId=0;
810  tmpTreeUV->SetBranchAddress("Par", &par);
811  tmpTreeUV->SetBranchAddress("Id", &detId);
812 
813 
814  modules_accepted=0;
815  //loop on entries
816  for (int j = 0; j < tmpTreeUV->GetEntries(); j++){ //loop on j (modules -> entries in each tree)
817  tmpTreeUV->GetEntry(j);
818 
819  bool passSubdetCut = true;
820  if (subDet > 0){ if (GetSubDet(detId) != subDet) passSubdetCut = false; }
821  if (passSubdetCut){
822  hsurvey[j]->SetBinContent(iter, double(par[currentPar]));
823 
824  modules_accepted++;
825  }
826 
827  }//end loop on j - alignables
828  cout << "alignables accepted at iteration " << iter << " = " << modules_accepted << endl;
829  delete tmpTreeUV;
830  }//end loop on iterations
831 
832  //Save
833  fout->Write();
834  delete fout;
835  delete fv;
836  cout << "Finished extractAlignableChiSquare" << endl;
837 }//end extractAlignableChiSquare
838 
839 
840 
841 
842 
843 
844 void HIPplots::plotAlignableChiSquare(char* plotName, float minChi2n){
845 
846  int i = 0;
847  TFile* f = new TFile(_outFile, "READ");
848  TCanvas* c_alichi2n=new TCanvas("can_alichi2n", "CAN_ALIGNABLECHI2N", 900, 900);
849  c_alichi2n->cd();
850  TDirectory* chi_d1=(TDirectory*)f->Get("AlignablesChi2n");
851  const int maxIteration= GetNIterations(chi_d1, "Chi2n_iter", 1)-1;
852  cout << "N iterations " << maxIteration << endl;
853  //take the histos prepared with extractAlignableChiSquare
854  gDirectory->cd("AlignablesChi2n");
855  TDirectory* chi_d=(TDirectory*)gDirectory->Get("AlignablewiseChi2n");
856  const int ndets= GetNIterations(chi_d, "Chi2n_");
857 
858  gDirectory->cd("AlignablewiseChi2n");
859  TH1D* hchi2n[ndets];
860  char histoname[64];
861  int sampling_ratio=1;
862  int ndets_plotted=(int)ndets/sampling_ratio;
863  cout << "Sampling " << ndets_plotted << " detectors over a total of " << ndets << endl;
864  double histomax=0.1, histomin=-0.1;
865  bool firstplotted=false;
866  int firstplottedindex=0;
867  int totalplotted=0;
868  bool plothisto=true;
869  while ((i<ndets_plotted) && (i*sampling_ratio<ndets)){
870  sprintf(histoname, "Chi2n_%d", i);
871  hchi2n[i]=(TH1D*)gDirectory->Get(histoname);
872 
873  //if required, check that the chi2n exceeds at any point the threshold set as input
874  bool chi2ncut=false;
875  if (minChi2n>0.0){
876  for (int bin=1; bin<=hchi2n[i]->GetNbinsX(); bin++){
877  if (hchi2n[i]->GetBinContent(bin)>minChi2n){ chi2ncut=true; break; }
878  }
879  }
880  else chi2ncut=true;
881 
882  //if required check that the chi2 is increasing three iterations in a row and plot only those histos
883  if (plotbadchi2)plothisto=CheckHistoRising(hchi2n[i]);
884  else plothisto=true;
885 
886  if (chi2ncut&&plothisto){
887  double tmpmax=hchi2n[i]->GetBinContent(hchi2n[i]->GetMaximumBin());
888  double tmpmin=hchi2n[i]->GetBinContent(hchi2n[i]->GetMinimumBin());
889  histomax=2.0;
890  histomin=0.0;
891  if (tmpmax>histomax)histomax=tmpmax;
892  if (tmpmin<histomin)histomin=tmpmin;
893 
894  hchi2n[i]->SetMaximum(histomax);
895  hchi2n[i]->SetMinimum(histomin);
896  hchi2n[i]->SetStats(0);
897 
898 
899  if (!firstplotted){
900  hchi2n[i]->SetXTitle("Iteration");
901  hchi2n[i]->SetYTitle("#chi^{2} / # dof");
902  // hchi2n[i]->SetYTitle("#chi^{2}");
903  hchi2n[i]->SetTitle("Reduced #chi^{2} for alignables");
904  hchi2n[i]->Draw("PL");
905  firstplotted=true;
906  firstplottedindex=i;
907  }
908  else hchi2n[i]->Draw("PLsame");
909  totalplotted++;
910  }//END IF plot it because it passed the chi2n cut
911  i++;
912  }//end while loop on alignables
913  hchi2n[firstplottedindex]->SetMaximum(histomax*1.2);
914  hchi2n[firstplottedindex]->SetMinimum(histomin*0.8);
915  //override auto-resize of axis scale
916  hchi2n[firstplottedindex]->SetMaximum(histomax);
917  hchi2n[firstplottedindex]->SetMinimum(histomin);
918 
919  cout << "Plotted " << totalplotted << " alignables over an initial sample of " << ndets_plotted << endl;
920  TText* txtchi2n_1=new TText();
921  txtchi2n_1->SetTextFont(63);
922  txtchi2n_1->SetTextSize(22);
923  char strchi2n_1[128];
924  sprintf(strchi2n_1, "Plotted alignables (Chi2n > %.3f): %d / %d", minChi2n, totalplotted, ndets_plotted);
925  txtchi2n_1->DrawText(1.2, 0.0, strchi2n_1);
926  char finplotname[192];
927  sprintf(finplotname, "%s.png", plotName);
928  c_alichi2n->SaveAs(finplotname);
929  delete c_alichi2n;
930  //delete hchi2n;
931 
932 
933  cout << "Doing distrib" << endl;
934  gDirectory->cd("../");
935  TCanvas* c_chi2ndist=new TCanvas("can_chi2ndistr", "CAN_CHI2N_DISTRIBUTION", 900, 900);
936  c_chi2ndist->cd();
937  TH1D* hiter[maxIteration];
938  TLegend* l=new TLegend(0.7, 0.7, 0.9, 0.9);
939  l->SetFillColor(0);
940  l->SetBorderSize(0);
941  int colors[10]={ 1, 2, 8, 4, 6, 7, 94, 52, 41, 45 };
942  int taken=0;
943  float tmpmax=1.0;
944  float newmax=0.0;
945  for (i=0; i<maxIteration; i++){
946  sprintf(histoname, "Chi2n_iter%d", i+1);
947  hiter[i]=(TH1D*)gDirectory->Get(histoname);
948  hiter[i]->SetXTitle("#chi^{2} / dof");
949  hiter[i]->SetYTitle("Alignables");
950  hiter[i]->SetTitle("Normalised #chi^{2} of Alignables");
951  hiter[i]->Rebin(5);
952  hiter[i]->GetXaxis()->SetRangeUser(0.0, 3.0);
953  hiter[i]->SetLineColor(i);
954 
955  char legend_entry[64];
956  float histmean=hiter[i]->GetMean();
957  if (i==0){
958  hiter[i]->SetLineColor(colors[taken]);
959  taken++;
960  sprintf(legend_entry, "Norm #chi^{2}; Iter %d; %.4f", i, histmean);
961  l->AddEntry(hiter[i], legend_entry, "L");
962  tmpmax=hiter[i]->GetBinContent(hiter[i]->GetMaximumBin());
963  newmax=tmpmax;
964  //hiter[i]->Draw("");
965  }
966  else if ((i+1)%5==0){
967  hiter[i]->SetLineColor(colors[taken]);
968  taken++;
969  sprintf(legend_entry, "Norm #chi^{2}; Iter %d; %.4f", i+1, histmean);
970  l->AddEntry(hiter[i], legend_entry, "L");
971  tmpmax=hiter[i]->GetBinContent(hiter[i]->GetMaximumBin());
972  if (tmpmax>newmax)newmax=tmpmax;
973  //hiter[i]->Draw("same");
974  }
975  }
976  cout << "NewMax after 1st loop -> " << newmax << endl;
977 
978  for (i=0; i<maxIteration; i++){
979  hiter[i]->SetMaximum(newmax*1.1);
980  if (i==1) hiter[i]->Draw("");
981  else if ((i+1)%5==0) hiter[i]->Draw("same");
982  }
983  l->Draw();
984 
985  sprintf(finplotname, "%s_distr.png", plotName);
986  cout << finplotname << endl;
987  c_chi2ndist->SaveAs(finplotname);
988  c_chi2ndist->SetLogy();
989  sprintf(finplotname, "%s_distrlog.png", plotName);
990  cout << finplotname << endl;
991  c_chi2ndist->SaveAs(finplotname);
992 
993  delete c_chi2ndist;
994  delete f;
995 }//end plotAlignableChiSquare
996 
997 
998 
999 //-----------------------------------------------------------------
1000 //private classes
1001 int HIPplots::GetNIterations(TDirectory* f, char* tag, int startingcounter){
1002  int fin_iter=0, i=startingcounter;
1003  bool obj_exist=kTRUE;
1004  while (obj_exist){
1005  char objname[32];
1006  sprintf(objname, "%s%d", tag, i);
1007  if (!f->FindObjectAny(objname))obj_exist=kFALSE;
1008  fin_iter=i;
1009  i++;
1010  }
1011  cout << "Max Iterations is " << fin_iter << endl;
1012 
1013  return fin_iter;
1014 }
1015 
1016 int HIPplots::GetSubDet(unsigned int id){
1017 
1018  // TO BE UPDATED FOR PHASE-1
1019  const int reserved_subdetectorstartbit=25;
1020  const int reserved_subdetectorfinalbit=27;
1021 
1022  unsigned int detID = id;
1023 
1024  int shift = 31-reserved_subdetectorfinalbit;
1025  detID = detID << (shift);
1026  shift = reserved_subdetectorstartbit + shift;
1027  detID = detID>>(shift);
1028 
1029  return detID;
1030 
1031 }
1032 
1033 int HIPplots::GetBarrelLayer(unsigned int id){
1034 
1035  const int reserved_layerstartbit=14;
1036  const int reserved_layermask=0x7;
1037 
1038  return int((id>>reserved_layerstartbit) & reserved_layermask);
1039 
1040 }
1041 
1042 void HIPplots::SetMinMax(TH1* h){
1043 
1044  Double_t xmin = 10000;
1045  Double_t xmax = -10000;
1046 
1047 
1048  for (int i = 1; i <= h->GetNbinsX(); ++i) {
1049  if ((h->GetBinContent(i) > 0)&&(h->GetBinCenter(i)>xmax)) xmax = h->GetBinCenter(i);
1050  }
1051  for (int i = h->GetNbinsX(); i >= 1; --i) {
1052  if ((h->GetBinContent(i) > 0)&&(h->GetBinCenter(i)<xmin)) xmin = h->GetBinCenter(i);
1053  }
1054 
1055  h->SetAxisRange((xmin-xmin*0.1), (xmax+xmax*0.1), "X");
1056  // std::cout << "xmin: " << xmin << ", xmax: " << xmax << std::endl;
1057 
1058 }
1059 
1060 
1061 
1062 void HIPplots::plotHitMap(char* outpath, int subDet, int minHits){
1063  cout << "Starting plotHitMap" << flush;
1064 
1065  TFile* falignable=new TFile(_inFile_HIPalign, "READ");
1066  cout << "\tLoaded file" << flush;
1067  //take the alignablewise tree and address useful branches
1068  TTree* talignable=(TTree*)falignable->Get("T2");
1069  cout << "\t Loaded tree" << endl;
1070 
1071  float eta=-999.0, phi=-55.0, xpos=-999.0, ypos=+999.0, zpos=-11111.0;
1072  int layer=-1, type=-1, nhit=-11111;
1073  int id=0;
1074  talignable->SetBranchAddress("Id", &id);
1075  talignable->SetBranchAddress("Layer", &layer);
1076  talignable->SetBranchAddress("Type", &type);
1077  talignable->SetBranchAddress("Nhit", &nhit);
1078  talignable->SetBranchAddress("Ypos", &ypos);
1079  talignable->SetBranchAddress("Eta", &eta);
1080  talignable->SetBranchAddress("Phi", &phi);
1081  talignable->SetBranchAddress("Ypos", &ypos);
1082  talignable->SetBranchAddress("Xpos", &xpos);
1083  talignable->SetBranchAddress("Zpos", &zpos);
1084 
1085  //loop on subdets
1086  char typetag[16];
1087  int maxLayers=0;
1088  if (subDet == TPBid){ sprintf(typetag, "TPB"); maxLayers=3; }
1089  else if (subDet == TPEid){ sprintf(typetag, "TPE"); maxLayers=2; }
1090  else if (subDet == TIBid){ sprintf(typetag, "TIB"); maxLayers=4; }
1091  else if (subDet == TIDid){ sprintf(typetag, "TID"); maxLayers=3; }
1092  else if (subDet == TOBid){ sprintf(typetag, "TOB"); maxLayers=6; }
1093  else if (subDet == TECid){ sprintf(typetag, "TEC"); maxLayers=9; }
1094  else { sprintf(typetag, "UNKNOWN"); }
1095  cout << "Starting to plot Hit Distributions for " << typetag << endl;
1096 
1097  bool printbinning=true;
1098  char psname[600];
1099  char ovtag[16];
1100  sprintf(psname, "%s/Hits_%s_Layers_Skimmed.eps", outpath, typetag);
1101 
1102  char binfilename[64];
1103  sprintf(binfilename, "./BinningHitMaps_%s.txt", typetag);
1104  ofstream binfile(binfilename, ios::out);
1105 
1106  if (printbinning){
1107  binfile << "******** Binning for Subdet " << typetag << "* ********" << endl << endl;
1108  }
1109 
1110  for (int layerindex=1; layerindex<=maxLayers; layerindex++){ //loop on layers
1111 
1112  cout << "\n\n*** Layer # " << layerindex << "* **" << endl;
1113  //activate only useful branches
1114  talignable->SetBranchStatus("*", 0);
1115  talignable->SetBranchStatus("Id", 1);
1116  talignable->SetBranchStatus("Type", 1);
1117  talignable->SetBranchStatus("Layer", 1);
1118  talignable->SetBranchStatus("Nhit", 1);
1119  talignable->SetBranchStatus("Eta", 1);
1120  talignable->SetBranchStatus("Phi", 1);
1121  talignable->SetBranchStatus("Ypos", 1);
1122  talignable->SetBranchStatus("Zpos", 1);
1123  talignable->SetBranchStatus("Xpos", 1);
1124 
1125  TCut selA, selECpos, selECneg;
1126  char selA_str[196], selECneg_str[196], selECpos_str[196];
1127  char varA_str[64], varB_str[64];
1128  char commonsense_str[128]="Xpos>-150.0&&Xpos<150.0&&Ypos>-150.0&&Ypos<150.0&&Zpos>-400.0&&Zpos<400.0";
1129  TCut commonsense=commonsense_str;
1130 
1131  sprintf(selECneg_str, "(Type==%d && Layer==%d && Zpos<0.0 && Nhit>=%d && sqrt(pow(Xpos,2)+pow(Ypos,2) )<125.0 )", subDet, layerindex, minHits);
1132  sprintf(selECpos_str, "(Type==%d && Layer==%d && Zpos>=0.0 && Nhit>=%d && sqrt(pow(Xpos,2)+pow(Ypos,2) )<125.0 )", subDet, layerindex, minHits);
1133  //sprintf(selB_str,"(Eta>0.0)&&(Eta<0.2)");
1134  sprintf(selA_str, "Type==%d && Layer==%d", subDet, layerindex);
1135 
1136  sprintf(varA_str, "Eta>>hvarx");
1137  sprintf(varB_str, "Phi>>hvary");
1138 
1139  selECneg=selECneg_str;
1140  selECpos=selECpos_str;
1141  selA=selA_str;
1142  cout << "Cuts defined as " << selA << endl;
1146 
1147  //--------- (2) START bin definition -----
1148  //cout << "Selection is " << sel << endl;
1149  //char sel[96];
1150 
1151  int nzentries= talignable->Draw("Zpos>>hZ(360,-270.0,270.0)", commonsense&&selA, "goff");
1152  TH1F* hZ=(TH1F*)gDirectory->Get("hZ");
1153  if (subDet == TOBid) SetPeakThreshold(8.0);
1154  else SetPeakThreshold(5.0);
1155  float Zpeaks[120];
1156  int bin=0;
1157  for (bin=0; bin<120; bin++){
1158  Zpeaks[bin]=-444.0;
1159  }
1160  const int nZpeaks=FindPeaks(hZ, Zpeaks, 99);
1161  const int nZbinlims=nZpeaks+1;
1162  float Zwidth=(Zpeaks[nZpeaks-1]-Zpeaks[0])/ (nZpeaks-1);
1163  float Zmin=Zpeaks[0]- Zwidth/2.0;
1164  float Zmax=Zpeaks[nZpeaks-1] + Zwidth/2.0;//((Zpeaks[nZbinlims-1]-Zpeaks[nZbinlims-2])/2.0) ;
1165  cout << "--> Zmin= " << Zmin << " - Zmax= " << Zmax << " Zwidth= " << Zwidth << " ; found " << nZpeaks << " Zpeaks" << endl;
1166  cout << "Zpeaks[0] is " << Zpeaks[0] << endl;
1167 
1168 
1169  float Phipeaks[120];
1170  if ((subDet==TIBid||subDet==TOBid)&&layerindex<=2) sprintf(selA_str, "%s&&Zpos>%f&&Zpos<%f", selA_str, Zpeaks[0]-2.0, Zpeaks[0]+2.0);//DS modules
1171  else sprintf(selA_str, "%s&&Zpos>%f&&Zpos<%f", selA_str, Zpeaks[0]-2.0, Zpeaks[0]+2.0);
1172  int nphientries=talignable->Draw("Phi>>hPhi", selA_str, "goff");
1173  cout << "N phi entries " << nphientries << " from sel " << selA_str << endl;
1174  TH1F* hPhi=(TH1F*)gDirectory->Get("hPhi");
1175  //nPhibins=FindPeaks(hPhi,Phipeaks,99);
1176  if (subDet==TPBid&&layerindex==1)nphientries=nphientries-1;
1177  const int nPhibins=nphientries;
1178  cout << "+ It would have found " << nPhibins << " phi bins" << endl;
1179 
1180  //fill Z array with binning. special case for DS layers of TIB with non equidistant z bins
1181  float phibin[nPhibins+1];
1182  float zbin[nZbinlims];
1183  float Phiwidth=6.28/nPhibins;
1184 
1185  if ((subDet==TIBid||subDet==TOBid)&&layerindex<=2){//DS modules in TIB and TOB
1186 
1188  cout << "Gonna loop over " << nZpeaks << " peaks / " << nZbinlims << " bin limits" << endl;
1189  for (bin=0; bin<nZbinlims-1; bin++){
1190  float zup=0.0;
1191  float zdown=0.0;
1192 
1193  if (bin==0){ //don't go underflow!
1194  zup=(Zpeaks[bin+1]-Zpeaks[bin])/2.0;
1195  zdown=zup;
1196  }
1197  else if (bin==nZbinlims-2){//don't go overflow!
1198  cout << "Don't go overflow !" << endl;
1199  zdown=(Zpeaks[bin]-Zpeaks[bin-1])/2.0;
1200  if (layerindex==1) zup=(Zpeaks[bin-1]-Zpeaks[bin-2])/2.0;
1201  else zup=zdown;
1202  }
1203  else{
1204  zup=(Zpeaks[bin+1]-Zpeaks[bin])/2.0;
1205  zdown=(Zpeaks[bin]-Zpeaks[bin-1])/2.0;
1206  }
1207  zbin[bin] = Zpeaks[bin]-zdown;
1208  zbin[bin+1]= Zpeaks[bin]+zup;
1209 
1210  }//end loop on z bin
1212 
1213 
1214  for (int bin=0; bin<=nPhibins; ++bin){
1215  if (bin==0)phibin[bin]=-3.14+bin*(Phiwidth)+Phiwidth/4;
1216  else if (bin==nPhibins-1)phibin[bin]=-3.14+bin*(Phiwidth)-Phiwidth/4;
1217  //else phibin[bin]=-3.14+bin*(Phiwidth);
1218  else phibin[bin]=phibin[bin-1]+ Phiwidth;
1219  }//end loop on phi bin
1220 
1221  }//end if DS layers of TIB/TOB
1222  else{
1223  for (int bin=0; bin<nZbinlims; ++bin){
1224  zbin[bin]=Zmin+bin*(Zwidth);
1225  }//end loop on z bin
1226 
1227  for (int bin=0; bin<=nPhibins; ++bin){
1228  phibin[bin]=-3.14+bin*(Phiwidth);
1229  }//end loop on phi bin
1230  }
1231 
1232 
1233  float Phimin=Phipeaks[0]- ((Phipeaks[1]-Phipeaks[0])/2.0);
1234  float Phimax=Phipeaks[nPhibins-1] + ((Phipeaks[nPhibins-1]-Phipeaks[nPhibins-2])/2.0);
1235 
1236  cout << "N Z bin LIMITs = " << nZbinlims << " N Phi bin LIMITS = " << nPhibins << endl;
1237  //--------- END bin definition -----
1238 
1239 
1240 
1241  char histoname[64];
1242  sprintf(histoname, "%s_Layer%d", typetag, layerindex);
1243  TH2F* hetaphi=new TH2F(histoname, histoname, nZpeaks, zbin, nPhibins, phibin);
1244 
1245 
1246  cout << "Starting to loop on entries" << flush;
1247  int nmods=0;
1248  int nlowentrycells=0;
1249 
1250  for (int j=0; j<talignable->GetEntries(); j++){ //loop on entries (-> alignables)
1251  if (j%1000==0)cout << "." << flush;
1252  talignable->GetEntry(j);
1253  if (type==subDet&&layer==layerindex){
1254  if (nhit>=minHits){
1255  //find the right eta bin
1256  hetaphi->Fill(zpos, phi, nhit);
1257  nmods++;
1258  }
1259  else{
1260  hetaphi->Fill(zpos, phi, -99);
1261  nlowentrycells++;
1262  }//end else if nhits < minhits
1263 
1264  }//end if the type and layer are the desired ones
1265  }//end loop on entries(->alignables)
1266 
1267  //Let's start with the funny things: fancy graphics!
1268  hetaphi->SetXTitle("Z [cm]");
1269  hetaphi->SetYTitle("#phi (rad)");
1270 
1271  int Nxbins=hetaphi->GetXaxis()->GetNbins();
1272  int Nybins=hetaphi->GetYaxis()->GetNbins();
1273  cout << "On x-axis there are " << Nxbins << " bins " << endl;
1274  cout << "On y-axis there are " << Nybins << " bins " << endl;
1275 
1276 
1277  bool smooth_etaphi=false;
1278  if (smooth_etaphi){
1279  for (int i=1; i<=Nxbins; i++){
1280  for (int j=1; j<=Nybins; j++){
1281  float bincont=hetaphi->GetBinContent(i, j);
1282  if (bincont==0){//average with the surrounding bins
1283  float binup1=hetaphi->GetBinContent(i, j+1);
1284  float bindown1=hetaphi->GetBinContent(i, j-1);
1285  float binlx1=hetaphi->GetBinContent(i-1, j);
1286  float binrx1=hetaphi->GetBinContent(i+1, j);
1287  if (i==1)binlx1=binrx1;//at the edges you don't have lx or rx bins; set a def value
1288  if (i==Nxbins)binrx1=binlx1;
1289  if (j==1)bindown1=binup1;
1290  if (j==Nybins)binup1=bindown1;
1291  int adiacentbins=0;
1292  if (binup1>0.0)adiacentbins++;
1293  if (bindown1>0.0)adiacentbins++;
1294  if (binlx1>0.0)adiacentbins++;
1295  if (binrx1>0.0)adiacentbins++;
1296  if (adiacentbins>=2){
1297  float avg=(binup1+bindown1+binlx1+binrx1)/adiacentbins;
1298  hetaphi->SetBinContent(i, j, avg);
1299  }
1300  }
1301  }
1302  }
1303  }//end if smooth_etaphi
1304 
1305 
1306  //for debugging purposes
1307  TGraph* gretaphi;
1308  bool plotAlignablePos=false;
1309  if (plotAlignablePos){
1310  const int ngrpoints=nmods;
1311  float etagr[ngrpoints], phigr[ngrpoints];
1312  nmods=0;
1313  for (int j=0; j<talignable->GetEntries(); j++){ //loop on entries (-> alignables)
1314  if (j%1000==0)cout << "." << flush;
1315  talignable->GetEntry(j);
1316  if (type==subDet&&layer==layerindex){
1317  etagr[nmods]=zpos;
1318  phigr[nmods]=phi;
1319  nmods++;
1320  }
1321  }
1322 
1323  gretaphi=new TGraph(ngrpoints, etagr, phigr);
1324  gretaphi->SetMarkerStyle(20);
1325  gretaphi->SetMarkerColor(1);
1326  gretaphi->SetMarkerSize(0.75);
1327  }
1328 
1329 
1330 
1332  //cout << "Layer #" << i << endl;
1333  float Zcellgr[512];
1334  float Phicellgr[512];
1335  int nemptycells=0;
1336  for (int zcells=1; zcells<=hetaphi->GetNbinsX(); zcells++){
1337  for (int phicells=1; phicells<=hetaphi->GetNbinsY(); phicells++){
1338  if (hetaphi->GetBinContent(zcells, phicells)==-99){
1339  hetaphi->SetBinContent(zcells, phicells, 0);
1340  Zcellgr[nemptycells]= float(hetaphi->GetXaxis()->GetBinCenter(zcells));
1341  Phicellgr[nemptycells]= float(hetaphi->GetYaxis()->GetBinCenter(phicells));
1342  nemptycells++;
1343  }
1344  //if(a==2)cout << "Finished Z value " << hetaphi->GetXaxis()->GetBinCenter(zcells) << " the emptycells are " << nemptycells << endl;
1345  }//end loop on Phi bins
1346  }//end loop on Z bins
1347  TGraph* gr_empty=new TGraph(nlowentrycells, Zcellgr, Phicellgr);
1348  sprintf(histoname, "gr_emptycells_subdet%d_layer%d", subDet, layerindex);
1349  //cout << "Graph name: " << histoname << " with " << gr_empty->GetN() << endl;
1350  gr_empty->SetName(histoname);
1351  gr_empty->SetMarkerStyle(5);
1353 
1354 
1355  cout << " Done! Used " << nmods << " alignables. Starting to plot !" << endl;
1356 
1357  //plot them and save the canvas appending it to a ps file
1358  gStyle->SetPalette(1, 0);
1359  TCanvas* c_barrels=new TCanvas("canvas_hits_barrels", "CANVAS_HITS_BARRELS", 1600, 1600);
1360  TCanvas* c_endcaps=new TCanvas("canvas_hits_endcaps", "CANVAS_HITS_ENDCAPS", 3200, 1600);
1361  TPad* pleft=new TPad("left_panel", "Left Panel", 0.0, 0.0, 0.49, 0.99);
1362  TPad* pcent=new TPad("central_up_panel", "Central Panel", 0.01, 0.00, 0.99, 0.99);
1363  TPad* pright=new TPad("right_panel", "Right Panel", 0.51, 0.0, 0.99, 0.99);
1364 
1365 
1366  if (subDet==1 ||subDet==3 ||subDet==5){
1367  c_barrels->cd();
1368  pcent->Draw();
1369  pcent->cd();
1370  gPad->SetRightMargin(0.15);
1371 
1372  //hbarrel->SetMarkerSize(2.0);
1373  //hbarrel->SetMarkerSize(21);
1374  hetaphi->SetStats(0);
1375  if (subDet==TPBid||subDet==TIBid||subDet==TOBid)hetaphi->Draw("COLZtext");//COLZtext
1376  if (plotAlignablePos) gretaphi->Draw("P");
1377  gr_empty->Draw("P");
1378 
1379  if (printbinning){
1380  binfile << "--> Layer #" << layerindex << endl;
1381  binfile << "Eta Binning: " << flush;
1382  for (int h=1; h<=hetaphi->GetNbinsX(); h++){
1383  binfile << hetaphi->GetXaxis()->GetBinLowEdge(h) << "\t";
1384  if (h==hetaphi->GetNbinsX()) binfile << hetaphi->GetXaxis()->GetBinLowEdge(h)+hetaphi->GetXaxis()->GetBinWidth(h);
1385  }
1386  binfile << endl;
1387  binfile << "Phi Binning: " << flush;
1388  for (int h=1; h<=hetaphi->GetNbinsY(); h++){
1389  binfile << hetaphi->GetYaxis()->GetBinLowEdge(h) << "\t";
1390  if (h==hetaphi->GetNbinsX()) binfile << hetaphi->GetYaxis()->GetBinLowEdge(h)+hetaphi->GetYaxis()->GetBinWidth(h);
1391  }
1392  binfile << endl;
1393  }
1394 
1395  }
1396  else{
1397  c_endcaps->cd();
1398  pleft->Draw();
1399  pright->Draw();
1400 
1401 
1402  pleft->cd();
1403  gPad->SetRightMargin(0.15);
1404  char selEC_str[192], varEC_str[128];
1405  float radlimit=0.0;
1406  if (subDet == TPBid){ radlimit=45.0; }
1407  else if (subDet == TPEid){ radlimit=45.0; }
1408  else if (subDet == TIBid){ radlimit=70.0; }
1409  else if (subDet == TIDid){ radlimit=70.0; }
1410  else if (subDet == TOBid){ radlimit=130.0; }
1411  else if (subDet == TECid){ radlimit=130.0; }
1412  else { radlimit=0.0; }
1413 
1414  sprintf(varEC_str, "Ypos:Xpos>>hxy_negz(30,%f,%f)", -radlimit, radlimit);
1415  sprintf(selEC_str, "Nhit*(%s&&%s)", selECneg_str, commonsense_str);
1416  cout << selEC_str << endl;
1417  int selentriesECneg=talignable->Draw(varEC_str, selEC_str, "goff"); //total number of alignables for thys type and layer
1418  if (selentriesECneg>0){
1419  TH2F* hxy_negz=(TH2F*)gDirectory->Get("hxy_negz");
1420  hxy_negz->GetXaxis()->SetRangeUser(-radlimit, radlimit);
1421  hxy_negz->GetYaxis()->SetRangeUser(-radlimit, radlimit);
1422  char histoname_negz[32];
1423  sprintf(histoname_negz, "%s (-Z)", histoname);
1424  hxy_negz->SetStats(0);
1425  hxy_negz->SetTitle(histoname_negz);
1426  hxy_negz->SetXTitle("X (cm)");
1427  hxy_negz->SetYTitle("Y (cm)");
1428  hxy_negz->Draw("COLZ");
1429  }
1430  else{
1431  cout << "WARNING !!!! No hits on this layer ! not plotting (-Z) !" << endl;
1432  }
1433 
1434 
1435 
1436  cout << "PAD 3" << endl;
1437  pright->cd();
1438  gPad->SetRightMargin(0.15);
1439  sprintf(selEC_str, "Nhit*(%s&&%s)", selECpos_str, commonsense_str);
1440  //selEC=selEC_str&&commonsense;
1441  cout << "(2)" << selEC_str << endl;
1442  sprintf(varEC_str, "Ypos:Xpos>>hxy_posz(30,%f,%f)", -radlimit, radlimit);
1443  int selentriesECpos=talignable->Draw(varEC_str, selEC_str, "goff"); //total number of alignables for thys type and layer
1444  if (selentriesECpos>0){
1445  TH2F* hxy_posz=(TH2F*)gDirectory->Get("hxy_posz");
1446  char histoname_posz[32];
1447  hxy_posz->GetXaxis()->SetRangeUser(-radlimit, radlimit);
1448  hxy_posz->GetYaxis()->SetRangeUser(-radlimit, radlimit);
1449  sprintf(histoname_posz, "%s (+Z)", histoname);
1450  hxy_posz->SetStats(0);
1451  hxy_posz->SetTitle(histoname_posz);
1452  hxy_posz->SetXTitle("X (cm)");
1453  hxy_posz->SetYTitle("Y (cm)");
1454  hxy_posz->Draw("COLZ");
1455  }
1456  else{
1457  cout << "WARNING !!!! No hits on this layer ! not plotting (+Z) !" << endl;
1458  }
1459 
1460 
1461  }
1462 
1463  //save in a ps file
1464  cout << "******* " << typetag << endl;
1465  char psnamefinal[600];
1466 
1467  if (layerindex==1) sprintf(psnamefinal, "%s(", psname);
1468  else if (layerindex==maxLayers) sprintf(psnamefinal, "%s)", psname);
1469  else sprintf(psnamefinal, "%s", psname);
1470 
1471  cout << "Saving in " << psnamefinal << endl;
1472  if (subDet==1 ||subDet==3 ||subDet==5)c_barrels->SaveAs(psnamefinal);
1473  else c_endcaps->SaveAs(psnamefinal);
1474 
1475  if (subDet==1 ||subDet==3 ||subDet==5) delete c_barrels;
1476  else delete c_endcaps;
1477 
1478 
1479  //ctmp->cd();
1480  //hbarrel->Draw("scat");
1481 
1482  }//end loop on layers
1483  cout << "Finished " << maxLayers << " of the " << typetag << endl;
1484 
1485  delete talignable;
1486  delete falignable;
1487 }//end PlotHitDistribution
1488 
1489 
1490 void HIPplots::dumpAlignedModules(int nhits){
1491 
1493 
1494 }//end dumpAlignedModules
1495 
1496 
1497 
1498 int HIPplots::FindPeaks(TH1F* h1, float* peaklist, const int maxNpeaks, int startbin, int endbin){
1499 
1500  int Npeaks=0;
1501  if (startbin<0)startbin=1;
1502  if (endbin<0)endbin=h1->GetNbinsX();
1503  int bin=startbin;
1504  int prevbin=startbin;
1505  bool rising=true;
1506  while (bin<=endbin){
1507 
1508  float bincont=h1->GetBinContent(bin);
1509  float prevbincont=h1->GetBinContent(prevbin);
1510 
1511  if (prevbincont>peakthreshold||bincont>1.0){
1512  if (bincont>=prevbincont){//we are rising, keep going until we don't go down
1513  rising=true;
1514  }
1515  else{//ehi! we are decreasing. But check if we were decreasing already at the step before
1516  if (!rising){//we were already decreasing, this is not a maximum
1517  rising=true;
1518  }
1519  else{//we found a maximum
1520  rising=false;
1521  peaklist[Npeaks]=h1->GetBinCenter(prevbin);
1522  Npeaks++;
1523  }
1524  }
1525  }//end if bincont>1.0
1526  if (Npeaks>=maxNpeaks){//max number reached
1527  bin=endbin;
1528  }
1529  bin++;
1530  prevbin=bin-1;
1531  }//end while loop on bins
1532 
1533 
1534  return Npeaks;
1535 }//End FindPeaks
1536 
1537 
1538 void HIPplots::SetPeakThreshold(float newpeakthreshold){
1539  peakthreshold=newpeakthreshold;
1540 }
1541 
1542 bool HIPplots::CheckFileExistence(TString filename){
1543  bool flag = false;
1544  fstream fin;
1545  fin.open(filename.Data(), ios::in);
1546  if (fin.is_open()) flag=true;
1547  fin.close();
1548  return flag;
1549 }
1550 
1551 void HIPplots::CheckFiles(int &ierr){
1552 
1553  if (!CheckFileExistence(_inFile_uservars)){
1554  cout << "Missing file " << _inFile_uservars << endl;
1555  ierr++;
1556  }
1557 
1558  if (!CheckFileExistence(_inFile_HIPalign)){
1559  cout << "Missing file " << _inFile_HIPalign << endl;
1560  ierr++;
1561  }
1562 
1563  if (!CheckFileExistence(_inFile_alipos)){
1564  cout << "Missing file " << _inFile_alipos << endl;
1565  ierr++;
1566  }
1567 
1568  if (CheckFileExistence(_outFile)){
1569  cout << "Output file already exists!" << endl;
1570  ierr=-1*(ierr+1);
1571  }
1572 
1573 }
1574 
1575 
1576 bool HIPplots::CheckHistoRising(TH1D* h){
1577 
1578  bool rise1=false, rise2=false, rise3=false;
1579  int totbins=h->GetNbinsX();
1580  //for(int i=1;i<=totbins;i++){
1581  int i=1;
1582 
1583 
1584  for (i=1; i<=totbins-3; i++){
1585  double cont1 = h->GetBinContent(i);
1586  if (h->GetBinContent(i+1)>cont1)rise1=true;
1587  if (rise1){
1588  if (h->GetBinContent(i+2)>h->GetBinContent(i+1))rise2=true;
1589  if (rise2){
1590  if (h->GetBinContent(i+3)>h->GetBinContent(i+2)){
1591  rise3=true;
1592  break;
1593  }
1594  }
1595  }
1596 
1597  }//emnd loop on bins
1598 
1599  return rise3;
1600 
1601 }//end CheckHistoRising
1602 
HIPplots::TPBid
Definition: HIPplots.h:42
histoParameters_cff.maxLayers
maxLayers
Definition: histoParameters_cff.py:65
HIPplots::_inFile_uservars
TString _inFile_uservars
Definition: HIPplots.h:36
HIPplots::_inFile_truepos
TString _inFile_truepos
Definition: HIPplots.h:37
mps_fire.i
i
Definition: mps_fire.py:428
HIPplots::TPEid
Definition: HIPplots.h:42
dqmMemoryStats.float
float
Definition: dqmMemoryStats.py:127
Phase2TrackerMonitorDigi_cff.Nxbins
Nxbins
Definition: Phase2TrackerMonitorDigi_cff.py:20
colors
Definition: colors.py:1
COLOR_CODE
#define COLOR_CODE(icolor)
HIPplots::GetSubDet
int GetSubDet(unsigned int id)
Definition: HIPplots.cc:1016
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
HIPplots::_outFile
TString _outFile
Definition: HIPplots.h:34
HIPplots::_inFile_mispos
TString _inFile_mispos
Definition: HIPplots.h:39
testProducerWithPsetDescEmpty_cfi.x2
x2
Definition: testProducerWithPsetDescEmpty_cfi.py:28
HIPplots::SetPeakThreshold
void SetPeakThreshold(float newpeakthreshold)
Definition: HIPplots.cc:1538
cuy.col
col
Definition: cuy.py:1010
gather_cfg.cout
cout
Definition: gather_cfg.py:144
HIPplots::GetBarrelLayer
int GetBarrelLayer(unsigned int id)
Definition: HIPplots.cc:1033
HIPplots::plotAlignParams
void plotAlignParams(string ShiftsOrParams, char *plotName="test.png")
Definition: HIPplots.cc:394
HIPplots::plotAlignableChiSquare
void plotAlignableChiSquare(char *plotName="testchi2.png", float minChi2n=-1.0)
Definition: HIPplots.cc:844
HIPplots::_inFile_surveys
TString _inFile_surveys
Definition: HIPplots.h:41
HIPplots::TECid
Definition: HIPplots.h:42
testProducerWithPsetDescEmpty_cfi.x1
x1
Definition: testProducerWithPsetDescEmpty_cfi.py:33
testProducerWithPsetDescEmpty_cfi.y1
y1
Definition: testProducerWithPsetDescEmpty_cfi.py:29
HIPplots::CheckHistoRising
bool CheckHistoRising(TH1D *h)
Definition: HIPplots.cc:1576
Gflash::Zmax
const double Zmax[kNumberCalorimeter]
Definition: GflashNameSpace.h:29
groupFilesInBlocks.fin
fin
Definition: groupFilesInBlocks.py:94
GlobalPosition_Frontier_DevDB_cff.tag
tag
Definition: GlobalPosition_Frontier_DevDB_cff.py:11
HIPplots::dumpAlignedModules
void dumpAlignedModules(int nhits=0)
Definition: HIPplots.cc:1490
HIPplots::extractAlignableChiSquare
void extractAlignableChiSquare(int minHits=0, int subDet=0, int doubleSided=0)
Definition: HIPplots.cc:636
HIPplots::CheckFileExistence
bool CheckFileExistence(TString filename)
Definition: HIPplots.cc:1542
PVValHelper::eta
Definition: PVValidationHelpers.h:69
HIPplots::_inFile_HIPalign
TString _inFile_HIPalign
Definition: HIPplots.h:40
L1TEGammaDiff_cfi.plotName
plotName
Definition: L1TEGammaDiff_cfi.py:21
corrVsCorr.filename
filename
Definition: corrVsCorr.py:123
h
HIPplots::CheckFiles
void CheckFiles(int &ierr)
Definition: HIPplots.cc:1551
nhits
Definition: HIMultiTrackSelector.h:42
L1TdeCSCTF_cfi.outFile
outFile
Definition: L1TdeCSCTF_cfi.py:5
HIPplots::_path
TString _path
Definition: HIPplots.h:33
testProducerWithPsetDescEmpty_cfi.y2
y2
Definition: testProducerWithPsetDescEmpty_cfi.py:30
HIPplots::HIPplots
HIPplots(int IOV, char *path, char *outFile)
Definition: HIPplots.cc:20
a
double a
Definition: hdecay.h:119
type
type
Definition: SiPixelVCal_PayloadInspector.cc:37
recoMuon::in
Definition: RecoMuonEnumerators.h:6
listHistos.legend
legend
Definition: listHistos.py:41
groupFilesInBlocks.fout
fout
Definition: groupFilesInBlocks.py:162
HIPplots::SetMinMax
void SetMinMax(TH1 *h)
Definition: HIPplots.cc:1042
createfilelist.int
int
Definition: createfilelist.py:10
HIPplots::FindPeaks
int FindPeaks(TH1F *h1, float *peaklist, const int maxNpeaks, int startbin=-1, int endbin=-1)
Definition: HIPplots.cc:1498
Gflash::Zmin
const double Zmin[kNumberCalorimeter]
Definition: GflashNameSpace.h:28
HIPplots::extractSurveyResiduals
void extractSurveyResiduals(int currentPar, int subDet=0)
Definition: HIPplots.cc:754
HIPAlignmentAlgorithm_cfi.outpath
outpath
Definition: HIPAlignmentAlgorithm_cfi.py:75
cmsLHEtoEOSManager.l
l
Definition: cmsLHEtoEOSManager.py:204
newFWLiteAna.bin
bin
Definition: newFWLiteAna.py:161
DDAxes::phi
HIPplots::GetNIterations
int GetNIterations(TDirectory *f, char *tag, int startingcounter=0)
Definition: HIPplots.cc:1001
HIPplots::_inFile_alipos
TString _inFile_alipos
Definition: HIPplots.h:38
HIPplots::plotbadchi2
bool plotbadchi2
Definition: HIPplots.h:54
HIPplots::peakthreshold
float peakthreshold
Definition: HIPplots.h:53
edm::shift
static unsigned const int shift
Definition: LuminosityBlockID.cc:7
HIPplots::plotAlignParamsAtIter
void plotAlignParamsAtIter(int iter, string ShiftsOrParams, char *plotName="test.png")
Definition: HIPplots.cc:539
triggerObjects_cff.id
id
Definition: triggerObjects_cff.py:31
HIPplots::TOBid
Definition: HIPplots.h:42
TrackerOfflineValidation_Dqm_cff.xmax
xmax
Definition: TrackerOfflineValidation_Dqm_cff.py:11
HIPplots::TIBid
Definition: HIPplots.h:42
Phase2TrackerMonitorDigi_cff.Nybins
Nybins
Definition: Phase2TrackerMonitorDigi_cff.py:23
HIPplots::_inFile_params
TString _inFile_params
Definition: HIPplots.h:35
benchmark_cfg.fa
fa
Definition: benchmark_cfg.py:13
MillePedeFileConverter_cfg.out
out
Definition: MillePedeFileConverter_cfg.py:31
HIPplots::extractAlignShifts
void extractAlignShifts(int i, int minHits=0, int subDet=0)
Definition: HIPplots.cc:199
dataAnalyzerFineBiningParameters_cff.Phimax
Phimax
Definition: dataAnalyzerFineBiningParameters_cff.py:7
ztail.d
d
Definition: ztail.py:151
castor_dqm_sourceclient_file_cfg.path
path
Definition: castor_dqm_sourceclient_file_cfg.py:37
TrackerOfflineValidation_Dqm_cff.xmin
xmin
Definition: TrackerOfflineValidation_Dqm_cff.py:10
dqmiolumiharvest.j
j
Definition: dqmiolumiharvest.py:66
dataAnalyzerFineBiningParameters_cff.Phimin
Phimin
Definition: dataAnalyzerFineBiningParameters_cff.py:8
HIPplots::extractAlignParams
void extractAlignParams(int i, int minHits=0, int subDet=0, int doubleSided=0)
Definition: HIPplots.cc:81
HIPplots.h
TtFullHadEvtBuilder_cfi.prob
prob
Definition: TtFullHadEvtBuilder_cfi.py:33
HIPplots::MakeLegend
TLegend * MakeLegend(double x1=0.1, double y1=0.1, double x2=0.1, double y2=0.1)
Definition: HIPplots.cc:46
HIPplots::plotHitMap
void plotHitMap(char *outpath, int subDet, int minHits=0)
Definition: HIPplots.cc:1062
reco_skim_cfg_mod.minHits
minHits
Definition: reco_skim_cfg_mod.py:177
g
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e g
Definition: Activities.doc:4
RemoveAddSevLevel.flag
flag
Definition: RemoveAddSevLevel.py:116
HIPplots::TIDid
Definition: HIPplots.h:42
Generated for CMSSW Reference Manual by doxygen 1.8.16