C++ TString::Replace方法代码示例

///
/// Scale the error of a given observable.
/// Both stat and syst errors are being scaled by the same factor.
/// In order to become effective, the PDF needs to be rebuild by
/// calling buildCov() and buildPdf().
///
/// \param obsname	- observable name. It may or may not include a unique ID string, both works.
/// \param scale	- the scale factor the current error is being multiplied with
/// \return		- true if successful
///
bool PDF_Abs::ScaleError(TString obsname, float scale)
{
	// remove unique ID if necessary
	TString UID = "UID";
	if ( obsname.Contains(UID) ){
		obsname.Replace(obsname.Index(UID), obsname.Length(), ""); // delete the unique ID. That should leave just the observable name.
	}
	// find the index of the observable - if it exists at all!
	if ( !hasObservable(obsname) ){
		cout << "PDF_Abs::ScaleError() : ERROR : observable '" << obsname << "' not found." << endl;
		return false;
	}
	int index = -1;
	for ( int i=0; i<getNobs(); i++ ){
		if ( observables->at(i)->GetName()==obsname ){
			index = i;
			break;
		}
	}
	if ( index==-1 ){
		// this should never happen...
		cout << "PDF_Abs::ScaleError() : ERROR : internal self inconsistency discovered. Exit." << endl;
		assert(0);
	}
	// scale error
	StatErr[index] *= scale;
	SystErr[index] *= scale;
	// update error source string
	obsErrSource += " (PDF_Abs::ScaleError(): scaled error of " + obsname + ")";
	return true;
}

TString nPart(Int_t part,TString string)
{
    if (part <= 0) return "";
    for (int i = 1; i < part; i++)    //part-1 times
    {
        if (string.Index(";") < 0) return "";
        string.Replace(0,string.Index(";")+1,"",0);
    }
    if (string.Index(";") >= 0)
        string.Remove(string.Index(";"));
    return string;
}

///
/// Compute the file name of the parameter file defining the Asimov
/// point where the Asimov toy is generated at.
/// The combiner name will be followed by the Asimov addition (getAsimovCombinerNameAddition()),
/// but without the number denoting the Asimov point in the file, as the file contains all points.
/// Format of returned filename:
///
/// plots/par/basename_combinernameAsimov[_+N][_-N]_var1[_var2]_genpoints.dat
///
/// \param c - Combiner object
/// \return - filename
///
TString FileNameBuilder::getFileNameAsimovPar(const Combiner *c)
{
	TString name = "plots/par/";
	name += getFileBaseName(c);
	// remove any string after the "Asimov" token and the first "_" after that
	// e.g.: "combinerAsimov3_" -> "combinerAsimov_"
	int startOfToken = name.Index(m_asimov);
	int startOfFirstUnderscore = name.Index("_",startOfToken);
	int length = startOfFirstUnderscore-(startOfToken+m_asimov.Sizeof())+1;
	name.Replace(startOfToken+m_asimov.Sizeof()-1, length, "");
	name += "_genpoints.dat";
	return name;
}

TVirtualPad* ParameterEvolutionPlotter::selectNewPad()
{
	TCanvas* c1 = m_canvases[m_canvases.size()-1];
	if ( m_padId>=6 ){
		// Create a new canvas that has the old title "foo 3" but with
		// incremented number: "foo 4". Only one-digit numbers are supported.
		TString title = c1->GetTitle();
		int oldNumber = TString(title[title.Sizeof()-2]).Atoi(); // get last character and turn into integer
		title.Replace(title.Sizeof()-2,1,Form("%i",++oldNumber)); // replace last character with incremented integer
		c1 = selectNewCanvas(title);
	}
	m_padId+=1;
	return c1->cd(m_padId);
}

void
RooStats::HypoTestInvTool::AnalyzeResult( HypoTestInverterResult * r,
                                          int calculatorType,
                                          int testStatType, 
                                          bool useCLs,  
                                          int npoints,
                                          const char * fileNameBase ){

   // analyze result produced by the inverter, optionally save it in a file 
   
  
   double lowerLimit = 0;
   double llError = 0;
#if defined ROOT_SVN_VERSION &&  ROOT_SVN_VERSION >= 44126
   if (r->IsTwoSided()) {
      lowerLimit = r->LowerLimit();
      llError = r->LowerLimitEstimatedError();
   }
#else
   lowerLimit = r->LowerLimit();
   llError = r->LowerLimitEstimatedError();
#endif

   double upperLimit = r->UpperLimit();
   double ulError = r->UpperLimitEstimatedError();

   //std::cout << "DEBUG : [ " << lowerLimit << " , " << upperLimit << "  ] " << std::endl;
      
   if (lowerLimit < upperLimit*(1.- 1.E-4) && lowerLimit != 0) 
      std::cout << "The computed lower limit is: " << lowerLimit << " +/- " << llError << std::endl;
   std::cout << "The computed upper limit is: " << upperLimit << " +/- " << ulError << std::endl;
  

   // compute expected limit
   std::cout << "Expected upper limits, using the B (alternate) model : " << std::endl;
   std::cout << " expected limit (median) " << r->GetExpectedUpperLimit(0) << std::endl;
   std::cout << " expected limit (-1 sig) " << r->GetExpectedUpperLimit(-1) << std::endl;
   std::cout << " expected limit (+1 sig) " << r->GetExpectedUpperLimit(1) << std::endl;
   std::cout << " expected limit (-2 sig) " << r->GetExpectedUpperLimit(-2) << std::endl;
   std::cout << " expected limit (+2 sig) " << r->GetExpectedUpperLimit(2) << std::endl;
  

   printf("\n\n  owen: (-2s,-1s,m,1s,2s)  %5.2f   %5.2f   %5.2f   %5.2f   %5.2f\n\n",
      r->GetExpectedUpperLimit(-2),
      r->GetExpectedUpperLimit(-1),
      r->GetExpectedUpperLimit(0),
      r->GetExpectedUpperLimit(1),
      r->GetExpectedUpperLimit(2) ) ;
  
   // write result in a file 
   if (r != NULL && mWriteResult) {
    
      // write to a file the results
      const char *  calcType = (calculatorType == 0) ? "Freq" : (calculatorType == 1) ? "Hybr" : "Asym";
      const char *  limitType = (useCLs) ? "CLs" : "Cls+b";
      const char * scanType = (npoints < 0) ? "auto" : "grid";
      if (mResultFileName.IsNull()) {
         mResultFileName = TString::Format("%s_%s_%s_ts%d_",calcType,limitType,scanType,testStatType);      
         //strip the / from the filename
         if (mMassValue.size()>0) {
            mResultFileName += mMassValue.c_str();
            mResultFileName += "_";
         }
    
         TString name = fileNameBase; 
         name.Replace(0, name.Last('/')+1, "");
         mResultFileName += name;
      }

      TFile * fileOut = new TFile(mResultFileName,"RECREATE");
      r->Write();
      fileOut->Close();                                                                     
   }   
  
  
   // plot the result ( p values vs scan points) 
   std::string typeName = "";
   if (calculatorType == 0 )
      typeName = "Frequentist";
   if (calculatorType == 1 )
      typeName = "Hybrid";   
   else if (calculatorType == 2 || calculatorType == 3) { 
      typeName = "Asymptotic";
      mPlotHypoTestResult = false; 
   }
  
   const char * resultName = r->GetName();
   TString plotTitle = TString::Format("%s CL Scan for workspace %s",typeName.c_str(),resultName);
   HypoTestInverterPlot *plot = new HypoTestInverterPlot("HTI_Result_Plot",plotTitle,r);

   // plot in a new canvas with style
   TString c1Name = TString::Format("%s_Scan",typeName.c_str());
   TCanvas * c1 = new TCanvas(c1Name); 
   c1->SetLogy(false);

   plot->Draw("CLb 2CL");  // plot all and Clb

   // if (useCLs) 
   //    plot->Draw("CLb 2CL");  // plot all and Clb
   // else 
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  grP->Draw("AP");
  //grC->Draw("same P");
  gr2->Draw("same P");
  grF->Draw("same");
  grFmax->Draw("same");
  grFmin->Draw("same");
  gr0->Draw("same P");
 // grPr->Draw("same P");
  //grPPr->Draw("same P");

  TLegend* leg2 = new TLegend(0.15,0.68,0.48,0.85);
  leg2->AddEntry(gr0,"STAR Run 12 - Low p_{T} Analysis","pe");
  leg2->AddEntry(gr2,"STAR Run 12 - High p_{T} Analysis","pe");
  //leg2->AddEntry(grC,"Combined Trigs","pe");
  leg2->AddEntry(grP,"STAR Run 6 Analysis (Stat. Uncertainty)","pe");
//  leg2->AddEntry(grPr,"Run 12 Data, Run 5/6 Templates)","pe");
  //leg2->AddEntry(grPPr,"Run 5/6 Refit (prev Template)","pe");
  leg2->AddEntry(grF,"FONLL Calculation","l");
  leg2->Draw("same");

  // Write to Root File if open
  if(makeROOT){
    file3->Close();
    file4->Close();
    file->cd();
    grP->Write("PreviousData");
    //grC->Write("same P");
    gr2->Write("HT2");
    grF->Write("FONLL");
    grFmax->Write("FONLLmax");
    grFmin->Write("FONLLmin");
    gr0->Write("HT0");
    // grPr->Write("PrevTempMyData");
    //grPPr->Write("PrevTempPreData");
  }

  // Make PDF with output canvases
  if(makePDF)
  {
    //Set front page
    TCanvas* fp = new TCanvas("fp","Front Page",100,0,1000,900);
    fp->cd();
    TBox *bLabel = new TBox(0.01, 0.88, 0.99, 0.99);
    bLabel->SetFillColor(38);
    bLabel->Draw();
    TLatex tl;
    tl.SetNDC();
    tl.SetTextColor(kWhite);
    tl.SetTextSize(0.033);
    char tlName[100];
    char tlName2[100];

    TString titlename = FileName;
    int found = titlename.Last('/');
    if(found >= 0){
      titlename.Replace(0, found+1, "");
    } 
    sprintf(tlName, "RUN 12 NPE-h   #Delta#phi Correlations");
    tl.SetTextSize(0.05);
    tl.SetTextColor(kWhite);
    tl.DrawLatex(0.05, 0.92,tlName);

    TBox *bFoot = new TBox(0.01, 0.01, 0.99, 0.12);
    bFoot->SetFillColor(38);
    bFoot->Draw();
    tl.SetTextColor(kWhite);
    tl.SetTextSize(0.05);
    tl.DrawLatex(0.05, 0.05, (new TDatime())->AsString());
    tl.SetTextColor(kBlack);
    tl.SetTextSize(0.03);
    tl.DrawLatex(0.1, 0.14, titlename);
    sprintf(tlName,"TEST");
    tl.DrawLatex(0.1, 0.8,tlName);

    // Place canvases in order
    TCanvas* temp = new TCanvas();
    sprintf(name, "FFOutput/%s.pdf[", FileName);
    temp->Print(name);
    sprintf(name, "FFOutput/%s.pdf", FileName);

    temp = deltaPhi; 
    temp->Print(name);
    temp = fitResult0;
    temp->Print(name);
    temp = fitResult2;
    temp->Print(name);
    // temp = fitResultC;
    // temp->Print(name);
    temp = c1;
    temp->Print(name);

    sprintf(name, "FFOutput/%s.pdf]", FileName);
    temp->Print(name);
  }
  if(makeROOT)
  {
    file->Write();
    file->Close();
  }
}

//___________________________________________________________________
Int_t studyBinZero(TString pathData, TString fileNameData, TString listName = "",
                   Double_t etaLow = -0.9, Double_t etaUp = 0.9,
                   Double_t lowerCentrality = -2, Double_t upperCentrality = -2)
{
  PrintSettingsAxisRangeForMultiplicityAxisForMB();
  
  TString pathNameData = Form("%s/%s", pathData.Data(), fileNameData.Data());
  
  if (listName == "") {
    listName = pathNameData;
    listName.Replace(0, listName.Last('/') + 1, "");
    listName.ReplaceAll(".root", "");
  }
  
  TObjArray* histList = 0x0;
  
  TFile* f = TFile::Open(pathNameData.Data());
  if (!f)  {
    std::cout << std::endl;
    std::cout << "Failed to open file \"" << pathNameData.Data() << "\"!" << std::endl;
    return -1;
  }
  
  histList = (TObjArray*)(f->Get(listName.Data()));
  if (!histList) {
    std::cout << std::endl;
    std::cout << "Failed to load list \"" << listName.Data() << "\"!" << std::endl;
    return -1;
  }
  
  // Extract the data histograms
  TH1* hNumEventsTriggerSel = dynamic_cast<TH1*>(histList->FindObject("fhEventsTriggerSel"));
  TH1* hNumEventsTriggerSelVtxCut = dynamic_cast<TH1*>(histList->FindObject("fhEventsTriggerSelVtxCut"));
  TH1* hNumEventsTriggerSelVtxCutZ= dynamic_cast<TH1*>(histList->FindObject("fhEventsProcessedNoPileUpRejection"));
  TH1* hNumEventsTriggerSelVtxCutZPileUpRej= dynamic_cast<TH1*>(histList->FindObject("fhEventsProcessed"));
  
  THnSparse* hDataTriggerSel = dynamic_cast<THnSparse*>(histList->FindObject("fChargedGenPrimariesTriggerSel"));
  THnSparse* hDataTriggerSelVtxCut = dynamic_cast<THnSparse*>(histList->FindObject("fChargedGenPrimariesTriggerSelVtxCut"));
  THnSparse* hDataTriggerSelVtxCutZ = dynamic_cast<THnSparse*>(histList->FindObject("fChargedGenPrimariesTriggerSelVtxCutZ"));
  THnSparse* hDataTriggerSelVtxCutZPileUpRej = dynamic_cast<THnSparse*>(histList->FindObject("fChargedGenPrimariesTriggerSelVtxCutZPileUpRej"));
  
  setSparseErrors(hDataTriggerSel);
  setSparseErrors(hDataTriggerSelVtxCut);
  setSparseErrors(hDataTriggerSelVtxCutZ);
  setSparseErrors(hDataTriggerSelVtxCutZPileUpRej);
  
  
  
  // Set multiplicity range, if desired. Note that mult and pT axes are the same for all histos
  Int_t lowerCentralityBinLimit = -1;
  Int_t upperCentralityBinLimit = -2;
  Bool_t restrictCentralityAxis = kFALSE;
  Double_t actualLowerCentrality = -1.;
  Double_t actualUpperCentrality = -1.;
  
  if (lowerCentrality >= -1 && upperCentrality >= -1) {
    // Add subtract a very small number to avoid problems with values right on the border between to bins
    lowerCentralityBinLimit = hNumEventsTriggerSel->GetXaxis()->FindFixBin(lowerCentrality + 0.001);
    upperCentralityBinLimit = hNumEventsTriggerSel->GetXaxis()->FindFixBin(upperCentrality - 0.001);
    
    // Check if the values look reasonable
    if (lowerCentralityBinLimit <= upperCentralityBinLimit && lowerCentralityBinLimit >= 1
        && upperCentralityBinLimit <= hNumEventsTriggerSel->GetXaxis()->GetNbins()) {
      actualLowerCentrality = hNumEventsTriggerSel->GetXaxis()->GetBinLowEdge(lowerCentralityBinLimit);
      actualUpperCentrality = hNumEventsTriggerSel->GetXaxis()->GetBinUpEdge(upperCentralityBinLimit);

      restrictCentralityAxis = kTRUE;
    }
    else {
      std::cout << std::endl;
      std::cout << "Requested centrality range out of limits or upper and lower limit are switched!" << std::endl;
      return -1;
    }
  }
  
  if (!restrictCentralityAxis)
    GetAxisRangeForMultiplicityAxisForMB(hNumEventsTriggerSel->GetXaxis(), lowerCentralityBinLimit, upperCentralityBinLimit);
  
  hNumEventsTriggerSel->GetXaxis()->SetRange(lowerCentralityBinLimit, upperCentralityBinLimit);
  
  actualLowerCentrality = hNumEventsTriggerSel->GetXaxis()->GetBinLowEdge(hNumEventsTriggerSel->GetXaxis()->GetFirst());
  actualUpperCentrality = hNumEventsTriggerSel->GetXaxis()->GetBinUpEdge(hNumEventsTriggerSel->GetXaxis()->GetLast());
  
  std::cout << "centrality: ";
  if (restrictCentralityAxis) {
    std::cout << actualLowerCentrality << " - " << actualUpperCentrality << std::endl;
    std::cout << "WARNING: Does it really make sense to restrict the centrality? Usually, centrality estimators imply centrality >= 0!"
              << std::endl;
  }
  else
    std::cout << "MB (" << actualLowerCentrality << " - " << actualUpperCentrality << ")" << std::endl;
    
  if (restrictCentralityAxis) {
    hDataTriggerSel->GetAxis(iMult)->SetRange(lowerCentralityBinLimit, upperCentralityBinLimit);
    hDataTriggerSelVtxCut->GetAxis(iMult)->SetRange(lowerCentralityBinLimit, upperCentralityBinLimit);
    hDataTriggerSelVtxCutZ->GetAxis(iMult)->SetRange(lowerCentralityBinLimit, upperCentralityBinLimit);
    hDataTriggerSelVtxCutZPileUpRej->GetAxis(iMult)->SetRange(lowerCentralityBinLimit, upperCentralityBinLimit);
  }
  
//.........这里部分代码省略.........

//.........这里部分代码省略.........
        // Get Fit information and store to use in corrections
        TF1 *fitResult = projZDCxHad[ptbin][trig]->GetFunction("pol1");
        pu[0][ptbin][trig] = fitResult->GetParameter(0);
        pu[1][ptbin][trig] = fitResult->GetParameter(1);
           	cout << trig << " " << ptbin << ": " << pu[0][ptbin][trig] << " " << pu[1][ptbin][trig] << endl;
        }*/
    }

    // Draw on "SinglePlot" canvas for saving single plots from grid
    TPad* pNew = (TPad*)result[2]->GetPad(4)->Clone();
    singlePlot->cd();
    pNew->ResizePad();
    pNew->Draw();

    // Make PDF with output canvases
    if(makePDF)
    {
        //Set front page
        TCanvas* fp = new TCanvas("fp","Front Page",100,0,1000,900);
        fp->cd();
        TBox *bLabel = new TBox(0.01, 0.88, 0.99, 0.99);
        bLabel->SetFillColor(38);
        bLabel->Draw();
        TLatex tl;
        tl.SetNDC();
        tl.SetTextColor(kWhite);
        tl.SetTextSize(0.033);
        char tlName[100];
        char tlName2[100];

        TString titlename = FileName;
        int found = titlename.Last('/');
        if(found >= 0) {
            titlename.Replace(0, found+1, "");
        }
        sprintf(tlName, "RUN 12 pp 200 GeV NPE-h    #Delta#phi Analysis");
        tl.SetTextSize(0.05);
        tl.SetTextColor(kWhite);
        tl.DrawLatex(0.05, 0.92,tlName);

        TBox *bFoot = new TBox(0.01, 0.01, 0.99, 0.12);
        bFoot->SetFillColor(38);
        bFoot->Draw();
        tl.SetTextColor(kWhite);
        tl.SetTextSize(0.05);
        tl.DrawLatex(0.05, 0.05, (new TDatime())->AsString());
        tl.SetTextColor(kBlack);
        tl.SetTextSize(0.03);
        tl.DrawLatex(0.1, 0.14, titlename);
        sprintf(tlName,"eID: -1 < n  #sigma_{e TPC} < 3;  #left|gDCA #right| < 1 cm; 0.3 < p/E < 1.5;");
        tl.DrawLatex(0.1, 0.8,tlName);
        sprintf(tlName,"       nHitsFit > 20; nHits   #frac{dE}{dx} > 15; nHitFit/Max > 0.52;    #left|#eta#right| < 0.7;");
        tl.DrawLatex(0.1, 0.75,tlName);
        sprintf(tlName,"       n #phi > 1; n #eta > 1;  #left|dZ#right| < 3 cm;  #left|d#phi#right| < 0.015;");
        tl.DrawLatex(0.1, 0.7,tlName);
        sprintf(tlName,"hID: p_{T} > 0.5;  #left|#eta#right| < 1; nHitsFit > 15; nHits   #frac{dE}{dx} > 10; DCA < 1 cm;");
        tl.DrawLatex(0.1, 0.6,tlName);
        sprintf(tlName,"Event:  #left|V_{z}#right| < 35 cm;");
        tl.DrawLatex(0.1, 0.5,tlName);
        sprintf(tlName,"Triggers:  BHT0; BHT2;");
        tl.DrawLatex(0.1, 0.4,tlName);


        // Place canvases in order
        TCanvas* temp = new TCanvas();
        sprintf(name, "%s.pdf[", FileName);

//___________________________________________________________________
Int_t extractPtResolution(TString pathNameData, TString listName,
                          Int_t chargeMode /*kNegCharge = -1, kAllCharged = 0, kPosCharge = 1*/,
                          Double_t lowerCentrality /*= -2*/, Double_t upperCentrality /*= -2*/,
                          Double_t lowerJetPt /*= -1*/ , Double_t upperJetPt/* = -1*/)
{
  if (listName == "") {
    listName = pathNameData;
    listName.Replace(0, listName.Last('/') + 1, "");
    listName.ReplaceAll(".root", "");
  }
  
  TString pathData = pathNameData;
  pathData.Replace(pathData.Last('/'), pathData.Length(), "");
  
  TH1D* hPtResolutionFit[AliPID::kSPECIES] = {0x0, };
  TH2D* hPtResolution[AliPID::kSPECIES] = {0x0, };
  THnSparse* hPtResolutionRaw[AliPID::kSPECIES] = {0x0, };

  TFile* fileData = TFile::Open(pathNameData.Data());
  if (!fileData) {
    printf("Failed to open data file \"%s\"\n", pathNameData.Data());
    return -1;
  }
  
  TObjArray* histList = (TObjArray*)(fileData->Get(listName.Data()));
  
  if (!histList) {
    printf("Failed to load list!\n");
    return -1;
  }
  
  Double_t actualLowerCentrality = -2;
  Double_t actualUpperCentrality = -2;
  
  Double_t actualLowerJetPt = -1.;
  Double_t actualUpperJetPt = -1.;
  
  Bool_t restrictJetPtAxis = (lowerJetPt >= 0 && upperJetPt >= 0);
  const Bool_t restrictCentrality = ((lowerCentrality >= -1) && (upperCentrality >= -1));
  
  for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
    const TString sparseName = Form("fPtResolution_%s", AliPID::ParticleShortName(species));
    hPtResolutionRaw[species] = (THnSparse*)histList->FindObject(sparseName.Data());
    
    if (!hPtResolutionRaw[species]) {
      printf("Failed to load THnSparse for %s: %s!\n", AliPID::ParticleShortName(species), sparseName.Data());
      return -1;
    }
    
    // Set proper errors, if not yet calculated
    if (!hPtResolutionRaw[species]->GetCalculateErrors()) {
      std::cout << "Re-calculating errors of " << hPtResolutionRaw[species]->GetName() << "..." << std::endl;
      
      hPtResolutionRaw[species]->Sumw2();
      
      Long64_t nBinsPtResolutionRaw = hPtResolutionRaw[species]->GetNbins();
      Double_t binContent = 0;
      for (Long64_t bin = 0; bin < nBinsPtResolutionRaw; bin++) {
        binContent = hPtResolutionRaw[species]->GetBinContent(bin);
        hPtResolutionRaw[species]->SetBinError(bin, TMath::Sqrt(binContent));
      }
    }
    
    // Integral(lowerCentBinLimit, uppCentBinLimit) will not be restricted if these values are kept
    const Int_t lowerCentralityBinLimit = restrictCentrality ? hPtResolutionRaw[species]->GetAxis(kPtResCentrality)->FindBin(lowerCentrality + 0.001) 
                                                                    : -1;
    const Int_t upperCentralityBinLimit = restrictCentrality ? hPtResolutionRaw[species]->GetAxis(kPtResCentrality)->FindBin(upperCentrality - 0.001) 
                                                                    : -2;
    
    
    if (restrictCentrality) {
      actualLowerCentrality = hPtResolutionRaw[species]->GetAxis(kPtResCentrality)->GetBinLowEdge(lowerCentralityBinLimit);
      actualUpperCentrality = hPtResolutionRaw[species]->GetAxis(kPtResCentrality)->GetBinLowEdge(upperCentralityBinLimit);
      
      hPtResolutionRaw[species]->GetAxis(kPtResCentrality)->SetRange(lowerCentralityBinLimit, upperCentralityBinLimit);
    }
    
    const Bool_t restrictCharge = (chargeMode != kAllCharged);
    
    Int_t lowerChargeBinLimit = -1;
    Int_t upperChargeBinLimit = -2;
      
    if (restrictCharge) {
      // Add subtract a very small number to avoid problems with values right on the border between to bins
      if (chargeMode == kNegCharge) {
        lowerChargeBinLimit = hPtResolutionRaw[species]->GetAxis(kPtResCharge)->FindBin(-1. + 0.001);
        upperChargeBinLimit = hPtResolutionRaw[species]->GetAxis(kPtResCharge)->FindBin(0. - 0.001);
      }
      else if (chargeMode == kPosCharge) {
        lowerChargeBinLimit = hPtResolutionRaw[species]->GetAxis(kPtResCharge)->FindBin(0. + 0.001);
        upperChargeBinLimit = hPtResolutionRaw[species]->GetAxis(kPtResCharge)->FindBin(1. - 0.001);
      }
      
      // Check if the values look reasonable
      if (lowerChargeBinLimit <= upperChargeBinLimit && lowerChargeBinLimit >= 1
          && upperChargeBinLimit <= hPtResolutionRaw[species]->GetAxis(kPtResCharge)->GetNbins()) {
        // OK
      }
      else {
//.........这里部分代码省略.........

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

  gr0->SetTitle("Bottom Contribution");
  gr0->GetXaxis()->SetTitle("p_{T,e}");
  gr0->GetYaxis()->SetTitle("#frac{r_{B}}{(r_{B}+r_{C})}");
  gr0->SetMarkerStyle(20);
  gr0->SetMarkerSize(1);
  gr0->SetLineColor(kBlue);
  gr0->SetMarkerColor(kBlue);
  gr2->SetMarkerStyle(22);
  gr2->SetMarkerSize(1);
  gr2->SetLineColor(kRed);
  gr2->SetMarkerColor(kRed);
  grC->SetMarkerStyle(21);
  grC->SetMarkerSize(1);
  grC->SetLineColor(kRed);
  grC->SetMarkerColor(kRed);
  gr0->GetXaxis()->SetLimits(1,14);
  gr0->GetYaxis()->SetRangeUser(0,1);
  grF->SetLineStyle(1);
  grFmax->SetLineStyle(2);
  grFmin->SetLineStyle(2);
  grP->SetMarkerStyle(33);
  grP->SetMarkerColor(kBlack);
  
  
  gr0->Draw("AP");
  // grC->Draw("same P");
  gr2->Draw("same P");
  grF->Draw("same");
  grFmax->Draw("same");
  grFmin->Draw("same");
  grP->Draw("same P");

  TLegend* leg2 = new TLegend(0.15,0.68,0.4,0.85);
  leg2->AddEntry(gr0,"High Tower 0 Trigs","pe");
  leg2->AddEntry(gr2,"High Tower 2 Trigs","pe");
  // leg2->AddEntry(grC,"Combined Trigs","pe");
  leg2->AddEntry(grF,"FONLL (Uncertainty: Scale Only)","l");
  leg2->AddEntry(grP,"Run 5/6 Analysis (Stat Uncertainty)","pe");
  leg2->Draw("same");
  
   // Make PDF with output canvases
  if(makePDF)
    {
      //Set front page
      TCanvas* fp = new TCanvas("fp","Front Page",100,0,1000,900);
      fp->cd();
      TBox *bLabel = new TBox(0.01, 0.88, 0.99, 0.99);
      bLabel->SetFillColor(38);
      bLabel->Draw();
      TLatex tl;
      tl.SetNDC();
      tl.SetTextColor(kWhite);
      tl.SetTextSize(0.033);
      char tlName[100];
      char tlName2[100];
      
      TString titlename = FileName;
      int found = titlename.Last('/');
      if(found >= 0){
	titlename.Replace(0, found+1, "");
      } 
      sprintf(tlName, "RUN 12 NPE-h   #Delta#phi Correlations");
      tl.SetTextSize(0.05);
      tl.SetTextColor(kWhite);
      tl.DrawLatex(0.05, 0.92,tlName);
      
      TBox *bFoot = new TBox(0.01, 0.01, 0.99, 0.12);
      bFoot->SetFillColor(38);
      bFoot->Draw();
      tl.SetTextColor(kWhite);
      tl.SetTextSize(0.05);
      tl.DrawLatex(0.05, 0.05, (new TDatime())->AsString());
      tl.SetTextColor(kBlack);
      tl.SetTextSize(0.03);
      tl.DrawLatex(0.1, 0.14, titlename);
      sprintf(tlName,"TEST");
      tl.DrawLatex(0.1, 0.8,tlName);
      
      // Place canvases in order
      TCanvas* temp = new TCanvas();
      sprintf(name, "FFOutput/%s.pdf[", FileName);
      temp->Print(name);
      sprintf(name, "FFOutput/%s.pdf", FileName);

      temp = deltaPhi; 
      temp->Print(name);
      temp = fitResult0;
      temp->Print(name);
      temp = fitResult2;
      temp->Print(name);
      temp = fitResultC;
      temp->Print(name);
      temp = c1;
      temp->Print(name);
      
      sprintf(name, "FFOutput/%s.pdf]", FileName);
      temp->Print(name);
    }
}

void
RooStats::HypoTestInvTool::AnalyzeResult( HypoTestInverterResult * r,
                                          int calculatorType,
                                          int testStatType, 
                                          bool useCLs,  
                                          int npoints,
                                          const char * fileNameBase ){

   // analyize result produced by the inverter, optionally save it in a file 
   
   double upperLimit = r->UpperLimit();
   double ulError = r->UpperLimitEstimatedError();
  
   std::cout << "The computed upper limit is: " << upperLimit << " +/- " << ulError << std::endl;
  
   // compute expected limit
   std::cout << " expected limit (median) " << r->GetExpectedUpperLimit(0) << std::endl;
   std::cout << " expected limit (-1 sig) " << r->GetExpectedUpperLimit(-1) << std::endl;
   std::cout << " expected limit (+1 sig) " << r->GetExpectedUpperLimit(1) << std::endl;
   std::cout << " expected limit (-2 sig) " << r->GetExpectedUpperLimit(-2) << std::endl;
   std::cout << " expected limit (+2 sig) " << r->GetExpectedUpperLimit(2) << std::endl;
  
  
   // write result in a file 
   if (r != NULL && mWriteResult) {
    
      // write to a file the results
      const char *  calcType = (calculatorType == 0) ? "Freq" : (calculatorType == 1) ? "Hybr" : "Asym";
      const char *  limitType = (useCLs) ? "CLs" : "Cls+b";
      const char * scanType = (npoints < 0) ? "auto" : "grid";
      TString resultFileName = TString::Format("%s_%s_%s_ts%d_",calcType,limitType,scanType,testStatType);      
      //strip the / from the filename
      if (mMassValue.size()>0) {
         resultFileName += mMassValue.c_str();
         resultFileName += "_";
      }
    
      TString name = fileNameBase; 
      name.Replace(0, name.Last('/')+1, "");
      resultFileName += name;
    
      TFile * fileOut = new TFile(resultFileName,"RECREATE");
      r->Write();
      fileOut->Close();                                                                     
   }   
  
  
   // plot the result ( p values vs scan points) 
   std::string typeName = "";
   if (calculatorType == 0 )
      typeName = "Frequentist";
   if (calculatorType == 1 )
      typeName = "Hybrid";   
   else if (calculatorType == 2 ) { 
      typeName = "Asymptotic";
      mPlotHypoTestResult = false; 
   }
  
   const char * resultName = r->GetName();
   TString plotTitle = TString::Format("%s CL Scan for workspace %s",typeName.c_str(),resultName);
   HypoTestInverterPlot *plot = new HypoTestInverterPlot("HTI_Result_Plot",plotTitle,r);
   plot->Draw("CLb 2CL");  // plot all and Clb
  
   const int nEntries = r->ArraySize();
  
   // plot test statistics distributions for the two hypothesis 
   if (mPlotHypoTestResult) { 
      TCanvas * c2 = new TCanvas();
      if (nEntries > 1) { 
         int ny = TMath::CeilNint( sqrt(nEntries) );
         int nx = TMath::CeilNint(double(nEntries)/ny);
         c2->Divide( nx,ny);
      }
      for (int i=0; i<nEntries; i++) {
         if (nEntries > 1) c2->cd(i+1);
         SamplingDistPlot * pl = plot->MakeTestStatPlot(i);
         pl->SetLogYaxis(true);
         pl->Draw();
      }
   }
}

void write_ist_pednoise() {
	gSystem->Setenv("DB_ACCESS_MODE", "write");
  	gROOT->Macro("LoadLogger.C");
  	gSystem->Load("St_base.so");
  	gSystem->Load("libStDb_Tables.so");
  	gSystem->Load("StDbLib.so");

	std::ifstream fPDlist("/star/u/ypwang/disk01/offline_hft/DB/pedNoiseTable4DB/pedNoiseFiles_Corr/pedRunList_Corr.txt");
        if(!fPDlist.is_open())
        {
                std::cout << " There is no corrected pedRunList file! " << endl;
                exit(0);
        }

        while(!fPDlist.eof()) {
                TString  sPDfile;
                fPDlist >> sPDfile;

                int found = sPDfile.Last('_');
                if(found >= 0)
                        sPDfile.Replace(0, found + 1, "");

                found = sPDfile.Last('.');
                if(found >= 0)
                        sPDfile.Replace(found, sPDfile.Length(), "");

                int runnumb = sPDfile.Atoi();
                if(runnumb)
                        pedRunNumVec.push_back(runnumb);
        }
        fPDlist.close();
	const int nPedRuns = pedRunNumVec.size();
	cout << nPedRuns << " good pedestal runs recorded! " << endl;
        if(debug) {
                for(int i=0; i<nPedRuns; i++)
                        cout << i << "th pedestal run: " << pedRunNumVec[i] << endl;
        }


  	std::ifstream fCosmicRunList("/star/u/ypwang/disk01/offline_hft/DB/pedNoiseTable4DB/cmNoise_cosmic/cosmicRunlist.txt");
        if(!fCosmicRunList.is_open())
        {
                std::cout << " There is no cosmicRunList file! " << endl;
                exit(0);
        }

        while(!fCosmicRunList.eof()) {
                TString  sCMNfile;
                fCosmicRunList >> sCMNfile;

                int found = sCMNfile.First('.');
                if(found >= 0)
                        sCMNfile.Replace(found, sCMNfile.Length(), "");

                int runnumb = sCMNfile.Atoi();
                if(runnumb)
                        cosmicRunNumVec.push_back(runnumb);
        }
        fCosmicRunList.close();
	const int nCosmicRuns = cosmicRunNumVec.size();
	cout << nCosmicRuns << " good cosmic runs recorded! " << endl;
        if(debug) {
                for(int i=0; i<nCosmicRuns; i++)
                        cout << i << "th cosmic run: " << cosmicRunNumVec[i] << endl;
        }

	Int_t   isBadChannel[110592];
        Int_t   isBadChip[864];
	//loop pedestal runs
	for(int iRun=0; iRun<nPedRuns; iRun++)
        {
                int runNumber = pedRunNumVec[iRun];
		//below only for testing
		//if(runNumber<15065020) continue;

                int TbinPD = PDTBin(runNumber);
		cout << "Reading pedestal run " << runNumber << " with " << TbinPD << " time bins! " << endl;
		
		//retrieve time stamp information
		char cmd[256];
            	sprintf(cmd, "/usr/bin/mysql -h dbbak.starp.bnl.gov --port=3413 -e \"SELECT FROM_UNIXTIME(startRunTime) as beginTime_human FROM RunLog.runDescriptor WHERE runNumber = %d\" >timeStamp_tmp", runNumber);
            	system(cmd);

            	TString mTimeStamp;
            	FILE *inR = fopen("timeStamp_tmp","r");
            	while(!feof(inR)) {
                	char buff[64];
                	if(fgets(buff,sizeof(buff),inR) == 0) continue ;
                	switch(buff[0]) {
                    		case 'b' : //beginTime_human
                        		continue ;
                	}
                	mTimeStamp = Form(buff);
            	}
            	fclose(inR);

            	sprintf(cmd, "/bin/rm timeStamp_tmp");
            	system(cmd);

		int found = 19;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  gPad-> SetLogy();
  ptHat->GetXaxis()->SetTitle("pT-Hat (GeV/c)");
  ptHat->SetTitle("Raw pT Hat");
  ptHat->Draw();
  ptHatC->cd(2);
  gPad-> SetLogy();
  ptHatCorr->GetXaxis()->SetTitle("pT-Hat (GeV/c)");
  ptHatCorr->SetTitle("Weighted pT Hat");
  ptHatCorr->Draw();
  
  for(Int_t ptbin=0; ptbin<numPtBins; ptbin++)
    {
      // Init necessary plotting tools
      lbl[ptbin] = new TPaveText(.2,.8,.5,.85,Form("NB NDC%i",ptbin));
      sprintf(textLabel,"%.1f < P_{T,e} < %.1f",lowpt[ptbin],highpt[ptbin]);
      lbl[ptbin]->AddText(textLabel);
      lbl[ptbin]->SetFillColor(kWhite);

      deltaPhi->cd(ptbin+1);
      delPhi[ptbin]->GetXaxis()->SetTitle("#Delta#phi_{eh}");
      // delPhi[ptbin]->Sumw2();
      //cout << totalNorm[ptbin] << endl;
      //delPhi[ptbin]->Scale(wt);
      delPhi[ptbin]->GetYaxis()->SetTitle("1/N_{NPE} #upoint dN/d(#Delta)#phi");
      delPhi[ptbin]->GetXaxis()->SetRangeUser(-3.5,3.5);
      if(ptbin == 0)
	{
	  if(mode == 1)
	    delPhi[ptbin]->SetTitle("Pythia NPE-had #Delta#phi - c/#bar{c}");
	  if(mode == 2)
	    delPhi[ptbin]->SetTitle("Pythia NPE-had #Delta#phi - b/#bar{b}");
	}
      else
	delPhi[ptbin]->SetTitle("");
      if(ptbin < 13){
	delPhi[ptbin]->Draw("E");
	lbl[ptbin]->Draw("same");
      }
    }
  
      
  // Make PDF with output canvases
  if(makePDF)
    {
      //Set front page
      TCanvas* fp = new TCanvas("fp","Front Page",100,0,1000,900);
      fp->cd();
      TBox *bLabel = new TBox(0.01, 0.88, 0.99, 0.99);
      bLabel->SetFillColor(38);
      bLabel->Draw();
      TLatex tl;
      tl.SetNDC();
      tl.SetTextColor(kWhite);
      tl.SetTextSize(0.033);
      char tlName[100];
      char tlName2[100];
      
      TString titlename = FileName;
      int found = titlename.Last('/');
      if(found >= 0){
	titlename.Replace(0, found+1, "");
      } 
      sprintf(tlName, "RUN 12 NPE-h   #Delta#phi Pythia Templates");
      tl.SetTextSize(0.05);
      tl.SetTextColor(kWhite);
      tl.DrawLatex(0.05, 0.92,tlName);
      
      TBox *bFoot = new TBox(0.01, 0.01, 0.99, 0.12);
      bFoot->SetFillColor(38);
      bFoot->Draw();
      tl.SetTextColor(kWhite);
      tl.SetTextSize(0.05);
      tl.DrawLatex(0.05, 0.05, (new TDatime())->AsString());
      tl.SetTextColor(kBlack);
      tl.SetTextSize(0.03);
      tl.DrawLatex(0.1, 0.14, titlename);
      sprintf(tlName,"TEST");
      tl.DrawLatex(0.1, 0.8,tlName);
      
      // Place canvases in order
      TCanvas* temp = new TCanvas();
      sprintf(name, "%s.pdf[", FileName);
      temp->Print(name);
      sprintf(name, "%s.pdf", FileName);
      temp = fp; // print front page
      temp->Print(name);
      temp = ptHatC;
      temp->Print(name);
      temp = deltaPhi;
      temp->Print(name);
      sprintf(name, "%s.pdf]", FileName);
      temp->Print(name);
    }

  if(makeROOT)
    {
      file->Write();
      file->Close();
    }
}

void TSelectionClass::saveHistos()
{
	cout<<"save Histo: "<<hFiducialCutSilicon->GetTitle()<<endl;
	TString name = hFiducialCutSilicon->GetName();
	name.Insert(0,"c");
	TCanvas *c1= fiducialCuts->getAllFiducialCutsCanvas(hFiducialCutSilicon);
	c1->SetName(name);
	histSaver->SaveCanvas(c1);
	delete c1;
	c1 = 0;
    delete hFiducialCutSilicon;

	name = hFiducialCutSiliconDiamondHit->GetName();
	name.Insert(0,"c");
	c1 = fiducialCuts->getAllFiducialCutsCanvas(hFiducialCutSiliconDiamondHit,true);
	c1->SetName(name);
	histSaver->SaveCanvas(c1);
	delete c1;
	c1=0;
    delete hFiducialCutSiliconDiamondHit;

	name = "c";
	name.Append(hFiducialCutSiliconOneAndOnlyOneDiamondHit->GetName());
	c1= fiducialCuts->getAllFiducialCutsCanvas(hFiducialCutSiliconOneAndOnlyOneDiamondHit,true);
	c1->SetName(name);
	histSaver->SaveCanvas(c1);
    delete c1;
    c1=0;
    delete hFiducialCutSiliconOneAndOnlyOneDiamondHit;

	name = "c";
	name.Append(hSelectedEvents->GetName());
	c1 = fiducialCuts->getAllFiducialCutsCanvas(hSelectedEvents,true);
	c1->SetName(name);
	histSaver->SaveCanvas(c1);
    delete c1;
    c1=0;
    delete hSelectedEvents;

    map<Int_t,TH2F*>::iterator it;
    for (it = mapFiducialCutSiliconDiamondHitSamePattern.begin(); it!=mapFiducialCutSiliconDiamondHitSamePattern.end(); it++){
        TH2F* histo = (*it).second;
        name = histo->GetName();
        name.Replace(0,1,"c");
        c1 = fiducialCuts->getAllFiducialCutsCanvas(histo,true);
        c1->SetName(name);
        histSaver->SaveCanvas(c1);
        delete c1;
        c1=0;
        delete histo;
    }
	hAnalysisFraction->Scale(.1);
	hAnalysisFraction->SetStats(false);
	//	hAnalysisFraction->GetYaxis()->SetRangeUser(0,100);
	histSaver->SaveHistogram(hAnalysisFraction);
	delete hAnalysisFraction;


    histSaver->SaveHistogram(hDiamondPatternFiducialPattern);

    histSaver->SaveHistogram(hDiamondPatternFiducialPatternNoMapping);
    name = "stackPatternMapping";
    hDiamondPatternFiducialPattern->SetLineColor(kGreen);
    hDiamondPatternFiducialPatternNoMapping->SetLineColor(kRed);
    THStack *stack = new THStack(name,name);
    stack->Add(hDiamondPatternFiducialPattern);
    stack->Add(hDiamondPatternFiducialPatternNoMapping);
    stack->Draw();
    if(stack->GetXaxis()) stack->GetXaxis()->SetTitle("pattern no.");
    if(stack->GetYaxis()) stack->GetYaxis()->SetTitle("number of entries #");
    histSaver->SaveStack(stack,"hist",true);
    if(stack) delete stack;
    histSaver->SaveHistogram(pDiamondPatternFiducialPatternProfile);
    if(hDiamondPatternFiducialPatternNoMapping) delete hDiamondPatternFiducialPatternNoMapping;
    if(hDiamondPatternFiducialPattern) delete hDiamondPatternFiducialPattern;
    if (pDiamondPatternFiducialPatternProfile) delete pDiamondPatternFiducialPatternProfile;
}