CMSSW Reference Manual

00001 #include "comparisonPlots.h"
00002 #include <string>
00003 #include <sstream>
00004 
00005 #include "TProfile.h"
00006 #include "TList.h"
00007 #include "TNtuple.h"
00008 #include "TString.h"
00009 #include <iostream>
00010 #include "TStyle.h"
00011 
00012 comparisonPlots::comparisonPlots(std::string filename, std::string outputDir, std::string outputFilename)
00013 {
00014         
00015         _outputDir = outputDir;
00016         
00017         fin = new TFile(filename.c_str());
00018         fin->cd();
00019         
00020         output = new TFile((outputDir+outputFilename).c_str(),"recreate");
00021         output->cd();
00022         
00023         readTree();
00024         
00025 }
00026 
00027 void comparisonPlots::readTree(){
00028         
00029         data = (TTree*)fin->Get("alignTree");
00030         data->SetBranchAddress("id",&id_);
00031         data->SetBranchAddress("mid",&mid_);
00032         data->SetBranchAddress("level",&level_);
00033         data->SetBranchAddress("mlevel",&mlevel_);
00034         data->SetBranchAddress("sublevel",&sublevel_);
00035         data->SetBranchAddress("x",&x_);
00036         data->SetBranchAddress("y",&y_);
00037         data->SetBranchAddress("z",&z_);
00038         data->SetBranchAddress("alpha",&alpha_);
00039         data->SetBranchAddress("beta",&beta_);
00040         data->SetBranchAddress("gamma",&gamma_);
00041         data->SetBranchAddress("phi",&phi_);
00042         data->SetBranchAddress("eta",&eta_);
00043         data->SetBranchAddress("r",&r_);
00044         data->SetBranchAddress("dx",&dx_);
00045         data->SetBranchAddress("dy",&dy_);
00046         data->SetBranchAddress("dz",&dz_);
00047         data->SetBranchAddress("dphi",&dphi_);
00048         data->SetBranchAddress("dr",&dr_);      
00049         data->SetBranchAddress("dalpha",&dalpha_);
00050         data->SetBranchAddress("dbeta",&dbeta_);
00051         data->SetBranchAddress("dgamma",&dgamma_);
00052         data->SetBranchAddress("useDetId",&useDetId_);
00053         data->SetBranchAddress("detDim",&detDim_);
00054 }
00055 
00056 void comparisonPlots::plot3x5(TCut Cut, char* dirName, bool savePlot, std::string plotName, bool autolimits){
00057         
00058         // ---------  create directory for histograms ---------
00059         //const char* dirName = Cut;
00060         TDirectory* plotDir = output->mkdir( dirName );
00061         
00062         // ---------  get right limits for histogram ---------
00063         double minimumR, maximumR;
00064         double minimumZ, maximumZ;
00065         double minimumPhi, maximumPhi;
00066         double minimumDR, maximumDR;
00067         double minimumDZ, maximumDZ;
00068         double minimumRDPhi, maximumRDPhi;
00069         double minimumDX, maximumDX;
00070         double minimumDY, maximumDY;
00071         if (autolimits){
00072         // ---------  get right limits for histogram ---------
00073         TH1F* phr = new TH1F("phr", "phr", 200, 0, 200);
00074         TH1F* phz = new TH1F("phz", "phz", 400, -200, 200);
00075         TH1F* phphi = new TH1F("phphi", "phphi", 200, -3.15, 3.15);
00076         TH1F* phdr = new TH1F("phdr", "phdr", 2000, -10, 10);
00077         TH1F* phdz = new TH1F("phdz", "phdz", 2000, -10, 10);
00078         TH1F* phrdphi = new TH1F("phrdphi", "phrdphi", 200, -10, 10);
00079         TH1F* phdx = new TH1F("phdx", "phy", 2000, -10, 10);
00080         TH1F* phdy = new TH1F("phdy", "phy", 2000, -10, 10);
00081         data->Project("phr","r",Cut);
00082         data->Project("phz","z",Cut);
00083         data->Project("phphi","phi",Cut);
00084         data->Project("phdr","dr",Cut);
00085         data->Project("phdz","dz",Cut);
00086         data->Project("phrdphi","r*dphi",Cut);
00087         data->Project("phdx","dx",Cut);
00088         data->Project("phdy","dy",Cut);
00089         getHistMaxMin(phr, maximumR, minimumR, 0);
00090         getHistMaxMin(phz, maximumZ, minimumZ, 0);
00091         getHistMaxMin(phphi, maximumPhi, minimumPhi, 0);
00092         getHistMaxMin(phdr, maximumDR, minimumDR, 1);
00093         getHistMaxMin(phdz, maximumDZ, minimumDZ, 1);
00094         getHistMaxMin(phrdphi, maximumRDPhi, minimumRDPhi, 1);
00095         getHistMaxMin(phdx, maximumDX, minimumDX, 1);
00096         getHistMaxMin(phdy, maximumDY, minimumDY, 1);
00097         }
00098         else{
00099         minimumR = 0., maximumR = 200.; 
00100         minimumZ = -200., maximumZ = 200.; 
00101         minimumPhi = -3.15, maximumPhi = 3.15;
00102         minimumDR = -1, maximumDR = 1;
00103         minimumDZ = -1, maximumDZ = 1;
00104         minimumRDPhi = -1, maximumRDPhi = 1;
00105         minimumDX = -1, maximumDX = 1;
00106         minimumDY = -1, maximumDY = 1;
00107         }
00108                 
00109         
00110         // ---------  declare histograms ---------
00111         TH1F* h_dr = new TH1F("h_dr", "#Delta r", 2000, minimumDR, maximumDR);
00112         TH1F* h_dz = new TH1F("h_dz", "#Delta z", 2000, minimumDZ, maximumDZ);
00113         TH1F* h_rdphi = new TH1F("h_rdphi", "r* #Delta #phi", 2000, minimumRDPhi, maximumRDPhi);
00114         TH2F* h_drVr = new TH2F("h_drVr","#Delta r vs. r",200,minimumR,maximumR,200,minimumDR,maximumDR);
00115         TH2F* h_dzVr = new TH2F("h_dzVr","#Delta z vs. r",200,minimumR,maximumR,200,minimumDZ,maximumDZ);
00116         TH2F* h_rdphiVr = new TH2F("h_rdphiVr","r#Delta #phi vs. r",200,minimumR,maximumR,200,minimumRDPhi,maximumRDPhi);
00117         TH2F* h_dxVr = new TH2F("h_dxVr","#Delta x vs. r", 200,minimumR,maximumR, 200,minimumDX,maximumDX);
00118         TH2F* h_dyVr = new TH2F("h_dyVr","#Delta y vs. r", 200,minimumR,maximumR, 200,minimumDY,maximumDY);
00119         TH2F* h_drVz = new TH2F("h_drVz","#Delta r vs. z", 200,minimumZ,maximumZ, 200,minimumDR,maximumDR);
00120         TH2F* h_dzVz = new TH2F("h_dzVz","#Delta z vs. z", 200,minimumZ,maximumZ, 200,minimumDZ,maximumDZ);
00121         TH2F* h_rdphiVz = new TH2F("h_rdphiVz","r#Delta #phi vs. z", 200,minimumZ,maximumZ, 200,minimumRDPhi,maximumRDPhi);
00122         TH2F* h_dxVz = new TH2F("h_dxVz","#Delta x vs. z", 200,minimumZ,maximumZ, 200,minimumDX,maximumDX);
00123         TH2F* h_dyVz = new TH2F("h_dyVz","#Delta y vs. z", 200,minimumZ,maximumZ, 200,minimumDY,maximumDY);
00124         TH2F* h_drVphi = new TH2F("h_drVphi","#Delta r vs. #phi", 200,minimumPhi,maximumPhi,200,minimumDR,maximumDR);
00125         TH2F* h_dzVphi = new TH2F("h_dzVphi","#Delta z vs. #phi", 200,minimumPhi,maximumPhi, 200,minimumDZ,maximumDZ);
00126         TH2F* h_rdphiVphi = new TH2F("h_rdphiVphi","r#Delta #phi vs. #phi", 200,minimumPhi,maximumPhi,200,minimumRDPhi,maximumRDPhi);
00127         TH2F* h_dxVphi = new TH2F("h_dxVphi","#Delta x vs. #phi", 200,minimumPhi,maximumPhi, 200,minimumDX,maximumDX);
00128         TH2F* h_dyVphi = new TH2F("h_dyVphi","#Delta y vs. #phi", 200,minimumPhi,maximumPhi, 200,minimumDY,maximumDY);
00129         
00130         // ---------  project tree onto histograms ---------
00131         data->Project("h_dr","dr",Cut);
00132         data->Project("h_dz","dz",Cut);
00133         data->Project("h_rdphi","r*dphi",Cut);
00134         data->Project("h_drVr", "dr:r",Cut);
00135         data->Project("h_dzVr", "dz:r",Cut);
00136         data->Project("h_rdphiVr", "r*dphi:r",Cut);
00137         data->Project("h_dxVr", "dx:r",Cut);
00138         data->Project("h_dyVr", "dy:r",Cut);
00139         data->Project("h_drVz", "dr:z",Cut);
00140         data->Project("h_dzVz", "dz:z",Cut);
00141         data->Project("h_rdphiVz", "r*dphi:z",Cut);
00142         data->Project("h_dxVz", "dx:z",Cut);
00143         data->Project("h_dyVz", "dy:z",Cut);
00144         data->Project("h_drVphi", "dr:phi",Cut);
00145         data->Project("h_dzVphi", "dz:phi",Cut);
00146         data->Project("h_rdphiVphi", "r*dphi:phi",Cut);
00147         data->Project("h_dxVphi", "dx:phi",Cut);
00148         data->Project("h_dyVphi", "dy:phi",Cut);
00149          
00150         
00151         // ---------  draw histograms ---------
00152         TCanvas* c0 = new TCanvas("c0", "c0", 200, 10, 900, 300);
00153         c0->SetFillColor(0);
00154         c0->Divide(3,1);
00155         c0->cd(1);
00156         h_dr->Draw();
00157         c0->cd(2);
00158         h_dz->Draw();
00159         c0->cd(3);
00160         h_rdphi->Draw();
00161         if (savePlot) c0->Print((_outputDir+"plot3x1_"+plotName).c_str());
00162 
00163         
00164         // ---------  draw histograms ---------
00165         TCanvas* c = new TCanvas("c", "c", 200, 10, 1200, 700);
00166         c->SetFillColor(0);
00167         data->SetMarkerSize(0.5);
00168         data->SetMarkerStyle(6);
00169         c->Divide(5,3);
00170         c->cd(1);
00171         h_drVr->Draw();
00172         c->cd(2);
00173         h_dzVr->Draw();
00174         c->cd(3);
00175         h_rdphiVr->Draw();
00176         c->cd(4);
00177         h_dxVr->Draw();
00178         c->cd(5);
00179         h_dyVr->Draw();
00180         c->cd(6);
00181         h_drVz->Draw();
00182         c->cd(7);
00183         h_dzVz->Draw();
00184         c->cd(8);
00185         h_rdphiVz->Draw();
00186         c->cd(9);
00187         h_dxVz->Draw();
00188         c->cd(10);
00189         h_dyVz->Draw();
00190         c->cd(11);
00191         h_drVphi->Draw();
00192         c->cd(12);
00193         h_dzVphi->Draw();
00194         c->cd(13);
00195         h_rdphiVphi->Draw();
00196         c->cd(14);
00197         h_dxVphi->Draw();
00198         c->cd(15);
00199         h_dyVphi->Draw();
00200         
00201         // ---------  set output directory for histograms ---------
00202         plotDir->cd();
00203         h_dr->Write(); h_dz->Write(); h_rdphi->Write();
00204         h_drVr->Write(); h_dzVr->Write(); h_rdphiVr->Write(); h_dxVr->Write(); h_dyVr->Write();
00205         h_drVz->Write(); h_dzVz->Write(); h_rdphiVz->Write(); h_dxVz->Write(); h_dyVz->Write();
00206         h_drVphi->Write(); h_dzVphi->Write(); h_rdphiVphi->Write(); h_dxVphi->Write(); h_dyVphi->Write();
00207         
00208         if (savePlot) c->Print((_outputDir+"plot3x5_"+plotName).c_str());
00209         
00210 }
00211 
00212 void comparisonPlots::plot3x5Profile(TCut Cut, char* dirName, int nBins, bool savePlot, std::string plotName, bool autolimits){
00213         
00214         // ---------  create directory for histograms ---------
00215         //const char* dirName = Cut;
00216         string s;// = "profile";
00217         s = s + dirName;
00218         s.append("_profile");
00219         TDirectory* plotDir = output->mkdir( s.data() );
00220 
00221 
00222         double minimumR, maximumR;
00223         double minimumZ, maximumZ;
00224         double minimumPhi, maximumPhi;
00225         double minimumDR, maximumDR;
00226         double minimumDZ, maximumDZ;
00227         double minimumRDPhi, maximumRDPhi;
00228         double minimumDX, maximumDX;
00229         double minimumDY, maximumDY;
00230         if (autolimits){
00231         // ---------  get right limits for histogram ---------
00232         TH1F* phr = new TH1F("phr", "phr", 200, 0, 200);
00233         TH1F* phz = new TH1F("phz", "phz", 400, -200, 200);
00234         TH1F* phphi = new TH1F("phphi", "phphi", 200, -3.15, 3.15);
00235         TH1F* phdr = new TH1F("phdr", "phdr", 2000, -10, 10);
00236         TH1F* phdz = new TH1F("phdz", "phdz", 2000, -10, 10);
00237         TH1F* phrdphi = new TH1F("phrdphi", "phrdphi", 200, -10, 10);
00238         TH1F* phdx = new TH1F("phdx", "phy", 2000, -10, 10);
00239         TH1F* phdy = new TH1F("phdy", "phy", 2000, -10, 10);
00240         data->Project("phr","r",Cut);
00241         data->Project("phz","z",Cut);
00242         data->Project("phphi","phi",Cut);
00243         data->Project("phdr","dr",Cut);
00244         data->Project("phdz","dz",Cut);
00245         data->Project("phrdphi","r*dphi",Cut);
00246         data->Project("phdx","dx",Cut);
00247         data->Project("phdy","dy",Cut);
00248         getHistMaxMin(phr, maximumR, minimumR, 0);
00249         getHistMaxMin(phz, maximumZ, minimumZ, 0);
00250         getHistMaxMin(phphi, maximumPhi, minimumPhi, 0);
00251         getHistMaxMin(phdr, maximumDR, minimumDR, 1);
00252         getHistMaxMin(phdz, maximumDZ, minimumDZ, 1);
00253         getHistMaxMin(phrdphi, maximumRDPhi, minimumRDPhi, 1);
00254         getHistMaxMin(phdx, maximumDX, minimumDX, 1);
00255         getHistMaxMin(phdy, maximumDY, minimumDY, 1);
00256         }
00257         else{
00258         minimumR = 0., maximumR = 200.; 
00259         minimumZ = -200., maximumZ = 200.; 
00260         minimumPhi = -3.15, maximumPhi = 3.15;
00261         minimumDR = -1, maximumDR = 1;
00262         minimumDZ = -1, maximumDZ = 1;
00263         minimumRDPhi = -1, maximumRDPhi = 1;
00264         minimumDX = -1, maximumDX = 1;
00265         minimumDY = -1, maximumDY = 1;
00266         }
00267         
00268         // ---------  declare histograms ---------
00269         TProfile* hprof_drVr = new TProfile("hprof_drVr","#Delta r vs. r",nBins,minimumR,maximumR,minimumDR,maximumDR);
00270         TProfile* hprof_dzVr = new TProfile("hprof_dzVr","#Delta z vs. r",nBins,minimumR,maximumR,minimumDZ,maximumDZ);
00271         TProfile* hprof_rdphiVr = new TProfile("hprof_rdphiVr","r#Delta #phi vs. r",nBins,minimumR,maximumR,minimumRDPhi,maximumRDPhi);
00272         TProfile* hprof_dxVr = new TProfile("hprof_dxVr","#Delta x vs. r", nBins,minimumR,maximumR,minimumDX,maximumDX);
00273         TProfile* hprof_dyVr = new TProfile("hprof_dyVr","#Delta y vs. r", nBins,minimumR,maximumR,minimumDY,maximumDY);
00274         TProfile* hprof_drVz = new TProfile("hprof_drVz","#Delta r vs. z", nBins,minimumZ,maximumZ,minimumDR,maximumDR);
00275         TProfile* hprof_dzVz = new TProfile("hprof_dzVz","#Delta z vs. z", nBins,minimumZ,maximumZ,minimumDZ,maximumDZ);
00276         TProfile* hprof_rdphiVz = new TProfile("hprof_rdphiVz","r#Delta #phi vs. z", nBins,minimumZ,maximumZ,minimumRDPhi,maximumRDPhi);
00277         TProfile* hprof_dxVz = new TProfile("hprof_dxVz","#Delta x vs. z", nBins,minimumZ,maximumZ,minimumDX,maximumDX);
00278         TProfile* hprof_dyVz = new TProfile("hprof_dyVz","#Delta y vs. z", nBins,minimumZ,maximumZ,minimumDY,maximumDY);
00279         TProfile* hprof_drVphi = new TProfile("hprof_drVphi","#Delta r vs. #phi", nBins,minimumPhi,maximumPhi,minimumDR,maximumDR);
00280         TProfile* hprof_dzVphi = new TProfile("hprof_dzVphi","#Delta z vs. #phi", nBins,minimumPhi,maximumPhi,minimumDZ,maximumDZ);
00281         TProfile* hprof_rdphiVphi = new TProfile("hprof_rdphiVphi","r#Delta #phi vs. #phi", nBins,minimumPhi,maximumPhi,minimumRDPhi,maximumRDPhi);
00282         TProfile* hprof_dxVphi = new TProfile("hprof_dxVphi","#Delta x vs. #phi", nBins,minimumPhi,maximumPhi,minimumDX,maximumDX);
00283         TProfile* hprof_dyVphi = new TProfile("hprof_dyVphi","#Delta y vs. #phi", nBins,minimumPhi,maximumPhi,minimumDY,maximumDY);
00284         
00285         // ---------  project tree onto histograms ---------
00286         data->Project("hprof_drVr", "dr:r",Cut,"prof");
00287         data->Project("hprof_dzVr", "dz:r",Cut,"prof");
00288         data->Project("hprof_rdphiVr", "r*dphi:r",Cut,"prof");
00289         data->Project("hprof_dxVr", "dx:r",Cut,"prof");
00290         data->Project("hprof_dyVr", "dy:r",Cut,"prof");
00291         data->Project("hprof_drVz", "dr:z",Cut,"prof");
00292         data->Project("hprof_dzVz", "dz:z",Cut,"prof");
00293         data->Project("hprof_rdphiVz", "r*dphi:z",Cut,"prof");
00294         data->Project("hprof_dxVz", "dx:z",Cut,"prof");
00295         data->Project("hprof_dyVz", "dy:z",Cut,"prof");
00296         data->Project("hprof_drVphi", "dr:phi",Cut,"prof");
00297         data->Project("hprof_dzVphi", "dz:phi",Cut,"prof");
00298         data->Project("hprof_rdphiVphi", "r*dphi:phi",Cut,"prof");
00299         data->Project("hprof_dxVphi", "dx:phi",Cut,"prof");
00300         data->Project("hprof_dyVphi", "dy:phi",Cut,"prof");
00301         
00302         // ---------  draw histograms ---------
00303         TCanvas* cp = new TCanvas("cp", "cp", 200, 10, 1200, 700);
00304         cp->SetFillColor(0);
00305         data->SetMarkerSize(0.5);
00306         data->SetMarkerStyle(6);
00307         cp->Divide(5,3);
00308         cp->cd(1);
00309         hprof_drVr->Draw();
00310         cp->cd(2);
00311         hprof_dzVr->Draw();
00312         cp->cd(3);
00313         hprof_rdphiVr->Draw();
00314         cp->cd(4);
00315         hprof_dxVr->Draw();
00316         cp->cd(5);
00317         hprof_dyVr->Draw();
00318         cp->cd(6);
00319         hprof_drVz->Draw();
00320         cp->cd(7);
00321         hprof_dzVz->Draw();
00322         cp->cd(8);
00323         hprof_rdphiVz->Draw();
00324         cp->cd(9);
00325         hprof_dxVz->Draw();
00326         cp->cd(10);
00327         hprof_dyVz->Draw();
00328         cp->cd(11);
00329         hprof_drVphi->Draw();
00330         cp->cd(12);
00331         hprof_dzVphi->Draw();
00332         cp->cd(13);
00333         hprof_rdphiVphi->Draw();
00334         cp->cd(14);
00335         hprof_dxVphi->Draw();
00336         cp->cd(15);
00337         hprof_dyVphi->Draw();
00338         
00339         // ---------  set output directory for histograms ---------
00340         plotDir->cd();
00341         hprof_drVr->Write(); hprof_dzVr->Write(); hprof_rdphiVr->Write(); hprof_dxVr->Write(); hprof_dyVr->Write();
00342         hprof_drVz->Write(); hprof_dzVz->Write(); hprof_rdphiVz->Write(); hprof_dxVz->Write(); hprof_dyVz->Write();
00343         hprof_drVphi->Write(); hprof_dzVphi->Write(); hprof_rdphiVphi->Write(); hprof_dxVphi->Write(); hprof_dyVphi->Write();
00344         
00345         if (savePlot) cp->Print((_outputDir+"plot3x5Prof_"+plotName).c_str());
00346 }
00347 
00348 
00349 void comparisonPlots::getMaxMin(){
00350         
00351         data->GetEntry(0);
00352         minR = r_; maxR = r_;
00353         minZ = z_; maxZ = z_;
00354         minPhi = phi_; maxPhi = phi_;
00355         minDR = dr_; maxDR = dr_;
00356         minDZ = dz_; maxDZ = dz_;
00357         minRDPhi = r_*dphi_; maxRDPhi = r_*dphi_;
00358         minDX = dx_; maxDX = dx_;
00359         minDY = dy_; maxDY = dy_;
00360         
00361         int nEntries = data->GetEntries();
00362         for (int i = 1; i < nEntries; ++i){
00363                 data->GetEntry(i);
00364                 
00365                 if (r_ < minR) minR = r_;
00366                 if (r_ > maxR) maxR = r_;
00367                 if (z_ < minZ) minZ = z_;
00368                 if (z_ > maxZ) maxZ = z_;
00369                 if (phi_ < minPhi) minPhi = phi_;
00370                 if (phi_ > maxPhi) maxPhi = phi_;
00371                 if (dr_ < minDR) minDR = dr_;
00372                 if (dr_ > maxDR) maxDR = dr_;
00373                 if (dz_ < minDZ) minDZ = dz_;
00374                 if (dz_ > maxDZ) maxDZ = dz_;
00375                 if (r_*dphi_ < minRDPhi) minRDPhi = r_*dphi_;
00376                 if (r_*dphi_ > maxRDPhi) maxRDPhi = r_*dphi_;
00377                 if (dx_ < minDX) minDX = dx_;
00378                 if (dx_ > maxDX) maxDX = dx_;
00379                 if (dy_ < minDY) minDY = dy_;
00380                 if (dy_ > maxDY) maxDY = dy_;
00381         }
00382 }
00383         
00384 
00385 void comparisonPlots::getHistMaxMin(TH1* hist, double &max, double &min, int flag){
00386         
00387         int nBins = hist->GetNbinsX();
00388         for (int i = 0; i < nBins; ++i){
00389                 double binContent = hist->GetBinContent(i);
00390                 if (binContent > 0){
00391                         //double binWidth = hist->GetBinLowEdge(i) - hist->GetBinLowEdge(i-1);
00392                         //std::cout << "bin width1: " << hist->GetBinWidth(i) << ", bin width2: " << binWidth << std::endl;
00393                         if (flag == 0) max = hist->GetBinLowEdge(i) + 2.*hist->GetBinWidth(i);
00394                         if (flag == 1) max = hist->GetBinLowEdge(i) + hist->GetBinWidth(i);
00395                 }
00396         }
00397         for (int i = (nBins-1); i >= 0; i--){
00398                 double binContent = hist->GetBinContent(i);
00399                 if (binContent > 0) min = hist->GetBinLowEdge(i);
00400         }
00401         //std::cout << "max: " << max << ", min: " << min << std::endl;
00402 }
00403 
00404 void comparisonPlots::Write()
00405 {
00406         output->Write();
00407 }
00408 
00409
comparisonPlots.cc
