SiPixelTools/rootscripts/make_scurves.C

0001 #include <iostream>
0002 #include <vector>
0003 #include "TKey.h"
0004 #include "TMath.h"
0005
0006 Double_t mycurvefunc(Double_t *x, Double_t *par){
0007   double xval = x[0]-par[0];
0008   xval /=sqrt(2.)*par[1];
0009   double val=0.5+0.5*TMath::Erf(xval);
0010   return val;
0011 }
0012 void make_scurves(){
0013
0014   TFile *file = new TFile("SCurve_DQM_Calibrationsmall.root");
0015   file->cd();
0016   // make a loop over all plots
0017   TList *list = file->GetListOfKeys();
0018   gStyle->SetOptStat(0);
0019   Int_t nkeys = file->GetNkeys();
0020   TDirectory *dir = file->GetDirectory("DQMData");
0021   if(dir==0)
0022     return;
0023
0024
0025   TLatex CMSprelim(0.3,0.3,"CMS Preliminary");
0026   CMSprelim.SetTextColor(2);
0027   CMSprelim.SetTextSize(CMSprelim.GetTextSize()*1.2);
0028
0029   TString comparestring = "Module";
0030   TString curvestring = "row";
0031   std::vector<TString> keylist;
0032   std::vector<TString> hist1list;
0033   std::vector<TString> hist2list;
0034   std::vector<TString> dirlist;
0035   std::vector<TString> notdonelist;
0036   std::vector<int> nsubdirs;
0037   TDirectory *dirsav = dir;
0038   list = dir->GetListOfKeys();
0039   int ikey=0;
0040   int localkey=0;
0041   int ntimes=0;
0042
0043   TCanvas *curvecanvas = new TCanvas();
0044   for(ikey=0;ikey<list->GetEntries();  ikey++){
0045     TKey *thekey = (TKey*)list->At(ikey);
0046     if(thekey==0)
0047       continue;
0048     TString keyname=thekey->GetName();
0049     //    keyname.ReplaceAll(" ","");
0050     TString keytype=thekey->GetClassName();
0051     keytype.ReplaceAll(" ","");
0052     if(keyname=="EventInfo")
0053       continue;
0054     //    std::cout <<  keytype << " " << keyname << std::endl;
0055     if(keytype=="TDirectoryFile"){
0056       TString dirname=dir->GetPath();
0057       dirname+="/";
0058       dirname+=keyname;
0059       //      std::cout << dirname << std::endl;
0060       dir=file->GetDirectory(dirname);
0061
0062       list=dir->GetListOfKeys();
0063       if(dirname.Contains(comparestring))
0064     dirlist.push_back(dirname);
0065       else{
0066     notdonelist.push_back(dirname);
0067     nsubdirs.push_back(-1);
0068       }
0069     }
0070   }
0071   int nempty=0;
0072   while(nempty!=notdonelist.size()){
0073     for(int idir=0; idir<notdonelist.size(); ++idir){
0074       if(nsubdirs[idir]==0)
0075     continue;
0076       //      std::cout << "now examining " << notdonelist[idir]<< " " << nsubdirs[idir] <<  std::endl;
0077       dir = file->GetDirectory(notdonelist[idir]);
0078       //      std::cout << dir->GetName() << std::endl;
0079       list= dir->GetListOfKeys();
0080       //      std::cout << list->GetEntries() << std::endl;
0081       int ndirectories=0;
0082       for(ikey=0;ikey<list->GetEntries();  ikey++){
0083     TKey *thekey = (TKey*)list->At(ikey);
0084     if(thekey==0)
0085       continue;
0086     TString keyname=thekey->GetName();
0087     //  keyname.ReplaceAll(" ","");
0088     TString keytype=thekey->GetClassName();
0089     keytype.ReplaceAll(" ","");
0090     if(keytype=="TDirectoryFile"){
0091       TString dirname=dir->GetPath();
0092       dirname+="/";
0093       dirname+=keyname;
0094       //      std::cout << dirname << std::endl;
0095       ndirectories++;
0096       if(dirname.Contains(comparestring))
0097         dirlist.push_back(dirname);
0098       else{
0099         notdonelist.push_back(dirname);
0100         nsubdirs.push_back(-1);
0101       }
0102     }
0103       }
0104       nsubdirs[idir]=ndirectories;
0105       // std::cout << "now done examining " << notdonelist[idir]<< " " << nsubdirs[idir] <<  std::endl;
0106     }
0107     // count number of done dirs;
0108     nempty=0;
0109     for(int i=0; i<nsubdirs.size(); i++){
0110       if(nsubdirs[i]!=-1)
0111     nempty++;
0112     }
0113   }
0114   gStyle->SetOptStat(0);
0115   int ncurves=0;
0116   for(int i=0; i<dirlist.size() ; ++i){
0117     //    std::cout << "good dir "  << dirlist[i] << std::endl;
0118     // now count histograms:
0119     dir = file->GetDirectory(dirlist[i]);
0120     list= dir->GetListOfKeys();
0121     //      std::cout << list->GetEntries() << std::endl;
0122     for(ikey=0;ikey<list->GetEntries();  ikey++){
0123       //      std::cout << ikey << std::endl;
0124       TKey *thekey = (TKey*)list->At(ikey);
0125       if(thekey==0)
0126     continue;
0127       TString keyname=thekey->GetName();
0128       keyname.ReplaceAll(" ","");
0129       TString keytype=thekey->GetClassName();
0130       keytype.ReplaceAll(" ","");
0131       //      std::cout << ikey << " " <<  keyname << std::endl;
0132
0133       if(keytype=="TH1F"){
0134         std::cout << keyname << " " << keytype << std::endl;
0135     if(keyname.Contains(curvestring)){
0136       std::cout << keyname << " " << keytype << std::endl;
0137       dir=file->GetDirectory(dirlist[i]);
0138       TH1F* temp = (TH1F*) dir->Get(keyname);
0139       ncurves++;
0140       std::cout << keyname << std::endl;
0141       curvecanvas->cd();
0142       curvecanvas->Clear();
0143       dir=file->GetDirectory(dirlist[i]);
0144       dir->cd();
0145       gStyle->SetOptStat(0);
0146       gStyle->SetOptFit(0);
0147       temp->SetMarkerSize(2.0);
0148       temp->GetXaxis()->SetTitle("VCAL");
0149       temp->GetYaxis()->SetTitle("Efficiency");
0150       temp->DrawClone();
0151       float gainandped[2];
0152       getgainresult(keyname,file, dirlist[i],gainandped);
0153       std::cout << gainandped[0] << " " << gainandped[1] << std::endl;
0154       if(gainandped[0]>0. && gainandped[1]>0.){
0155         TString funcname = keyname;
0156         funcname+="func";
0157         TF1 *func = new TF1(funcname,mycurvefunc,0,256,2);
0158         func->SetLineWidth(3);
0159         func->SetParameter(0,gainandped[1]);
0160         func->SetParameter(1,gainandped[0]);
0161         func->DrawClone("lsame");
0162         CMSprelim.DrawText(gainandped[1]*1.05,temp->GetMinimum()+0.01,"CMS Pixel");
0163         TLegend *legend = new TLegend(0.45,0.2,0.89,0.28);
0164         //      legend->SetFillStyle(0);
0165         legend->SetFillColor(0);
0166         legend->SetBorderSize(0);
0167         TString fitstring = "y=0.5+0.5 Erf(#frac{x-";
0168         fitstring+= gainandped[1];
0169         fitstring+="}{" ;
0170         fitstring+= gainandped[0];
0171         fitstring+=" #sqrt{2}})";
0172         legend->AddEntry(func,fitstring,"l");
0173         legend->DrawClone("same");
0174
0175       }
0176       curvecanvas->Update();
0177       curvecanvas->Print(keyname+".jpg");
0178       curvecanvas->Print(keyname+".eps");
0179     }
0180       }
0181     }
0182   }
0183   //  for(int i=0; i<notdonelist.size();++i){
0184   //    std::cout << notdonelist[i] << std::endl;
0185   //  }
0186
0187
0188 }
0189
0190
0191 float getgainresult(TString keyname, TFile *file, TString dirname, float *gainandped){
0192
0193   //  std::cout << "examining histo " << keyname << std::endl;
0194   std::string strversion = keyname;
0195   std::string cutstring[20];
0196   int foundunderscore=0;
0197   for(int i=0; i<20 && foundunderscore<strversion.size(); i++){
0198     int found=strversion.find("_",foundunderscore);
0199     std::string thesubstr = strversion.substr(foundunderscore,found-foundunderscore);
0200     //std::cout << thesubstr << " " << found << std::endl;
0201     cutstring[i]=thesubstr;
0202     foundunderscore=found+1;
0203   }
0204
0205   int row = atoi(cutstring[2].c_str());
0206   int col = atoi(cutstring[4].c_str());
0207   int detid = atoi(cutstring[6].c_str());
0208   //  std::cout << row << " " << col << " " << detid << std::endl;
0209   TString gainhisto = dirname;
0210   gainhisto+="/ScurveSigmas_siPixelCalibDigis_";
0211   gainhisto+=detid;
0212   TString pedhisto =  dirname;
0213   pedhisto+="/ScurveThresholds_siPixelCalibDigis_";
0214   pedhisto+=detid;
0215   //  std::cout << gainhisto << " " << pedhisto << std::endl;
0216   file->cd();
0217   TH2F *gain2d = file->Get(gainhisto);
0218   float gainval=0;
0219   float pedval=0;
0220   if(gain2d!=0)
0221     gainval = gain2d->GetBinContent(col,row);
0222   else return 0;
0223   TH2F *ped2d = file->Get(pedhisto);
0224   if(ped2d)
0225     pedval=ped2d->GetBinContent(col,row);
0226   else return 0;
0227   //std::cout << ped2d->GetNbinsX() << " " << ped2d->GetNbinsY() << std::endl;
0228
0229   //std::cout << pedhisto << " " << gainhisto << " " << col+1 << " " << row+1 << " " <<  ped2d->GetBinContent(col+1,row+1) << " " << gain2d->GetBinContent(col+1,row+1) <<  std::endl;
0230   gainval = gain2d->GetBinContent(col,row);
0231   pedval = ped2d->GetBinContent(col,row);
0232   //  std::cout << gainhisto << " " << row << " " << col << " " << gainval << " " << pedval << std::endl;
0233   gainandped[0]=gainval;
0234   gainandped[1]=pedval;
0235   return 1;
0236 }
0237