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

	NBD(TString datastr, TString Glauberstr,TString histonamev){
		dataname.SetString(datastr.Data());
		Glaubername.SetString(Glauberstr.Data());
		histoname.SetString(histonamev.Data());
	        mubest.ResizeTo(1);  kbest.ResizeTo(1);
	        chis.ResizeTo(1);  Ndf.ResizeTo(1);
	};

// ------------ method called once each job just before starting event loop  ------------
void
LYZ::beginJob()
{
    double Vmax[nbin], eps[nbin];
    for(int ibin=0; ibin<nbin ; ibin++) {
        Vmax[ibin]=0.065*(trkbin[ibin]+30);
        eps[ibin]=0.00025*(trkbin[ibin]+30);
    }

    Nevent.ResizeTo(nbin);
    totmultall.ResizeTo(nbin), tottrk.ResizeTo(nbin);;
    Nevent.Zero();
    totmultall.Zero(),	tottrk.Zero();
    for(int ibin=0; ibin<nbin; ibin++) {
        for(int iptbin=0; iptbin<nptV; iptbin++) {
            Qx2[ibin][iptbin]=0;
            Qy2[ibin][iptbin]=0;
        }
    }

    for(int itheta=0; itheta<ntheta; itheta++) {
        theta[itheta]=itheta*TMath::Pi()/ntheta/nn;
        for(int ibin=0; ibin<nbin; ibin++) {
            for(int iptbin=0; iptbin<nptV; iptbin++) {
                GRe[ibin][iptbin][itheta].ResizeTo(nstepr);
                GRe[ibin][iptbin][itheta].Zero();
                GIm[ibin][iptbin][itheta].ResizeTo(nstepr);
                GIm[ibin][iptbin][itheta].Zero();
            }
        }
    }

    for(int ibin=0; ibin<nbin; ibin++) {
        for(int ir=0; ir<nstepr; ir++)
            if(isSimple==0)	r[ibin][ir]=j01/(Vmax[ibin]-eps[ibin]*ir);
            else	r[ibin][ir]=0.00025*20*(ir+1);
        Qx1[ibin].ResizeTo(nptV);
        Qy1[ibin].ResizeTo(nptV);
        Q2[ibin].ResizeTo(nptV);
        Qx1[ibin].Zero();
        Qy1[ibin].Zero();
        Q2[ibin].Zero();
        totmult[ibin].ResizeTo(nptv);
        totmult[ibin].Zero();
        totpt[ibin].ResizeTo(nptv);
        totpt[ibin].Zero();
        dDRe[ibin].ResizeTo(ntheta);
        dDRe[ibin].Zero();
        dDIm[ibin].ResizeTo(ntheta);
        dDIm[ibin].Zero();
        for(int itheta=0; itheta<ntheta; itheta++) {
            dNRe[ibin][itheta].ResizeTo(nptv);
            dNRe[ibin][itheta].Zero();
            dNIm[ibin][itheta].ResizeTo(nptv);
            dNIm[ibin][itheta].Zero();

        }
    }
}

void merge(){
	TVectorD Nevent;	Nevent.ResizeTo(nbin);  Nevent.Zero();
     //   TVectorD totmultall;	totmultall.ResizeTo(nbin);      totmultall.Zero();
     //   TVectorD avgmultall;	avgmultall.ResizeTo(nbin);      avgmultall.Zero();
        TVectorD tottrk;	tottrk.ResizeTo(nbin);      tottrk.Zero();
//	TVectorD totptall;      totptall.ResizeTo(nbin);    totptall.Zero();
//	TVectorD totetaall;      totetaall.ResizeTo(nbin);    totetaall.Zero();
        TVectorD avgtrk;	avgtrk.ResizeTo(nbin);      avgtrk.Zero();
//	TVectorD avgmult;       avgmult.ResizeTo(nbin);
            
        TH2F* s[nbin];
        TH2F* b[nbin];
        for(int ibin=0;ibin<nbin;ibin++){
            s[ibin] = new TH2F(Form("s_%d",ibin),Form("signal",ibin),detastep,detamin,detamax,dphistep,dphimin,dphimax);
            s[ibin]->Sumw2();
            b[ibin]  = new TH2F(Form("b_%d",ibin), "background",detastep,detamin,detamax,dphistep,dphimin,dphimax);
            b[ibin]->Sumw2();
        }
        TFile *fout = new TFile(Form("Anav3_merged.root"),"Recreate");
        TFile *f[nFileAll];
        for(int ifile=; ifile<63; ifile++){
                f[ifile] = TFile::Open(Form("%s/Anav3_%d.root",outdir.Data(),ifile));
                if(!f[ifile]) continue;
		TVectorD* Nevent_t =  (TVectorD*)f[ifile]->Get(Form("Nevent"));
	//	TVectorD* totmultall_t =  (TVectorD*)f[ifile]->Get(Form("totmultall"));
		TVectorD* tottrk_t =  (TVectorD*)f[ifile]->Get(Form("tottrk"));
	//	TVectorD* totptall_t =  (TVectorD*)f[ifile]->Get(Form("totptall"));
//		TVectorD* totetaall_t =  (TVectorD*)f[ifile]->Get(Form("totetaall"));
		for(int ibin=0;ibin<nbin;ibin++){
			//	totptall[ibin] += (*totptall_t)[ibin];
			//	totetaall[ibin] += (*totetaall_t)[ibin];
			        Nevent[ibin] += (*Nevent_t)[ibin];
			  //      totmultall[ibin] += (*totmultall_t)[ibin];	
			        tottrk[ibin] += (*tottrk_t)[ibin];	
                                TH1D* s_t = (TH1D*)f[ifile]->Get(Form("signal_%d",ibin));
                                TH1D* b_t = (TH1D*)f[ifile]->Get(Form("background_%d",ibin));
                                s[ibin]->Add(s_t);
                                b[ibin]->Add(b_t);
		}
		f[ifile]->Close();
        }
    for(int ibin=0;ibin<nbin;ibin++){
    avgtrk[ibin] = tottrk[ibin]/Nevent[ibin];
    TH1D* hr = (TH1D*)s[ibin]->Clone("hr");
    hr->Divide(b[ibin]);
    hr->Scale(b[ibin]->GetBinContent(b[ibin]->FindBin(0,0)));

    fout->cd();
    Nevent.Write("Nevent");
    avgtrk.Write("avgtrk");
    TDirectory *dir0 = (TDirectory*)fout->mkdir(Form("D_%d",ibin));
    dir0->cd();
    s[ibin]->Write("s");
    b[ibin]->Write("b");
    hr->Write();
    }
}

void loadVectorsFromFile(const char *filename, 
			 const char *scanfmt,
			 TVectorD&   v1,
			 TVectorD&   v2,
			 TVectorD&   v3)
{
  char linein[LINELEN];
  vector<double> v;

  if (TString(scanfmt).Contains("%f")) {
    cerr << "Must use %lf format for doubles, sorry " << endl;
    exit(-1);
  }

  FILE *fp = fopen(filename, "r");

  if (!fp) {
    cerr << "File failed to open, " << filename << endl;
    exit(-1);
  }

  if (gl_verbose)
    cout << "Loading vectors from file " << filename
	 << " with scan format \"" << scanfmt << "\"" << endl;

  while (!feof(fp) && fgets(linein,LINELEN,fp)) {
    double x1, x2, x3;
    if( linein[0]=='#' ) continue;                // comments are welcome
    if( sscanf(linein, scanfmt, &x1, &x2, &x3) != 3 ) {
      cerr << "scan failed, file " << filename << ", line = " << linein;
      cerr << ", scanfmt = " << scanfmt;
      cerr << ", skipping" << endl;
      //return;
      continue;
    }
    else {
      v.push_back(x1); v.push_back(x2); v.push_back(x3);
    }
  }

  int vecsize = v.size()/3;
  v1.ResizeTo(vecsize);
  v2.ResizeTo(vecsize);
  v3.ResizeTo(vecsize);

  if (gl_verbose) cout << "; read " << vecsize << " lines" << endl;

  for (int i=0; i<vecsize; i++) {
    v1[i] = v[3*i];
    v2[i] = v[3*i+1];
    v3[i] = v[3*i+2];
  }

  fclose(fp);
}                                                 // loadVectorsFromFile

void loadVectorsFromFile(const char *filename, 
			 const char *scanfmt,
			 TVectorD&   vx,
			 TVectorD&   vy,
			 char        xheader[],
			 char        yheader[])
{
  char linein[LINELEN];
  vector<double> v;

  FILE *fp = fopen(filename, "r");

  if (!fp) {
    cerr << "File failed to open, " << filename << endl;
    return;
  }

  if (gl_verbose)
    cout << "Loading vectors from file " << filename;

  while (!feof(fp) && fgets(linein,LINELEN,fp)) {
    double x, y;
    if( linein[0]=='#' ) {
      if (!v.size()) {           // first line, try to read headers
	TString hscanfmt(scanfmt);
	hscanfmt.ReplaceAll("lf","s");
	if( sscanf(&linein[1],hscanfmt.Data(),xheader,yheader) != 2 ) {
	  cerr << "failed to read in column headers" << endl;
	}
      }
      continue;                // comments are welcome
    }

    if( sscanf(linein, scanfmt, &x, &y) != 2 ) {
      cerr << "scan failed, file " << filename << ", line = " << linein << endl;
      return;
    }
    else {
      v.push_back(x); v.push_back(y);
    }
  }

  int vecsize = v.size()/2;
  vx.ResizeTo(vecsize);
  vy.ResizeTo(vecsize);

  if (gl_verbose) cout << "; read " << vecsize << " lines" << endl;

  for (int i=0; i<vecsize; i++) {
    vx[i] = v[2*i];
    vy[i] = v[2*i+1];
  }
}                                                 // loadVectorsFromFile

void FillVar(){
  ///

  npoints =tree->Draw("side",cutA+cut1+cut2);
  vside.ResizeTo(npoints);  vside.SetElements(tree->GetV1());
  //
  npoints =tree->Draw("z",cutA+cut1+cut2);
  vz.ResizeTo(npoints);  vz.SetElements(tree->GetV1());
  //
  npoints =tree->Draw("dr",cutA+cut1+cut2);
  vdr.ResizeTo(npoints);  vdr.SetElements(tree->GetV1());
  //
  npoints =tree->Draw("tl",cutA+cut1+cut2);
  vtl.ResizeTo(npoints);  vtl.SetElements(tree->GetV1());
  //
  npoints =tree->Draw("sa",cutA+cut1+cut2);
  vsa.ResizeTo(npoints);  vsa.SetElements(tree->GetV1());
  //
  npoints =tree->Draw("ca",cutA+cut1+cut2);
  vca.ResizeTo(npoints);  vca.SetElements(tree->GetV1());
  //
  npoints =tree->Draw("dz",cutA+cut1+cut2);
  vdz.ResizeTo(npoints);  vdz.SetElements(tree->GetV1());
  //
  npoints =tree->Draw("dtl",cutA+cut1+cut2);
  vdtl.ResizeTo(npoints);  vdtl.SetElements(tree->GetV1());
  //
  npoints =tree->Draw("ez",cutA+cut1+cut2);
  vez.ResizeTo(npoints);  vez.SetElements(tree->GetV1());
  //
  npoints =tree->Draw("etl*0.5",cutA+cut1+cut2);
  vetl.ResizeTo(npoints);  vetl.SetElements(tree->GetV1()); 
  //
  npoints = tree->Draw("run","run>0");
  Int_t *np=new Int_t[npoints]; for (Int_t i=0;i<npoints;i++) np[i]=tree->GetV1()[i];
  Int_t *nsort=new Int_t[npoints];
  Int_t nruns = toolkit.Freq(npoints,np,nsort,kTRUE); 
}

void DumpRecenterParv2(){
    TString pro = "Pro104";
    int taxi = 8583;
    TFile *fin;
    int nrun = GetTotalRun();
    if(nrun<0) exit("Empty run list file!");

    TFile *fout = new TFile(Form("Recentering.root"),"Recreate");
    TVectorD vecmean;
    TVectorD vecrms;
    vecmean.ResizeTo(ncent*nbbcz*nhar*nsub*nxy);
    vecrms.ResizeTo(ncent*nbbcz*nhar*nsub*nxy);
    for(int irun=0;irun<nrun;irun++){
     cout<<irun<<" of total "<<nrun<<" runs"<<endl;
     int RunNumber=GetRun(irun);
     fin = TFile::Open(Form("Run15pAu200MinBias/output_fvtxwithcntrecenter_%d.root",RunNumber));
     //ofstream fout(Form("Calibration/%s/Recentering_%d.dat",dataset.Data(),GetRun(irun)));
     for(int icent=0;icent<ncent;icent++){
      for(int ibbcz=0;ibbcz<nbbcz;ibbcz++){
       for(int ihar=0;ihar<nhar;ihar++){
        for(int isub=0;isub<nsub;isub++){
         for(int ixy=0;ixy<nxy;ixy++){
             TH1F* q = (TH1F*)fin->Get(Form("q_%d_%d_%d_%d_%d",icent,ibbcz,ihar,isub,ixy));
             float mean = q->GetMean();
             float rms = q->GetRMS();
             vecmean[icent*nbbcz*nhar*nsub*nxy+ibbcz*nhar*nsub*nxy+ihar*nsub*nxy+isub*nxy+ixy] = mean;
             vecrms[icent*nbbcz*nhar*nsub*nxy+ibbcz*nhar*nsub*nxy+ihar*nsub*nxy+isub*nxy+ixy] = rms;
         }
        }
       }
      }
     }
    fout->cd();
    vecmean.Write(Form("mean_%d",RunNumber));
    vecrms.Write(Form("rms_%d",RunNumber));
    fin->Close();
    }
}

        void initN(double *binv, int Nv, double methodv){
                centbin.ResizeTo(Nv);   kpoint.ResizeTo(Nv);
                centbin_.ResizeTo(Nv);  kpoint_.ResizeTo(Nv);
		method.ResizeTo(1);
                if(methodv==0){
                        centbin=TVectorD(Nv,binv);
                        centbin_=TVectorD(Nv,binv);
                }
                else{
                        kpoint=TVectorD(Nv,binv);
                        kpoint_=TVectorD(Nv,binv);
                }
                N=Nv-1;
                method[0]=methodv;
                NpartAver.ResizeTo(N); NcollAver.ResizeTo(N);BAver.ResizeTo(N);
                Npartdis = new TObjArray(0);
        };

void MultiGaus(const TVectorD& parMeans, const TMatrixDSym& covMatrix, TVectorD& genPars)
{

  TRandom3 rnd(0);

  int nPars = parMeans.GetNrows();
  if(nPars <= 0) {
    Error("MultiGaus", "Must have >0 pars");
    return;
  }
  if(covMatrix.GetNrows() != nPars) {
    Error("MultiGaus", "parMeans.GetNrows() != covMatrix.GetNrows()");
    return;
  }
 
  // Check that covMatrix is symmetric
  for(int iRow = 0; iRow < nPars; iRow++) {
    for(int iCol = iRow; iCol < nPars; iCol++) {
      if(covMatrix(iRow, iCol) != covMatrix(iCol, iRow)) {
        Error("MultiGaus", "malformed cov matrix at row %d, col %d", iRow, iCol);
        return;
      }
    }
  }

  genPars.ResizeTo(nPars);

  TMatrixDSymEigen eigenvariances(covMatrix);
  
  TMatrixD V = eigenvariances.GetEigenVectors();

  TVectorD rotParMeans = V * parMeans;

  for(int iPar = 0; iPar < nPars; iPar++) {
    double variance = eigenvariances.GetEigenValues()[iPar];
    // check for positive-definiteness of covMatrix
    if(variance < 0) {
      Error("MultiGaus", "Got a negative eigenvariance (%f) on iPar = %d", variance, iPar);
    }
    genPars[iPar] = rnd.Gaus(rotParMeans[iPar], sqrt(variance));
  }

  V.Invert();
  
  genPars = V * genPars;

}

// ------------ method called once each job just before starting event loop  ------------
void 
LYZ::beginJob(int ispt_)
{
  double Vmax[nbin], eps[nbin];
  for(int ibin=0; ibin<nbin ;ibin++){
        Vmax[ibin]=0.065*(trkbin[ibin]+30);
        eps[ibin]=0.00025*(trkbin[ibin]+30);
  }
    ispt = ispt_;
    if(ispt_){             nvv = nptv;       binv = ptbinv;}
    else{             nvv = netav;          binv = etabinv;}

	Nevent.ResizeTo(nbin);	totmultall.ResizeTo(nbin), tottrk.ResizeTo(nbin), totptall.ResizeTo(nbin), totetaall.ResizeTo(nbin); totwall.ResizeTo(nbin); 
		Qx1.ResizeTo(nbin);	Qy1.ResizeTo(nbin);	Q2.ResizeTo(nbin);
	Nevent.Zero();	totmultall.Zero(); totwall.Zero(); tottrk.Zero(); totptall.Zero(); totetaall.Zero();
		Qx1.Zero();	Qy1.Zero();	Q2.Zero();
	for(int ibin=0; ibin<nbin; ibin++){
        		Qx2[ibin]=0;  Qy2[ibin]=0;
	}

        	for(int itheta=0;itheta<ntheta;itheta++){
               		theta[itheta]=itheta*TMath::Pi()/ntheta/nn;
			for(int ibin=0; ibin<nbin; ibin++){
					GRe[ibin][itheta].ResizeTo(nstepr); 	GRe[ibin][itheta].Zero();
					GIm[ibin][itheta].ResizeTo(nstepr);	GIm[ibin][itheta].Zero();
				}
		}
	
	for(int ibin=0; ibin<nbin; ibin++){
        	    for(int ir=0; ir<nstepr; ir++)
                        if(isSimple==0)	r[ibin][ir]=j01/(Vmax[ibin]-eps[ibin]*ir);
                	else	r[ibin][ir]=0.00025*20*(ir+1);
                totmultv[ibin].ResizeTo(nvv);	totmultv[ibin].Zero();  
                totptv[ibin].ResizeTo(nvv);	totptv[ibin].Zero();  
                totetav[ibin].ResizeTo(nvv);	totetav[ibin].Zero();  
		dDRe[ibin].ResizeTo(ntheta);	dDRe[ibin].Zero();
		dDIm[ibin].ResizeTo(ntheta);	dDIm[ibin].Zero();
                for(int itheta=0;itheta<ntheta;itheta++){
                    dNRe[ibin][itheta].ResizeTo(nvv);	dNRe[ibin][itheta].Zero();
                    dNIm[ibin][itheta].ResizeTo(nvv);	dNIm[ibin][itheta].Zero();
                }
        }
}

// ------------ method called once each job just before starting event loop  ------------
void 
qfit::beginJob(int ispt_)
{
    ispt = ispt_;
    if(ispt_){             nvv = nptv;       binv = ptbinv;}
    else{             nvv = netav;          binv = etabinv;}

	Nevent.ResizeTo(nbin);	totmultall.ResizeTo(nbin), tottrk.ResizeTo(nbin), totptall.ResizeTo(nbin), totptall2.ResizeTo(nbin); totetaall.ResizeTo(nbin);
	Nevent.Zero();	totmultall.Zero(),	tottrk.Zero(); totptall.Zero(); totptall2.Zero(); totetaall.Zero();
        q22.ResizeTo(nbin);q24.ResizeTo(nbin);
        q22.Zero();q24.Zero();

        	for(int itheta=0;itheta<ntheta;itheta++){
               		theta[itheta]=itheta*TMath::Pi()/ntheta/nn;
		}
	
	for(int ibin=0; ibin<nbin; ibin++){
                totmultv[ibin].ResizeTo(nvv);	totmultv[ibin].Zero();  
                totptv[ibin].ResizeTo(nvv);	totptv[ibin].Zero();  
                totetav[ibin].ResizeTo(nvv);	totetav[ibin].Zero();  
                for(int itheta=0;itheta<ntheta;itheta++){
                    hq[ibin][itheta] = new TH1D(Form("hq_%d_%d",ibin,itheta),Form("hq_%d_%d",ibin,itheta),1000,0,10);
                    hq[ibin][itheta]->Sumw2();
                }
                hqx[ibin] = new TH1D(Form("hqx_%d",ibin),Form("hqx_%d",ibin),1000,0,10);
                hqy[ibin] = new TH1D(Form("hqy_%d",ibin),Form("hqy_%d",ibin),1000,0,10);
                hq2[ibin] = new TH1D(Form("hq2_%d",ibin),Form("hq2_%d",ibin),1000,0,10);
                hvq2[ibin] = new TH1D(Form("hvq2_%d",ibin),Form("hvq2_%d",ibin),1000,0,1);
                hq2nonf[ibin] = new TH1D(Form("hq2nonf_%d",ibin),Form("hq2nonf_%d",ibin),1000,0,10);
                hqx[ibin]->Sumw2();
                hqy[ibin]->Sumw2();
                hq2[ibin]->Sumw2();
                hvq2[ibin]->Sumw2();
                hq2nonf[ibin]->Sumw2();
        }
}

	void initx(double xminv, double xmaxv){
		xmin.ResizeTo(1); xmax.ResizeTo(1);
		xmin[0]=xminv;	xmax[0]=xmaxv;
	};

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

   Double_t threshold = 0;
   Double_t Tthreshold = 1.;

   Double_t charge, Ccharge;
   Double_t Rm,Cm,Zcor,Zm;

   //Matrix for the charge map
   Double_t **padCharge=new Double_t*[numberOfRows];
   Double_t **padTime=new Double_t*[numberOfRows];
   Double_t **padHeight=new Double_t*[numberOfRows];
   for(Int_t j=0;j<numberOfRows;j++){
     padCharge[j]=new Double_t[numberOfColumns];
     padTime[j]=new Double_t[numberOfColumns];
   }
   
   //Double_t *padChargeTest=new Double_t[numberOfRows];
   //Double_t *padTimeTest=new Double_t[numberOfRows];

   //TMatrixD *padCharge=0;
   //TMatrixD *padTime=0;
   //TMatrixD *padHeight=0;

   //TMatrixD padCharge(32,64,2000);
   //TMatrixD padTime(32,64,2000);
   //TMatrixD padHeight(32,64,2000);

   //TMatrixD *padChargeTest=0;
   //TMatrixD *padTimeTest=0;
   TMatrixD padChargeTest;
   TMatrixD padTimeTest;

   TMatrixD a(32,64);
   padChargeTest.ResizeTo(a);
   padTimeTest.ResizeTo(a);

   TMatrixD *pointerCharge = &padChargeTest;
   TMatrixD *pointerTime = &padTimeTest;

   //Silicon charge
   //Double_t SilCharge[16];
   //Int_t SilID[16];

   TVectorD SilCharge;
   SilCharge.ResizeTo(16);
   TVectorD SilID;
   SilID.ResizeTo(16);

   TVectorD *pointerSilCharge=&SilCharge;
   TVectorD *pointerSilID=&SilID;
   //Double_t *SilChargePointer=&SilCharge;

   //Branching the out tree
   //out_tree->Branch("PadCharge","TMatrixD",&padCharge,64000,0);
   //out_tree->Branch("PadTime","TMatrixD",&padTime,64000,0);
   //out_tree->Branch("SilCharge",&SilCharge,64000,0);
   out_tree->Branch("PadCharge","TMatrixD",&pointerCharge,64000,0);
   out_tree->Branch("PadTime","TMatrixD",&pointerTime,64000,0);
   out_tree->Branch("SilCharge","TVectorD",&pointerSilCharge,64000,0);
   out_tree->Branch("SilID","TVectorD",&pointerSilID,64000,0);
   //out_tree->Branch("SilCharge",&SilCharge,"energy/D");


   // out_tree->Branch("Energy",&Qtot,"energy/D");
   // out_tree->Branch("HAngle",&HAngle,"angle/D");
   // out_tree->Branch("VAngle",&VAngle,"angle/D");

//.........这里部分代码省略.........
hqy->GetXaxis()->SetTitle("q_{y}");
hqy->GetYaxis()->SetTitle("#frac{dN}{dq_{y}}");
hqy->GetYaxis()->SetRangeUser(0,1);
hqy->SetMarkerStyle(24);
hqy->SetMarkerSize(0.5);
hqy->Fit(Form("f1fit"),"R","",0,10);
hqy->Draw("Psame");
c3->cd(2)->SetLogy();
TH1D* hqy_cp = (TH1D*)hqy->Clone("hqy_cp");
fixedFontHist(hqy_cp,1.0,2.0);
divideByBinCenter(hqy_cp);
hqy_cp->GetYaxis()->SetTitle("#frac{dN}{q_{y}dq_{y}}");
hqy_cp->GetYaxis()->SetRangeUser(1e-10,10);
hqy_cp->SetMarkerStyle(24);
hqy_cp->SetMarkerSize(0.5);
hqy_cp->Fit(Form("ffit"),"R","",0,10);
hqy_cp->Draw("Psame");
t->DrawLatex(0.5,0.2,Form("N_{trk}^{offline} = %.2f",avgtrk[ibin]));
*/
TCanvas *c4 = new TCanvas("c4","c4",1000,500);
c4->Divide(2,1);
c4->cd(1)->SetLeftMargin(0.18);
fixedFontHist(hq2,1.0,2.0);
hq2->SetTitle("");
hq2->GetXaxis()->SetTitle("q_{2}");
hq2->GetYaxis()->SetTitle("#frac{dN}{dq_{2}}");
hq2->GetYaxis()->SetRangeUser(0,1);
hq2->SetMarkerStyle(24);
hq2->SetMarkerColor(4);
hq2->SetLineColor(4);
hq2->SetMarkerSize(0.5);
hq2->Fit(Form("f1fit"),"R","",0,10);
TVectorD vecr;
vecr.ResizeTo(8);
vecr[0]=f1fit->GetParameter(0);
vecr[1]=f1fit->GetParError(0);
vecr[2]=f1fit->GetParameter(1);
vecr[3]=f1fit->GetParError(1);
vecr[4]=avgmult[ibin];
vecr[5]=avgtrk[ibin];
hq2->Draw("Psame");
c4->cd(2)->SetLogy();
c4->cd(2)->SetLeftMargin(0.2);
TH1D* hq2_cp = (TH1D*)hq2->Clone("hq2_cp");
fixedFontHist(hq2_cp,1.0,2.0);
divideByBinCenter(hq2_cp);
hq2_cp->GetYaxis()->SetTitle("#frac{dN}{q_{2}dq_{2}}");
hq2_cp->GetYaxis()->SetRangeUser(1e-10,10);
hq2_cp->SetMarkerStyle(24);
hq2_cp->SetMarkerColor(4);
hq2_cp->SetLineColor(4);
hq2_cp->SetMarkerSize(0.5);
hq2_cp->Fit(Form("ffit"),"R","",0,10);
hq2_cp->Draw("Psame");
t->DrawLatex(0.5,0.2,Form("N_{trk}^{offline} = %.2f", avgtrk[ibin]));

TCanvas *c5 = new TCanvas("c5","c5",1000,500);
c5->Divide(2,1);
c5->cd(1);
fixedFontHist(hq2nonf,1.0,2.0);
hq2nonf->SetTitle("");
hq2nonf->GetXaxis()->SetTitle("q_{2}");
hq2nonf->GetYaxis()->SetTitle("#frac{dN}{dq_{2}}");
hq2nonf->GetYaxis()->SetRangeUser(0,1);
hq2nonf->SetMarkerStyle(24);
hq2nonf->SetMarkerColor(4);

void getResv(){

	double theta[ntheta];
	TFile *f[nFileAlltrack];
		
	for(int itheta=0;itheta<ntheta;itheta++)
        	theta[itheta]=itheta*TMath::Pi()/ntheta/nn;


	string SumorProd = getenv("SUMORPROD");
	string mdir = getenv("DIR");
	ofstream fstrv;
	if(SumorProd=="Sum")fstrv.open("v.txt");
	else fstrv.open("v_2.txt");
	TVectorD totmult[nbin], totpt[nbin];	TVectorD Nevent, totmultall, totmultall_;
	TVectorD totmulthisto[nbin], totmulthistocorr[nbin];
	TVectorD V_int, V_interr;
	TVectorD V_intcorr, V_intcorrerr;
	TVectorD* V_mean;
	TVectorD* deltaV_mean;
	TVectorD* V[nbin]; TVectorD* chi[nbin];
	TVectorD* dDRe[nbin]; TVectorD* dDIm[nbin];
	TVectorD* dNRe[nbin][ntheta]; TVectorD* dNIm[nbin][ntheta];
	TComplex dD[nbin][ntheta], dN[nbin][ntheta][nptv];
	TVectorD avgmult[nbin];	TVectorD avgmultall;
	TVectorD avgpt[nbin];
	TVectorD v[nbin][ntheta],vmean[nbin],deltav[nbin][ntheta],deltavmean[nbin];
	if(SumorProd=="Sum")	TFile *infile = TFile::Open("mergedV_Sum.root");
	else	TFile *infile = TFile::Open("mergedV_Prod.root");
	Nevent.ResizeTo(nbin);	Nevent.Zero();
	totmultall.ResizeTo(nbin);	totmultall.Zero();
	totmultall_.ResizeTo(nbin);	totmultall_.Zero();
	avgmultall.ResizeTo(nbin);
	V_int.ResizeTo(nbin);	V_int.Zero();
	V_interr.ResizeTo(nbin);	V_interr.Zero();
	V_intcorr.ResizeTo(nbin);	V_intcorr.Zero();
	V_intcorrerr.ResizeTo(nbin);	V_intcorrerr.Zero();
	for(int ibin=0;ibin<nbin;ibin++){
		totmult[ibin].ResizeTo(nptv);	totmult[ibin].Zero();
		totmulthisto[ibin].ResizeTo(nptv);	
		totmulthistocorr[ibin].ResizeTo(nptv);
		avgmult[ibin].ResizeTo(nptv);	avgpt[ibin].ResizeTo(nptv);
		totpt[ibin].ResizeTo(nptv);	totpt[ibin].Zero();
		vmean[ibin].ResizeTo(nptv);	deltavmean[ibin].ResizeTo(nptv);
		V[ibin] = (TVectorD*) infile->Get(Form("D_%d/D_0/V",ibin));
		chi[ibin] = (TVectorD*) infile->Get(Form("D_%d/chi",ibin));
		V_mean = (TVectorD*) infile->Get(Form("D_%d/Vmean",ibin));
		deltaV_mean = (TVectorD*) infile->Get(Form("D_%d/deltaVmean",ibin));
		for(int itheta=0;itheta<ntheta;itheta++){
			v[ibin][itheta].ResizeTo(nptv);	deltav[ibin][itheta].ResizeTo(nptv);
			dD[ibin][itheta]=0;
			for(int iptbin=0;iptbin<nptv;iptbin++){
				dN[ibin][itheta][iptbin]=0;
			}
		}
	}

        for(int ifile=0; ifile<nFileAlltrack; ifile++){
	        if(SumorProd=="Sum") f[ifile] = TFile::Open(Form("/lio/lfs/cms/store/user/qixu/flow/pbsjoboutput/PFcandpt01to10tracknormcpt03to6/%s/Anav_Prod_%d.root",mdir.c_str(),ifile));
	        else f[ifile] = TFile::Open(Form("/lio/lfs/cms/store/user/qixu/flow/pbsjoboutput/PFcandpt01to10tracknormcpt03to6/%s/Anav_Prod2_%d.root",mdir.c_str(),ifile));
		TVectorD* Nevent_t = (TVectorD*)f[ifile]->Get("Nevent");	
		TVectorD* totmultall_t = (TVectorD*)f[ifile]->Get("totmultall");
		for(int ibin=0;ibin<nbin;ibin++){
			TVectorD* totmult_t = (TVectorD*)f[ifile]->Get(Form("totmult_%d",ibin));
			TVectorD* totpt_t = (TVectorD*)f[ifile]->Get(Form("totpt_%d",ibin));
			dDRe[ibin] = (TVectorD*)f[ifile]->Get(Form("dDRe_%d",ibin));
			dDIm[ibin] = (TVectorD*)f[ifile]->Get(Form("dDIm_%d",ibin));
			Nevent[ibin]+=(*Nevent_t)[ibin];
			totmultall[ibin]+=(*totmultall_t)[ibin];
			for(int iptbin=0;iptbin<nptv;iptbin++){
				totmult[ibin][iptbin]+=(*totmult_t)[iptbin];
				totpt[ibin][iptbin]+=(*totpt_t)[iptbin];
			}
			for(int itheta=0;itheta<ntheta;itheta++){
				dNRe[ibin][itheta] = (TVectorD*)f[ifile]->Get(Form("dNRe_%d_%d",ibin,itheta));
				dNIm[ibin][itheta] = (TVectorD*)f[ifile]->Get(Form("dNIm_%d_%d",ibin,itheta));
				dD[ibin][itheta]+=TComplex((*dDRe[ibin])[itheta],(*dDIm[ibin])[itheta]);
				for(int iptbin=0;iptbin<nptv;iptbin++)
					dN[ibin][itheta][iptbin]+=TComplex((*dNRe[ibin][itheta])[iptbin],(*dNIm[ibin][itheta])[iptbin]);
			}
		}
	f[ifile]->Close();
	}
	
	TH1D* hpt[nbin];
	TH1D* hpteffcorr[nbin];
	TFile *fhisto = TFile::Open("histomerged.root");
	TFile *feff = TFile::Open("/home/xuq7/HI/flow/LYZ/v2/TrackCorrections_HIJING_538_OFFICIAL_Mar24.root");
	TH2F* heff = (TH2F*)feff->Get("rTotalEff3D");
      	TH1D* hpteff = (TH1D*)heff->ProjectionY("hpteff",heff->GetXaxis()->FindBin(-2.4),heff->GetXaxis()->FindBin(2.4)-1,"o");
	TArrayD *ptarr = (TArrayD*)heff->GetYaxis()->GetXbins();
	double *ptbinhisto = ptarr->GetArray();
	int NbinX = heff->GetXaxis()->GetNbins();
	int NbinY = heff->GetYaxis()->GetNbins();
	hpteff->Scale(1.0/NbinX);
        
	for(int ibin=0;ibin<nbin;ibin++)
		for(int itheta=0;itheta<ntheta;itheta++)
			dD[ibin][itheta]*=TComplex::Power(TComplex::I(),mm-1)/(Double_t)Nevent[ibin];
	
//.........这里部分代码省略.........