Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:14:29

0001 #define analyse_residuals_cxx
0002 #include "analyse_residuals.h"
0003 #include <TH1F.h>
0004 #include <TH2.h>
0005 #include <TStyle.h>
0006 #include <TCanvas.h>
0007 #include <iostream>
0008 #include <TF1.h>
0009 #include <TMath.h>
0010 #include "TFile.h"
0011 
0012 double ng(Double_t* xx, Double_t* par)
0013 {
0014   Double_t x = xx[0];
0015   Double_t norm = par[0];
0016   Double_t mean = par[1];
0017   Double_t sigma = par[2];
0018 
0019   return norm * 1.0/sqrt(2.0*TMath::Pi())/sigma * exp(-0.5*(x-mean)*(x-mean)/sigma/sigma);
0020 }
0021 
0022 const double a_min = 1.37078;
0023 const double a_max = 1.77078;
0024 const double a_bin = 0.10000;
0025 
0026 double aa_min[3] = {1.50, 1.45, 1.40};
0027 double aa_max[3] = {1.75, 1.70, 1.65};
0028 
0029 double ys_bl[6];
0030 double ys_bh[6];
0031 
0032 TCanvas* can_y_barrel_sizey_alpha[6][4];
0033 TCanvas* can_y_barrel_sizey[6];
0034 TH1F* h_yres_npix_alpha_beta[6][4][10]; 
0035 TH1F* h_yres_npix_alpha[6][4]; 
0036 TH1F* h_yres_npix_alpha_rms[6][4]; 
0037 
0038 TCanvas* can_x_barrel_sizex_beta[3][4];
0039 TCanvas* can_x_barrel_sizex[3];
0040 TH1F* h_xres_npix_beta_alpha[3][4][10];
0041 TH1F* h_xres_npix_beta[3][4];
0042 TH1F* h_xres_npix_beta_rms[3][4];
0043 
0044 
0045 /*
0046 TCanvas* can_x_barrel_sizex_beta_flipy[3][4];
0047 TCanvas* can_x_barrel_sizex_flipy[3];
0048 TH1F* h_xres_npix_beta_alpha_flipy[3][4][10];
0049 TH1F* h_xres_npix_beta_flipy[3][4];
0050 TH1F* h_xres_npix_beta_rms_flipy[3][4];
0051 
0052 TCanvas* can_x_barrel_sizex_beta_flipn[3][4];
0053 TCanvas* can_x_barrel_sizex_flipn[3];
0054 TH1F* h_xres_npix_beta_alpha_flipn[3][4][10];
0055 TH1F* h_xres_npix_beta_flipn[3][4];
0056 TH1F* h_xres_npix_beta_rms_flipn[3][4];
0057 */
0058 
0059 
0060 TCanvas* can_y_forward_sizey[2];
0061 TCanvas* can_y_forward;
0062 TH1F* h_forward_yres_npix_beta[2][10]; 
0063 TH1F* h_forward_yres_npix[2]; 
0064 TH1F* h_forward_yres_npix_rms[2]; 
0065 
0066 TCanvas* can_x_forward_sizex[2];
0067 TCanvas* can_x_forward;
0068 TH1F* h_forward_xres_npix_alpha[2][10]; 
0069 TH1F* h_forward_xres_npix[2];
0070 TH1F* h_forward_xres_npix_rms[2];
0071 
0072 const char* fname;
0073 
0074 void analyse_residuals::Loop()
0075 {
0076   if (fChain == 0) return;
0077  
0078   Long64_t nentries = fChain->GetEntries();
0079   cout << "nentries = " << nentries << endl;
0080 
0081   bool do_residuals = true;
0082   bool do_plots = true;
0083  
0084   bool do_yb = true;
0085   bool do_xb = true;
0086   bool do_yf = true;
0087   bool do_xf = true;
0088 
0089   bool do_fix = true;
0090 
0091   if ( do_residuals )
0092     fname = "residuals.dat";
0093   else
0094     fname = "pulls.dat";
0095   
0096   printf("%s \n", fname);
0097   FILE* datfile; 
0098   if ( (datfile = fopen(fname, "w"))==NULL )
0099     {      cout << "could not open the output file." << endl;
0100       exit(-1);
0101     }
0102 
0103   char hname[100];
0104 
0105   ys_bl[0] = 0.05;
0106   ys_bh[0] = 0.50;
0107 
0108   ys_bl[1] = 0.15; 
0109   ys_bh[1] = 0.90;
0110 
0111   ys_bl[2] = 0.70; 
0112   ys_bh[2] = 1.05;
0113 
0114   ys_bl[3] = 0.95; 
0115   ys_bh[3] = 1.15;
0116 
0117   ys_bl[4] = 1.15; 
0118   ys_bh[4] = 1.20;
0119 
0120   ys_bl[5] = 1.20; 
0121   ys_bh[5] = 1.40;
0122 
0123   if ( do_yb )
0124     {
0125       for (int i=0; i<6; ++i)
0126         {
0127           for (int j=0; j<4; ++j)
0128             {
0129               sprintf(hname, "can_y_barrel_sizey_alpha_%i_%i", i, j);
0130               can_y_barrel_sizey_alpha[i][j] = new TCanvas(hname, hname, 1200, 500);
0131               can_y_barrel_sizey_alpha[i][j]->Divide(5,2);
0132             }
0133          
0134           sprintf(hname, "can_y_barrel_sizey_alpha_%i", i);
0135           can_y_barrel_sizey[i] = new TCanvas(hname, hname, 1200, 350);
0136           can_y_barrel_sizey[i]->Divide(4);
0137         }
0138     }
0139   
0140   if ( do_xb )
0141     {
0142       for (int i=0; i<3; ++i)
0143         {
0144           for (int j=0; j<4; ++j)
0145             {
0146               sprintf(hname, "can_x_barrel_sizex_beta_%i_%i", i, j);
0147               can_x_barrel_sizex_beta[i][j] = new TCanvas(hname, hname, 1200, 500);
0148               can_x_barrel_sizex_beta[i][j]->Divide(5,2);
0149             }
0150           
0151           sprintf(hname, "can_x_barrel_sizex_beta_%i", i);
0152           can_x_barrel_sizex[i] = new TCanvas(hname, hname, 1200, 350);
0153           can_x_barrel_sizex[i]->Divide(4);
0154         }
0155     }
0156   
0157   if ( do_yf )
0158     {
0159       for (int i=0; i<2; ++i)
0160         {
0161           sprintf(hname, "can_y_forward_sizey_%i", i);
0162           can_y_forward_sizey[i] = new TCanvas(hname, hname, 1200, 500);
0163           can_y_forward_sizey[i]->Divide(5,2);
0164         } 
0165       
0166       sprintf(hname, "can_y_forward");
0167       can_y_forward = new TCanvas(hname, hname, 800, 400);
0168       can_y_forward->Divide(2);
0169     }
0170 
0171   if ( do_xf )
0172     {
0173       for (int i=0; i<2; ++i)
0174         {
0175           sprintf(hname, "can_x_forward_sizex_%i", i);
0176           can_x_forward_sizex[i] = new TCanvas(hname, hname, 1200, 500);
0177           can_x_forward_sizex[i]->Divide(5,2);
0178         } 
0179       
0180       sprintf(hname, "can_x_forward");
0181       can_x_forward = new TCanvas(hname, hname, 800, 400);
0182       can_x_forward->Divide(2);
0183     }
0184 
0185   // barrel histograms --------------------------------------------------------------
0186   for (int i=0; i<6; ++i) // loop over size_y
0187     for (int j=0; j<4; ++j) // loop over alpha bins
0188       for (int k=0; k<10; ++k) // loop over be5a bins
0189         {
0190           sprintf(hname, "h_yres_npix_alpha_beta_%i_%i_%i", i, j, k );
0191           if ( do_residuals )
0192             h_yres_npix_alpha_beta[i][j][k]= new TH1F(hname, hname, 100, -0.02, 0.02);
0193           else // do pulls
0194             h_yres_npix_alpha_beta[i][j][k]= new TH1F(hname, hname, 100, -10.0, 10.0);
0195         }
0196 
0197   for (int i=0; i<6; ++i) // loop over size_y
0198     for (int j=0; j<4; ++j) // loop over alpha bins
0199       {
0200         sprintf(hname, "h_yres_npix_alpha_%i_%i", i, j );
0201         h_yres_npix_alpha[i][j] = new TH1F(hname, hname, 10, ys_bl[i], ys_bh[i]);
0202         
0203         sprintf(hname, "h_yres_npix_alpha_rms_%i_%i", i, j );
0204         h_yres_npix_alpha_rms[i][j] = new TH1F(hname, hname, 10, ys_bl[i], ys_bh[i]);
0205         h_yres_npix_alpha_rms[i][j]->SetLineColor(kRed);
0206         if ( do_residuals )
0207           {
0208             h_yres_npix_alpha_rms[i][j]->SetMinimum(0.0);
0209             h_yres_npix_alpha_rms[i][j]->SetMaximum(0.0060);
0210           }
0211         else
0212           {
0213             h_yres_npix_alpha_rms[i][j]->SetMinimum(0.0);
0214             h_yres_npix_alpha_rms[i][j]->SetMaximum(2.0);
0215           }
0216       }
0217  
0218   for (int i=0; i<3; ++i) // loop over size_x
0219     for (int j=0; j<4; ++j) // loop over beta bins
0220       for (int k=0; k<10; ++k) // loop over alpha bins
0221         {
0222           sprintf(hname, "h_xres_npix_beta_alpha_%i_%i_%i", i, j, k );
0223           if ( do_residuals )
0224             h_xres_npix_beta_alpha[i][j][k] = new TH1F(hname, hname, 100, -0.01, 0.01);
0225           else // do pulls 
0226             h_xres_npix_beta_alpha[i][j][k] = new TH1F(hname, hname, 100, -10.0, 10.0);
0227           
0228           /*
0229           sprintf(hname, "h_xres_npix_beta_alpha_%i_%i_%i_flipy", i, j, k );
0230           if ( do_residuals )
0231             h_xres_npix_beta_alpha_flipy[i][j][k] = new TH1F(hname, hname, 100, -0.01, 0.01);
0232           else // do pulls 
0233             h_xres_npix_beta_alpha_flipy[i][j][k] = new TH1F(hname, hname, 100, -10.0, 10.0);
0234         
0235           sprintf(hname, "h_xres_npix_beta_alpha_%i_%i_%i_flipn", i, j, k );
0236           if ( do_residuals )
0237             h_xres_npix_beta_alpha_flipn[i][j][k] = new TH1F(hname, hname, 100, -0.01, 0.01);
0238           else // do pulls 
0239             h_xres_npix_beta_alpha_flipn[i][j][k] = new TH1F(hname, hname, 100, -10.0, 10.0); 
0240           */
0241 
0242         }
0243 
0244   for (int i=0; i<3; ++i) // loop over size_x
0245     for (int j=0; j<4; ++j) // loop over beta bins
0246       {
0247         sprintf(hname, "h_xres_npix_beta_%i_%i", i, j );
0248         h_xres_npix_beta[i][j]= new TH1F(hname, hname, 10, a_min, a_max);
0249       
0250         sprintf(hname, "h_xres_npix_beta_rms_%i_%i", i, j );
0251         h_xres_npix_beta_rms[i][j]= new TH1F(hname, hname, 10, a_min, a_max);
0252         h_xres_npix_beta_rms[i][j]->SetLineColor(kRed);
0253         if ( do_residuals )
0254           {
0255             h_xres_npix_beta_rms[i][j]->SetMinimum(0.0);
0256             h_xres_npix_beta_rms[i][j]->SetMaximum(0.0060);
0257           }
0258         else
0259           {
0260             h_xres_npix_beta_rms[i][j]->SetMinimum(0.0);
0261             h_xres_npix_beta_rms[i][j]->SetMaximum(2.0);
0262           }
0263 
0264         /*
0265         sprintf(hname, "h_xres_npix_beta_%i_%i_flipy", i, j );
0266         h_xres_npix_beta_flipy[i][j]= new TH1F(hname, hname, 10, a_min, a_max);
0267       
0268         sprintf(hname, "h_xres_npix_beta_rms_%i_%i_flipy", i, j );
0269         h_xres_npix_beta_rms_flipy[i][j]= new TH1F(hname, hname, 10, a_min, a_max);
0270         h_xres_npix_beta_rms_flipy[i][j]->SetLineColor(kRed);
0271         if ( do_residuals )
0272           {
0273             h_xres_npix_beta_rms_flipy[i][j]->SetMinimum(0.0);
0274             h_xres_npix_beta_rms_flipy[i][j]->SetMaximum(0.0060);
0275           }
0276         else
0277           {
0278             h_xres_npix_beta_rms_flipy[i][j]->SetMinimum(0.0);
0279             h_xres_npix_beta_rms_flipy[i][j]->SetMaximum(2.0);
0280           }
0281         
0282         sprintf(hname, "h_xres_npix_beta_%i_%i_flipn", i, j );
0283         h_xres_npix_beta_flipn[i][j]= new TH1F(hname, hname, 10, a_min, a_max);
0284       
0285         sprintf(hname, "h_xres_npix_beta_rms_%i_%i_flipn", i, j );
0286         h_xres_npix_beta_rms_flipn[i][j]= new TH1F(hname, hname, 10, a_min, a_max);
0287         h_xres_npix_beta_rms_flipn[i][j]->SetLineColor(kRed);
0288         if ( do_residuals )
0289           {
0290             h_xres_npix_beta_rms_flipn[i][j]->SetMinimum(0.0);
0291             h_xres_npix_beta_rms_flipn[i][j]->SetMaximum(0.0060);
0292           }
0293         else
0294           {
0295             h_xres_npix_beta_rms_flipn[i][j]->SetMinimum(0.0);
0296             h_xres_npix_beta_rms_flipn[i][j]->SetMaximum(2.0);
0297           }
0298         */
0299 
0300       }
0301 
0302   // forward histograms --------------------------------------------------------------
0303   for (int i=0; i<2; ++i) // loop ove sizey bins
0304     for (int j=0; j<10; ++j) // loop over beta bins
0305       {
0306       sprintf(hname, "h_forward_yres_npix_beta_%i_%i", i, j );
0307       if ( do_residuals )
0308         h_forward_yres_npix_beta[i][j] = new TH1F(hname, hname, 100, -0.01, 0.01);
0309       else // do pulls
0310         h_forward_yres_npix_beta[i][j] = new TH1F(hname, hname, 100, -10.0, 10.0);
0311       }
0312 
0313   for (int i=0; i<2; ++i) // loop ove sizey bins
0314     {
0315       sprintf(hname, "h_forward_yres_npix_%i", i );
0316       h_forward_yres_npix[i] = new TH1F(hname, hname, 10, 0.3, 0.4);
0317     
0318       sprintf(hname, "h_forward_yres_npix_rms_%i", i );
0319       h_forward_yres_npix_rms[i] = new TH1F(hname, hname, 10, 0.3, 0.4);
0320       h_forward_yres_npix_rms[i]->SetLineColor(kRed);
0321       if ( do_residuals )
0322         {
0323           h_forward_yres_npix_rms[i]->SetMinimum(0.0);
0324           h_forward_yres_npix_rms[i]->SetMaximum(0.0040);
0325         }
0326       else
0327         {
0328           h_forward_yres_npix_rms[i]->SetMinimum(0.0);
0329           h_forward_yres_npix_rms[i]->SetMaximum(2.0);
0330         }
0331     }
0332 
0333   for (int i=0; i<2; ++i) // loop ove sizex bins
0334     for (int j=0; j<10; ++j) // loop over alpha bins
0335       {
0336         sprintf(hname, "h_forward_xres_npix_alpha_%i_%i", i, j );
0337         if ( do_residuals )
0338           h_forward_xres_npix_alpha[i][j] = new TH1F(hname, hname, 100, -0.01, 0.01);
0339         else // do pulls
0340           h_forward_xres_npix_alpha[i][j] = new TH1F(hname, hname, 100, -10.0, 10.0);
0341       }
0342   
0343   sprintf(hname, "h_forward_xres_npix_%i", 0 );
0344   h_forward_xres_npix[0] = new TH1F(hname, hname, 10, 0.15, 0.3);
0345   sprintf(hname, "h_forward_xres_npix_%i", 1 );
0346   h_forward_xres_npix[1] = new TH1F(hname, hname, 10, 0.15, 0.5);
0347 
0348   sprintf(hname, "h_forward_xres_npix_rms_%i", 0 );
0349   h_forward_xres_npix_rms[0] = new TH1F(hname, hname, 10, 0.15, 0.3);
0350   h_forward_xres_npix_rms[0]->SetLineColor(kRed);
0351   sprintf(hname, "h_forward_xres_npix_rms_%i", 1 );
0352   h_forward_xres_npix_rms[1] = new TH1F(hname, hname, 10, 0.15, 0.5);
0353   h_forward_xres_npix_rms[1]->SetLineColor(kRed);
0354 
0355   if ( do_residuals )
0356     {
0357       h_forward_xres_npix_rms[0]->SetMinimum(0.0);
0358       h_forward_xres_npix_rms[0]->SetMaximum(0.0040);
0359       h_forward_xres_npix_rms[1]->SetMinimum(0.0);
0360       h_forward_xres_npix_rms[1]->SetMaximum(0.0040);
0361     }
0362   else
0363     {
0364       h_forward_xres_npix_rms[0]->SetMinimum(0.0);
0365       h_forward_xres_npix_rms[0]->SetMaximum(2.0);
0366       h_forward_xres_npix_rms[1]->SetMinimum(0.0);
0367       h_forward_xres_npix_rms[1]->SetMaximum(2.0);
0368     }
0369 
0370   for (Long64_t jentry=0; jentry<nentries; jentry++) 
0371     {
0372       Long64_t ientry = LoadTree(jentry);
0373       if (ientry < 0) break;
0374       fChain->GetEntry(jentry);
0375       
0376       if ( !do_residuals )
0377         {
0378           rechitresx = rechitpullx;
0379           rechitresy = rechitpully;
0380         }
0381       
0382       if ( do_fix )
0383         {
0384           alpha = trk_alpha;
0385           beta  = trk_beta ;
0386         }
0387 
0388       double alpha_rad = fabs(alpha)/180.0*TMath::Pi();
0389       double beta_rad  = fabs(beta) /180.0*TMath::Pi();
0390       double betap_rad = fabs( TMath::Pi()/2.0 - beta_rad );
0391       double alphap_rad = fabs( TMath::Pi()/2.0 - alpha_rad );
0392  
0393       if ( subdetId == 1 )
0394         {
0395           // y residuals----------------------------------------------------------------
0396           int sizey = nypix;
0397           
0398           // skip ( sizey == 1 && bigy == 1 ) clusters; the associated error is pitch_y/sqrt(12.0) 
0399           if ( !( sizey == 1 && bigy == 1  ) ) 
0400             {
0401               if ( sizey > 6 ) sizey = 6;
0402               
0403               int ind_sizey = sizey - 1;
0404               int ind_alpha = -9999;
0405               int ind_beta  = -9999; 
0406               
0407               if      ( alpha_rad <= a_min ) ind_alpha = 0;
0408               else if ( alpha_rad >= a_max ) ind_alpha = 3;
0409               else if ( alpha_rad > a_min && 
0410                         alpha_rad < a_max ) 
0411                 {
0412                   double binw = ( a_max - a_min ) / 4.0;
0413                   ind_alpha = (int)( ( alpha_rad - a_min ) / binw );
0414                 }               
0415               else cout << " Wrong alpha_rad = " << alpha_rad << endl << endl;
0416 
0417               if      ( betap_rad <= ys_bl[sizey-1] ) ind_beta = 0;
0418               else if ( betap_rad >= ys_bh[sizey-1] ) ind_beta = 9;
0419               else if ( betap_rad >  ys_bl[sizey-1] && 
0420                         betap_rad <  ys_bh[sizey-1] ) 
0421                 {
0422                   double binw = ( ys_bh[sizey-1] - ys_bl[sizey-1] ) / 8.0;
0423                   ind_beta = 1 + (int)( ( betap_rad - ys_bl[sizey-1] ) / binw );
0424                 }               
0425               else cout << " Wrong betap_rad = " << betap_rad << endl << endl;
0426                               
0427               h_yres_npix_alpha_beta[ind_sizey][ind_alpha][ind_beta]->Fill( rechitresy ); 
0428               
0429             } // if ( !( sizey == 1 && bigy == 1 ) )
0430           
0431           // x residuals----------------------------------------------------------------
0432           int sizex = nxpix;
0433           // skip ( sizex == 1 && bigx == 1 ) clusters; the associated error is pitch_x/sqrt(12.0) 
0434           if ( !( sizex == 1 && bigx == 1 ) ) 
0435             {
0436               if ( sizex > 3 ) sizex = 3;
0437               
0438               int ind_sizex = sizex - 1;
0439               int ind_beta  = -9999;
0440               int ind_alpha = -9999;
0441               
0442               if      (                     betap_rad <= 0.7 ) ind_beta = 0;
0443               else if ( 0.7 <  betap_rad && betap_rad <= 1.0 ) ind_beta = 1;
0444               else if ( 1.0 <  betap_rad && betap_rad <= 1.2 ) ind_beta = 2;
0445               else if ( 1.2 <= betap_rad                     ) ind_beta = 3;
0446               else cout << " Wrong betap_rad = " << betap_rad << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! " << endl;
0447                       
0448               if      ( alpha_rad <= aa_min[ind_sizex] ) ind_alpha = 0;
0449               else if ( alpha_rad >= aa_max[ind_sizex] ) ind_alpha = 9;
0450               else
0451                 ind_alpha = (int) ( ( alpha_rad - aa_min[ind_sizex] ) / ( ( aa_max[ind_sizex] - aa_min[ind_sizex] ) / 10.0 ) );  
0452               
0453               h_xres_npix_beta_alpha[ind_sizex][ind_beta][ind_alpha]->Fill( rechitresx ); 
0454               
0455               /*
0456               if ( flipped )
0457                 h_xres_npix_beta_alpha_flipy[ind_sizex][ind_beta][ind_alpha]->Fill( rechitresx );
0458               else
0459                 h_xres_npix_beta_alpha_flipn[ind_sizex][ind_beta][ind_alpha]->Fill( rechitresx );
0460               */
0461 
0462  
0463             } //  if ( !( sizex == 1 && bigx == 1 ) )
0464           
0465         } // if ( subdetId == 1 )
0466       else if ( subdetId == 2 )
0467         {
0468           // forward y residuals----------------------------------------------------------------
0469           int sizey = nypix;
0470           // skip ( sizex == y && bigx == y ) clusters; the associated error is pitch_y/sqrt(12.0)
0471           if ( !( sizey == 1 && bigy == 1 ) )  
0472             {
0473               if ( sizey > 2 ) sizey = 2;
0474               
0475               int ind_sizey = sizey - 1;
0476               int ind_beta  = -9999; 
0477               
0478               if      ( betap_rad < 0.3 ) ind_beta = 0;
0479               else if ( betap_rad > 0.4 ) ind_beta = 9;
0480               else 
0481                 ind_beta = (int) ( ( betap_rad - 0.3 ) / ( ( 0.4 - 0.3 ) / 10.0 ) );  
0482               
0483               h_forward_yres_npix_beta[ind_sizey][ind_beta]->Fill( rechitresy ); 
0484             
0485             } // if ( !( sizey == 1 && bigy == 1 ) )
0486           
0487           // forward x residuals----------------------------------------------------------------
0488           int sizex = nxpix;
0489           // skip ( sizex == 1 && bigx == 1 ) clusters; the associated error is pitch_x/sqrt(12.0)
0490           if ( !( sizex == 1 && bigx == 1 ) )  
0491             {
0492               if ( sizex > 2 ) sizex = 2;
0493               
0494               int ind_sizex = sizex - 1;
0495               int ind_alpha  = -9999; 
0496               
0497               if ( sizex == 1 )
0498                 {
0499                   if      ( alphap_rad < 0.15 ) ind_alpha = 0;
0500                   else if ( alphap_rad > 0.30 ) ind_alpha = 9;
0501                   else 
0502                     ind_alpha = (int) ( ( alphap_rad - 0.15 ) / ( ( 0.3 - 0.15 ) / 10.0 ) );  
0503                 }
0504               if ( sizex > 1 )
0505                 {
0506                   if      ( alphap_rad < 0.15 ) ind_alpha = 0;
0507                   else if ( alphap_rad > 0.50 ) ind_alpha = 9;
0508                   else 
0509                     ind_alpha = (int) ( ( alphap_rad - 0.15 ) / ( ( 0.5 - 0.15 ) / 10.0 ) );  
0510                 }
0511             
0512               h_forward_xres_npix_alpha[ind_sizex][ind_alpha]->Fill( rechitresx ); 
0513             
0514             } // if ( !( sizex == 1 && bigx == 1 ) )
0515           
0516         } // else if ( subdetId == 2 )
0517       else
0518         cout << " Wrong Detector ID !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! " << endl;
0519       
0520     } //  for (Long64_t jentry=0; jentry<nentries; jentry++) 
0521   
0522   float low = -999.9;
0523   float high = -999.9;
0524   if ( do_residuals )
0525     {
0526       low = -0.02;
0527       high = 0.02;
0528     }
0529   else
0530     {
0531       low = -10.0;
0532       high = 10.0;
0533     }
0534   TF1* myfunc = new TF1("myfunc", ng, low, high, 3);
0535   myfunc->SetParNames("norm", "mean", "sigma");
0536   myfunc->SetParameter(0, 100.0);
0537   myfunc->FixParameter(1, 0.0);
0538   myfunc->SetParameter(2, 0.00020);
0539   myfunc->SetLineColor(kRed);
0540   
0541   Double_t sigma  = -99999.9;
0542   Double_t ssigma = -99999.9;
0543 
0544   // ----------------------------------------- barrel Y ------------------------------------------------------------
0545 
0546   if ( do_yb )
0547     {
0548       for (int k=0; k<6; ++k)
0549         {
0550           for (int i=0; i<4; ++i)
0551             for (int j=0; j<10; ++j)
0552               {
0553                 can_y_barrel_sizey_alpha[k][i]->cd(j + 1);
0554                 h_yres_npix_alpha_beta[k][i][j]->Draw();
0555                 
0556                 double n = h_yres_npix_alpha_beta[k][i][j]->GetEntries();
0557                 double rms = h_yres_npix_alpha_beta[k][i][j]->GetRMS();
0558                 double binw = h_yres_npix_alpha_beta[k][i][j]->GetBinWidth(1);
0559               
0560                 if ( n != 0.0 )
0561                   {
0562                     double norm = n*binw;
0563                     
0564                     myfunc->SetParameter(0, norm);
0565                     myfunc->FixParameter(1, 0.0);
0566                     myfunc->SetParameter(2, rms/3.0);
0567                     
0568                     
0569                     //if ( do_residuals && k==0 && j<3 )
0570                     if ( k==0 && j<3 )
0571                       h_yres_npix_alpha_beta[k][i][j]->Fit("myfunc", "LQR");
0572                     else
0573                       h_yres_npix_alpha_beta[k][i][j]->Fit("myfunc", "QR");
0574                     
0575                     sigma = myfunc->GetParameter(2);
0576                     ssigma = myfunc->GetParError(2);
0577                     h_yres_npix_alpha[k][i]->SetBinContent(j+1, sigma);
0578                     h_yres_npix_alpha[k][i]->SetBinError(j+1, ssigma);
0579                   
0580                     h_yres_npix_alpha_rms[k][i]->SetBinContent(j+1, rms);
0581                   }
0582 
0583         if ( sigma < 0.0 )
0584           {
0585             sigma = rms;
0586             cout << "Bad error, check fit convergence !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
0587           }
0588 
0589                 fprintf( datfile,
0590                          "%d %d %d %d %f %f \n", 
0591                          1, k, i, j, sigma, rms );
0592                 
0593               } //  for (int j=0; j<3; ++j)
0594          
0595           for (int i=0; i<4; ++i)
0596             {
0597               if ( !do_residuals )
0598                 {
0599                   h_yres_npix_alpha[k][i]->SetMinimum(0.0);
0600                   h_yres_npix_alpha[k][i]->SetMinimum(2.0);
0601                 }
0602               can_y_barrel_sizey[k]->cd(i+1);
0603               h_yres_npix_alpha_rms[k][i]->Draw("");
0604               h_yres_npix_alpha[k][i]->Draw("same");
0605             }
0606         }
0607 
0608       if ( do_plots )
0609         {
0610           for (int i=0; i<6; ++i)
0611             {
0612               for (int j=0; j<4; ++j)
0613                 {
0614                   if ( do_residuals )
0615                     sprintf(hname, "res_y_barrel_sizey_alpha_%i_%i.eps", i, j);
0616                   else
0617                     sprintf(hname, "pull_y_barrel_sizey_alpha_%i_%i.eps", i, j);
0618                 
0619                   can_y_barrel_sizey_alpha[i][j]->SaveAs(hname);
0620                 }
0621 
0622               if ( do_residuals )
0623                 sprintf(hname, "res_y_barrel_sizey_%i.eps", i);
0624               else
0625                 sprintf(hname, "pull_y_barrel_sizey_%i.eps", i);
0626               
0627               can_y_barrel_sizey[i]->SaveAs(hname);
0628             }
0629         }
0630     } // if ( do_yb )
0631 
0632 
0633   // ------------------------------------------- barrel X --------------------------------------------------------
0634 
0635 
0636   if ( do_xb )
0637     {
0638       for (int k=0; k<3; ++k)
0639         {
0640           for (int i=0; i<4; ++i)
0641             for (int j=0; j<10; ++j)
0642               {
0643                 Double_t sigma = -99999.9;
0644                 can_x_barrel_sizex_beta[k][i]->cd(j + 1);
0645                 h_xres_npix_beta_alpha[k][i][j]->Draw();
0646                 
0647                 double n = h_xres_npix_beta_alpha[k][i][j]->GetEntries();
0648                 double rms = h_xres_npix_beta_alpha[k][i][j]->GetRMS();
0649                 double binw = h_xres_npix_beta_alpha[k][i][j]->GetBinWidth(1);
0650                 
0651                 if ( n != 0.0 )
0652                   {
0653                     double norm = n*binw;
0654                     
0655                     myfunc->SetParameter(0, norm);
0656                     myfunc->FixParameter(1, 0.0);
0657                     myfunc->SetParameter(2, rms/3.0);
0658                     
0659                     h_xres_npix_beta_alpha[k][i][j]->Fit("myfunc", "QR");
0660                     sigma = myfunc->GetParameter(2);
0661                     ssigma = myfunc->GetParError(2);
0662                     h_xres_npix_beta[k][i]->SetBinContent(j+1, sigma);
0663                     h_xres_npix_beta[k][i]->SetBinError(j+1, ssigma);
0664 
0665                     h_xres_npix_beta_rms[k][i]->SetBinContent(j+1, rms);
0666                   }
0667 
0668         if ( sigma < 0.0 )
0669           {
0670             sigma = rms;
0671             cout << "Bad error, check fit convergence !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
0672           }
0673         
0674                 fprintf( datfile,
0675                          "%d %d %d %d %f %f \n",
0676                          2, k, i, j, sigma, rms );
0677                 
0678               } //  for (int j=0; j<3; ++j)
0679           
0680           for (int i=0; i<4; ++i)
0681             {
0682               if ( !do_residuals )
0683                 {
0684                   h_xres_npix_beta[k][i]->SetMinimum(0.0);
0685                   h_xres_npix_beta[k][i]->SetMinimum(2.0);
0686                 }
0687               can_x_barrel_sizex[k]->cd(i+1);
0688               h_xres_npix_beta_rms[k][i]->Draw("");
0689               h_xres_npix_beta[k][i]->Draw("same");
0690             }
0691         }
0692 
0693       if ( do_plots )
0694         {
0695           for (int i=0; i<3; ++i)
0696             {
0697               for (int j=0; j<4; ++j)
0698                 {
0699                   if ( do_residuals )
0700                     sprintf(hname, "res_x_barrel_sizex_beta_%i_%i.eps", i, j);
0701                   else
0702                     sprintf(hname, "pull_x_barrel_sizex_beta_%i_%i.eps", i, j);
0703 
0704                   can_x_barrel_sizex_beta[i][j]->SaveAs(hname);
0705                 }
0706               
0707               if ( do_residuals )
0708                 sprintf(hname, "res_x_barrel_sizex_%i.eps", i);
0709               else
0710                 sprintf(hname, "pull_x_barrel_sizex_%i.eps", i);
0711 
0712               can_x_barrel_sizex[i]->SaveAs(hname);
0713             }
0714         }
0715       
0716     } //  if ( do_xb )
0717 
0718 
0719   // ---------------------------------------------- forward Y -------------------------------------------------
0720 
0721   if ( do_yf )
0722     {
0723       for ( int k=0; k<2; ++k )
0724         for (int i=0; i<10; ++i)
0725           {
0726             can_y_forward_sizey[k]->cd(i+1);
0727             h_forward_yres_npix_beta[k][i]->Draw();
0728             
0729             double n = h_forward_yres_npix_beta[k][i]->GetEntries();
0730             double rms = h_forward_yres_npix_beta[k][i]->GetRMS();
0731             double binw = h_forward_yres_npix_beta[k][i]->GetBinWidth(1);
0732             
0733             if ( n != 0.0 )
0734               {
0735                 double norm = n*binw;
0736                 
0737                 myfunc->SetParameter(0, norm);
0738                 myfunc->FixParameter(1, 0.0);
0739                 myfunc->SetParameter(2, rms/3.0);
0740                 
0741                 h_forward_yres_npix_beta[k][i]->Fit("myfunc", "QR");
0742                 sigma = myfunc->GetParameter(2);
0743                 ssigma = myfunc->GetParError(2);
0744                 h_forward_yres_npix[k]->SetBinContent(i+1, sigma);
0745                 h_forward_yres_npix[k]->SetBinError(i+1, ssigma);
0746                 
0747                 h_forward_yres_npix_rms[k]->SetBinContent(i+1, rms);
0748               }
0749 
0750         if ( sigma < 0.0 )
0751           {
0752         sigma = rms;
0753         cout << "Bad error, check fit convergence !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
0754           }
0755         
0756             fprintf( datfile,
0757                      "%d %d %d %d %f %f \n",
0758                      3, 0, k, i, sigma, rms );
0759             
0760           }   
0761       
0762       for (int i=0; i<2; ++i)
0763         {
0764           if ( !do_residuals )
0765             {
0766               h_forward_yres_npix[i]->SetMinimum(0.0);
0767               h_forward_yres_npix[i]->SetMinimum(2.0);
0768             }
0769           can_y_forward->cd(i+1);
0770           h_forward_yres_npix_rms[i]->Draw("");
0771           h_forward_yres_npix[i]->Draw("same");
0772         }
0773       
0774       if ( do_plots )
0775         {
0776           for (int i=0; i<2; ++i)
0777             {
0778               if ( do_residuals )
0779                 sprintf(hname, "res_y_forward_sizey_%i.eps", i);
0780               else
0781                 sprintf(hname, "pull_y_forward_sizey_%i.eps", i);
0782 
0783               can_y_forward_sizey[i]->SaveAs(hname);
0784             }
0785           
0786           if ( do_residuals )
0787             sprintf(hname, "res_y_forward.eps");
0788           else
0789             sprintf(hname, "pull_y_forward.eps");
0790 
0791           can_y_forward->SaveAs(hname);
0792         }
0793     } // if ( do_yf )
0794   
0795 
0796   // ----------------------------------------------------- forward X -----------------------------------------
0797 
0798   if ( do_xf )
0799     {
0800       for (int k=0; k<2; ++k)
0801         for (int i=0; i<10; ++i)
0802           {
0803             can_x_forward_sizex[k]->cd(i+1);
0804             h_forward_xres_npix_alpha[k][i]->Draw();
0805             
0806             double n = h_forward_xres_npix_alpha[k][i]->GetEntries();
0807             double rms = h_forward_xres_npix_alpha[k][i]->GetRMS();
0808             double binw = h_forward_xres_npix_alpha[k][i]->GetBinWidth(1);
0809             
0810             if ( n != 0.0 )
0811               {
0812                 double norm = n*binw;
0813                 
0814                 myfunc->SetParameter(0, norm);
0815                 myfunc->FixParameter(1, 0.0);
0816                 myfunc->SetParameter(2, rms/3.0);
0817                 
0818                 h_forward_xres_npix_alpha[k][i]->Fit("myfunc", "QR");
0819                 sigma = myfunc->GetParameter(2);
0820                 ssigma = myfunc->GetParError(2);
0821                 h_forward_xres_npix[k]->SetBinContent(i+1, sigma);
0822                 h_forward_xres_npix[k]->SetBinError(i+1, ssigma);
0823         
0824                 h_forward_xres_npix_rms[k]->SetBinContent(i+1, rms);
0825                 
0826               }
0827 
0828         if ( sigma < 0.0 )
0829           {
0830         sigma = rms;
0831         cout << "Bad error, check fit convergence !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
0832           }
0833         
0834             fprintf( datfile,
0835                      "%d %d %d %d %f %f \n",
0836                      4, 0, k, i, sigma, rms );
0837             
0838           }
0839       
0840       for (int i=0; i<2; ++i)
0841         {
0842           if ( !do_residuals )
0843             {
0844               h_forward_xres_npix[i]->SetMinimum(0.0);
0845               h_forward_xres_npix[i]->SetMinimum(2.0);
0846             }
0847           can_x_forward->cd(i+1);
0848           h_forward_xres_npix_rms[i]->Draw("");
0849           h_forward_xres_npix[i]->Draw("same");
0850         }
0851       
0852       if ( do_plots )
0853         {
0854           for (int i=0; i<2; ++i)
0855             {
0856               if ( do_residuals )
0857                 sprintf(hname, "res_x_forward_sizex_%i.eps", i);
0858               else
0859                 sprintf(hname, "pull_x_forward_sizex_%i.eps", i);
0860 
0861               can_x_forward_sizex[i]->SaveAs(hname);
0862             }
0863 
0864           if ( do_residuals )
0865             sprintf(hname, "res_x_forward.eps");
0866           else
0867             sprintf(hname, "pull_x_forward.eps");
0868 
0869           can_x_forward->SaveAs(hname);
0870         }
0871 
0872     } // if ( do_yf )
0873 
0874 }