C++ TFile::WriteObject方法代码示例

void regressCovU2(std::string iName) {
  TString  lName       = iName; TFile *lInput = TFile::Open(lName); 
  //TString  lName       = "../Data/s12-zjets_v29b.root"; TFile *lInput = TFile::Open(lName); 
  //TString  lName       = "../Data/r11-dimu_jetv3_v29b.root"; TFile *lInput = TFile::Open(lName); 
  //TString  lName       = "../Data/f11-zjets_jetv4_v11.root"; TFile *lInput = TFile::Open(lName); 
  TTree   *lRegress    = (TTree*)lInput    ->Get("Flat");
  
  GBRTrainer *train = new GBRTrainer;
  train->AddTree(lRegress);
  train->SetTrainingCut("(nbtag == 0)*(pt_z > 0)");//*PtWeight");
  train->SetMinEvents(250);
  train->SetShrinkage(0.3);
  train->SetMinCutSignificance(3.0);  
  
  train->SetTargetVar( "abs(pt_z*sin(phi_z-uphix_mva))*sqrt(TMath::Pi()/2.)/ux_mva");
  
  std::vector<std::string> *lVec = new std::vector<std::string>;
  lVec->push_back( "sumet"      );
  lVec->push_back( "npv"        );
  lVec->push_back( "pfu"        );
  lVec->push_back( "pfuphi"     );
  lVec->push_back( "tksumet"    );
  lVec->push_back( "tku"        );
  lVec->push_back( "tkuphi"     );
  lVec->push_back( "nopusumet"  );
  lVec->push_back( "nopuu"      );
  lVec->push_back( "nopuuphi"   );
  lVec->push_back( "pusumet"    );
  lVec->push_back( "pumet"      );
  lVec->push_back( "pumetphi"   );
  lVec->push_back( "pucsumet"   );
  lVec->push_back( "pucu"       );
  lVec->push_back( "pucuphi"    );
  lVec->push_back( "jetpt1"     );
  lVec->push_back( "jeteta1"    );
  lVec->push_back( "jetphi1"    );
  lVec->push_back( "jetpt2"     );
  lVec->push_back( "jeteta2"    );
  lVec->push_back( "jetphi2"    );
  lVec->push_back( "nalljet"    );
  lVec->push_back( "njet"       );
  lVec->push_back( "uphix_mva"  );
  lVec->push_back( "ux_mva"     );

  for (int i=0; i<int(lVec->size()); ++i) {
     train->AddInputVar(lVec->at(i));
  }
  
  ROOT::Cintex::Cintex::Enable();   
  const GBRForest *forest = train->TrainForest(300);
  
  TFile *fout = new TFile("gbru2cov.root","RECREATE");    
  fout->WriteObject(forest,"CovU2");
  fout->WriteObject(lVec, "varlist");
  fout->Close();
}

void test()
{

    TFile* fFakeFactorW  = TFile::Open("/afs/cern.ch/user/j/jsauvan/workspace/Projects/Htautau_Run2/Studies/FakeRate/ComputeFakeRates/plots/FakeFactors_Data_HighMT_2D/FakeFactors_Data_HighMT_2D.root");
    TFile* fFakeFactorQCD  = TFile::Open("/afs/cern.ch/user/j/jsauvan/workspace/Projects/Htautau_Run2/Studies/FakeRate/ComputeFakeRates/plots/FakeFactors_Data_QCDSS_2D/FakeFactors_Data_QCDSS_2D.root");
    TFile* fMtCorrection = TFile::Open("/afs/cern.ch/user/j/jsauvan/workspace/Projects/Htautau_Run2/Studies/FakeRate/ComputeMTCorrection/results/mtCorrections.root");
    TFile* fFractions    = TFile::Open("/afs/cern.ch/user/j/jsauvan/workspace/Projects/Htautau_Run2/Studies/FakeRate/ComputeBackgroundFractions/results/backgroundFraction_Iso_Medium_mvis_vs_mt.root");

    TH2F* fakeFactorW   = (TH2F*)fFakeFactorW->Get("FakeFactors_Data_HighMT_2D_Iso_Medium_InvertIso_Medium_tau_pt_vs_decayMode");
    TH2F* fakeFactorQCD = (TH2F*)fFakeFactorQCD->Get("FakeFactors_Data_QCDSS_2D_Iso_Medium_InvertIso_Medium_tau_pt_vs_decayMode");
    TGraph* mtCorrection = (TGraph*)fMtCorrection->Get("mt_correction");
    TH2F* fractionW = (TH2F*)fFractions->Get("h_backgroundFraction_Iso_Medium_mvis_vs_mt_W_Nom");
    TH2F* fractionQCD = (TH2F*)fFractions->Get("h_backgroundFraction_Iso_Medium_mvis_vs_mt_QCD_Nom");
    TH2F* fractionTT = (TH2F*)fFractions->Get("h_backgroundFraction_Iso_Medium_mvis_vs_mt_TT_Nom");
    TH2F* fractionVV = (TH2F*)fFractions->Get("h_backgroundFraction_Iso_Medium_mvis_vs_mt_VV_Nom");
    TH2F* fractionZJ = (TH2F*)fFractions->Get("h_backgroundFraction_Iso_Medium_mvis_vs_mt_ZJ_Nom");


    //wrappers
    WrapperTH2F* wFakeFactorW    = new WrapperTH2F(*fakeFactorW, "FF_W");
    WrapperTH2F* wFakeFactorQCD  = new WrapperTH2F(*fakeFactorQCD, "FF_QCD");
    WrapperTGraph* wMtCorrection = new WrapperTGraph(*mtCorrection, "MT_Corr");
    WrapperTH2F* wFractionW      = new WrapperTH2F(*fractionW, "f_W");
    WrapperTH2F* wFractionQCD    = new WrapperTH2F(*fractionQCD, "f_QCD");
    WrapperTH2F* wFractionTT     = new WrapperTH2F(*fractionTT, "f_TT");
    WrapperTH2F* wFractionVV     = new WrapperTH2F(*fractionVV, "f_VV");
    WrapperTH2F* wFractionZJ     = new WrapperTH2F(*fractionZJ, "f_ZJ");

    // formulas
    TFormula mtCorr("mtCorr", "x[0]*x[1]");
    WrapperTFormula* wFakeFactorWCorr = new WrapperTFormula(mtCorr, "FF_WCorr");

    TFormula combination("combination", "x[0]*x[2]+x[1]*(x[3]+x[4]+x[5]+x[6])");
    WrapperTFormula* wFakeFactorComb = new WrapperTFormula(combination, "FF_Comb");

    // fake factor
    // tau_pt = 0
    // tau_decay = 1
    // mt = 2
    // mvis = 3
    FakeFactor* factor = new FakeFactor();
    factor->addNode(wFakeFactorW, {}, {0,1});
    factor->addNode(wFakeFactorQCD, {}, {0,1});
    factor->addNode(wMtCorrection, {}, {2});
    factor->addNode(wFractionW, {}, {3,2});
    factor->addNode(wFractionQCD, {}, {3,2});
    factor->addNode(wFractionTT, {}, {3,2});
    factor->addNode(wFractionVV, {}, {3,2});
    factor->addNode(wFractionZJ, {}, {3,2});
    factor->addNode(wFakeFactorWCorr, {0,2}, {});
    factor->addNode(wFakeFactorComb, {1,8,4,3,5,6,7}, {});

    TFile* file = TFile::Open("test.root", "recreate");
    file->WriteObject(factor, "ff");
    file->Close();
    std::cout<<"Done\n";

}

void makeclassification() {
  
  Float_t *vars = new Float_t[10];
  
  //initialize TMVA Reader (example here is diphoton mva from higgs->gamma gamma mva analysis)
  TMVA::Reader* tmva = new TMVA::Reader();
  tmva->AddVariable("masserrsmeared/mass",            &vars[0]);
  tmva->AddVariable("masserrsmearedwrongvtx/mass",    &vars[1]);
  tmva->AddVariable("vtxprob",                        &vars[2]);
  tmva->AddVariable("ph1.pt/mass",                    &vars[3]);
  tmva->AddVariable("ph2.pt/mass",                    &vars[4]);
  tmva->AddVariable("ph1.eta",                        &vars[5]);
  tmva->AddVariable("ph2.eta",                        &vars[6]);
  tmva->AddVariable("TMath::Cos(ph1.phi-ph2.phi)"   , &vars[7]);
  tmva->AddVariable("ph1.idmva",                      &vars[8]);
  tmva->AddVariable("ph2.idmva",                      &vars[9]);
  
  tmva->BookMVA("BDTG","/afs/cern.ch/user/b/bendavid/cmspublic/diphotonmvaApr1/HggBambu_SMDipho_Jan16_BDTG.weights.xml");
  //tmva->BookMVA("BDTG","/scratch/bendavid/root/HggBambu_SMDipho_Jan16_BDTG.weights.xml");
  
  TMVA::MethodBDT *bdt = dynamic_cast<TMVA::MethodBDT*>(tmva->FindMVA("BDTG"));

 
  //enable root i/o for objects with reflex dictionaries in standalone root mode
  ROOT::Cintex::Cintex::Enable();   

  
  //open output root file
  TFile *fout = new TFile("gbrtest.root","RECREATE");
  
  //create GBRForest from tmva object
  GBRForest *gbr = new GBRForest(bdt);  
  
  //write to file
  fout->WriteObject(gbr,"gbrtest");

  fout->Close();
  
  
}

TObjArray* PrepSummaryKMC(const char* sumf, int cls, const char* pref, const char* outF)
{
  if (!gROOT->GetClass("KMCDetector")) gROOT->LoadMacro("KMCDetector.cxx+");
  TFile* fl = TFile::Open(sumf);
  if (!fl) {printf("No file %s\n",sumf); return 0;}
  TObjArray* arrs = (TObjArray*)gDirectory->Get("trSum");
  if (!arrs) {printf("No summary in file %s\n",sumf); return 0;}
  //
  TObjArray* sums =  ProcessSummary(arrs,cls,pref);
  if (!sums) return 0;
  //
  if (outF) {
    TFile* flOut = TFile::Open(outF,"update");
    if (!flOut) {printf("Failed to open output file %s\n",outF);}
    else {
      flOut->WriteObject(sums,sums->GetName(),"kSingleKey");
      flOut->Close();
      delete flOut;
      printf("Stored array %s in %s\n",sums->GetName(), outF);
    }
  }
  //
  return sums;
}

//.........这里部分代码省略.........
  for (int i=1; i<= runFilter.getLenght(); i++){
    valuesSilvia.push_back(atof(runFilter.getElementFromPropertyContainer(i,6).c_str()));
  }
  
  
  Run2011.LoadSelectionByDO(&runFilter, filterCol);
  
  //Run2011.PrintValidRuns();
  
  TFile* f = new TFile("Outputs/OutputFile.root", "RECREATE");
  TH2F* plot = new TH2F("RateVsLumi", "Rate Vs Lumi", 2000, 0, 2000, 100, 0, 10);
  
  plot->SetMarkerStyle(21);
  contFile->Close();
  
  int n=0; 
  int nSkip = 100;

  
  while (Run2011.hasNext()){
    
    n++;
    double NormalizedRate = 0;
    
    Run* currentRun = Run2011.next();
    //check on Runs
    //    if (n <= nSkip) continue;
    
    //otherwise no directories
    contFile = TFile::Open("Container20062011.root");
    currentRun->FillFromContainer(contFile, areaMap);
    contFile->Close();
    
    if (currentRun->getID() > 163869) break; //163869
    if (!Run2011.isValid(currentRun)) continue;
    
    //    cout << "continue on: " << currentRun->getID() << endl;
    double TotRate = 0;
    double TotArea = 0;
    int totAmount = 0;
    int k=0;    
    if (currentRun->isFilled() == true){
      
      const map<string,ExRoll* > barrelMap = currentRun->GetBarrel();
      
      for (map<string,ExRoll* >::const_iterator it = barrelMap.begin(); it != barrelMap.end(); it++){
	
	/*
	//Typical selection for Endcap
	//like RE+3_2_10
       	if (atoi(it->second->getRollID().c_str()) != 10) continue;
	if (it->second->getSector() != 2) continue;
	if (it->second->getWheel_Disk_Number() != 3) continue;
	*/

	
	//Typical selection for Barrel
	//like tot_W+1_RB4_04
	if (it->second->getWheel_Disk_Number() != 2) continue;
	if (it->second->getRollID() != "RB2in") continue;
	if (it->second->getSector() != 4) continue;
//	if ( it->second->getOptionalRollID().compare("F") ) continue;
//	cout << "SEC " << it->second->getSector()  << endl;
	cout << "OPT: " << it->second-> getOptionalRollID() << endl;
	cout << "PASS: " << it->second->getFullOnlineRollID() << endl;

	//Roll Rate
        for (int k=1; k<97; k++){
          ExStrip * strip = it->second->getStrip(k);
          TotRate+= strip->getRate();
	  if (TotRate == 0)
	    delete strip;
       	}
	
      	//Roll Area
        double rollArea = it->second->getRollArea();
        TotArea += rollArea;
	cout << rollArea << endl;
      }
      double NormalizedRate = TotRate / TotArea;
      
        if (TotRate != 0 && TotArea != 0 && currentRun->getDelivered() != 0 && currentRun->isFilled()){
		
	cout << currentRun->getID() << " "
	     << currentRun->getLumi() << " "
	     <<  "; rate: " << TotRate << "; area: "
	     << TotArea << "; norm: " 
	     << NormalizedRate << endl;

 	plot->Fill(currentRun->getLumi(), NormalizedRate);
	k++;
      }
    }
  }
  
    
  f->WriteObject(plot, "Plot");
  
  return 0;
}

//.........这里部分代码省略.........
    RollIDs["2"] = new TGraph();
    RollIDs["3"] = new TGraph();
    nSector = 36;
    
    Rolls = currentRun->GetEndCapMinus();
  }
  
  //Filling                                                                                                         
  contFile = TFile::Open( ContainerName.c_str() );
  
  cout << "Begin Loop" << endl;  
  for (map<string, TGraph*>::const_iterator itRoll = RollIDs.begin(); itRoll != RollIDs.end(); itRoll++){
    int nPoint = 0;
  bool setSize = false;    
    string name;
    
    if (subDet == "Barrel")
      name = "Barrel_" + itRoll->first;
    else if (subDet == "EndCapPlus" || subDet == "EndCapMinus")
      name = "EndCap_" + itRoll->first;
    
    itRoll->second->SetName(name.c_str());


 
  for (int iSector = minSector; iSector <= nSector; iSector++){
    
    double TotRate = 0;
    double TotArea = 0;
    
    
    for (map<string,ExRoll* >::const_iterator it = Rolls.begin(); it != Rolls.end(); it++){

      string th1path = it->second->getTH1Path();
      
      if (subDet == "Barrel"){
	string opt1 = "F";
	string opt2 = "+F";
	string opt3 = "-F";

	if (it->second->getSector() != iSector) continue;     
	
	if (it->second->getRollID() != itRoll->first) continue;      
	

	if (WheelDiskNo != -9)
	  if (it->second->getWheel_Disk_Number() != WheelDiskNo) continue;

	if (it->second->getOptionalRollID().find(opt1) == string::npos && it->second->getOptionalRollID().find(opt2) == string::npos) {
		 continue;
	}	
	if (it->second->getOptionalRollID().find(opt3) != string::npos){
		continue;
	}	
	
	
      }
      
      else if (subDet == "EndCapPlus" || subDet == "EndCapMinus"){
	
	if (atoi(it->second->getRollID().c_str()) != iSector) continue;
        
	if (it->second->getSector() != atoi(itRoll->first.c_str())) continue;
	
	if (WheelDiskNo != -9)
	  if (it->second->getWheel_Disk_Number() != WheelDiskNo) continue;
      }
      int nClones = it->second->getClones();
      emptyStrip += it->second->getEmptyStrip();

      TH1F* histo = (TH1F*)(contFile->Get(th1path.c_str()));

      int cNumber = partitionNumber;
      if (partitionNumber > nClones-1) continue;
      cout << "This roll has " << nClones << endl;
      double TotRate = 0;
      double TotArea = 0;
      int cEffectiveNumber =  cNumber + 1;
      
      NoiseAnalyzer nAnalyzer;
      TH1F* cHisto = it->second->getRateHistoForClone(histo, cEffectiveNumber);
      TotRate = nAnalyzer.HealClone(it->second, cHisto, cEffectiveNumber);
      double cloneArea = it->second->getCloneArea(cNumber);
      TotArea += cloneArea;
     
      double NormalizedRate = (TotRate / TotArea);
      itRoll->second->SetPoint(nPoint, (iSector-1)*6.28/(double)nSector, NormalizedRate/2);
      cout << "\t\tPOINT: " << nPoint << " " << (iSector-1)*6.28/(double)nSector << " " << NormalizedRate/2 << endl;
      nPoint++;
    }
  }  
  }


  for (map<string, TGraph*>::const_iterator it = RollIDs.begin(); it != RollIDs.end(); it++)
    f->WriteObject( it->second, it->second->GetName() );
  
  f->Close();
  return 0;
}

void eregmerge(bool doele) {
  
  TString dirname = "/afs/cern.ch/user/b/bendavid/CMSSWhgg/CMSSW_5_3_11_patch5/src/HiggsAnalysis/GBRLikelihoodEGTools/data/"; 
  gSystem->mkdir(dirname,true);
  gSystem->cd(dirname);    
  
  TString fnameeb;
  TString fnameee;
  if (doele) { 
    fnameeb = "wereg_ele_eb.root";
    fnameee = "wereg_ele_ee.root";
  }
  else if (!doele) {
    fnameeb = "wereg_ph_eb.root";
    fnameee = "wereg_ph_ee.root";
  }    
    
   
  TString infileeb = TString::Format("/afs/cern.ch/work/b/bendavid/bare/eregAug10RCalphafix/%s",fnameeb.Data());
  TString infileee = TString::Format("/afs/cern.ch/work/b/bendavid/bare/eregAug10RCalphafix/%s",fnameee.Data());
  
  TFile *fwseb = TFile::Open(infileeb); 
  TFile *fwsee = TFile::Open(infileee); 
  
  RooWorkspace *wseb = (RooWorkspace*)fwseb->Get("wereg");
  RooWorkspace *wsee = (RooWorkspace*)fwsee->Get("wereg");
  
  RooAbsPdf *sigpdfeborig = wseb->pdf("sigpdf");
  RooAbsPdf *sigpdfeeorig = wsee->pdf("sigpdf");
  
  RooAbsPdf *sigpdfeb = static_cast<RooAbsPdf*>(cloneRecursiveRename(sigpdfeborig,"EB"));
  RooAbsPdf *sigpdfee = static_cast<RooAbsPdf*>(cloneRecursiveRename(sigpdfeeorig,"EE"));
    
  RooWorkspace *wsout = new RooWorkspace("EGRegressionWorkspace");
  wsout->import(*sigpdfeb);
  wsout->import(*sigpdfee);
  
  TString outname;
  if (doele) outname = "regweights_v4_ele.root";
  else outname = "regweights_v4_ph.root";
  
  wsout->writeToFile(outname);
  
  RooArgList pdfeblist;
  RooArgSet *pdfebcomps = sigpdfeb->getComponents();
  RooArgSet *pdfebvars = sigpdfeb->getVariables();
  pdfeblist.add(*pdfebcomps);
  pdfeblist.add(*pdfebvars);
  delete pdfebcomps;
  delete pdfebvars;
  
  
  RooArgList pdfeelist;
  RooArgSet *pdfeecomps = sigpdfee->getComponents();
  RooArgSet *pdfeevars = sigpdfee->getVariables();
  pdfeelist.add(*pdfeecomps);
  pdfeelist.add(*pdfeevars);
  delete pdfeecomps;
  delete pdfeevars;  
  
  
//   RooArgList components(ws->components());
//   for (int iarg=0; iarg<components.getSize(); ++iarg) {
//     components.at(iarg)->SetName(TString::Format("%s_1",components.at(iarg)->GetName()));
//   }
  
  RooGBRFunction *funceb = static_cast<RooGBRFunction*>(pdfeblist.find("func_EB"));
  RooGBRFunction *funcee = static_cast<RooGBRFunction*>(pdfeelist.find("func_EE"));
  
//   funceb->Vars().Print("V");
//   funcee->Vars().Print("V");

  for (int ivar=0; ivar<funceb->Vars().getSize(); ++ivar) {
    printf("%i: %s, %s\n",ivar,funceb->Vars().at(ivar)->GetName(),funceb->Vars().at(ivar)->GetTitle());
  }
  
  for (int ivar=0; ivar<funcee->Vars().getSize(); ++ivar) {
    printf("%i: %s, %s\n",ivar,funcee->Vars().at(ivar)->GetName(),funcee->Vars().at(ivar)->GetTitle());
  }
  
  TString outnameforest;
  if (doele) outnameforest = "regweights_v4_forest_ele.root";
  else outnameforest = "regweights_v4_forest_ph.root";  
  
  TFile *fforest = new TFile(outnameforest,"RECREATE");
  fforest->WriteObject(funceb->Forest(),"EGRegressionForest_EB");
  fforest->WriteObject(funcee->Forest(),"EGRegressionForest_EE");
  fforest->Close();
  
}

//.........这里部分代码省略.........
 
                    
          int label = cl->GetLabel(0);
          TParticle* part = 0;
          if (label>=0 && (part=stack->Particle(label)) ) {
            cSum.pdg = part->GetPdgCode();
            cSum.eta = part->Eta();
            cSum.pt  = part->Pt();
            cSum.phi = part->Phi();
            cSum.prim = stack->IsPhysicalPrimary(label);
          } 
          cSum.ntr = 0;
          for (int ilb=0;ilb<3;ilb++) if (cl->GetLabel(ilb)>=0) cSum.ntr++;
          for (int i=0;i<3;i++) cSum.xyz[i] = xyzClGloF[i];
          //
          trOut->Fill();
          /*
          if (clsize==5) {
            printf("\nL%d(%c) Mod%d, Cl:%d | %+5.1f %+5.1f (%d/%d)|H:%e %e %e | C:%e %e %e\n",ilr,cl->TestBit(kSplit) ? 'S':'N',
             modID,icl,(txyzH[0]-xyzClTr[0])*1e4,(txyzH[2]-xyzClTr[2])*1e4, row,col,
             gxyzH[0],gxyzH[1],gxyzH[2],xyzClGlo[0],xyzClGlo[1],xyzClGlo[2]);
            cl->Print();
            pHit->Print();
            //
            double a0,b0,c0,a1,b1,c1,e0;
            pHit->GetPositionL0(a0,b0,c0,e0);
            pHit->GetPositionL(a1,b1,c1);
            float cloc[3];
            cl->GetLocalXYZ(cloc);
            printf("LocH: %e %e %e | %e %e %e\n",a0,b0,c0,a1,b1,c1);
            printf("LocC: %e %e %e | %e %e %e\n",cloc[0],cloc[1],cloc[2],xyzClTr[0],xyzClTr[1],xyzClTr[2]);
          }
          */
          //
        }
      }
    }
    
    //    layerClus.Clear();
    //
    arrMCTracks.Delete();
  }//event loop
  //
  DB.EndAndSort();
  DB.SetThresholdCumulative(0.95);
  cout << "Over threshold: : "<< DB.GetOverThr()<<endl;
  DB.Grouping(10,10);
  DB.PrintDB("Database1.txt"); 
  flOut->cd();
  trOut->Write();
  delete trOut;
  flOut->Close();
  flOut->Delete();
  DrawReport("clinfo.ps",&histoArr);
  TFile* flDB = TFile::Open("TopologyDatabase.root", "recreate");
  flDB->WriteObject(&DB,"DB","kSingleKey");
  flDB->Close();
  delete flDB;

  TCanvas* cnv123 = new TCanvas("cnv123","cnv123");
  cnv123->Divide(1,2);
  cnv123->Print("anglesdistr.pdf[");
  cnv123->cd(1);
  hL0A->Draw();
  cnv123->cd(2);
  hL0B->Draw();
  cnv123->Print("anglesdistr.pdf");
  cnv123->cd(1);
  hL1A->Draw();
  cnv123->cd(2);
  hL1B->Draw();
  cnv123->Print("anglesdistr.pdf");
  cnv123->cd(1);
  hL2A->Draw();
  cnv123->cd(2);
  hL2B->Draw();
  cnv123->Print("anglesdistr.pdf");
  cnv123->cd(1);
  hL3A->Draw();
  cnv123->cd(2);
  hL3B->Draw();
  cnv123->Print("anglesdistr.pdf");
  cnv123->cd(1);
  hL4A->Draw();
  cnv123->cd(2);
  hL4B->Draw();
  cnv123->Print("anglesdistr.pdf");
  cnv123->cd(1);
  hL5A->Draw();
  cnv123->cd(2);
  hL5B->Draw();
  cnv123->Print("anglesdistr.pdf");
  cnv123->cd(1);
  hL6A->Draw();
  cnv123->cd(2);
  hL6B->Draw();
  cnv123->Print("anglesdistr.pdf");
  cnv123->Print("anglesdistr.pdf]");
  //
}

void MakeSnapshot(Int_t run, const char* defStorage, TMap* specStorages, const char* snapshotFileName)
{
    AliCDBManager *cdb = AliCDBManager::Instance();
    cdb->SetDefaultStorage(defStorage);
    cdb->SetRun(run);

    TIter iter(specStorages->GetTable());
    TPair *pair = 0;
    while((pair = dynamic_cast<TPair*> (iter.Next()))){
	TObjString* caltype = dynamic_cast<TObjString*> (pair->Key());
	TObjString* specstor= dynamic_cast<TObjString*> (pair->Value());
	if (caltype && specstor)
	    //TString calType = caltype->GetString();
	    //TString specStor = specstor->GetString();
	    //cdb->SetSpecificStorage(calType.Data(),specStor.Data());
	    cdb->SetSpecificStorage(caltype->GetString().Data(),specstor->GetString().Data());
	else
	    //AliFatal("Error reading info for specific storage")
	    Printf("Error reading info for specific storage");
    }

    // ********************************** GRP ******************************************
    cdb->Get("GRP/CTP/Config");
    cdb->Get("GRP/Calib/LHCClockPhase");
    cdb->Get("GRP/GRP/Data");
    cdb->Get("GRP/Align/Data");
    cdb->Get("GRP/Calib/MeanVertexSPD");
    cdb->Get("GRP/Calib/MeanVertex");
    cdb->Get("GRP/Calib/MeanVertexTPC");
    cdb->Get("GRP/Calib/CosmicTriggers");
    cdb->Get("GRP/CTP/Scalers");
    cdb->Get("GRP/CTP/CTPtiming");
    cdb->Get("GRP/CTP/TimeAlign");
    cdb->Get("GRP/GRP/LHCData");
    cdb->Get("GRP/Calib/RecoParam");

    // ********************************** ALL ******************************************
    TString detStr = ("ITS TPC TRD TOF PHOS HMPID EMCAL MUON ZDC PMD T0 VZERO");
    //TString detStr = ("ITS MUON TPC");
    TObjArray *arr = detStr.Tokenize(' ');
    for (Int_t iDet=0; iDet<arr->GetEntries(); iDet++) {
	TObjString *detOStr = dynamic_cast<TObjString*>(arr->At(iDet));
	AliCDBManager::Instance()->GetAll(Form("%s/Calib/*",detOStr->GetString().Data()));
	AliCDBManager::Instance()->Get(Form("%s/Align/Data",detOStr->GetString().Data()));
    }

    // ******************************** TRIGGER ****************************************
    // Temporary fix - one has to define the correct policy in order
    // to load the trigger OCDB entries only for the detectors that
    // in the trigger or that are needed in order to put correct
    // information in ESD
    AliCDBManager::Instance()->GetAll("TRIGGER/*/*");

    // ********************************** HLT ******************************************
    // cdb->Get("HLT/ConfigHLT/esdLayout");
    // cdb->Get("HLT/Calib/StreamerInfo");

    TMap* entriesMap = const_cast<TMap*>(cdb->GetEntryCache());
    Printf("\nentriesMap has %d entries!\n", entriesMap->GetEntries());

    TList* entriesList = const_cast<TList*>(cdb->GetRetrievedIds());
    Printf("\nentriesList has %d entries!\n", entriesList->GetEntries());

    //TString filename(TString::Format("CDBsnapshot_Run%d.root",run));
    TString filename(snapshotFileName);
    TFile *f = new TFile(filename.Data(),"recreate");
    f->cd();
    f->WriteObject(entriesMap,"entriesMap");
    f->WriteObject(entriesList,"entriesList");
    f->Close();
    entriesMap->SetOwnerKeyValue(kFALSE,kFALSE);
    entriesList->SetOwner(kFALSE);
}

void plotRatio() {


  // set proper style for plots
  gROOT->ProcessLine(".L mystyle.C");
  setTDRStyle();
  tdrStyle->SetErrorX(0.5);
  tdrStyle->SetPadLeftMargin(0.18);
  tdrStyle->SetPadRightMargin(0.10);
  tdrStyle->SetPadBottomMargin(0.16);
  tdrStyle->SetLegendBorderSize(0);
  tdrStyle->SetTitleYOffset(1.5);


  TFile* f1 = new TFile("histograms_7TeV.root","read");
  TFile* f2 = new TFile("histograms_10TeV.root","read");
  gROOT->cd();
  TH1D ratioHist[nHists];
  for (int i=0; i<nHists; i++) {
    TH1D* numHist = (TH1D* ) f1->Get( histNames[i] );
    TH1D* denHist = (TH1D* ) f2->Get( histNames[i] );
    ratioHist[i] = (TH1D*) numHist->Clone(histNames[i]);
    ratioHist[i]->Divide( denHist );
    ratioHist[i]->GetYaxis()->SetTitle("Events @ 7 TeV / @ 10 TeV");
  }

  TFile* file = new TFile( "histograms_ratio.root", "RECREATE" );
  for (int i=0; i<nHists; i++) {
    file->WriteObject(ratioHist[i], histNames[i]);
  }

  TH1D* ZPt_Gen = (TH1D* ) file->Get("ZPt_Gen");
  TH1D* ZPt_Reco = (TH1D* ) file->Get("ZPt_Reco");
  TH1D* GenJetPt = (TH1D* ) file->Get("GenJetPt");
  TH1D* CaloJetPt = (TH1D* ) file->Get("CaloJetPt");
  TH1D* PFJetPt = (TH1D* ) file->Get("PFJetPt");

  TH1D* GenJetEta = (TH1D* ) file->Get("GenJetEta");
  TH1D* CaloJetEta = (TH1D* ) file->Get("CaloJetEta");
  TH1D* PFJetEta = (TH1D* ) file->Get("PFJetEta");
  
  TH1D* GenJetPt2 = (TH1D* ) file->Get("GenJetPt2");
  TH1D* CaloJetPt2 = (TH1D* ) file->Get("CaloJetPt2");
  TH1D* PFJetPt2 = (TH1D* ) file->Get("PFJetPt2");

  TH1D* GenJetPt3 = (TH1D* ) file->Get("GenJetPt3");
  TH1D* CaloJetPt3 = (TH1D* ) file->Get("CaloJetPt3");
  TH1D* PFJetPt3 = (TH1D* ) file->Get("PFJetPt3");

  TH1D* CaloJetPtRatio2over1 = (TH1D* ) file->Get("CaloJetPtRatio2over1");
  TH1D* CaloJetPtRatio3over2 = (TH1D* ) file->Get("CaloJetPtRatio3over2");
  TH1D* PFJetPtRatio2over1 = (TH1D* ) file->Get("PFJetPtRatio2over1");
  TH1D* PFJetPtRatio3over2 = (TH1D* ) file->Get("PFJetPtRatio3over2");
  TH1D* GenJetPtRatio2over1 = (TH1D* ) file->Get("GenJetPtRatio2over1");
  TH1D* GenJetPtRatio3over2 = (TH1D* ) file->Get("GenJetPtRatio3over2");


  TH1D* nJetsCalo = (TH1D* ) file->Get("nJetsCalo");
  TH1D* nJetsPF = (TH1D* ) file->Get("nJetsPF");
  TH1D* nJetsGen = (TH1D* ) file->Get("nJetsGen");


  makeplotTwo(*ZPt_Reco, *ZPt_Gen, "ZPt_spectrum", 1);
  makeplotThree(*CaloJetPt, *GenJetPt, *PFJetPt, "Jet_spectrum", 1);
  makeplotThree(*CaloJetPt2, *GenJetPt2, *PFJetPt2, "Jet_spectrum2", 1);
  makeplotThree(*CaloJetPt3, *GenJetPt3, *PFJetPt3, "Jet_spectrum3", 1);
  makeplotThree(*CaloJetEta, *GenJetEta, *PFJetEta, "Jet_eta_spectrum", 2);

  makeplotThree(*CaloJetPtRatio2over1, *GenJetPtRatio2over1, *PFJetPtRatio2over1, "JetPtRatio2over1", 0);
  makeplotThree(*CaloJetPtRatio3over2, *GenJetPtRatio3over2, *PFJetPtRatio3over2, "JetPtRatio3over2", 0);
  makeplotThree( *nJetsCalo, *nJetsGen, *nJetsPF, "JetMultiplicity", 1);

}

//.........这里部分代码省略.........

    currentRun->FillFromContainer(contFile, areaMap);

        
    double TotRate = 0;
    double TotArea = 0;
    int totAmount = 0;
    if (currentRun->isFilled() == false){
	cout << "Not filled: " << currentRun->getID() << endl;
    }
    if (currentRun->isFilled() == true){

      map<string,ExRoll* > barrelMap;

      if (zone == "Barrel")
	barrelMap = currentRun->GetBarrel();

      else if (zone == "EndCapPlus")
	barrelMap = currentRun->GetEndCapPlus();

      else if (zone == "EndCapMinus")
	barrelMap = currentRun->GetEndCapMinus();

      for (map<string,ExRoll* >::const_iterator it = barrelMap.begin(); it != barrelMap.end(); it++){

	if (RollID != "-9")
	  if (it->second->getRollID().c_str() != RollID) continue;

	if (Sector != -9)
	  if (it->second->getSector() != Sector) continue;

	if (WheelDiskNo != -9)
	  if (it->second->getWheel_Disk_Number() != WheelDiskNo) continue;

	//When we will decide to plot partition by partition
	if (optionalID != "-9"){
//	cout << optionalID << endl;
//	cout << it->second->getOptionalRollID() << endl;
                if ( it->second->getOptionalRollID().compare(optionalID) ) continue;
	}

      
        emptyStrip += it->second->getEmptyStrip();
//	if ( it->second->getOptionalRollID().compare("B") ) continue;
//	cout << "SEC " << it->second->getSector()  << endl;
/*	cout << "OPT: " << it->second-> getOptionalRollID() << endl;
	cout << "RAW: " << it->second->getRawIDofClone(1) << endl;*/
//	cout << "PASS: " << it->second->getFullOnlineRollID() << endl;
/*	TotRate += it->second->CalculateRateFromTH1(corrected);
	string th1path = it->second->getTH1Path();
	TH1F* histo = (TH1F*)(contFile->Get(th1path.c_str()));*/


        NoiseAnalyzer nAnalyzer;
        TotRate += nAnalyzer.HealCloneByClone(it->second);
//	cout << "VALUE: " << nAnalyzer.HealCloneByClone(it->second) << endl;
        double rollArea = it->second->getRollArea();
        TotArea += rollArea;	
      }

      //Free RAM
      for (map<string,ExRoll* >::const_iterator it = barrelMap.begin(); it != barrelMap.end(); it++){
	delete it->second;
      }

      double NormalizedRate = TotRate / TotArea;
//       cout << NormalizedRate << endl;	      
  string zone = argv[1];                                                                                                                                                  
  string RollID = argv[4]; 
  int Sector = atoi(argv[3]);                                                                                                                                          
  int WheelDiskNo = atoi(argv[2]);     


        if (TotRate != 0 && TotArea != 0 && currentRun->getDelivered() != 0 && currentRun->isFilled()){
//	currentRun->PrintStatistics();	
             
             cout << currentRun->getID() <<
	     "\t" << currentRun->getLumiSection() << "\t" << currentRun->getDelivered() 
	     << "\t" << currentRun->getLumi() <<
	     "\t" << TotRate << "\t" << TotArea << "\t" << NormalizedRate/2 << endl;
		
// 	plot->SetBinContent(plot->FindBin(currentRun->getLumi()), NormalizedRate);
	cout << "POINT: " << k << " "  << currentRun->getLumi()<< " " << NormalizedRate / 2 << endl;; 
	plot->SetPoint(k , currentRun->getLumi(), NormalizedRate / 2);
	plotRun->SetPoint(k, currentRun->getID(), NormalizedRate / 2);
	k++;
	}
    }
    delete currentRun;
    contFile->Close();  
    contFile->Delete();
  }

  plot->SetDrawOption("P");
  f->WriteObject(plotRun, plotRun->GetName());
  f->WriteObject(plot, plot->GetName());

  
  return 0;
}

//.........这里部分代码省略.........
		if ( bHf->Integral() <=0 )  fitPullf = false;
		if (fitPullf) {bHf->Fit("gaus"); bHf->GetFunction("gaus")->SetLineColor(2);}

		c->Clear();
		
 		TPad pad1("t1","",0.01,0.01,0.66,0.95);
 		TPad pad2("t2","",0.70,0.20,0.98,0.80);

		pad1.SetNumber(1); pad2.SetNumber(2);

                if ( isFitted ) {pad2.Draw();}

		pad1.Draw();
		pad1.SetGrid(true);


		TLatex *titletext = new TLatex();titletext->SetNDC();

		
		c->cd(1); bH->Draw(); bHf->Draw("same");
		TLegend *legend = new TLegend(0.6,0.8,0.9,0.89);
		legend->SetFillColor(0);
		legend->AddEntry(bH,"All Toys","L");
		legend->AddEntry(bHf,selectString.c_str(),"L");
		legend->Draw();

		if (fitPull){
			c->cd(2);
			double gap;
			TLatex *tlatex = new TLatex(); tlatex->SetNDC(); 
			if (fitPullf) {tlatex->SetTextSize(0.09); gap=0.12;}
			else  {tlatex->SetTextSize(0.11);gap=0.14;}

			tlatex->SetTextColor(4);
			tlatex->DrawLatex(0.11,0.80,Form("Mean    : %.3f #pm %.3f",bH->GetFunction("gaus")->GetParameter(1),bH->GetFunction("gaus")->GetParError(1)));
			tlatex->DrawLatex(0.11,0.80-gap,Form("Sigma   : %.3f #pm %.3f",bH->GetFunction("gaus")->GetParameter(2),bH->GetFunction("gaus")->GetParError(2)));

			if (fitPullf){ 
				tlatex->SetTextColor(2);
				tlatex->DrawLatex(0.11,0.60,Form("Mean    : %.3f #pm %.3f",bHf->GetFunction("gaus")->GetParameter(1),bHf->GetFunction("gaus")->GetParError(1)));
				tlatex->DrawLatex(0.11,0.60-gap,Form("Sigma   : %.3f #pm %.3f",bHf->GetFunction("gaus")->GetParameter(2),bHf->GetFunction("gaus")->GetParError(2)));
			}

			tlatex->SetTextSize(0.10);
			tlatex->SetTextColor(1);
				
                        tlatex->DrawLatex(0.11,0.33,Form("Pre-fit: %.3f",p_mean));
			
			pullSummaryMap[name]=std::make_pair<double,double>(bH->GetFunction("gaus")->GetParameter(1),bH->GetFunction("gaus")->GetParameter(2));
			nPulls++;

		}

		// double titleSize = isFitted ? 0.1 : 0.028;
		//titletext->SetTextSize(titleSize);titletext->SetTextAlign(21); titletext->DrawLatex(0.55,0.92,name);
		c->SaveAs(Form("%s_normresiduals_%s.pdf",treename.c_str(),(do7TeV ? "7TeV": "8TeV")));
                c->SaveAs(Form("mlfit/%s_residual_%s_%s.pdf",name,treename.c_str(),(do7TeV ? "7TeV": "8TeV")));
		fOut->WriteObject(c,Form("%s_%s",treename.c_str(),name));
	}
	
	if (nPulls>0){
	  
	    std::cout << "Generating Pull Summaries" <<std::endl; 
	    int nRemainingPulls = nPulls;
	    TCanvas *hc = new TCanvas("hc","",3000,2000); hc->SetGrid(0);
	    std::map < const char*, std::pair <double,double> >::iterator pull_it = pullSummaryMap.begin();
	    std::map < const char*, std::pair <double,double> >::iterator pull_end = pullSummaryMap.end();

	    int pullPlots = 1;
	    while (nRemainingPulls > 0){

		int nThisPulls = min(maxPullsPerPlot,nRemainingPulls);

		TH1F pullSummaryHist("pullSummary","",nThisPulls,0,nThisPulls);
		for (int pi=1;pull_it!=pull_end && pi<=nThisPulls ;pull_it++,pi++){
			pullSummaryHist.GetXaxis()->SetBinLabel(pi,(*pull_it).first);
			pullSummaryHist.SetBinContent(pi,((*pull_it).second).first);
			pullSummaryHist.SetBinError(pi,((*pull_it).second).second);
			nRemainingPulls--;
		}		

		pullSummaryHist.SetMarkerStyle(21);pullSummaryHist.SetMarkerSize(1.5);pullSummaryHist.SetMarkerColor(2);pullSummaryHist.SetLabelSize(pullLabelSize);
		pullSummaryHist.GetYaxis()->SetRangeUser(-1,1);pullSummaryHist.GetYaxis()->SetTitle("residual summary (relative)");pullSummaryHist.Draw("E1");
		hc->SaveAs(Form("%s_normresiduals_%s.pdf",treename.c_str(),(do7TeV ? "7TeV": "8TeV")));
                hc->SaveAs(Form("mlfit/residual_summary_%d_%s_%s.pdf",pullPlots,treename.c_str(),(do7TeV ? "7TeV": "8TeV")));
		fOut->WriteObject(hc,Form("comb_pulls_%s_%d",treename.c_str(),pullPlots));
	//	hc->SaveAs(Form("comb_pulls_%s_%d.pdf",treename.c_str(),pullPlots));
		pullPlots++;
	   }

	    delete hc;
	}

	c->SaveAs(Form("%s_normresiduals_%s.pdf]",treename.c_str(),(do7TeV ? "7TeV": "8TeV")));
	fOut->Close();
	delete c;
	return;


}

//.........这里部分代码省略.........
			//h2d_corr->SetMarkerColor(4);
			//h2d_corr->SetTitle("");
			//h2d_corr->GetXaxis()->SetTitle(Form("%s_In",name));
			//h2d_corr->GetYaxis()->SetTitle(Form("%s",name));
			titletext->SetTextAlign(11);
			titletext->SetTextSize(0.05);
			titletext->DrawLatex(0.05,0.02,Form("%s_In",name));
			titletext->SetTextAngle(90);
			titletext->DrawLatex(0.04,0.06,Form("%s",name));
			titletext->SetTextAngle(0);
		}

		
		c->cd(1); bH->Draw(); bHf->Draw("same");
		TLegend *legend = new TLegend(0.6,0.8,0.9,0.89);
		legend->SetFillColor(0);
		legend->AddEntry(bH,"All Toys","L");
		legend->AddEntry(bHf,selectString.c_str(),"L");
		legend->Draw();

		if (doPull && plotLH) {
			c->cd(2); gr->Draw("ALP");
		}

		if (fitPull){
			c->cd(3);
			double gap;
			TLatex *tlatex = new TLatex(); tlatex->SetNDC(); 
			if (fitPullf) {tlatex->SetTextSize(0.09); gap=0.12;}
			else  {tlatex->SetTextSize(0.11);gap=0.14;}

			tlatex->SetTextColor(4);
			tlatex->DrawLatex(0.11,0.80,Form("Mean    : %.3f #pm %.3f",bH->GetFunction("gaus")->GetParameter(1),bH->GetFunction("gaus")->GetParError(1)));
			tlatex->DrawLatex(0.11,0.80-gap,Form("Sigma   : %.3f #pm %.3f",bH->GetFunction("gaus")->GetParameter(2),bH->GetFunction("gaus")->GetParError(2)));

			if (fitPullf){ 
				tlatex->SetTextColor(2);
				tlatex->DrawLatex(0.11,0.60,Form("Mean    : %.3f #pm %.3f",bHf->GetFunction("gaus")->GetParameter(1),bHf->GetFunction("gaus")->GetParError(1)));
				tlatex->DrawLatex(0.11,0.60-gap,Form("Sigma   : %.3f #pm %.3f",bHf->GetFunction("gaus")->GetParameter(2),bHf->GetFunction("gaus")->GetParError(2)));
			}

			tlatex->SetTextSize(0.10);
			tlatex->SetTextColor(1);
				
                        if(namestr.find("n_exp")!=string::npos) tlatex->DrawLatex(0.11,0.33,Form("Pre-fit: %.3f",prevals_[name].first));
			else tlatex->DrawLatex(0.11,0.33,Form("Pre-fit #pm #sigma_{#theta}: %.3f #pm %.3f",prevals_[name].first, p_err));
			tlatex->DrawLatex(0.11,0.18,Form("Best-fit (#theta_{B})  : %.3f ",p_mean));
			tlatex->DrawLatex(0.11,0.03,Form("Best-fit (#theta_{S+B}): %.3f ",bfvals_sb_[name].first));
			
			pullSummaryMap[name]=std::make_pair<double,double>(bH->GetFunction("gaus")->GetParameter(1),bH->GetFunction("gaus")->GetParameter(2));
			nPulls++;

		}

		double titleSize = isFitted ? 0.1 : 0.028;
		titletext->SetTextSize(titleSize);titletext->SetTextAlign(21); titletext->DrawLatex(0.55,0.92,name);
		c->SaveAs(Form("%s.pdf",treename.c_str()));
		fOut->WriteObject(c,Form("%s_%s",treename.c_str(),name));
		//c->SaveAs(Form("%s_%s.pdf",treename.c_str(),name));
	}
	
	if (doPull && nPulls>0){
	  
	    std::cout << "Generating Pull Summaries" <<std::endl; 
	    int nRemainingPulls = nPulls;
	    TCanvas *hc = new TCanvas("hc","",3000,2000); hc->SetGrid(0);
	    std::map < const char*, std::pair <double,double> >::iterator pull_it = pullSummaryMap.begin();
	    std::map < const char*, std::pair <double,double> >::iterator pull_end = pullSummaryMap.end();

	    int pullPlots = 1;
	    while (nRemainingPulls > 0){

		int nThisPulls = min(maxPullsPerPlot,nRemainingPulls);

		TH1F pullSummaryHist("pullSummary","",nThisPulls,0,nThisPulls);
		for (int pi=1;pull_it!=pull_end && pi<=nThisPulls ;pull_it++,pi++){
			pullSummaryHist.GetXaxis()->SetBinLabel(pi,(*pull_it).first);
			pullSummaryHist.SetBinContent(pi,((*pull_it).second).first);
			pullSummaryHist.SetBinError(pi,((*pull_it).second).second);
			nRemainingPulls--;
		}		

		pullSummaryHist.SetMarkerStyle(21);pullSummaryHist.SetMarkerSize(1.5);pullSummaryHist.SetMarkerColor(2);pullSummaryHist.SetLabelSize(pullLabelSize);
		pullSummaryHist.GetYaxis()->SetRangeUser(-3,3);pullSummaryHist.GetYaxis()->SetTitle("pull summary (n#sigma)");pullSummaryHist.Draw("E1");
		hc->SaveAs(Form("%s.pdf",treename.c_str()));
		fOut->WriteObject(hc,Form("comb_pulls_%s_%d",treename.c_str(),pullPlots));
	//	hc->SaveAs(Form("comb_pulls_%s_%d.pdf",treename.c_str(),pullPlots));
		pullPlots++;
	   }

	    delete hc;
	}

	c->SaveAs(Form("%s.pdf]",treename.c_str()));
	fOut->Close();
	delete c;
	return;


}

int main( int argc ,char *argv[] ){
  
  DataObject areas;
  areas.getDataFromFile("resources/static_resources/area_noise_CMSSW_IDS.txt",2);

  map<string, double> areaMap;
  
  for (int i=1; i<= areas.getLenght(); i++){
    double area = atof(areas.getElementFromPropertyContainer(i,2).c_str());
    string roll = areas.getElementFromPropertyContainer(i,1);
    areaMap[roll] = area;
  }

  //Create a new Acquisition Registry
  TFile* contFile = TFile::Open("Container.root");
  Acquisition Run2011(contFile);
  TFile* f = new TFile("OutputFile.root", "RECREATE");
  TH2F* plot = new TH2F("RateVsLumi", "Rate Vs Lumi", 2000, 0, 2000, 100, 0, 100);
  contFile->Close();
  while (Run2011.hasNext()){

    double NormalizedRate = 0;

    Run* currentRun = Run2011.next();
    //otherwise no directories
    contFile = TFile::Open("Container.root");
    currentRun->FillFromContainer(contFile, areaMap);
    contFile->Close();

    double TotRate = 0;
    double TotArea = 0;
    int totAmount = 0;

    if (currentRun->isFilled() == true){

      const map<string,ExRoll* > barrelMap = currentRun->GetBarrel();

      for (map<string,ExRoll* >::const_iterator it = barrelMap.begin(); it != barrelMap.end(); it++){

       //Roll Rate
        for (int k=1; k<=96; k++){
          ExStrip * strip = it->second->getStrip(k);
          TotRate = TotRate + strip->getRate();
       	}

      	//Roll Area
        double rollArea = it->second->getRollArea();
        TotArea = TotArea + rollArea;
      }
      double NormalizedRate = TotRate / TotArea;

      if (TotRate != 0 && currentRun->getDelivered() != 0 ){
        cout << currentRun->getID() << " " << currentRun->getLumi() << " " << TotRate << " " << TotArea << " " << NormalizedRate <<endl;
 	plot->Fill(currentRun->getLumi(), NormalizedRate);
	cout  << currentRun->getLumi() << " " << NormalizedRate << endl; 
  }
    }
  }
  f->WriteObject(plot, "Plot");
  return 0;
}

void TrigWeightFit1D(){
  TFile *f = new TFile("../data/scale_factors/HLTEffParkedABCD_1D.root");
  TFile *out = new TFile("HLTEffParkedABCD_1DFit.root","RECREATE");
  
  unsigned nVars=4;
  std::vector<std::string> histstofit[4];
  std::vector<std::string> mjjhists;
  mjjhists.push_back("hData_MJJ_1D_A");
  mjjhists.push_back("hData_MJJ_1D_BC");
  mjjhists.push_back("hData_MJJ_1D_D");
  std::vector<std::string> metl1hists;
  metl1hists.push_back("hData_MET_L1_A");
  metl1hists.push_back("hData_MET_L1_BC");
  metl1hists.push_back("hData_MET_L1_D");
  std::vector<std::string> methists;
  methists.push_back("hData_MET_1D_A");
  methists.push_back("hData_MET_1D_BC");
  methists.push_back("hData_MET_1D_D");
  std::vector<std::string> jet2hists;
  jet2hists.push_back("hData_JET2_1D_A");
  jet2hists.push_back("hData_JET2_1D_BC");
  jet2hists.push_back("hData_JET2_1D_D");

  histstofit[0]=mjjhists;
  histstofit[1]=metl1hists;
  histstofit[2]=methists;
  histstofit[3]=jet2hists;

  //create a function with 3 parameters in the range [-3,3]
  for(unsigned iVar=0;iVar<nVars;iVar++){
    for(unsigned iTurnon=0;iTurnon<histstofit[iVar].size();iTurnon++){
      TH1F *hpx = (TH1F*)f->Get(histstofit[iVar][iTurnon].c_str()); //LOOP OVER THE TRIG TURN ONS
      double rangelow,rangehigh;
      if(iVar==0){
	rangelow=600;
	rangehigh=2200;
      }
      else if(iVar==1){
	rangelow=60;
	rangehigh=400;
      }
      else if(iVar==2){
	rangelow=0;
	rangehigh=400;
      }
      else if(iVar==3){
	rangelow=0;
	rangehigh=100;
      }
      TF1 *func = new TF1("erf","0.5*(1+TMath::Erf((x-[0])/(sqrt([1]))))",rangelow,rangehigh);
      func->SetParameters(100,10000);
      func->SetParNames("Turn On","Sigma");
      func->SetParLimits(0,-10000,10000);
      func->SetParLimits(1,0,100000);
      hpx->Fit("erf","R");
      TCanvas *c1=new TCanvas();
      hpx->Draw();
      c1->SaveAs((histstofit[iVar][iTurnon]+".pdf").c_str());
      TF1 *fitresult=hpx->GetFunction("erf");
      //Get name to save as
      std::string outname=histstofit[iVar][iTurnon].substr(1);
      out->WriteObject(fitresult,("f"+outname).c_str());
    }
  }
}