Back to home page

Project CMSSW displayed by LXR

 
 

    


Warning, /Validation/DTRecHits/test/plotValidation.r is written in an unsupported language. File is not indexed.

0001 //------------------------------
0002 //
0003 // A macro to decorate validation plots produced with "local" option.
0004 //
0005 // Usage:
0006 //
0007 // .x plotValidation.r("file",wheel,station)
0008 //
0009 // Some drawing options are set in the code (see below).
0010 //
0011 //
0012 // Author: N. Amapane
0013 //
0014 //------------------------------
0015 
0016 
0017 void plotValidation(){
0018   cout << endl << "Usage: .x plotValidation.r(\"inputFile.root\",<wheel>,<station>)" << endl << endl;
0019 }
0020 
0021 
0022 void plotValidation(TString filename, int wheel, int station) {
0023 
0024   if (! TString(gSystem->GetLibraries()).Contains("Histograms_h")) {
0025     gROOT->LoadMacro("$CMSSW_BASE/src/Validation/DTRecHits/interface/Histograms.h+");
0026     gROOT->LoadMacro("macros.C");
0027   }
0028 
0029   //----------------------------------------------------------------------
0030   //  Configurable options
0031    addProfile = false;
0032    addSlice = true;
0033 
0034   int rbx =2; // rebin x in scatter plots
0035   int rby =1; // rebin y in scatter plots
0036   int rbp = 1; // rebin profiles
0037   float nsigma = 2; // interval for the fit of residual distributions
0038 
0039   // Canvases to plot
0040   bool doPhiAndThetaS3 =true;
0041   bool doHitPull = true;
0042   bool doEff = true;
0043   bool doT0= true;
0044   bool doSegRes=true;
0045   bool doSegPull = true;
0046   bool doAngularDeps = true;
0047   bool doEff4D = true;
0048   bool doNSeg = true;
0049 
0050 
0051   //----------------------------------------------------------------------
0052  
0053   TStyle * style = getStyle("tdr");
0054   style->cd();  
0055   setPalette();
0056   gStyle->SetTitleSize(0.05,"XYZ"); // Set larger axis titles
0057   gStyle->SetTitleOffset(1.3,"Y");
0058   //  gStyle->SetOptTitle(0); // remove histogram titles
0059 
0060   //  gStyle->SetOptStat(111);
0061 
0062   float cmToMicron = 10000.;
0063   float vdrift = 54.3;
0064 
0065 
0066   TFile *file = new TFile(filename);
0067   
0068   HRes1DHit *hResPhi1  = new HRes1DHit(file, wheel, station, 1, "S3");
0069   HRes1DHit *hResTheta = new HRes1DHit(file, wheel, station, 2, "S3");
0070   //  HRes1DHit *hResPhi2  = new HRes1DHit(file, wheel, station, 4, "S3");
0071   HRes1DHit *hResPhi = hResPhi1;
0072 
0073   HEff1DHit* hEffS1RPhi= new HEff1DHit(file, wheel, station, 1, "S1");
0074   HEff1DHit* hEffS3RPhi= new HEff1DHit(file, wheel, station, 1, "S3");
0075   HEff1DHit* hEffS1RZ=0;
0076   HEff1DHit* hEffS3RZ=0;
0077   if (station!=4) {
0078     hEffS1RZ=   new HEff1DHit(file, wheel, station, 2, "S1");
0079     hEffS3RZ=   new HEff1DHit(file, wheel, station, 2, "S3");
0080   }
0081   
0082 
0083   HRes4DHit* hRes4D= new HRes4DHit(file, wheel, station, 0);
0084   HEff4DHit* hEff4D = new HEff4DHit(file, wheel, station, 0);
0085 
0086 
0087   // Result of fits
0088   float m_phi = 0.;
0089   float s_phi = 0.;
0090   float m_theta = 0.;
0091   float s_theta = 0.;
0092   float m_phi1 = 0.;
0093   float s_phi1 = 0.;
0094   float m_phi2 = 0.;
0095   float s_phi2 = 0.;
0096   float m_phiS1 = 0.;
0097   float s_phiS1 = 0.;
0098   float m_phiS2 = 0.;
0099   float s_phiS2 = 0.;
0100   float m_thetaS1 = 0.;
0101   float s_thetaS1 = 0.;
0102   float m_thetaS2 = 0.;
0103   float s_thetaS2 = 0.;
0104   float t0phi   =0.;
0105   float t0theta =0.;
0106 
0107   
0108   TString canvbasename = filename;
0109   canvbasename = canvbasename.Replace(canvbasename.Length()-5,5,"") + TString("_W") + (long) wheel + "_St" + (long) station ;
0110 
0111 
0112   //-------------------- Hit Residuals at step 3 in phi and theta (full distrib and vs distance from wire)
0113   if (doPhiAndThetaS3) {
0114     TCanvas* c1= new TCanvas;
0115     c1->SetTitle(canvbasename+"_ResPhiTheta"); 
0116     c1->SetName(canvbasename+"_ResPhiTheta");
0117     c1->Divide(2,2);
0118     
0119     c1->cd(1);
0120 
0121     TH1F* hRes;
0122     hRes=hResPhi->hRes;
0123 
0124     hResPhi->hRes->Rebin(2);
0125     TF1* fphi=drawGFit(hRes, nsigma, -0.4, 0.4);
0126 
0127     c1->cd(2);
0128 
0129     plotAndProfileX(hResPhi->hResVsPos,rbx,rby,rbp,-.1, .1, 0, 2.1);
0130 
0131 
0132     m_phi = fphi->GetParameter("Mean")*cmToMicron;
0133     s_phi = fphi->GetParameter("Sigma")*cmToMicron;
0134 
0135     if (hResTheta->hRes) {
0136 
0137       c1->cd(3);
0138       hRes=hResTheta->hRes;
0139 
0140       hResTheta->hRes->Rebin(2);
0141       TF1* ftheta=drawGFit(hRes, nsigma, -0.4, 0.4);
0142 
0143       c1->cd(4);  
0144       plotAndProfileX(hResTheta->hResVsPos,rbx,rby,rbp,-.1, .1, 0, 2.1);  
0145     
0146       m_theta = ftheta->GetParameter("Mean")*cmToMicron;
0147       s_theta = ftheta->GetParameter("Sigma")*cmToMicron;  
0148     }
0149     
0150 
0151     cout << canvbasename << "  Step3 W" << wheel << " St" << station << endl
0152          << "   Res:          Phi: M= " << int(floor(m_phi+0.5))
0153          << " S= "      << int(s_phi+0.5)
0154          << "; Theta: M= " << int(floor(m_theta+0.5))
0155          << " S= "  << int(s_theta+0.5) << endl;
0156   }
0157 
0158 
0159   //-------------------- Hit pulls
0160   if (doHitPull){
0161     TCanvas* c1= new TCanvas;
0162     c1->SetTitle(canvbasename+"_PullPhiTheta"); 
0163     c1->SetName(canvbasename+"_PullPhiTheta");
0164     c1->Divide(2,2);
0165     
0166     c1->cd(1);
0167 
0168     TH1F* hPull;
0169     hPull=hResPhi->hPull;
0170 
0171     //    hResPhi->hPull->Rebin(2);
0172     TF1* fphi=drawGFit(hPull, nsigma, -5, 5);
0173 
0174     c1->cd(2);
0175 
0176     plotAndProfileX(hResPhi->hPullVsPos,rbx,rby,rbp,-5, 5, 0, 2.1);
0177 
0178 
0179     m_phi = fphi->GetParameter("Mean")*cmToMicron;
0180     s_phi = fphi->GetParameter("Sigma")*cmToMicron;
0181 
0182     if (hResTheta->hPull) {
0183 
0184       c1->cd(3);
0185       hPull=hResTheta->hPull;
0186 
0187       //      hResTheta->hPull->Rebin(2);
0188       TF1* ftheta=drawGFit(hPull, nsigma, -5, 5);
0189 
0190       c1->cd(4);  
0191       plotAndProfileX(hResTheta->hPullVsPos,rbx,rby,rbp,-5, 5, 0, 2.1);  
0192     
0193       m_theta = ftheta->GetParameter("Mean");
0194       s_theta = ftheta->GetParameter("Sigma");  
0195     }
0196     
0197 
0198     cout << canvbasename << "  Step3 W" << wheel << " St" << station << endl
0199          << "   Pulls:        Phi: M= " << int(floor(m_phi+0.5))
0200          << " S= "      << int(s_phi+0.5)
0201          << "; Theta: M= " << int(floor(m_theta+0.5))
0202          << " S= "  << int(s_theta+0.5) << endl;
0203   }
0204 
0205 
0206   //-------------------- Hit efficiencies as a function of distance from wire
0207   if (doEff) {
0208     TCanvas* c1= new TCanvas;
0209     c1->SetTitle(canvbasename+"_EffPhiTheta");  
0210     c1->SetName(canvbasename+"_EffPhiTheta");
0211     c1->SetWindowSize(325,750);
0212     c1->Divide(1,2);
0213     c1->cd(1);
0214     plotEff(hEffS1RPhi->hEffVsDist, hEffS3RPhi->hEffVsDist);
0215     c1->cd(2);
0216     if (station!=4) plotEff(hEffS1RZ->hEffVsDist, hEffS3RZ->hEffVsDist);
0217 //     c1->cd(2);
0218 //     //plotEff(hEffS1RPhi->hEffVsPhi, hEffS3RPhi->hEffVsPhi);
0219 //     plotEff(hEffS1RPhi->hEffVsEta, hEffS3RPhi->hEffVsEta);
0220 //     c1->cd(4);
0221 //     //    if (station!=4) plotEff(hEffS1RZ->hEffVsPhi, hEffS3RZ->hEffVsPhi);
0222 //     if (station!=4) plotEff(hEffS1RZ->hEffVsEta, hEffS3RZ->hEffVsEta);
0223 
0224   }
0225 
0226 
0227   //-------------------- #hits, t0s
0228   if (doT0) {
0229     TCanvas* c1= new TCanvas;
0230     c1->SetTitle(canvbasename+"_NHitsT0");   
0231     c1->SetName(canvbasename+"_NHitsT0");
0232     c1->Divide(2,2);
0233     c1->cd(1);
0234     
0235     TH2F* hNh =  hRes4D->hHitMult;
0236  
0237     hNh->SetXTitle("#phi hits");
0238     hNh->SetYTitle("#theta hits");
0239     hNh->Draw("BOX");
0240     c1->cd(2);
0241     hNh->ProjectionY()->Draw();
0242     c1->cd(3);
0243     hNh->ProjectionX()->Draw();
0244     c1->cd(4);
0245    
0246     TH2F* ht =  hRes4D->ht0;
0247     ht->SetXTitle("t0 #phi");
0248     ht->SetYTitle("t0 #theta");
0249     ht->Draw("BOX");
0250 
0251 
0252   }
0253 
0254   //-------------------- Segment x, y, alpha, beta resolutions
0255   if (doSegRes){
0256     TCanvas* c1= new TCanvas;
0257     c1->SetTitle(canvbasename+"_ResSeg"); 
0258     c1->SetName(canvbasename+"_ResSeg");
0259     c1->Divide(2,2);
0260     c1->cd(1);
0261     //    hRes4D->hResX->Rebin(2);
0262     drawGFit(hRes4D->hResX, nsigma, -0.1, 0.1);
0263     c1->cd(2);
0264     //    hRes4D->hResAlpha->Rebin(2);
0265     drawGFit(hRes4D->hResAlpha, nsigma, -0.01, 0.01);
0266     c1->cd(3);
0267     hRes4D->hResYRZ->Rebin(2);
0268     drawGFit(hRes4D->hResYRZ, nsigma, -0.4, 0.4);
0269     c1->cd(4);
0270     hRes4D->hResBeta->Rebin(2);
0271     drawGFit(hRes4D->hResBeta, nsigma, -0.4, 0.4);
0272   }
0273 
0274 
0275  //-------------------- Angular dependencies
0276   if (doAngularDeps) {
0277     TCanvas* c1= new TCanvas;
0278     
0279     float min;
0280     float max;
0281     c1->SetTitle(canvbasename+"_MeanVsAngles");
0282     c1->SetName(canvbasename+"_MeanVsAngles");
0283     c1->SetTitle(canvbasename+" Angles");
0284     
0285     c1->Divide(2,2);
0286 
0287     c1->cd(1);
0288     hRes4D->hSimAlpha->SetLineColor(kGreen);
0289     hRes4D->hSimAlpha->Draw();
0290     hRes4D->hRecAlpha->Draw("same");
0291 
0292     c1->cd(2);  
0293     plotAndProfileX(hResPhi->hResVsAngle,1,1,1,-.04, .04, -1., 1.);
0294 
0295     c1->cd(3);
0296     hRes4D->hSimBetaRZ->SetLineColor(kGreen);
0297     hRes4D->hSimBetaRZ->Draw();
0298     hRes4D->hRecBetaRZ->Draw("same");
0299        
0300     c1->cd(4);
0301     plotAndProfileX(hResTheta->hResVsAngle,1,1,1,-.04, .04, -1.2.,1.2);
0302    
0303 
0304   }
0305 
0306 
0307 
0308   //------------------- Efficiencies Vs X, Y. alpha, beta
0309 
0310   if(doEff4D){
0311     // FIXME: should rebin histograms.
0312     
0313     TH1F* hEffX;
0314     TH1F* hEffY;
0315     TH1F* hEffalpha;
0316     TH1F* hEffbeta;
0317 
0318     TCanvas* c1= new TCanvas;
0319     c1->SetTitle(canvbasename+"_Efficiencies"); 
0320     c1->SetName(canvbasename+"_Efficiencies");
0321     c1->SetTitle(canvbasename+"_Efficiencies");
0322     
0323     c1->Divide(2,2);
0324     
0325     c1->cd(1);
0326     hEffX = getEffPlot(hEff4D->hXRecHit, hEff4D->hXSimSegm,2);
0327     hEffX->Draw();
0328 
0329     c1->cd(2);
0330     hEffalpha = getEffPlot(hEff4D->hAlphaRecHit, hEff4D->hAlphaSimSegm,2);
0331     hEffalpha->Draw();
0332 
0333     c1->cd(3);
0334     hEffY = getEffPlot(hEff4D->hYRecHit, hEff4D->hYSimSegm,2);
0335     hEffY->Draw();
0336 
0337     c1->cd(4);
0338     hEffbeta = getEffPlot(hEff4D->hBetaRecHit, hEff4D->hBetaSimSegm,2);
0339     hEffbeta->Draw();
0340 
0341 
0342 
0343     }
0344 
0345  //-------------------- Segment x, y, alpha, beta pull
0346   if (doSegPull){
0347     TCanvas* c1= new TCanvas;
0348     c1->SetTitle(canvbasename+"_PullSeg"); 
0349     c1->SetName(canvbasename+"_PullSeg");
0350     c1->Divide(2,2);
0351     c1->cd(1);
0352 
0353     hRes4D->hPullX->Rebin(2);
0354     drawGFit(hRes4D->hPullX, nsigma, -10.,10.);
0355     c1->cd(2);
0356     hRes4D->hPullAlpha->Rebin(2);
0357     drawGFit(hRes4D->hPullAlpha, nsigma, -10., 10.);
0358     c1->cd(3);
0359     hRes4D->hPullYRZ->Rebin(2);
0360     drawGFit(hRes4D->hPullYRZ, nsigma, -10., 10.);
0361     c1->cd(4);
0362     hRes4D->hPullBetaRZ->Rebin(2);
0363     drawGFit(hRes4D->hPullBetaRZ, nsigma, -10.,10.);
0364 
0365     //Fixme: Move these to another canvas. Note that the error used for hPullY is not computed correctly.q
0366 //     c1->cd(3);
0367 //     hRes4D->hPullY->Rebin(2);
0368 //     drawGFit(hRes4D->hPullY, nsigma, -10., 10.);
0369 //     c1->cd(4);
0370 //     hRes4D->hPullBeta->Rebin(2);
0371 //     drawGFit(hRes4D->hPullBeta, nsigma, -10.,10.);  
0372 
0373 
0374   }
0375 
0376   //-------------------- #segments
0377   if (doNSeg) {
0378     TCanvas* c1= new TCanvas;
0379     c1->SetTitle(canvbasename+"_NSeg");   
0380     c1->SetName(canvbasename+"_NSeg");
0381     // c1->Divide(2,2);
0382 //     c1->cd(1);
0383     
0384     TH1F* hNs =  hEff4D->hNSeg;
0385  
0386     hNs->SetXTitle("#segments");
0387     hNs->Draw();
0388 
0389     double int_1 = hNs->Integral(2,21);
0390     double int_2 = hNs->Integral(3,21);
0391     
0392     double ratio = int_2/int_1;
0393     
0394     cout << "int_1: " << int_1 <<" int_2: " << int_2 << " ratio: " << ratio << endl;
0395     
0396   }
0397 
0398 
0399 }