HIPAlignmentAlgorithm/plots/HIPplots.cc

0001 #include "TROOT.h"
0002 #include "TAttFill.h"
0003 #include "TColor.h"
0004 #include "HIPplots.h"
0005 #include <string>
0006 #include <sstream>
0007 #include <cmath>
0008
0009 #include "TProfile.h"
0010 #include "TPaveStats.h"
0011 #include "TList.h"
0012 #include "TNtuple.h"
0013 #include "TString.h"
0014 #include "TTree.h"
0015
0016 #include "TMatrixD.h"
0017 #include "TVectorD.h"
0018 #include "TLegend.h"
0019
0020 HIPplots::HIPplots(int IOV, char* path, char* outFile):
0021 _IOV(IOV),
0022 _path(path),
0023 _outFile(outFile)
0024 {
0025   //Id,ObjId,Nhit,Jtvj,Jtve,AlignableChi2,AlignableNdof
0026   _inFile_uservars = Form("%s/IOUserVariables.root", _path.Data());
0027   if (!CheckFileExistence(_inFile_uservars)) _inFile_uservars = Form("%s/IOUserVariables_%d.root", _path.Data(), IOV);
0028
0029   //aligned absolute positions
0030   _inFile_alipos = Form("%s/IOAlignedPositions.root", _path.Data());
0031   if (!CheckFileExistence(_inFile_alipos)) _inFile_alipos = Form("%s/IOAlignedPositions_%d.root", _path.Data(), IOV);
0032
0033   _inFile_HIPalign = Form("%s/HIPAlignmentAlignables.root", _path.Data());
0034   if (!CheckFileExistence(_inFile_HIPalign)) _inFile_HIPalign = Form("%s/HIPAlignmentAlignables_%d.root", _path.Data(), IOV);
0035
0036   _inFile_params = Form("%s/IOAlignmentParameters.root", _path.Data()); //Id (raw id), ObjId, Par (alignable movements), CovariantMatrix
0037   _inFile_mispos = Form("%s/IOMisalignedPositions.root", _path.Data());
0038   _inFile_surveys = Form("%s/HIPSurveyResiduals.root", _path.Data());
0039   _inFile_truepos = Form("%s/IOTruePositions.root", _path.Data()); //redundant
0040
0041   SetPeakThreshold(8.0);
0042   plotbadchi2=true;
0043 }
0044
0045
0046 TLegend*  HIPplots::MakeLegend(double x1,
0047   double y1,
0048   double x2,
0049   double y2)
0050 {
0051   TLegend*  legend = new TLegend(x1, y1, x2, y2, "", "NBNDC");
0052   legend->SetNColumns(6);
0053   legend->SetFillColor(0);
0054   legend->SetBorderSize(0);
0055   int COLOR_CODE[6]={ 28, 2, 3, 4, 6, 7 };
0056
0057   // TO BE UPDATED FOR PHASE-1
0058   TString detNames[6];
0059   detNames[0] = TString("PXB");
0060   detNames[1] = TString("PXF");
0061   detNames[2] = TString("TIB");
0062   detNames[3] = TString("TID");
0063   detNames[4] = TString("TOB");
0064   detNames[5] = TString("TEC");
0065
0066   for (unsigned int isublevel = 0; isublevel < 6; isublevel++)
0067   {
0068     TGraph*  g = new TGraph(0);
0069     g->SetLineColor(COLOR_CODE[isublevel]);
0070     g->SetMarkerColor(COLOR_CODE[isublevel]);
0071     g->SetFillColor(COLOR_CODE[isublevel]);
0072     g->SetFillColor(1);
0073     g->SetMarkerStyle(8);
0074     g->SetMarkerSize(1);
0075     legend->AddEntry(g, detNames[isublevel], "lp");
0076   }
0077   return legend;
0078 }
0079
0080
0081 void HIPplots::extractAlignParams(int currentPar, int minHits, int subDet, int doubleSided){
0082
0083   cout << "--- extracting AlignParams ; Par " << currentPar << " ---" << endl << endl; //0-5: u,v,w,alpha,beta,gamma
0084   //void ExtractAlignPars(int pp, string SubDetString){
0085   //    const int par_now=pp;
0086   //load first tree
0087
0088   //    TFile* fp = new TFile( _inFile_params,"READ");
0089   //    const TList* keysp = fp->GetListOfKeys();
0090   //    const unsigned int maxIteration = keysp->GetSize() - 1;
0091   cout << "Loaded par file -. OK" << endl;
0092   TFile* fv = new TFile(_inFile_uservars, "READ");
0093   const TList* keysv = fv->GetListOfKeys();
0094   const unsigned int maxIteration = keysv->GetSize() - 1;
0095
0096   char fileaction[16];
0097   if (CheckFileExistence(_outFile))sprintf(fileaction, "UPDATE");
0098   else sprintf(fileaction, "NEW");
0099
0100   TFile* fout = new TFile(_outFile, fileaction);
0101
0102   TTree* tree0 = (TTree*)fv->Get(keysv->At(0)->GetName());
0103   //const char* tree0v_Name = keysv->At(0)->GetName();
0104   //tree0->AddFriend( tree0v_Name, fv );
0105
0106   int nHit;
0107   unsigned int detId;
0108   double par[6];
0109   unsigned int detId0;
0110
0111   tree0->SetBranchAddress("Id", &detId0);
0112
0113   const int ndets=tree0->GetEntries();
0114   TH1D* hpar1[ndets];
0115   char ppdirname[16];
0116   sprintf(ppdirname, "ShiftsPar%d", currentPar);
0117   fout->mkdir(ppdirname);
0118   gDirectory->cd(ppdirname);
0119   //delare histos
0120
0121
0122   TH1D* hpariter[maxIteration];
0123   for (unsigned int a = 0; a < maxIteration; a++){
0124     char histoname[32], histotitle[32];
0125     sprintf(histoname, "Par_%d_Iter_%d", currentPar, a);
0126     sprintf(histotitle, "Par %d for iteration #%d", currentPar, a);
0127     hpariter[a]=new TH1D(histoname, histoname, 1000, -1.0, 1.0);
0128   }
0129
0130   for (int i = 0; i < ndets; i++){
0131     char histoname[32], histotitle[32];
0132     sprintf(histoname, "Par_%d_SiDet_%d", currentPar, i);
0133     sprintf(histotitle, "Par %d for detector #%d", currentPar, i);
0134     hpar1[i]=new TH1D(histoname, histoname, maxIteration+1, -0.5, maxIteration+0.5);
0135   }
0136
0137   //file where to save aligned det id list
0138   ofstream flist("./List_aligned_dets.txt", ios::out);
0139
0140
0141   //Loop on trees
0142   int modules_accepted = 0;
0143   for (unsigned int iter = 0; iter <= maxIteration; iter++){//loop on i (HIP iterations -> trees in file)
0144
0145     //      TTree* tmpTree = (TTree*) fp->Get(keysp->At(iter)->GetName() );
0146     //      const char* tmpTreev = keysv->At(iter)->GetName();
0147     TTree* tmpTree = (TTree*)fv->Get(keysv->At(iter)->GetName());
0148     //      tmpTree->AddFriend( tmpTreev, fv ); //get access to both IOAlignmentParameters and IOUserVariables
0149     tmpTree->SetBranchAddress("Id", &detId);
0150     tmpTree->SetBranchAddress("Nhit", &nHit); //taken from IOUserVariables
0151     tmpTree->SetBranchAddress("Par", &par); //taken from IOAlignmentParameters
0152
0153     std::cout << "iteration: " << iter << "..." << std::endl;
0154
0155     modules_accepted=0;
0156     //loop on entries
0157     for (int j = 0; j < tmpTree->GetEntries(); j++){ //loop on j (modules -> entries in each tree)
0158
0159       tmpTree->GetEntry(j);
0160       bool passSubdetCut = true;
0161       if (subDet > 0){ if (GetSubDet(detId) != subDet) passSubdetCut = false; } //find the subDet from module detId
0162       if (doubleSided > 0){
0163         if (GetSubDet(detId)%2 == 1 && doubleSided == 1 && GetBarrelLayer(detId) < 3) passSubdetCut = false;
0164         if (GetSubDet(detId)%2 == 1 && doubleSided == 2 && GetBarrelLayer(detId) > 2) passSubdetCut = false;
0165       }
0166
0167       if ((nHit >= minHits)&&(passSubdetCut)){  //cut on min nhits per module
0168         hpar1[j]->SetBinContent(iter+1, par[currentPar]); //x-iteration, y-movement in this iteration
0169         //std::cout << "iteration: " << iter << ",par" << currentPar << "= " << par[currentPar] << std::endl;
0170
0171         int COLOR_CODE[6]={ 28, 2, 3, 4, 6, 7 };
0172         hpar1[j]->SetLineColor(COLOR_CODE[GetSubDet(detId)-1]);
0173         hpar1[j]->SetMarkerColor(COLOR_CODE[GetSubDet(detId)-1]);
0174
0175         //hpariter[iter]->Fill(par[currentPar]);
0176         //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;
0177
0178         if ((iter==maxIteration-1)&&(currentPar==0))flist << detId << endl;
0179         modules_accepted++;
0180
0181       }
0182     }//end loop on entries (i.e. modules)
0183     delete tmpTree;
0184   }//end loop on iteration
0185   std::cout << "Modules accepted: " << modules_accepted << std::endl;
0186   std::cout << "Writing..." << std::endl;
0187   //Save
0188   fout->Write();
0189   flist.close();
0190   //delete hpar1;
0191   std::cout << "Deleting..." << std::endl;
0192   delete fout;
0193   //    std::cout << "Deleting..." << std::endl;
0194   //    delete fp;
0195   //std::cout << "Deleting..." << std::endl;
0196   delete fv;
0197 }
0198
0199 void HIPplots::extractAlignShifts(int currentPar, int minHits, int subDet){
0200
0201   cout << "\n--- extracting AlignShifts ; Par " << currentPar << " ---" << endl << endl;
0202   TFile* fa = new TFile(_inFile_alipos, "READ");
0203   const TList* keysa = fa->GetListOfKeys();
0204   const unsigned int maxIteration = keysa->GetSize();
0205
0206   TFile* fv = new TFile(_inFile_uservars, "READ");
0207   const TList* keysv = fv->GetListOfKeys();
0208
0209   //    TFile* fm = new TFile( _inFile_mispos,"READ");
0210   //TFile* ft = new TFile( _inFile_truepos,"READ");
0211
0212   char fileaction[16];
0213   if (CheckFileExistence(_outFile))sprintf(fileaction, "UPDATE");
0214   else sprintf(fileaction, "NEW");
0215   TFile* fout = new TFile(_outFile, fileaction);
0216
0217   TTree* tree0 = (TTree*)fa->Get(keysa->At(0)->GetName());
0218   const char* tree0v_Name = keysv->At(0)->GetName();
0219   tree0->AddFriend(tree0v_Name, fv);  //get access to both IOAlignedPositions and IOUserVariables
0220
0221   unsigned int detId0;
0222
0223   tree0->SetBranchAddress("Id", &detId0);
0224
0225   const int ndets=tree0->GetEntries();
0226   TH1D* hshift[ndets]; //hshift[i] absolute position change for det i, at all iterations
0227   //    TH1D* hiter_mis[maxIteration]; //hiter_mis[i] misalignment position change at iteration i (defualt case : misalignment=0)
0228   TH1D* hiter_ali[maxIteration]; //hiter_ali[i] absolute position change at iteration i, for all detectors
0229   char ppdirname[16];
0230   sprintf(ppdirname, "Shifts%d", currentPar);
0231   fout->mkdir(ppdirname);
0232   gDirectory->cd(ppdirname);
0233   //delare histos
0234   for (unsigned int a = 0; a < maxIteration; a++){
0235     char histoname[64], histotitle[64];
0236     //      sprintf(histoname,"MisalignedShift_%d_Iter_%d",currentPar,a);
0237     //      sprintf(histotitle,"Misaligned Shift %d for iteration #%d",currentPar,a);
0238     //      hiter_mis[a]=new TH1D(histoname,histoname,1000,-1.0,1.0);
0239     sprintf(histoname, "AlignedShift_%d_Iter_%d", currentPar, a);
0240     sprintf(histotitle, "Aligned Shift %d for iteration #%d", currentPar, a);
0241     hiter_ali[a]=new TH1D(histoname, histoname, ndets, 0, ndets);
0242   }
0243   for (int i = 0; i < ndets; i++){
0244     char histoname[64], histotitle[64];
0245     sprintf(histoname, "Shift_%d_SiDet_%d", currentPar, i);
0246     sprintf(histotitle, "Shift %d for detector #%d", currentPar, i);
0247     hshift[i]=new TH1D(histoname, histoname, maxIteration, -0.5, maxIteration-0.5);
0248   }
0249
0250   //Loop on trees
0251   int modules_accepted = 0;
0252   int nHit;
0253   unsigned int detId;
0254   double tpos[3];
0255   double apos[3];
0256   //    double mpos[3];
0257   double trot[9];
0258   double arot[9];
0259   //    double mrot[9];
0260   double aerr[6];
0261
0262   TString detNames[6];
0263   detNames[0] = TString("PXB");
0264   detNames[1] = TString("PXF");
0265   detNames[2] = TString("TIB");
0266   detNames[3] = TString("TID");
0267   detNames[4] = TString("TOB");
0268   detNames[5] = TString("TEC");
0269
0270
0271   //    TTree* tmpTree_m = (TTree*)fm->Get("AlignablesAbsPos_1"); //misaligned position
0272   //TTree* tmpTree_true = (TTree*)ft->Get("AlignablesOrgPos_1");
0273   TTree* tmpTree_true = (TTree*)fa->Get("AlignablesAbsPos_0"); //starting position
0274
0275   if (currentPar < 3){
0276     //      tmpTree_m->SetBranchAddress("Pos", &mpos );
0277     tmpTree_true->SetBranchAddress("Pos", &tpos);
0278   }
0279   if (currentPar >= 3){
0280     //      tmpTree_m->SetBranchAddress("Rot", &mrot );
0281     tmpTree_true->SetBranchAddress("Rot", &trot);
0282   }
0283
0284   for (unsigned int iter = 1; iter < maxIteration; iter++){//loop on i (HIP iterations -> trees in file)
0285
0286     TTree* tmpTree = (TTree*)fa->Get(keysa->At(iter)->GetName());
0287     const char* tmpTreev = keysv->At(iter)->GetName();
0288     tmpTree->AddFriend(tmpTreev, fv);
0289     tmpTree->SetBranchAddress("Id", &detId);
0290     tmpTree->SetBranchAddress("Nhit", &nHit);
0291     tmpTree->SetBranchAddress("ParError", &aerr);
0292
0293     if (currentPar < 3){ tmpTree->SetBranchAddress("Pos", &apos); }
0294     tmpTree->SetBranchAddress("Rot", &arot);
0295
0296     //std::cout << "iteration: " << iter << "..." << std::endl;
0297     modules_accepted=0;
0298     //loop on entries
0299     for (int j = 0; j < tmpTree->GetEntries(); j++){ //loop on j (modules -> entries in each tree)
0300       tmpTree_true->GetEntry(j);
0301       tmpTree->GetEntry(j);
0302       //cout << "det" << j << "," << detId << "," << nHit << endl;
0303       //            tmpTree_m->GetEntry(j);
0304       //            tmpTree_true->GetEntry(j);
0305
0306
0307       bool passSubdetCut = true;
0308       int mysubDet=GetSubDet(detId);
0309       if (subDet > 0){ if (GetSubDet(detId) != subDet) passSubdetCut = false; }
0310
0311       int COLOR_CODE[6]={ 28, 2, 3, 4, 6, 7 };
0312       hshift[j]->SetLineColor(COLOR_CODE[GetSubDet(detId)-1]);
0313       hshift[j]->SetMarkerColor(COLOR_CODE[GetSubDet(detId)-1]);
0314
0315       if ((nHit >= minHits)&&(passSubdetCut)){
0316         //maths
0317         if (currentPar < 3){
0318           //                    TVectorD dr_mis(3, mpos);
0319           TVectorD dr_ali(3, apos);
0320           TVectorD r_true(3, tpos);
0321           TMatrixD R_true(3, 3, trot);
0322           //std::cout << "dr_ali 00 : " << dr_ali[currentPar] << std::endl;
0323           //                    dr_mis -= r_true;
0324           dr_ali -= r_true;
0325           //std::cout << "dr_ali 0 : " << dr_ali[currentPar] << std::endl;
0326           //to local
0327           //if (dr_mis != 0.) dr_mis = R_true*  dr_mis;
0328           //if (dr_ali != 0.) dr_ali = R_true*  dr_ali;
0329           //std::cout << "currentPar: " << currentPar << std::endl;
0330           //std::cout << "dr_mis: " << dr_mis[currentPar] << std::endl;
0331           //std::cout << "dr_ali: " << dr_ali[currentPar] << std::endl;
0332           //                    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;
0333           hshift[j]->SetBinContent(iter+1, dr_ali[currentPar]); //aligned position - start position
0334           //                                  hiter_mis[iter]->SetBinContent(j+1,dr_mis[currentPar]);
0335           //if(j=0&&currentPar==5) std::cout << "iter=" << iter << ",dr_ali: " << dr_ali[currentPar] << std::endl;
0336           //                                  cout << "iter=" << iter << "dr_ali: " << dr_ali[currentPar] << endl;
0337           hiter_ali[iter]->SetBinContent(j+1, dr_ali[currentPar]);
0338           //                                  cout << "bin: " << hiter_ali[iter]->GetBinContent(j+1) << endl;
0339           //                                  hiter_ali[iter]->SetBinError(j+1,5);
0340         }
0341         if (currentPar >= 3){
0342           //                    TMatrixD dR_mis(3, 3, mrot);
0343           TMatrixD dR_ali(3, 3, arot);
0344           TMatrixD R_true(3, 3, trot);
0345           //                    dR_mis = dR_mis*  TMatrixD(TMatrixD::kTransposed, R_true);
0346           dR_ali = dR_ali*  TMatrixD(TMatrixD::kTransposed, R_true);
0347           //-std::atan2(dR(2, 1), dR(2, 2)), std::asin(dR(2, 0)), -std::atan2(dR(1, 0), dR(0, 0)));
0348           //                    double dR_mis_euler = 0;
0349           double dR_ali_euler = 0;
0350           if (currentPar == 3){
0351             //                      dR_mis_euler = -std::atan2(dR_mis(2, 1), dR_mis(2, 2));
0352             dR_ali_euler = -std::atan2(dR_ali(2, 1), dR_ali(2, 2));
0353           }
0354           if (currentPar == 4){
0355             //                      dR_mis_euler = -std::asin(dR_mis(2, 0));
0356             dR_ali_euler = -std::asin(dR_ali(2, 0));
0357           }
0358           if (currentPar == 5){
0359             //                      dR_mis_euler = -std::atan2(dR_mis(1, 0), dR_mis(0, 0));
0360             dR_ali_euler = -std::atan2(dR_ali(1, 0), dR_ali(0, 0));
0361           }
0362           hshift[j]->SetBinContent(iter+1, dR_ali_euler);
0363           //                                  hiter_mis[iter]->SetBinContent(j+1,dR_mis_euler);
0364           hiter_ali[iter]->SetBinContent(j+1, dR_ali_euler);
0365         }
0366         modules_accepted++;
0367         hiter_ali[iter]->GetXaxis()->SetBinLabel(j+1, detNames[GetSubDet(detId)-1]);
0368         hiter_ali[iter]->SetBinError(j+1, aerr[currentPar]);
0369         //hiter_ali[iter]->SetBinError(j+1,0.001);
0370         hiter_ali[iter]->LabelsOption("u", "x");
0371       }
0372     }
0373     //      delete tmpTree;
0374   }
0375   delete tmpTree_true;
0376
0377   std::cout << "Modules accepted: " << modules_accepted << std::endl;
0378   std::cout << "Writing..." << std::endl;
0379   //Save
0380   fout->Write();
0381   //delete hpar1;
0382   std::cout << "Deleting..." << std::flush;
0383   delete fout;
0384   std::cout << "Deleting..." << std::flush;
0385   delete fa;
0386   //    delete ft;
0387   //    delete fm;
0388   std::cout << "Deleting..." << std::endl;
0389   delete fv;
0390
0391
0392 }
0393
0394 void HIPplots::plotAlignParams(string ShiftsOrParams, char* plotName){
0395
0396
0397   cout << "_|_| plotting AlignParams |_|_" << endl << "---> " << ShiftsOrParams  << endl;
0398   bool bParams = false;
0399   bool bShifts = false;
0400   if (ShiftsOrParams == "PARAMS") bParams = true;
0401   if (ShiftsOrParams == "SHIFTS") bShifts = true;
0402
0403   int i = 0;
0404
0405   TFile* f = new TFile(_outFile, "READ");
0406   //    f->ls();
0407
0408   TCanvas* c_params = new TCanvas("can_params", "CAN_PARAMS", 1200, 900);
0409   c_params->Divide(3, 2);
0410   //    cout << "(1) I am in " << gDirectory->GetPath() << endl;
0411   TDirectory* d;
0412   int ndets = 0;
0413   if (bParams){
0414     d = (TDirectory*)f->Get("ShiftsPar0");
0415     ndets = GetNIterations(d, "Par_0_SiDet_");
0416     //  std::cout << "INHERE!" << std::endl;
0417   }
0418   if (bShifts){
0419     d = (TDirectory*)f->Get("Shifts0");
0420     ndets = GetNIterations(d, "Shift_0_SiDet_");
0421     //  std::cout << "INHERE!" << std::endl;
0422   }
0423
0424   for (int iPar = 0; iPar < 6; iPar++){
0425
0426     c_params->cd(iPar+1);
0427     gPad->SetTopMargin(0.15);
0428     gPad->SetBottomMargin(0.15);
0429     char ppdirname[32];
0430     if (bParams) sprintf(ppdirname, "ShiftsPar%d", iPar);
0431     if (bShifts) sprintf(ppdirname, "Shifts%d", iPar);
0432     if (iPar > 0)gDirectory->cd("../");
0433     gDirectory->cd(ppdirname);
0434     //  cout << "(2) I am in " << gDirectory->GetPath() << endl;
0435
0436     TH1D* hpar1[ndets];
0437     char histoname[16];
0438     int sampling_ratio=1;
0439     int ndets_plotted=(int)ndets/sampling_ratio;
0440     cout << "Plotting " << ndets_plotted << " detectors over a total of " << ndets << endl;
0441     i=0;
0442     double histomax, histomin;
0443     if (iPar>=3){    //alpha, beta, gamma
0444       histomax=0.5;//in mrad
0445       histomin=-0.5;
0446     }
0447     else if (iPar>=2) {  //w
0448       if (bShifts){
0449         histomax=200.0;//in microns
0450         histomin=-200.0;
0451       }
0452       else{
0453         histomax=100.0;//in microns
0454         histomin=-100.0;
0455       }
0456     }
0457     else {  //u, v
0458       if (bShifts){
0459         histomax=200.0;//in microns
0460         histomin=-200.0;
0461       }
0462       else{
0463         histomax=100.0;//in microns
0464         histomin=-100.0;
0465       }
0466     }
0467     while ((i<ndets_plotted) && (i*sampling_ratio<ndets)){
0468       if (bParams) sprintf(histoname, "Par_%d_SiDet_%d", iPar, i*sampling_ratio);
0469       if (bShifts) sprintf(histoname, "Shift_%d_SiDet_%d", iPar, i*sampling_ratio);
0470       hpar1[i]=(TH1D*)gDirectory->Get(histoname);
0471
0472       if (iPar>=3)hpar1[i]->Scale(1000.0); //convert from rad to mrad
0473       else hpar1[i]->Scale(10000.0); //convert from cm to um
0474
0475       hpar1[i]->SetMarkerStyle(7);
0476       hpar1[i]->SetStats(0);
0477
0478       double tmpmax=hpar1[i]->GetBinContent(hpar1[i]->GetMaximumBin());
0479       double tmpmin=hpar1[i]->GetBinContent(hpar1[i]->GetMinimumBin());
0480
0481       //Auto-adjust axis range
0482       if (tmpmax>histomax)histomax=tmpmax*1.2;
0483       if (tmpmin<histomin)histomin=tmpmin*1.2;
0484
0485       //if(i%1300==0)cout << "Actual maximum is " << histomax << " Actual minimum is " << histomin << endl;
0486       if (i==0){
0487
0488         hpar1[i]->SetXTitle("Iteration");
0489
0490         if (bParams){
0491           if (iPar==0)hpar1[i]->SetYTitle("#delta u param (#mum)");
0492           else if (iPar==1)hpar1[i]->SetYTitle("#delta v param (#mum)");
0493           else if (iPar==2)hpar1[i]->SetYTitle("#delta w param (#mum)");
0494           else if (iPar==3)hpar1[i]->SetYTitle("#delta #alpha param (mrad)");
0495           else if (iPar==4)hpar1[i]->SetYTitle("#delta #beta param (mrad)");
0496           else if (iPar==5)hpar1[i]->SetYTitle("#delta #gamma param (mrad)");
0497           else hpar1[i]->SetYTitle("dunno");
0498         }
0499         else if (bShifts){
0500           if (iPar==0)hpar1[i]->SetYTitle("#delta u shift (#mum)");
0501           else if (iPar==1)hpar1[i]->SetYTitle("#delta v shift (#mum)");
0502           else if (iPar==2)hpar1[i]->SetYTitle("#delta w shift (#mum)");
0503           else if (iPar==3)hpar1[i]->SetYTitle("#delta #alpha shift (mrad)");
0504           else if (iPar==4)hpar1[i]->SetYTitle("#delta #beta shift (mrad)");
0505           else if (iPar==5)hpar1[i]->SetYTitle("#delta #gamma shift (mrad)");
0506           else hpar1[i]->SetYTitle("dunno");
0507         }
0508
0509         hpar1[i]->GetYaxis()->SetTitleOffset(1.5);
0510         //hpar1[i]->SetTitle("Par1: x shift");
0511         hpar1[i]->SetTitle("");
0512         hpar1[i]->SetMaximum(histomax);
0513         hpar1[i]->SetMinimum(histomin);
0514         hpar1[i]->Draw("PL");
0515       }
0516
0517       else hpar1[i]->Draw("PLsame");
0518       i++;
0519     }//end loop on NDets
0520     hpar1[0]->SetMaximum(histomax);
0521     hpar1[0]->SetMinimum(histomin);
0522     //hpar1[0]->SetLineColor(1);
0523     cout << "Plotted " << i << " aligned detectors" << endl;
0524   }//end loop on pars
0525   c_params->cd();
0526   TLegend*  legend = MakeLegend(.1, .93, .9, .98);
0527   legend->Draw();
0528   c_params->SaveAs(plotName);
0529   std::cout << "Deleting..." << std::flush;
0530   delete c_params;
0531   std::cout << "Deleting..." << std::flush;
0532   //delete f;
0533   std::cout << "Deleting..." << std::endl;
0534
0535 }//end PlotAlignPars
0536
0537
0538
0539 void HIPplots::plotAlignParamsAtIter(int iter, string ShiftsOrParams, char* plotName){
0540
0541   cout << "Welcome to  HIPplots::plotAlignParamsAtIter " << iter << endl;
0542   bool bParams = false;
0543   bool bShifts = false;
0544   if (ShiftsOrParams == "PARAMS") bParams = true;
0545   else if (ShiftsOrParams == "SHIFTS") bShifts = true;
0546   else { cout << "ERROR in plotAliParamsAtIter!!! Wrong input argument: " << ShiftsOrParams << " . Exiting" << endl; return; }
0547
0548   int i = 0;
0549   TFile* f = new TFile(_outFile, "READ");
0550   //f->ls();
0551
0552   TCanvas* c_params = new TCanvas("can_params", "CAN_PARAMS", 1200, 700);
0553   c_params->Divide(3, 2);
0554   //cout << "(1) I am in " << gDirectory->GetPath() << endl;
0555
0556   //TDirectory* d = (TDirectory*)f->Get("ShiftsPar0");
0557   //const int ndets = GetNIterations(d,"Par_0_SiDet_");
0558
0559   for (int iPar = 0; iPar < 6; iPar++){
0560
0561     c_params->cd(iPar+1);
0562
0563     gPad->SetTopMargin(0.15);
0564     gPad->SetBottomMargin(0.15);
0565     char ppdirname[32];
0566     if (bParams) sprintf(ppdirname, "ShiftsPar%d", iPar);
0567     if (bShifts) sprintf(ppdirname, "Shifts%d", iPar);
0568     if (iPar > 0)gDirectory->cd("../");
0569     gDirectory->cd(ppdirname);
0570     //cout << "(2) I am in " << gDirectory->GetPath() << endl;
0571
0572     TH1D* hiter;
0573     char histoname[16];
0574     if (bParams) sprintf(histoname, "Par_%d_Iter_%d", iPar, iter);
0575     if (bShifts) sprintf(histoname, "AlignedShift_%d_Iter_%d", iPar, iter);
0576     hiter = (TH1D*)gDirectory->Get(histoname);
0577
0578
0579     //hiter->SetMarkerStyle(1);
0580     hiter->SetStats(1);
0581
0582     hiter->SetXTitle("Sub-Det");
0583     if (iPar==0)hiter->SetYTitle("#delta u (#mum)");
0584     else if (iPar==1)hiter->SetYTitle("#delta v (#mum)");
0585     else if (iPar==2)hiter->SetYTitle("#delta w (#mum)");
0586     else if (iPar==3)hiter->SetYTitle("#delta #alpha (#murad)");
0587     else if (iPar==4)hiter->SetYTitle("#delta #beta (#murad)");
0588     else if (iPar==5)hiter->SetYTitle("#delta #gamma (#murad)");
0589     else hiter->SetXTitle("dunno");
0590     hiter->GetYaxis()->SetTitleOffset(1.5);
0591
0592     double histomax, histomin;
0593     if (iPar==2||iPar==5) { histomax=200; histomin=-200; }
0594     else { histomax=100; histomin=-100; }
0595
0596     if (iPar>=3)hiter->Scale(1000000.0); //convert from rad to micro rad
0597     else hiter->Scale(10000.0); //convert from cm to um
0598     double tmpmax=hiter->GetBinContent(hiter->GetMaximumBin());
0599     double tmpmin=hiter->GetBinContent(hiter->GetMinimumBin());
0600     if (tmpmax>histomax)histomax=tmpmax*1.2;
0601     if (tmpmin<histomin)histomin=tmpmin*1.2;
0602     hiter->SetMaximum(histomax);
0603     hiter->SetMinimum(histomin);
0604     //hiter->SetAxisRange(0, 6, "X");
0605     hiter->SetLineColor(1);
0606     TColor* col = gROOT->GetColor(kCyan+1);
0607     col->SetAlpha(0.3);
0608     hiter->SetFillColor(col->GetNumber());
0609     //          hiter->SetFillColor(kCyan);
0610     hiter->SetMarkerStyle(7);
0611     hiter->Draw();
0612     hiter->Draw("PE1");
0613     /*
0614         if (bShifts) {
0615         TH1D* hiter2;
0616         char histoname2[16];
0617         sprintf( histoname2, "MisalignedShift_%d_Iter_%d", iPar, iter );
0618         hiter2 = (TH1D*) gDirectory->Get(histoname2);
0619         hiter2->SetLineColor(1);
0620         hiter2->SetFillColor(kRed-9);
0621         hiter2->Draw("SAME");
0622         }
0623
0624         */
0625     //      cout << "Plotted " << i << " aligned detectors" << endl;
0626   }//end loop on pars
0627   c_params->SaveAs(plotName);
0628   delete c_params;
0629   std::cout << "Deleting..." << std::endl;
0630   delete f;
0631
0632 }//end PlotAlignParamsAtIter
0633
0634
0635
0636 void HIPplots::extractAlignableChiSquare(int minHits, int subDet, int doubleSided){
0637
0638   //TFile* fp = new TFile( _inFile_params,"READ");
0639   //    const TList* keysp = fp->GetListOfKeys();
0640   //    const unsigned int maxIteration = keysp->GetSize();
0641   cout << "\n--- Welcome to extractAlignableChiSquare ---" << endl;
0642   cout << "\nInput parameters:\n\tMinimum number of hits per alignbale = " << minHits << "\n\tSubdetetctor selection" << subDet << endl;
0643
0644   if (minHits<1){
0645     cout << "Warning ! Allowing to select modules with NO hits. Chi2 not defined for them. Setting automatically minNhits=1 !!!" << endl;
0646     minHits=1;
0647   }
0648
0649   TFile* fv = new TFile(_inFile_uservars, "READ");
0650   const TList* keysv = fv->GetListOfKeys();
0651   const unsigned int maxIteration = keysv->GetSize() - 1;
0652   cout << "MaxIteration is " << maxIteration  << endl;
0653
0654   char fileaction[16];
0655   if (CheckFileExistence(_outFile))sprintf(fileaction, "UPDATE");
0656   else sprintf(fileaction, "NEW");
0657   TFile* fout = new TFile(_outFile, fileaction);
0658
0659   TTree* tree0 = (TTree*)fv->Get(keysv->At(0)->GetName());
0660   unsigned int detId0;
0661   tree0->SetBranchAddress("Id", &detId0);
0662   const int ndets=tree0->GetEntries();
0663   TH1D* halichi2n[ndets];//norm chi2 for each module as a function of iteration
0664   TH1D* htotchi2n[maxIteration];//distrib of norm chi2 for all modules at a given iteration
0665   TH1D* hprobdist[maxIteration];
0666   char ppdirname[16];
0667   sprintf(ppdirname, "AlignablesChi2n");
0668   fout->mkdir(ppdirname);
0669   gDirectory->cd(ppdirname);
0670
0671   for (int iter = 1; iter <=(int)maxIteration; iter++){
0672     char histoname[64], histotitle[128];
0673     sprintf(histoname, "Chi2n_iter%d", iter);
0674     sprintf(histotitle, "Distribution of Normalised #chi^{2} for All alignables  at iter %d", iter);
0675     htotchi2n[iter-1]=new TH1D(histoname, histotitle, 1000, 0.0, 10.0);
0676     sprintf(histoname, "ProbChi2_%d", iter);
0677     sprintf(histotitle, "Distribution of  Prob(#chi^{2},ndof) at iter %d", iter);
0678     hprobdist[iter-1]=new TH1D(histoname, histotitle, 100, 0.0, 1.0);
0679   }
0680   gDirectory->mkdir("AlignablewiseChi2n");
0681   gDirectory->cd("AlignablewiseChi2n");
0682   //declare histos
0683   for (int i = 0; i < (int)ndets; i++){
0684     tree0->GetEntry(i);
0685     char histoname[64], histotitle[64];
0686     sprintf(histoname, "Chi2n_%d", i);
0687     sprintf(histotitle, "Normalised #chi^{2} for detector #%d", i);
0688     halichi2n[i]=new TH1D(histoname, histotitle, maxIteration, 0.5, maxIteration+0.5);
0689   }
0690
0691   //Loop on trees and fill histos
0692   int modules_accepted = 0;
0693   for (unsigned int iter = 1; iter <= maxIteration; iter++){//loop on i (HIP iterations -> trees in file)
0694
0695     TTree* tmpTreeUV = (TTree*)fv->Get(keysv->At(iter)->GetName()); //get the UserVariable tree at each iteration
0696     cout << "Taking tree " << keysv->At(iter)->GetName() << endl;
0697     //tmpTreeUV->GetListOfLeaves()->ls();
0698     tmpTreeUV->SetBranchStatus("*", 0);
0699     tmpTreeUV->SetBranchStatus("Id", 1);
0700     tmpTreeUV->SetBranchStatus("Nhit", 1);
0701     tmpTreeUV->SetBranchStatus("AlignableChi2", 1);
0702     tmpTreeUV->SetBranchStatus("AlignableNdof", 1);
0703     double alichi2=0.0;
0704     unsigned int alindof=0;
0705     int nHit=0;
0706     unsigned int detId=0;
0707     tmpTreeUV->SetBranchAddress("AlignableChi2", &alichi2);
0708     tmpTreeUV->SetBranchAddress("AlignableNdof", &alindof);
0709     tmpTreeUV->SetBranchAddress("Id", &detId);
0710     tmpTreeUV->SetBranchAddress("Nhit", &nHit);
0711
0712     modules_accepted=0;
0713     //loop on entries
0714     for (int j = 0; j < tmpTreeUV->GetEntries(); j++){ //loop on j (modules -> entries in each tree)
0715       tmpTreeUV->GetEntry(j);
0716
0717       bool passSubdetCut = true;
0718       if (subDet > 0){ if (GetSubDet(detId) != subDet) passSubdetCut = false; } //get subDet Id from module detId,and plot only the selected subDet
0719       if (doubleSided > 0){
0720         if (GetSubDet(detId)%2 == 1 && doubleSided == 1 && GetBarrelLayer(detId) < 3) passSubdetCut = false;
0721         if (GetSubDet(detId)%2 == 1 && doubleSided == 2 && GetBarrelLayer(detId) > 2) passSubdetCut = false;
0722       }
0723
0724       if ((nHit >= minHits)&&(passSubdetCut)){ //cut on nhits per module
0725         halichi2n[j]->SetBinContent(iter, double(alichi2 / alindof));
0726         //  halichi2n[j]->SetBinContent(iter,double(alichi2));
0727         double prob=TMath::Prob(double(alichi2), int(alindof));
0728
0729         htotchi2n[iter-1]->Fill(double(alichi2 / alindof));
0730         hprobdist[iter-1]->Fill(prob);
0731         modules_accepted++;
0732       }
0733
0734     }//end loop on j - alignables
0735     cout << "alignables accepted at iteration " << iter << " = " << modules_accepted << endl;
0736     delete tmpTreeUV;
0737   }//end loop on iterations
0738
0739   //Ma che e'???
0740   /* cout << "Prob for chi2=0,02 ndof=40 -> " << TMath::Prob(0.02,40) << endl;
0741   cout << "Prob for chi2=20, ndof=40 -> " << TMath::Prob(20.0,40) << endl;
0742   cout << "Prob for chi2=40, ndof=40 -> " << TMath::Prob(40.0,40) << endl;
0743   cout << "Prob for chi2=60, ndof=40 -> " << TMath::Prob(60.0,40) << endl;
0744   cout << "Prob for chi2=2000, ndof=40 -> " << TMath::Prob(2000.0,40) << endl; */
0745
0746   //Save
0747   fout->Write();
0748   delete fout;
0749   delete fv;
0750   cout << "Finished extractAlignableChiSquare" << endl;
0751 }//end extractAlignableChiSquare
0752
0753
0754 void HIPplots::extractSurveyResiduals(int currentPar, int subDet){
0755
0756   cout << "\n---  extractSurveyResiduals has been called ---" << endl;
0757
0758   TFile* fv = new TFile(_inFile_surveys, "READ");
0759   const TList* keysv = fv->GetListOfKeys();
0760   const unsigned int maxIteration = keysv->GetSize();
0761   cout << "MaxIteration is " << maxIteration  << endl;
0762
0763   char fileaction[16];
0764   if (CheckFileExistence(_outFile))sprintf(fileaction, "UPDATE");
0765   else sprintf(fileaction, "NEW");
0766   TFile* fout = new TFile(_outFile, fileaction);
0767
0768   TTree* tree0 = (TTree*)fv->Get(keysv->At(0)->GetName());
0769   unsigned int detId0;
0770   tree0->SetBranchAddress("Id", &detId0);
0771   const int ndets=tree0->GetEntries();
0772   TH1D* hsurvey[ndets];//norm chi2 for each module as a function of iteration
0773   TH1D* htotres[maxIteration];//distrib of norm chi2 for all modules at a given iteration
0774
0775   char ppdirname[16];
0776   sprintf(ppdirname, "SurveyResiduals");
0777   fout->mkdir(ppdirname);
0778   gDirectory->cd(ppdirname);
0779
0780   for (int iter = 1; iter <=(int)maxIteration; iter++){
0781     char histoname[64], histotitle[128];
0782     sprintf(histoname, "SurveyRes_Par%d_iter%d", currentPar, iter);
0783     sprintf(histotitle, "Distribution of survey Residuals for All alignables ; Par=%d ; iter %d", currentPar, iter);
0784     htotres[iter-1]=new TH1D(histoname, histotitle, 1000, 0.0, 10.0);
0785   }
0786   gDirectory->mkdir("SurveyResiduals_alignables");
0787   gDirectory->cd("SurveyResiduals_alignables");
0788   //declare histos
0789   for (int i = 0; i < (int)ndets; i++){
0790     tree0->GetEntry(i);
0791     char histoname[64], histotitle[64];
0792     sprintf(histoname, "SurveyRes_Par%d_%d", currentPar, i);
0793     sprintf(histotitle, "Survey residual for detector #%d - Par=%d", i, currentPar);
0794     hsurvey[i]=new TH1D(histoname, histotitle, maxIteration, 0.5, maxIteration+0.5);
0795   }
0796
0797   //Loop on trees and fill histos
0798   int modules_accepted = 0;
0799   for (unsigned int iter = 1; iter <= maxIteration; iter++){//loop on i (HIP iterations -> trees in file)
0800
0801     TTree* tmpTreeUV = (TTree*)fv->Get(keysv->At(iter)->GetName()); //get the UserVariable tree at each iteration
0802     cout << "Taking tree " << keysv->At(iter)->GetName() << endl;
0803     //tmpTreeUV->GetListOfLeaves()->ls();
0804     tmpTreeUV->SetBranchStatus("*", 0);
0805     tmpTreeUV->SetBranchStatus("Id", 1);
0806     tmpTreeUV->SetBranchStatus("Par", 1);
0807
0808     double par[6];
0809     unsigned int detId=0;
0810     tmpTreeUV->SetBranchAddress("Par", &par);
0811     tmpTreeUV->SetBranchAddress("Id", &detId);
0812
0813
0814     modules_accepted=0;
0815     //loop on entries
0816     for (int j = 0; j < tmpTreeUV->GetEntries(); j++){ //loop on j (modules -> entries in each tree)
0817       tmpTreeUV->GetEntry(j);
0818
0819       bool passSubdetCut = true;
0820       if (subDet > 0){ if (GetSubDet(detId) != subDet) passSubdetCut = false; }
0821       if (passSubdetCut){
0822         hsurvey[j]->SetBinContent(iter, double(par[currentPar]));
0823
0824         modules_accepted++;
0825       }
0826
0827     }//end loop on j - alignables
0828     cout << "alignables accepted at iteration " << iter << " = " << modules_accepted << endl;
0829     delete tmpTreeUV;
0830   }//end loop on iterations
0831
0832   //Save
0833   fout->Write();
0834   delete fout;
0835   delete fv;
0836   cout << "Finished extractAlignableChiSquare" << endl;
0837 }//end extractAlignableChiSquare
0838
0839
0840
0841
0842
0843
0844 void HIPplots::plotAlignableChiSquare(char* plotName, float minChi2n){
0845
0846   int i = 0;
0847   TFile* f = new TFile(_outFile, "READ");
0848   TCanvas* c_alichi2n=new TCanvas("can_alichi2n", "CAN_ALIGNABLECHI2N", 900, 900);
0849   c_alichi2n->cd();
0850   TDirectory* chi_d1=(TDirectory*)f->Get("AlignablesChi2n");
0851   const int maxIteration= GetNIterations(chi_d1, "Chi2n_iter", 1)-1;
0852   cout << "N iterations " << maxIteration << endl;
0853   //take the histos prepared with extractAlignableChiSquare
0854   gDirectory->cd("AlignablesChi2n");
0855   TDirectory* chi_d=(TDirectory*)gDirectory->Get("AlignablewiseChi2n");
0856   const int ndets= GetNIterations(chi_d, "Chi2n_");
0857
0858   gDirectory->cd("AlignablewiseChi2n");
0859   TH1D* hchi2n[ndets];
0860   char histoname[64];
0861   int sampling_ratio=1;
0862   int ndets_plotted=(int)ndets/sampling_ratio;
0863   cout << "Sampling " << ndets_plotted << " detectors over a total of " << ndets << endl;
0864   double histomax=0.1, histomin=-0.1;
0865   bool firstplotted=false;
0866   int firstplottedindex=0;
0867   int totalplotted=0;
0868   bool plothisto=true;
0869   while ((i<ndets_plotted) && (i*sampling_ratio<ndets)){
0870     sprintf(histoname, "Chi2n_%d", i);
0871     hchi2n[i]=(TH1D*)gDirectory->Get(histoname);
0872
0873     //if required, check that the chi2n exceeds at any point the threshold set as input
0874     bool chi2ncut=false;
0875     if (minChi2n>0.0){
0876       for (int bin=1; bin<=hchi2n[i]->GetNbinsX(); bin++){
0877         if (hchi2n[i]->GetBinContent(bin)>minChi2n){ chi2ncut=true; break; }
0878       }
0879     }
0880     else chi2ncut=true;
0881
0882     //if required check that the chi2 is increasing three iterations in a row and plot only those histos
0883     if (plotbadchi2)plothisto=CheckHistoRising(hchi2n[i]);
0884     else plothisto=true;
0885
0886     if (chi2ncut&&plothisto){
0887       double tmpmax=hchi2n[i]->GetBinContent(hchi2n[i]->GetMaximumBin());
0888       double tmpmin=hchi2n[i]->GetBinContent(hchi2n[i]->GetMinimumBin());
0889       histomax=2.0;
0890       histomin=0.0;
0891       if (tmpmax>histomax)histomax=tmpmax;
0892       if (tmpmin<histomin)histomin=tmpmin;
0893
0894       hchi2n[i]->SetMaximum(histomax);
0895       hchi2n[i]->SetMinimum(histomin);
0896       hchi2n[i]->SetStats(0);
0897
0898
0899       if (!firstplotted){
0900         hchi2n[i]->SetXTitle("Iteration");
0901         hchi2n[i]->SetYTitle("#chi^{2} / # dof");
0902         //  hchi2n[i]->SetYTitle("#chi^{2}");
0903         hchi2n[i]->SetTitle("Reduced #chi^{2} for alignables");
0904         hchi2n[i]->Draw("PL");
0905         firstplotted=true;
0906         firstplottedindex=i;
0907       }
0908       else hchi2n[i]->Draw("PLsame");
0909       totalplotted++;
0910     }//END IF plot it because it passed the chi2n cut
0911     i++;
0912   }//end while loop on alignables
0913   hchi2n[firstplottedindex]->SetMaximum(histomax*1.2);
0914   hchi2n[firstplottedindex]->SetMinimum(histomin*0.8);
0915   //override auto-resize of axis scale
0916   hchi2n[firstplottedindex]->SetMaximum(histomax);
0917   hchi2n[firstplottedindex]->SetMinimum(histomin);
0918
0919   cout << "Plotted " << totalplotted << " alignables over an initial sample of " << ndets_plotted << endl;
0920   TText* txtchi2n_1=new TText();
0921   txtchi2n_1->SetTextFont(63);
0922   txtchi2n_1->SetTextSize(22);
0923   char strchi2n_1[128];
0924   sprintf(strchi2n_1, "Plotted alignables (Chi2n > %.3f): %d / %d", minChi2n, totalplotted, ndets_plotted);
0925   txtchi2n_1->DrawText(1.2, 0.0, strchi2n_1);
0926   char finplotname[192];
0927   sprintf(finplotname, "%s.png", plotName);
0928   c_alichi2n->SaveAs(finplotname);
0929   delete c_alichi2n;
0930   //delete hchi2n;
0931
0932
0933   cout << "Doing distrib" << endl;
0934   gDirectory->cd("../");
0935   TCanvas* c_chi2ndist=new TCanvas("can_chi2ndistr", "CAN_CHI2N_DISTRIBUTION", 900, 900);
0936   c_chi2ndist->cd();
0937   TH1D* hiter[maxIteration];
0938   TLegend* l=new TLegend(0.7, 0.7, 0.9, 0.9);
0939   l->SetFillColor(0);
0940   l->SetBorderSize(0);
0941   int colors[10]={ 1, 2, 8, 4, 6, 7, 94, 52, 41, 45 };
0942   int taken=0;
0943   float tmpmax=1.0;
0944   float newmax=0.0;
0945   for (i=0; i<maxIteration; i++){
0946     sprintf(histoname, "Chi2n_iter%d", i+1);
0947     hiter[i]=(TH1D*)gDirectory->Get(histoname);
0948     hiter[i]->SetXTitle("#chi^{2} / dof");
0949     hiter[i]->SetYTitle("Alignables");
0950     hiter[i]->SetTitle("Normalised #chi^{2} of Alignables");
0951     hiter[i]->Rebin(5);
0952     hiter[i]->GetXaxis()->SetRangeUser(0.0, 3.0);
0953     hiter[i]->SetLineColor(i);
0954
0955     char legend_entry[64];
0956     float histmean=hiter[i]->GetMean();
0957     if (i==0){
0958       hiter[i]->SetLineColor(colors[taken]);
0959       taken++;
0960       sprintf(legend_entry, "Norm #chi^{2}; Iter %d; %.4f", i, histmean);
0961       l->AddEntry(hiter[i], legend_entry, "L");
0962       tmpmax=hiter[i]->GetBinContent(hiter[i]->GetMaximumBin());
0963       newmax=tmpmax;
0964       //hiter[i]->Draw("");
0965     }
0966     else if ((i+1)%5==0){
0967       hiter[i]->SetLineColor(colors[taken]);
0968       taken++;
0969       sprintf(legend_entry, "Norm #chi^{2}; Iter %d; %.4f", i+1, histmean);
0970       l->AddEntry(hiter[i], legend_entry, "L");
0971       tmpmax=hiter[i]->GetBinContent(hiter[i]->GetMaximumBin());
0972       if (tmpmax>newmax)newmax=tmpmax;
0973       //hiter[i]->Draw("same");
0974     }
0975   }
0976   cout << "NewMax after 1st loop -> " << newmax << endl;
0977
0978   for (i=0; i<maxIteration; i++){
0979     hiter[i]->SetMaximum(newmax*1.1);
0980     if (i==1)      hiter[i]->Draw("");
0981     else if ((i+1)%5==0) hiter[i]->Draw("same");
0982   }
0983   l->Draw();
0984
0985   sprintf(finplotname, "%s_distr.png", plotName);
0986   cout << finplotname << endl;
0987   c_chi2ndist->SaveAs(finplotname);
0988   c_chi2ndist->SetLogy();
0989   sprintf(finplotname, "%s_distrlog.png", plotName);
0990   cout << finplotname << endl;
0991   c_chi2ndist->SaveAs(finplotname);
0992
0993   delete  c_chi2ndist;
0994   delete f;
0995 }//end plotAlignableChiSquare
0996
0997
0998
0999 //-----------------------------------------------------------------
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;
1143     ///////////////////////////
1144     //////////////////////////
1145     //////////////////////////
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
1187       ///-************
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
1211       ///-************
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
1331     ///////////////////
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);
1352     /////////
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
1492   ///
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