C++ TVectorD::GetNoElements方法代码示例

Parabola2D::Parabola2D(TVectorD x, TVectorD y, TVectorD z) : Parabola() {

  // Prepare the members
  int n = x.GetNoElements();
  double* vx = x.GetMatrixArray();
  double* vy = y.GetMatrixArray();
  double* vz = z.GetMatrixArray();
  vector<double> entry_x;
  vector<double> entry_y;
  vector<double> entry_z;
  for(int i=0; i<n; i++){
    entry_x.push_back(vx[i]);
    entry_y.push_back(vy[i]);
    entry_z.push_back(vz[i]);
  }
  entries.push_back(entry_z); // _Has_ to be the first.
  entries.push_back(entry_x);
  entries.push_back(entry_y);

  // And do the magic
  FillMatrix();    // Do the fitting...
  FindMinimum();   // Find the Minimum...
  
  // Write out the information. The order that the coefficients appear has to be in increasing key,
  // as explained above. THIS IS IMPORTANT.
  FitFunction = new TF2("FitFnc","[5]*y*y+[4]*y*x+[3]*y+[2]*x*x+[1]*x+[0]");
  for (int i = 0; i<6; i++) FitFunction->SetParameter(i,fit_coef_vec[i]);
  FitGraph = new TGraph2D(n,vx,vy,vz);

}

void dumpElements(TVectorD& a)
{
  cout << endl << endl;
  const int nrows = a.GetNoElements();
  cout << "------------------------------------------------------------------"
       << endl;
  for(int i=0; i< nrows; i++)
    {
      cout << a(i) << ", \t";    
      cout << endl;
    }
  cout << "------------------------------------------------------------------"
       << endl;
  cout << endl << endl;

  return;
}

void fit_fy_vs_y_from_all_angles(TString tgt_pos){

  // histo parameters
  gStyle->SetTitleYOffset(1.2);
  gStyle->SetTitleXOffset(1.0);
  gStyle->SetLabelSize(0.06,"x");
  gStyle->SetLabelSize(0.06,"y");
  gStyle->SetTitleSize(0.06,"x");
  gStyle->SetTitleSize(0.06,"y");
  gStyle->SetTitleX(0.4);
  gStyle->SetTitleW(0.6);
  gStyle->SetTitleBorderSize(0);
  gStyle->SetTitleFillColor(0);
  gStyle->SetTitleFontSize(0.09);
  gStyle->SetStripDecimals(0);

 
  // Fit and stat parameters
  gStyle->SetOptFit(1111);
  gStyle->SetOptStat(0);
  gStyle->SetStatFont(132);
  gStyle->SetStatY(0.99);
  gStyle->SetStatX(0.99);
  gStyle->SetStatW(0.20);
  gStyle->SetStatH(0.20);

  //Pad parameters
  gStyle->SetPadColor(0); 
  gStyle->SetPadBorderMode(0);
  gStyle->SetFrameBorderMode(0);
  gStyle->SetPadBorderSize(0);
  gStyle->SetPadTopMargin(0.18);
  gStyle->SetPadBottomMargin(0.16);
  gStyle->SetPadRightMargin(0.05);
  gStyle->SetPadLeftMargin(0.15);
  
  //set target flag
  if(tgt_pos.Contains("us")) 
    upstream = true;
  else
    upstream = false;

  // Open a file to output the fy values
  fp = fopen(Form("%s_al_fy_for_all_angles.txt",tgt_pos.Data()), "w+");
  // fprintf(fp, This is testing for fprintf...\n");
  // fputs("This is testing for fputs...\n", fp);
  
  //Delete all the objects stored in the current directory memmory
  gDirectory->Delete("*");
  
  std::vector<std::string> files;
  Double_t angles[6] = {18,22,26,32,40,50};

  // vectors to store the fy and y from individual angles and their errors
  TVectorD *fyv = new TVectorD();
  TVectorD *yv= new TVectorD();
  TVectorD *yev= new TVectorD();
  TVectorD *fyev= new TVectorD();
  
  // vectors to store the fy_data/fy_model ratio and its error
  TVectorD *ratio = new TVectorD();
  TVectorD *eratio = new TVectorD();

  // vectors to store all the data together
  TVectorD *tfyv = new TVectorD();
  TVectorD *tyv= new TVectorD();
  TVectorD *tyev= new TVectorD();
  TVectorD *tfyev= new TVectorD();
  

 // Open canvas to draw the comparision between Fy vs Y using sig_exp and sig_exp corrected for sig_di// s
  TCanvas * c = new TCanvas("c", "",0,0,1000,1000);  
  c->Draw();
  c->Divide(1,2);
  
  Int_t pp = 0;

  // select files for a target with the given angle and get fy and y
  for(int i=0;i<6;i++){

    std::cout<<"reading angle = "<<angles[i]<<std::endl;
    files =  get_files(angles[i], tgt_pos);
    if(nodata){
      //skip. Do nothing
      nodata = false;
    }
    else{

      // get the fy and y from data with their errors
      get_fy_y(files, fyv, yv, yev, fyev);
      for(int k =0;k<fyv->GetNoElements();k++){


	// Fill the full arrays
	(*tfyv).ResizeTo(pp+1);
	(*tfyev).ResizeTo(pp+1);
	(*tyv).ResizeTo(pp+1);
	(*tyev).ResizeTo(pp+1);

	(*tfyv).operator()(pp) = (*fyv)[k];
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    eventTree = (TTree*)infile->Get("Events"); assert(eventTree);

    // Find weight for events for this file
    // The first file in the list comes with weight 1,
    // all subsequent ones are normalized to xsection and luminosity
    double xsec=xsecv[ifile];
    AdjustXSectionForSkim(infile,xsec,eventTree->GetEntries(),1);
    lumiv[ifile] = eventTree->GetEntries()/xsec;
    double scale = lumiv[0]/lumiv[ifile];
    cout << "       -> sample weight is " << scale << endl;

    // Set branch address to structures that will store the info  
    eventTree->SetBranchAddress("Info",&info);                TBranch *infoBr       = eventTree->GetBranch("Info");
    eventTree->SetBranchAddress("Gen",&gen);                  TBranch *genBr = eventTree->GetBranch("Gen");
    eventTree->SetBranchAddress("Dielectron",&dielectronArr); TBranch *dielectronBr = eventTree->GetBranch("Dielectron");
  
    // loop over events    
    for(UInt_t ientry=0; ientry<eventTree->GetEntries(); ientry++) {
      //if (ientry>100000) break;

      genBr->GetEntry(ientry);
      infoBr->GetEntry(ientry);

      double reweight;
      if (systematicsMode!=DYTools::FSR_STUDY) reweight=1.0;
      else if (((gen->mass)-(gen->vmass))>massLimit) reweight=1.0;
      else reweight=reweightFsr;

      if (ientry<20) std::cout<<"reweight="<<reweight<<std::endl;


      // Use post-FSR generator level mass in unfolding
      int ibinFsr = DYTools::findMassBin(gen->mass);
      if(ibinFsr != -1 && ibinFsr < yieldsMcFsr.GetNoElements()){
         yieldsMcFsr[ibinFsr] += reweight * scale * gen->weight;
      }


      /*
      // For EPS2011 for both data and MC (starting from Summer11 production)
      // we use an OR of the twi triggers below. Both are unpresecaled.
      ULong_t eventTriggerBit = kHLT_Ele17_CaloIdL_CaloIsoVL_Ele8_CaloIdL_CaloIsoVL 
	| kHLT_Ele17_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele8_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL;
      ULong_t leadingTriggerObjectBit = kHLT_Ele17_CaloIdL_CaloIsoVL_Ele8_CaloIdL_CaloIsoVL_Ele1Obj
	| kHLT_Ele17_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele8_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele1Obj;
      ULong_t trailingTriggerObjectBit = kHLT_Ele17_CaloIdL_CaloIsoVL_Ele8_CaloIdL_CaloIsoVL_Ele2Obj
	| kHLT_Ele17_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele8_CaloIdT_TrkIdVL_CaloIsoVL_TrkIsoVL_Ele2Obj;
      */

      const bool isData=kFALSE;
      TriggerConstantSet constantsSet = Full2011DatasetTriggers; // Enum from TriggerSelection.hh
      TriggerSelection requiredTriggers(constantsSet, isData, info->runNum);
      ULong_t eventTriggerBit = requiredTriggers.getEventTriggerBit(info->runNum);
      ULong_t leadingTriggerObjectBit = requiredTriggers.getLeadingTriggerObjectBit(info->runNum);
      ULong_t trailingTriggerObjectBit = requiredTriggers.getTrailingTriggerObjectBit(info->runNum);

      
      if(!(info->triggerBits & eventTriggerBit)) continue;  // no trigger accept? Skip to next event...                                   

      // loop through dielectrons
      dielectronArr->Clear();
      dielectronBr->GetEntry(ientry);    
      for(Int_t i=0; i<dielectronArr->GetEntriesFast(); i++) {

        const mithep::TDielectron *dielectron = (mithep::TDielectron*)((*dielectronArr)[i]);
	

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

      // The A_FSR starting point: gen level quantities
      // The FSR subscript means we work with post-FSR generator level quantities
      double et1 = gen->pt_1;
      double et2 = gen->pt_2;
      double eta1 = gen->eta_1;
      double eta2 = gen->eta_2;

      // Apply acceptance requirement
      if( ! ( ( et1 > DYTools::etMinLead  && et2 > DYTools::etMinTrail)
	      || ( et1 > DYTools::etMinTrail  && et2 > DYTools::etMinLead) )) continue;
      if((fabs(eta1) >= 2.5)      || (fabs(eta2) >= 2.5))       continue;
      if((fabs(eta1) >= kGAP_LOW) && (fabs(eta1) <= kGAP_HIGH)) continue;
      if((fabs(eta2) >= kGAP_LOW) && (fabs(eta2) <= kGAP_HIGH)) continue;

      // These events are in acceptance, use them for efficiency denominator
      Bool_t isBGen1 = (fabs(eta1)<kGAP_LOW);
      Bool_t isBGen2 = (fabs(eta2)<kGAP_LOW);         
      // determine number of good vertices
      pvBr->GetEntry(ientry);
      int iPUBin=-1;
      int nGoodPV=-1;
      if ((gen->mass>=60) && (gen->mass<=120)) {
	nGoodPV=0;
	for (Int_t ipv=0; ipv<pvArr->GetEntriesFast(); ipv++) {
	  const mithep::TVertex *pv = (mithep::TVertex*)((*pvArr)[ipv]);
	  if(pv->nTracksFit                        < 1)  continue;
	  if(pv->ndof                              < 4)  continue;
	  if(fabs(pv->z)                           > 24) continue;
	  if(sqrt((pv->x)*(pv->x)+(pv->y)*(pv->y)) > 2)  continue;
	  nGoodPV++;
	}
	if (nGoodPV>0) {
           if (nGoodPV<=nEventsZPeakPURaw.GetNoElements()) 
	        nEventsZPeakPURaw[nGoodPV] += scale * gen->weight;
	iPUBin=DYTools::findMassBin(double(nGoodPV),DYTools::nPVBinCount,DYTools::nPVLimits);
	//std::cout << "iPUBin=" << iPUBin << ", nGoodPV=" << nGoodPV << "\n";
	if ((iPUBin!=-1) && (iPUBin < nEventsZPeakPU.GetNoElements())) {
	  nEventsZPeakPU[iPUBin] += scale * gen->weight;
	}
	else {
	  std::cout << "error in PU bin indexing iPUBin=" << iPUBin << ", nGoodPV=" << nGoodPV << "\n";
	}
      } 
	//else std::cout << "nGoodPV=" << nGoodPV << "\n";
      }

      // Use post-FSR generator level mass for binning
      int ibinGen = DYTools::findMassBin(gen->mass);
      // Accumulate denominator for efficiency calculations
      if(ibinGen != -1 && ibinGen < nEventsv.GetNoElements()){
	nEventsv[ibinGen] += scale * gen->weight;
	// Split events barrel/endcap using matched supercluster or particle eta
	if(isBGen1 && isBGen2)                                  { nEventsBBv[ibinGen] += scale * gen->weight; } 
	else if(!isBGen1 && !isBGen2)                           { nEventsEEv[ibinGen] += scale * gen->weight; } 
	else if((isBGen1 && !isBGen2) || (!isBGen1 && isBGen2)) { nEventsBEv[ibinGen] += scale * gen->weight; }
      }else if(ibinGen >= nEventsv.GetNoElements())
	cout << "ERROR: binning problem" << endl;

      if(!(info->triggerBits & eventTriggerBit)) continue;  // no trigger accept? Skip to next event...                                   

      // loop through dielectrons
      dielectronArr->Clear();
      dielectronBr->GetEntry(ientry);    
      for(Int_t i=0; i<dielectronArr->GetEntriesFast(); i++) {
        const mithep::TDielectron *dielectron = (mithep::TDielectron*)((*dielectronArr)[i]);

//.........这里部分代码省略.........
                i2++;
            }
        }
        else if(i1 < length1 && i2 >= length2)
        {
            binLimits.push_back(binLimits1[i1]);
            std::vector<int> tmp;
            tmp.push_back(i1);
            tmp.push_back(-1);
            preU.push_back(tmp);
            i1++;
        }
        else
        {
            binLimits.push_back(binLimits2[i2]);
            std::vector<int> tmp;
            tmp.push_back(-1);
            tmp.push_back(i2);
            preU.push_back(tmp);
            i2++;
        }
    }

    TVectorD dv(length1 + length2);
    for(int i = 0; i < length1 + length2; i++)
    {
        dv[i] = (i < length1) ? dv1[i] : dv2[i - length1];
    }

    TMatrixT<double> cm(length1 + length2, length1 + length2), U(length1 + length2, preU.size());
    for(int i = 0; i < length1; i++)
    {
        for(int j = 0; j < length1; j++)
        {
            cm[i][j] = cm1[i][j];
        }
    }
    for(int i = length1; i < length1 + length2; i++)
    {
        for(int j = length1; j < length1 + length2; j++)
        {
            cm[i][j] = cm2[i - length1][j - length1];
        }
    }

    for(unsigned int i = 0; i < preU.size(); i++)
    {
        if(preU[i][0] >= 0) U[preU[i][0]][i] = 1;
        if(preU[i][1] >= 0) U[preU[i][1] + length1][i] = 1;
        if(ftr > 1 && preU[i][0] >= 0 && preU[i][1] >= 0)  cm[preU[i][0]][preU[i][1] + length1] = cm[preU[i][1] + length1][preU[i][0]] = ccCov1[preU[i][0]]*ccCov2[preU[i][1]];
    }
    
    //    cm.Print();

    TMatrixT<double> Ut(U);
    Ut.T();

    TMatrixT<double> cmInv(cm);
    cmInv.Invert();
    TMatrixT<double> step1 = Ut * cmInv * U;
    TMatrixT<double> step2 = Ut * cmInv;
    TMatrixT<double> lambda = step1.Invert() * step2;
    TVectorD bV = lambda*dv;
    TMatrixT<double> bcm = (Ut * cmInv * U).Invert();

    printf("Done with combination.\n");

    //write output
    FILE *file;
    char bVoutName[128], CMoutName[128];
    sprintf(bVoutName, "%s_data.txt", outputNameStub);

    file = fopen(bVoutName, "w");
    if(file)
    {
        fprintf(file, "#\n#%9s %9s %9s %15s %15s\n", "Bin", "Y_min", "Y_max", "Value", "Uncertainty");
        for(int i = 0; i < bV.GetNoElements(); i++)
        {
            fprintf(file, " %9i %9.2f %9.2f %15e %15e\n", i + 1, binLimits[i][0], binLimits[i][1], bV[i], sqrt(bcm[i][i]));
        }
        fclose(file);
    }

    sprintf(CMoutName, "%s_covMat.txt", outputNameStub);

    file = fopen(CMoutName, "w");
    if(file)
    {
        fprintf(file, "#\n#%9s %9s %15s\n", "Bin i", "Bin j", "Value");
        for(int i = 0; i < bcm.GetNrows(); i++)
        {
            for(int j = 0; j < bcm.GetNcols(); j++)
            {
                fprintf(file, " %9i %9i %15e\n", i + 1, j + 1, bcm[i][j]);
            }
        }
        fclose(file);
    }
    printf("Output complete.\n");
}