UsefullCode/PayloadFromASCIIfile/PredictedVsObserved.C

0001
0002
0003 #include "TROOT.h"
0004 #include "TFile.h"
0005 #include "TDirectory.h"
0006 #include "TChain.h"
0007 #include "TObject.h"
0008 #include "TTree.h"
0009 #include "TROOT.h"
0010 #include "TFile.h"
0011 #include "TDirectory.h"
0012 #include "TChain.h"
0013 #include "TObject.h"
0014 #include "TCanvas.h"
0015 #include "TMath.h"
0016 #include "TLegend.h"
0017 #include "TGraph.h"
0018 #include "TH1.h"
0019 #include "TH2.h"
0020 #include "TH3.h"
0021 #include "TTree.h"
0022 #include "TF1.h"
0023 #include "TPaveText.h"
0024
0025 #include<vector>
0026 #include<tdrstyle.C>
0027
0028 #if !defined(__CINT__) && !defined(__MAKECINT__)
0029 #include "DataFormats/FWLite/interface/Handle.h"
0030 #include "DataFormats/FWLite/interface/Event.h"
0031 #include "DataFormats/FWLite/interface/ChainEvent.h"
0032
0033 #include "CommonTools/TrackerMap/interface/TrackerMap.h"
0034
0035 using namespace fwlite;
0036 using namespace std;
0037 using namespace edm;
0038 #endif
0039
0040
0041
0042
0043 void tkMap(string TextToPrint_="CMS Preliminary 2015"){
0044
0045 }
0046
0047
0048 void PredictedVsObserved(){
0049
0050    //read gain modified from a txt file
0051    std::map<unsigned int, double> gainModifiers;
0052    FILE* inFile = fopen("GainRatio.txt","r");
0053    char line[4096];
0054    while(fgets(line, 4096, inFile)!=NULL){
0055       unsigned int detId;  double modifier;
0056       sscanf(line, "%u %lf", &detId, &modifier);
0057 //      printf("(%i , %f) \t", detId, modifier);
0058       gainModifiers[detId] = modifier;
0059    }
0060    fclose(inFile);
0061
0062
0063    string inputFileTree = "../../Data_Run_251252_to_251252_PCL/Gains_Tree.root";
0064    string moduleName = "SiStripCalib";
0065
0066    unsigned int  tree1_Index;
0067    unsigned int  tree1_DetId;
0068    unsigned char tree1_APVId;
0069    double        tree1_Gain;
0070    double        tree1_PrevGain;
0071    double        tree1_NEntries;
0072
0073    TFile* f1     = new TFile(inputFileTree.c_str());
0074    TTree *t1     = (TTree*)f1->Get((moduleName+"/APVGain").c_str());
0075
0076
0077    t1->SetBranchAddress("Index"             ,&tree1_Index      );
0078    t1->SetBranchAddress("DetId"             ,&tree1_DetId      );
0079    t1->SetBranchAddress("APVId"             ,&tree1_APVId      );
0080    t1->SetBranchAddress("Gain"              ,&tree1_Gain       );
0081    t1->SetBranchAddress("PrevGain"          ,&tree1_PrevGain   );
0082    t1->SetBranchAddress("NEntries"          ,&tree1_NEntries   );
0083
0084
0085
0086    string inputFileTree2 = "../../Data_Run_252126_to_252226_CalibTree/Gains_Tree.root";
0087
0088    unsigned int  tree2_Index;
0089    unsigned int  tree2_DetId;
0090    unsigned char tree2_APVId;
0091    double        tree2_Gain;
0092    double        tree2_PrevGain;
0093    double        tree2_NEntries;
0094
0095
0096    TFile* f2     = new TFile(inputFileTree2.c_str());
0097    TTree *t2     = (TTree*)f2->Get((moduleName+"/APVGain").c_str());
0098
0099
0100    t2->SetBranchAddress("Index"             ,&tree2_Index      );
0101    t2->SetBranchAddress("DetId"             ,&tree2_DetId      );
0102    t2->SetBranchAddress("APVId"             ,&tree2_APVId      );
0103    t2->SetBranchAddress("Gain"              ,&tree2_Gain       );
0104    t2->SetBranchAddress("PrevGain"          ,&tree2_PrevGain   );
0105    t2->SetBranchAddress("NEntries"          ,&tree2_NEntries   );
0106
0107
0108
0109    TH1D*  HRatio = new TH1D("Ratio", "Ratio", 150, 0.5, 1.5);
0110    TH1D*  HRatioB = new TH1D("RatioB", "RatioB", 150, 0.5, 1.5);
0111    TH1D*  HRatioC = new TH1D("RatioC", "RatioC", 150, 0.5, 1.5);
0112    TH2D*  HVS    = new TH2D("VS", "VS", 150, 0.5, 1.5, 150, 0.5, 1.5);
0113
0114
0115    TrackerMap* tkmap = new TrackerMap("  Gain : Measured(R252116+R252226) / Predicted(R251252*StrenghRatio from DelayScan)");
0116
0117    FILE* pFile = fopen("Gains_ASCIIObserved.txt","w");
0118
0119    int SAME=0;
0120    int DIFF=0;
0121
0122    int PreviousId = -1;
0123    printf("Progressing Bar              :0%%       20%%       40%%       60%%       80%%       100%%\n");
0124    printf("Looping on the Tree          :");
0125    int TreeStep = t1->GetEntries()/50;if(TreeStep==0)TreeStep=1;
0126    std::vector<double> gains;
0127    double MEAN = 0;  int NAPV=0;
0128    for (unsigned int ientry = 0; ientry < t1->GetEntries(); ientry++) {
0129       if(ientry%TreeStep==0){printf(".");fflush(stdout);}
0130       t1->GetEntry(ientry);
0131       t2->GetEntry(ientry);
0132
0133 //     if(tree1_DetId==369153064){
0134      if(tree1_DetId==402673324){
0135          printf("%i | %i : %f --> %f  (%f)\n", tree1_DetId, tree1_APVId, tree1_PrevGain, tree1_Gain, tree1_NEntries);
0136          printf("%i | %i : %f --> %f  (%f)\n", tree2_DetId, tree2_APVId, tree2_PrevGain, tree2_Gain, tree2_NEntries);
0137          printf("\n");
0138       }
0139
0140
0141        if(fabs(tree2_PrevGain - tree1_Gain) > 0.05*tree2_PrevGain){
0142 //         printf("%i | %i : %f != %f | %i | %i\n", tree1_DetId, tree1_APVId, tree1_Gain, tree2_PrevGain, tree2_APVId, tree2_DetId);
0143          DIFF++;
0144       }else{
0145          SAME++;
0146       }
0147
0148
0149       if((tree1_APVId==0 && PreviousId>0) || ientry==t1->GetEntries()-1){  //new module or very last APV
0150          fprintf(pFile,"%i ",PreviousId); for(unsigned int i=0;i<gains.size();i++){fprintf(pFile,"%f ", gains[i]);}   fprintf(pFile, "\n");
0151
0152           if(NAPV>0)tkmap->fill(tree1_DetId, MEAN/NAPV);
0153           MEAN = 0; NAPV=0;
0154          gains.clear();
0155       }
0156       PreviousId = tree1_DetId;
0157
0158       if(tree2_Gain==1.0000)continue;
0159
0160       double modifier = 1.0;
0161       if(gainModifiers.find(tree1_DetId)!=gainModifiers.end()){modifier = gainModifiers[tree1_DetId]; }
0162 //      gains.push_back(tree1_Gain * modifier);
0163       gains.push_back(tree2_Gain / (tree1_Gain * modifier) );
0164
0165
0166       HRatio->Fill(tree2_Gain / (tree1_Gain * modifier));
0167       if((fabs(tree2_PrevGain-tree1_Gain*modifier)/tree2_PrevGain)>0.03)  HRatioB->Fill(tree2_Gain / (tree1_Gain * modifier));
0168       if((fabs(tree2_PrevGain-tree1_Gain*modifier)/tree2_PrevGain)<0.03){
0169          HRatioC->Fill(tree2_Gain / (tree1_Gain * modifier));
0170          MEAN += tree2_Gain / (tree1_Gain * modifier);
0171          NAPV++;
0172          HVS   ->Fill((tree1_Gain * modifier), tree2_Gain);
0173       }
0174
0175    }printf("\n");
0176    fclose(pFile);
0177
0178    printf("%i (same) + %i (diff) = %i\n", SAME, DIFF, SAME+DIFF);
0179
0180 //   tkmap->setTitle("Gain: Predicted/Current    (white=unmodified)");
0181    tkmap->save(true, 0.95, 1.1, "tkMapPredictedVsObserved.png");
0182    tkmap->reset();
0183
0184
0185
0186
0187    gROOT->Reset();
0188    setTDRStyle();
0189    gStyle->SetPadTopMargin   (0.05);
0190    gStyle->SetPadBottomMargin(0.10);
0191    gStyle->SetPadRightMargin (0.18);
0192    gStyle->SetPadLeftMargin  (0.13);
0193    gStyle->SetTitleSize(0.04, "XYZ");
0194    gStyle->SetTitleXOffset(1.1);
0195    gStyle->SetTitleYOffset(1.35);
0196    gStyle->SetPalette(1);
0197    gStyle->SetCanvasColor(0);
0198    gStyle->SetBarOffset(0);
0199
0200
0201
0202
0203
0204
0205 if(true){
0206    TCanvas* c1 = new TCanvas("c1","c1,",900,600);
0207
0208    HRatio->SetLineWidth(2);
0209    HRatio->SetLineColor(1);
0210    HRatioB->SetLineWidth(2);
0211    HRatioB->SetLineColor(2);
0212
0213    TH1D* frame = HRatio;
0214    frame->GetXaxis()->SetNdivisions(505);
0215    frame->SetTitle("");
0216    frame->SetStats(kFALSE);
0217 //   frame->GetXaxis()->SetTitle("Integrated Lumi (fb^{-1})");
0218      frame->GetXaxis()->SetTitle("Gain : Measured / Predicted");
0219 //   frame->GetYaxis()->SetTitle("Average Gain");
0220 //   frame->SetMaximum(1.25);
0221 //   frame->SetMinimum(0.80);
0222    frame->GetYaxis()->SetTitle("#APV");
0223 //   frame->SetMaximum(1.10);
0224 //   frame->SetMinimum(0.85);
0225    frame->GetYaxis()->SetTitleOffset(1.50);
0226    frame->Draw("");
0227
0228    HRatioB->Draw("same");
0229 //   c1->SetLogy(true);
0230
0231
0232    TLegend* leg = new TLegend(0.15,0.93,0.80, 0.75);
0233    leg->SetFillStyle(0);
0234    leg->SetBorderSize(0);
0235    leg->SetTextFont(43);
0236    leg->SetTextSize(20);
0237    leg->SetNColumns(1);
0238    leg->Clear();
0239    leg->AddEntry(HRatio, "All APVs", "L");
0240    leg->AddEntry(HRatioB, "APVs with bugged predicted G2", "L");
0241    leg->Draw();
0242
0243
0244    c1->SaveAs("Ratio.png");
0245
0246    HRatioC->SaveAs("histo.root");
0247 }
0248
0249
0250
0251 if(true){
0252    TCanvas* c1 = new TCanvas("c1","c1,",900,600);
0253
0254    TH2D* frame = HVS;
0255    frame->GetXaxis()->SetNdivisions(505);
0256    frame->SetTitle("");
0257    frame->SetStats(kFALSE);
0258 //   frame->GetXaxis()->SetTitle("Integrated Lumi (fb^{-1})");
0259      frame->GetXaxis()->SetTitle("Gain: Predicted(R251252*StrenghRatio from DelayScan)");
0260 //   frame->GetYaxis()->SetTitle("Average Gain");
0261 //   frame->SetMaximum(1.25);
0262 //   frame->SetMinimum(0.80);
0263    frame->GetYaxis()->SetTitle("Gain : Measured(R252116+R252226)");
0264 //   frame->SetMaximum(1.10);
0265 //   frame->SetMinimum(0.85);
0266    frame->GetYaxis()->SetTitleOffset(1.50);
0267    frame->Draw("COLZ");
0268    c1->SetLogz(true);
0269    c1->SaveAs("SCATTER.png");
0270 }
0271
0272
0273 }