C++ TFitResultPtr::Get方法代码示例

void plotPedestalAnalysis(string inputFileName, string outputDIR, bool testDoubleGaussianChannels){

  system(("mkdir -p "+outputDIR).c_str());

  setTDRStyle();
  gROOT->SetBatch(kTRUE);

  TFile* inputFile = TFile::Open(inputFileName.c_str(),"READ");
  inputFile->cd();
  TTree* tree = (TTree*) inputFile->Get("pedestalFullNoise");

  uint32_t detid,fedKey;
  uint16_t fecCrate,fecSlot, fecRing, ccuAdd, ccuChan, lldChannel, fedId, fedCh, apvId, stripId;
  float    fitChi2Probab, kSProbab, jBProbab, aDProbab, fitChi2;
  float    noiseSkewness, noiseKurtosis;
  float    fitGausMean, fitGausSigma, fitGausNormalization;
  float    fitGausMeanError, fitGausSigmaError, fitGausNormalizationError;
  vector<float>* noiseDistribution = 0;
  vector<float>* noiseDistributionError = 0;
  float    nBin, xMin, xMax;

  tree->SetBranchStatus("*",kFALSE);
  tree->SetBranchStatus("detid",kTRUE);
  tree->SetBranchStatus("fedKey",kTRUE);
  tree->SetBranchStatus("fecCrate",kTRUE);
  tree->SetBranchStatus("fecSlot",kTRUE);
  tree->SetBranchStatus("fecRing",kTRUE);
  tree->SetBranchStatus("ccuAdd",kTRUE);
  tree->SetBranchStatus("ccuChan",kTRUE);
  tree->SetBranchStatus("lldChannel",kTRUE);
  tree->SetBranchStatus("fedId",kTRUE);
  tree->SetBranchStatus("fedCh",kTRUE);
  tree->SetBranchStatus("apvId",kTRUE);
  tree->SetBranchStatus("stripId",kTRUE);
  tree->SetBranchStatus("fitChi2",kTRUE);
  tree->SetBranchStatus("fitChi2Probab",kTRUE);
  tree->SetBranchStatus("kSProbab",kTRUE);
  tree->SetBranchStatus("jBProbab",kTRUE);
  tree->SetBranchStatus("aDProbab",kTRUE);
  tree->SetBranchStatus("fitGausNormalization",kTRUE);
  tree->SetBranchStatus("fitGausMean",kTRUE);
  tree->SetBranchStatus("fitGausSigma",kTRUE);
  tree->SetBranchStatus("fitGausNormalizationError",kTRUE);
  tree->SetBranchStatus("fitGausMeanError",kTRUE);
  tree->SetBranchStatus("fitGausSigmaError",kTRUE);
  tree->SetBranchStatus("noiseSkewness",kTRUE);
  tree->SetBranchStatus("noiseKurtosis",kTRUE);
  tree->SetBranchStatus("noiseDistribution",kTRUE);
  tree->SetBranchStatus("noiseDistributionError",kTRUE);
  tree->SetBranchStatus("nBin",kTRUE);
  tree->SetBranchStatus("xMin",kTRUE);
  tree->SetBranchStatus("xMax",kTRUE);

  tree->SetBranchAddress("detid",&detid);
  tree->SetBranchAddress("fedKey",&fedKey);
  tree->SetBranchAddress("fecCrate",&fecCrate);
  tree->SetBranchAddress("fecSlot",&fecSlot);
  tree->SetBranchAddress("fecRing",&fecRing);
  tree->SetBranchAddress("ccuAdd",&ccuAdd);
  tree->SetBranchAddress("ccuChan",&ccuChan);
  tree->SetBranchAddress("lldChannel",&lldChannel);
  tree->SetBranchAddress("fedId",&fedId);
  tree->SetBranchAddress("fedCh",&fedCh);
  tree->SetBranchAddress("apvId",&apvId);
  tree->SetBranchAddress("stripId",&stripId);
  tree->SetBranchAddress("fitGausNormalization",&fitGausNormalization);
  tree->SetBranchAddress("fitGausMean",&fitGausMean);
  tree->SetBranchAddress("fitGausSigma",&fitGausSigma);
  tree->SetBranchAddress("fitGausNormalizationError",&fitGausNormalizationError);
  tree->SetBranchAddress("fitGausMeanError",&fitGausMeanError);
  tree->SetBranchAddress("fitGausSigmaError",&fitGausSigmaError);
  tree->SetBranchAddress("fitChi2",&fitChi2);
  tree->SetBranchAddress("fitChi2Probab",&fitChi2Probab);
  tree->SetBranchAddress("noiseSkewness",&noiseSkewness);
  tree->SetBranchAddress("noiseKurtosis",&noiseKurtosis);
  tree->SetBranchAddress("kSProbab",&kSProbab);
  tree->SetBranchAddress("aDProbab",&aDProbab);
  tree->SetBranchAddress("jBProbab",&jBProbab);
  tree->SetBranchAddress("noiseDistribution",&noiseDistribution);
  tree->SetBranchAddress("noiseDistributionError",&noiseDistributionError);
  tree->SetBranchAddress("nBin",&nBin);
  tree->SetBranchAddress("xMin",&xMin);
  tree->SetBranchAddress("xMax",&xMax);

  TFile* badStripsNFilledBins = new TFile((outputDIR+"/badStripsNFilledBins.root").c_str(),"RECREATE");
  TFile* badKsTest = new TFile((outputDIR+"/badStripsKsTest.root").c_str(),"RECREATE");
  TFile* badjBTest = new TFile((outputDIR+"/badStripsjBTest.root").c_str(),"RECREATE");
  TFile* badChi2Test = new TFile((outputDIR+"/badStripsChi2Test.root").c_str(),"RECREATE");
  TFile* badaDTest = new TFile((outputDIR+"/badStripsaDTest.root").c_str(),"RECREATE");
  TFile* badCombinedTest = new TFile((outputDIR+"/badStripsCombined.root").c_str(),"RECREATE");
  TFile* badJBNotKSTest = new TFile((outputDIR+"/badStripsjBNotKS.root").c_str(),"RECREATE");
  TFile* badaDNotKSandjBTest = new TFile((outputDIR+"/badStripaDNotKSNotjB.root").c_str(),"RECREATE");
  TFile* badChi2NotKSandjBandaDTest = new TFile((outputDIR+"/badStripsChi2NotKsandjBandaD.root").c_str(),"RECREATE");

  long int nbadNFilledBins = 0;
  long int nbadKsTest = 0;
  long int nbadjBTest = 0;
  long int nbadaDTest = 0;
  long int nbadChi2Test = 0;
  long int nbadCombinedTest = 0;
//.........这里部分代码省略.........

void perform_smart_fit(TGraphErrors * gabscor, TF1 * fabscor) {

  int maxFitIter = 50;
  int fitIter = 0;
  vector<double> bestPars;
  double bestRChi2 = 0;
  do {
    //
    // do the fit, get the results and the parameters of the fitted function
    //

    TFitResultPtr fitResPtr = gabscor->Fit(fabscor,"RQ0S");
    vector<double> auxPars = fitResPtr.Get()->Parameters();

    //
    // compute the reduced chi2 of this fit and if it is the best fit so far
    // then save the parameters
    //
    double rchi2 = fitResPtr.Get()->Chi2()/ fitResPtr.Get()->Ndf();
    if (fitResPtr.Get()->Ndf() == 0) rchi2 = 0;
    if (rchi2 > 0 && (rchi2<bestRChi2 || bestRChi2==0)){
      bestRChi2 = rchi2;
      bestPars  = auxPars;
    }
    
    //
    // increment the counter
    //
    fitIter++;
  }while(( bestRChi2 > 2 || bestRChi2 == 0 ) && fitIter < maxFitIter);
  
  //
  // set the best parameters and chi2 to the fit function
  //

  TF1 * ffh = gabscor->GetFunction("fit");
  for (unsigned int np=0;np < bestPars.size() ; np++){
    ffh->SetParameter(np,bestPars[np]);
    fabscor->SetParameter(np,bestPars[np]);
  }
  fabscor->SetChisquare(bestRChi2 * fabscor->GetNDF());
  ffh->SetChisquare(bestRChi2 * fabscor->GetNDF());

  //
  // warn if the fit diverges at low pt
  //
  //if (fabscor->Integral(0,10) > 25)
  //cout << "\t***ERROR***, fit for histo " << gabscor->GetName() << " diverges at low pt (<10 GeV/c)" << endl;
  
  //
  // check for failed fits
  // a chi2 of zero is symptomatic of a failed fit.
  //
  
if (bestRChi2 < 0.001){
    cout<<"\t***ERROR***, FIT HAS FAILED for histo "<<gabscor->GetName()
        <<" which has a reduced chi2="<<bestRChi2
        <<" after "<<fitIter<<" iterations. "<<endl;
  }

  //
  // check for large reduced chi2's
  // above 10 is a plain error; between 5 and 10 is a warning
  //

 if (bestRChi2 > 5){
    if (bestRChi2 > 10)
      cout<<"\t***ERROR***,";
    else
      cout<<"\tWARNING,";

    cout<<" fit for histo "<<gabscor->GetName()
        <<" has a reduced chi2="<<bestRChi2
        <<" after "<<fitIter<<" iterations"<<endl;
  }
}

//.........这里部分代码省略.........
  outputTree->Branch("noiseIntegral3SigmaFromFit",&noiseIntegral3SigmaFromFit_,"noiseIntegral3SigmaFromFit/F");
  outputTree->Branch("noiseIntegral4Sigma",&noiseIntegral4Sigma_,"noiseIntegral4Sigma/F");
  outputTree->Branch("noiseIntegral4SigmaFromFit",&noiseIntegral4SigmaFromFit_,"noiseIntegral4SigmaFromFit/F");
  outputTree->Branch("noiseIntegral5Sigma",&noiseIntegral4Sigma_,"noiseIntegral5Sigma/F");
  outputTree->Branch("noiseIntegral5SigmaFromFit",&noiseIntegral4SigmaFromFit_,"noiseIntegral5SigmaFromFit/F");
  outputTree->Branch("kSValue",&kSValue_,"kSValue/F");
  outputTree->Branch("jBValue",&jBValue_,"jBValue/F");
  outputTree->Branch("aDValue",&aDValue_,"aDValue/F");
  outputTree->Branch("kSProbab",&kSProbab_,"kSProbab/F");
  outputTree->Branch("jBProbab",&jBProbab_,"jBProbab/F");
  outputTree->Branch("aDProbab",&aDProbab_,"aDProbab/F");
  outputTree->Branch("xMin",&xMin_,"xMin/F");
  outputTree->Branch("xMax",&xMax_,"xMax/F");
  outputTree->Branch("nBin",&nBin_,"nBin/F");

  bool histoBranches = false;

  // Loop on the file list to extract each histogram 2D DQM histo with full noise distribution  
  TH1F* histoNoiseStrip = NULL;
  TF1*  fitFunc = NULL;
  TH1F* randomHisto = NULL;
  TFitResultPtr result;
  for(auto file : fileList){
    cout<<"input file: "<<file<<endl;
    TFile* inputFile = TFile::Open(file.c_str(),"READ");
    inputFile->cd();
    // take into account that the DQM file structure for strips is always the same --> use cabling map to browse the histograms
    reader.SetEntry(0);
    TH2* histoNoise = NULL;
    long int iChannel = 0;
    int noFitResult = 0;
    while(reader.Next()){
      cout.flush();
      if(iChannel %10 == 0) cout<<"\r"<<"iChannel "<<100*double(iChannel)/(readoutMap->GetEntries()/reductionFactor)<<" % ";
      if(iChannel > double(readoutMap->GetEntries())/reductionFactor) break;
      iChannel++;
      TString objName;
      uint32_t fedKey =  SiStripFedKey(*fedId,SiStripFedKey::feUnit(*fedCh),SiStripFedKey::feChan(*fedCh)).key();
      std::stringstream stream;
      stream << std::hex << fedKey;
      string fedKeyStr = stream.str();
      if(fedKeyStr.size() == 4)
	objName = Form("DQMData/SiStrip/ControlView/FecCrate%d/FecSlot%d/FecRing%d/CcuAddr%d/CcuChan%d/ExpertHisto_PedsFullNoise_FedKey0x0000%s_LldChannel%d_Noise2D",*fecCrate,*fecSlot,*fecRing,*ccuAdd,*ccuChan,fedKeyStr.c_str(),*lldChannel);      
      else if(fedKeyStr.size() == 5)
	objName = Form("DQMData/SiStrip/ControlView/FecCrate%d/FecSlot%d/FecRing%d/CcuAddr%d/CcuChan%d/ExpertHisto_PedsFullNoise_FedKey0x000%s_LldChannel%d_Noise2D",*fecCrate,*fecSlot,*fecRing,*ccuAdd,*ccuChan,fedKeyStr.c_str(),*lldChannel);      
      else
	cerr<<"hex number to short "<<fedKeyStr<<" --> please check "<<endl;

      inputFile->GetObject(objName.Data(),histoNoise);
      // extract single strip noise histogram --> loop on the y-axis
      uint16_t apvID = 0;
      uint16_t stripID = 0;       
      if(histoNoiseStrip == 0 or histoNoiseStrip == NULL){
	histoNoiseStrip = new TH1F ("histoNoiseStrip","",histoNoise->GetNbinsX(),histoNoise->GetXaxis()->GetXmin(),histoNoise->GetXaxis()->GetXmax());
	histoNoiseStrip->Sumw2();
      }
      for(int iBinY = 0; iBinY < histoNoise->GetNbinsY(); iBinY++){
	histoNoiseStrip->Reset();
	histoNoiseStrip->SetDirectory(0);
	// two multiplexed APV per line
	if(iBinY < histoNoise->GetNbinsY()/2) apvID = 1;
	else apvID = 2;
	// strip id
	stripID++;
	if(stripID > 128) stripID = 1;
	// loop on x-axis bin