void corrhp(const Int_t year, const Int_t iweek, string *week, const Int_t iconfig, const Int_t nbinxbj, Double_t *x, Double_t *err_x, Double_t *A1p, Double_t *err_A1p, Double_t *CorrAV, Double_t *err_CorrAV, Int_t *HADP, Int_t *HADN, Int_t *pRun, Int_t *nRun, Int_t &RunMin, Int_t &RunMax){
  ofstream hout(Form("A1_contas_%i_%s_%i.out", year, week[iweek].c_str(),iconfig)); // output file for configuration i

  FILE *hadout;
  hadout = fopen("hadron.out","a"); // output file for configuration i
  FILE *A1out;
  A1out= fopen(Form("A1_histo_%i_%s_%i.out", year, week[iweek].c_str(),iconfig),"w");
  const int nbinxbj1=12;
  Float_t wxbj[nbinxbj1+1] = {0.003, 0.006, 0.010, 0.020, 0.030, 0.040, 0.060, 0.100, 0.150, 0.200, 0.300, 0.400, 0.700}; // x_Bj window

  //creates the histograms - array allocation 
  Double_t xbinsXBj[101];
  for(Int_t i=0; i<101; i++ ) xbinsXBj[i] = 1e-3 * pow(10,0.03*i);

  TH1D** histoUPP     = new TH1D*[nbinxbj]; // Sum of weights in target U, field positive, positive hadrons
  TH1D** histoDPP     = new TH1D*[nbinxbj]; // Sum of weights in target D, field positive, positive hadrons
  TH1D** histoUNP     = new TH1D*[nbinxbj]; // Sum of weights in target U, field negative, positive hadrons
  TH1D** histoDNP     = new TH1D*[nbinxbj]; // Sum of weights in target D, field negative, positive hadrons
  TH1D** histoUPN     = new TH1D*[nbinxbj]; // Sum of weights in target U, field positive, negative hadrons
  TH1D** histoDPN     = new TH1D*[nbinxbj]; // Sum of weights in target D, field positive, negative hadrons
  TH1D** histoUNN     = new TH1D*[nbinxbj]; // Sum of weights in target U, field negative, negative hadrons
  TH1D** histoDNN     = new TH1D*[nbinxbj]; // Sum of weights in target D, field negative, negative hadrons
  TH1D** histoBetaUPP  = new TH1D*[nbinxbj]; // Sum of weights^{2}*Beta in target U, field positive, positive hadrons
  TH1D** histoBetaDPP  = new TH1D*[nbinxbj]; // Sum of weights^{2}*Beta in target D, field positive, positive hadrons
  TH1D** histoBetaUNP  = new TH1D*[nbinxbj]; // Sum of weights^{2}*Beta in target U, field negative, positive hadrons
  TH1D** histoBetaDNP  = new TH1D*[nbinxbj]; // Sum of weights^{2}*Beta in target D, field negative, positive hadrons
  TH1D** histoBetaUPN  = new TH1D*[nbinxbj]; // Sum of weights^{2}*Beta in target U, field positive, negative hadrons
  TH1D** histoBetaDPN  = new TH1D*[nbinxbj]; // Sum of weights^{2}*Beta in target D, field positive, negative hadrons
  TH1D** histoBetaUNN  = new TH1D*[nbinxbj]; // Sum of weights^{2}*Beta in target U, field negative, negative hadrons
  TH1D** histoBetaDNN  = new TH1D*[nbinxbj]; // Sum of weights^{2}*Beta in target D, field negative, negative hadrons


  for(int i=0; i<nbinxbj; i++ ) {
    histoUPP[i]      = new TH1D(Form("histoUPP%d",i),"A1_{1,d} - Sum of weights in target U, B +",100, xbinsXBj);
    histoDPP[i]      = new TH1D(Form("histoDPP%d",i),"A1_{1,d} - Sum of weights in target D, B +",100, xbinsXBj);
    histoUNP[i]      = new TH1D(Form("histoUNP%d",i),"A1_{1,d} - Sum of weights in target U, B -",100, xbinsXBj);
    histoDNP[i]      = new TH1D(Form("histoDNP%d",i),"A1_{1,d} - Sum of weights in target D, B -",100, xbinsXBj);
    histoUPN[i]      = new TH1D(Form("histoUPN%d",i),"A1_{1,d} - Sum of weights in target U, B +",100, xbinsXBj);
    histoDPN[i]      = new TH1D(Form("histoDPN%d",i),"A1_{1,d} - Sum of weights in target D, B +",100, xbinsXBj);
    histoUNN[i]      = new TH1D(Form("histoUNN%d",i),"A1_{1,d} - Sum of weights in target U, B -",100, xbinsXBj);
    histoDNN[i]      = new TH1D(Form("histoDNN%d",i),"A1_{1,d} - Sum of weights in target D, B -",100, xbinsXBj);
    histoBetaUPP[i]   = new TH1D(Form("histoBetaUPP%d",i),"A1_{1,d} - Sum of weights*beta in target U, B +",100, xbinsXBj);
    histoBetaDPP[i]   = new TH1D(Form("histoBetaDPP%d",i),"A1_{1,d} - Sum of weights*beta in target D, B +",100, xbinsXBj);
    histoBetaUNP[i]   = new TH1D(Form("histoBetaUNP%d",i),"A1_{1,d} - Sum of weights*beta in target U, B -",100, xbinsXBj);
    histoBetaDNP[i]   = new TH1D(Form("histoBetaDNP%d",i),"A1_{1,d} - Sum of weights*beta in target D, B -",100, xbinsXBj);
    histoBetaUPN[i]   = new TH1D(Form("histoBetaUPN%d",i),"A1_{1,d} - Sum of weights*beta in target U, B +",100, xbinsXBj);
    histoBetaDPN[i]   = new TH1D(Form("histoBetaDPN%d",i),"A1_{1,d} - Sum of weights*beta in target D, B +",100, xbinsXBj);
    histoBetaUNN[i]   = new TH1D(Form("histoBetaUNN%d",i),"A1_{1,d} - Sum of weights*beta in target U, B -",100, xbinsXBj);
    histoBetaDNN[i]   = new TH1D(Form("histoBetaDNN%d",i),"A1_{1,d} - Sum of weights*beta in target D, B -",100, xbinsXBj);
  }

  //////////////////////////////////  Begining of field positive ////////////////////////////////////////////// 


//   FILE * pFile;FILE * nFile;
//   char pLine [10000], nLine [10000]; 
//   string stpor;
//   string stnpor;
//   Int_t runp[100], runn[100];
//   size_t posp, posn;
  // Int_t n = 0;

  TChain * t_p = new TChain  ("USRSI");
//   gSystem->Exec(Form("awk '/%i / {print substr($0,5,300)}' /comp02/comp/helena/DB/Grouping_%i_%s.list > mytree_%i_%s_%i.list",pRun[iconfig],year,week[iweek].c_str(),year,week[iweek].c_str(),pRun[iconfig]));
//   pFile = fopen (Form("mytree_%i_%s_%i.list",year,week[iweek].c_str(),pRun[iconfig]), "r"); 
//   while (!feof(pFile)) {
//     fgetc (pFile);
//     n++;
//   }
//   rewind (pFile);
//   stpor=fgets (pLine , 10000 , pFile);
//   posp = stpor.find(" ");
//   for (int i=0; i<n/6;i++){
//     sscanf (stpor.substr (posp,6).c_str(), "%d", &runp[i]);
//     posp = posp+6;
    t_p->Add("his_35259-2-7.root.03");
	     // }
    //  fclose (pFile);


  int Run; 
  int Nout;
  int Hcharge;
  long long EvN;
  int minRp = 0;
  int maxRp = 0;
  int minRn = 0;
  int maxRn = 0;
  float Xbj, Q2;
  Float_t like[6];
  Int_t TM;
  Double_t P_UP, P_CENTER, P_DOWN, weight, z;

  bool inTargetD;
  bool inTargetU;
  bool inTargetC;
  bool LastHadron, Muon, TI;

  int UPEP[nbinxbj], UPEN[nbinxbj], DPEP[nbinxbj], DPEN[nbinxbj];
  int UNEP[nbinxbj], UNEN[nbinxbj], DNEP[nbinxbj], DNEN[nbinxbj];
  double rfactor[nbinxbj], err_rfactor[nbinxbj]; 
  double EvUPP[nbinxbj],EvDPP[nbinxbj],EvUNP[nbinxbj],EvDNP[nbinxbj], EvUPN[nbinxbj],EvDPN[nbinxbj],EvUNN[nbinxbj],EvDNN[nbinxbj];// Ev:= event; UP,DP := Upstream, Downstream (field Positive); UN,DN := Upstream, Downstream (field Negative); last P or N:= hadron Positive or Negative 
  double CovUP[nbinxbj],CovDP[nbinxbj],CovUN[nbinxbj],CovDN[nbinxbj]; 
  double CorrUP[nbinxbj],CorrDP[nbinxbj],CorrUN[nbinxbj],CorrDN[nbinxbj]; 
  double GetAcceptance[nbinxbj];
  float aA = 0.967, bA = 0.256; //parameters of ratio of acceptances

  int hadp=0;
  int hadn=0;
  int hadp_z2=0;
  int hadn_z2=0;
  for (Int_t i = 0; i < nbinxbj; i++){ 
    HADP[i] = 0;
    HADN[i] = 0;
    EvUPP[i]=0; EvDPP[i]=0; EvUPN[i]=0; EvDPN[i]=0; CovUP[i]=0; CovDP[i]=0; CorrUP[i]=0; CorrDP[i]=0;
    UPEP[i]=0; UPEN[i]=0; DPEP[i]=0; DPEN[i]=0, rfactor[i]=0, err_rfactor[i]=0;
  }

  Int_t PrEvN = -1;

  t_p->SetBranchAddress("weight",&weight); 
  t_p->SetBranchAddress("Xbj",&Xbj); 
  t_p->SetBranchAddress("Q2",&Q2); 
  t_p->SetBranchAddress("P_UP",&P_UP);
  t_p->SetBranchAddress("P_DOWN",&P_DOWN);
  t_p->SetBranchAddress("P_CENTER",&P_CENTER);
  t_p->SetBranchAddress("inTargetU",&inTargetU); 
  t_p->SetBranchAddress("inTargetC",&inTargetC);  
  t_p->SetBranchAddress("inTargetD",&inTargetD);  
  t_p->SetBranchAddress("Run",&Run); 
  t_p->SetBranchAddress("EvN",&EvN); 
  t_p->SetBranchAddress("Nout",&Nout);
  t_p->SetBranchAddress("LastHadron",&LastHadron);
  t_p->SetBranchAddress("Muon",&Muon);
  t_p->SetBranchAddress("Hcharge",&Hcharge);
  t_p->SetBranchAddress("z",&z);
  t_p->SetBranchAddress("like",&like);
  t_p->SetBranchAddress("TM",&TM);  

  //reads all entries and fills the histograms
  Int_t nentries = (Int_t)t_p->GetEntries(); 
  for (Int_t j=0; j<nentries; j++) { 
  
    t_p->GetEntry(j); 
    if( (TM&8) || (TM&256) || (TM&1024) )    TI = true; // O, M, J/psi
    else 							  TI = false;
    cout<<z<<" "<<Nout<<endl;
    if(Hcharge > 0) {
      hadp++;
    }else if(Hcharge > 0) {
      hadn++;
    }
    //  if (Nout > 1){

    if  (z > 0.2 && z < 1 && !Muon){
      if(Hcharge > 0) {
	for (Int_t i=0; i < nbinxbj; i++){
	  if (Xbj >= wxbj[i] && Xbj < wxbj[i+1]){
	    if (Run < 41399) {
	      if(inTargetU){
		 
		histoUPP[i] ->Fill(Xbj,weight); 
		histoBetaUPP[i]->Fill(Xbj,P_UP*weight*weight);
	      } else if(inTargetD){
		histoDPP[i] ->Fill(Xbj,weight);
		histoBetaDPP[i]->Fill(Xbj,P_DOWN*weight*weight);
	      }      
	    } else {
	      if(inTargetU){
		histoUPP[i] ->Fill(Xbj,weight); 
		histoBetaUPP[i]->Fill(Xbj,P_UP*weight*weight);
	      }else if(inTargetD){
		histoUPP[i] ->Fill(Xbj,weight); 
		histoBetaUPP[i]->Fill(Xbj,P_DOWN*weight*weight);
	      }else if (inTargetC){
		histoDPP[i] ->Fill(Xbj,weight);
		histoBetaDPP[i]->Fill(Xbj,P_CENTER*weight*weight);
	      }
	    }
	  }
	}
      }	else {
	for (Int_t i=0; i < nbinxbj; i++){
	  if (Xbj >= wxbj[i] && Xbj < wxbj[i+1]){
	    if (Run < 41399) {
	      if(inTargetU){
		histoUPN[i] ->Fill(Xbj,weight); 
		histoBetaUPN[i]->Fill(Xbj,P_UP*weight*weight); 
	      }else if(inTargetD){
		histoDPN[i] ->Fill(Xbj,weight);
		histoBetaDPN[i]->Fill(Xbj,P_DOWN*weight*weight);
	      }
	    }	else {
	      if(inTargetU){
		histoUPN[i] ->Fill(Xbj,weight); 
		histoBetaUPN[i]->Fill(Xbj,P_UP*weight*weight);
	      }else if(inTargetD){
		histoUPN[i] ->Fill(Xbj,weight); 
		histoBetaUPN[i]->Fill(Xbj,P_DOWN*weight*weight);
	      }else if (inTargetC){
		histoDPN[i] ->Fill(Xbj,weight);
		histoBetaDPN[i]->Fill(Xbj,P_CENTER*weight*weight);
	      }
	    }
	  }
	}
      }
    }//   if (z > 0.2 && z < 1){
    //}//   if (Nout > 1){
    if (Run < 41399) {
      if(inTargetU){//for the correlation matrices
	for (Int_t i=0; i < nbinxbj; i++){
	  if(EvN != PrEvN) {
	    UPEP[i]=0; UPEN[i]=0;
	  }
	  if (Xbj >= wxbj[i] && Xbj < wxbj[i+1]){
	    if(Hcharge < 0 && z > 0.2 && z < 1 && !Muon) { //and not a muon
	      UPEP[i]++;
	    }else if((Muon && (TI || (!TI && Nout > 1)))){//else if is muon
	      UPEN[i]=1;
	    }
	    if(LastHadron){
	      EvUPP[i] += UPEP[i]*weight*UPEP[i]*weight;
	      EvUPN[i] += UPEN[i]*weight*UPEN[i]*weight;
	      CovUP[i] += UPEP[i]*weight*UPEN[i]*weight;
	    }
	  }
	}
      }else if(inTargetD){
	for (Int_t i=0; i < nbinxbj; i++){
	  if(EvN != PrEvN) {
	    DPEP[i]=0; DPEN[i]=0;
	  }
	  if (Xbj >= wxbj[i] && Xbj < wxbj[i+1]){
	    if(Hcharge < 0 && z > 0.2 && z < 1 && !Muon) {
	      DPEP[i]++;
	      //cout<<"h+"<<LastHadron<<"evn"<<EvN<<endl;
	      //	    }else if(Muon){
	    }else if((Muon && (TI || (!TI && Nout > 1)))){//else if is muon
	      //cout<<"h-"<<LastHadron<<"evn"<<EvN<<endl;
	      DPEN[i]=1;
	    }
	    if(LastHadron){
	      // cout<<LastHadron<<"evn"<<EvN<<endl;
	      EvDPP[i] += DPEP[i]*weight*DPEP[i]*weight;
	      EvDPN[i] += DPEN[i]*weight*DPEN[i]*weight;
	      CovDP[i] += DPEP[i]*weight*DPEN[i]*weight;
	    }
	  }	  
	}
      }
    }else {
      if(inTargetU){
	for (Int_t i=0; i < nbinxbj; i++){
	  if(EvN != PrEvN) {
	    UPEP[i]=0; UPEN[i]=0;
	  }
	  if (Xbj >= wxbj[i] && Xbj < wxbj[i+1]){
	    if(Hcharge < 0 && z > 0.2 && z < 1 && !Muon) {
	      hadp_z2++;
	      UPEP[i]++;
	    }else if((Muon && (TI || (!TI && Nout > 1)))){//else if is muon
	      //	    }else if(Muon){
	      hadn_z2++;
	      UPEN[i]=1;
	    }
	    if(LastHadron){
	      EvUPP[i] += UPEP[i]*weight*UPEP[i]*weight;
	      EvUPN[i] += UPEN[i]*weight*UPEN[i]*weight;
	      CovUP[i] += UPEP[i]*weight*UPEN[i]*weight;
	    }
	  }
	}
      }else if(inTargetD){
	for (Int_t i=0; i < nbinxbj; i++){
	  if(EvN != PrEvN) {
	    UPEP[i]=0; UPEN[i]=0;
	  }
	  if (Xbj >= wxbj[i] && Xbj < wxbj[i+1]){
	    if(Hcharge < 0 && z > 0.2 && z < 1 && !Muon) {
	      hadp_z2++;
	      UPEP[i]++;
	    }else if((Muon && (TI || (!TI && Nout > 1)))){//else if is muon
	      //	    }else if(Muon){
	      hadn_z2++;
	      UPEN[i]=1;
	    }
	    if(LastHadron){
	      EvUPP[i] += UPEP[i]*weight*UPEP[i]*weight;
	      EvUPN[i] += UPEN[i]*weight*UPEN[i]*weight;
	      CovUP[i] += UPEP[i]*weight*UPEN[i]*weight;
	    }
	  }	  
	}
      }else if(inTargetC){
	for (Int_t i=0; i < nbinxbj; i++){
	  if(EvN != PrEvN) {
	    DPEP[i]=0; DPEN[i]=0;
	  }
	  if (Xbj >= wxbj[i] && Xbj < wxbj[i+1]){
	    if(Hcharge < 0 && z > 0.2 && z < 1 && !Muon) {
	      hadp_z2++;
	      DPEP[i]++;
	    }else if((Muon && (TI || (!TI && Nout > 1)))){//else if is muon
	      //	    }else if(Muon){
	      hadn_z2++;
	      DPEN[i]=1;
	    }
	    if(LastHadron){
	      EvDPP[i] += DPEP[i]*weight*DPEP[i]*weight;
	      EvDPN[i] += DPEN[i]*weight*DPEN[i]*weight;
	      CovDP[i] += DPEP[i]*weight*DPEN[i]*weight;
	    }
	  }	  
	}
      }
    }

    // for next hadron
    PrEvN = EvN;

  
    ///////////////end of filling the correlation matrix
    if(j == 1) { // search min. run number
      minRp = Run;
    }else if(j==nentries-1){ // search max. run number
      maxRp = Run;
    }
  }//  for (Int_t j=0; j<nentries; j++) {

  fprintf  (A1out, "hadp = %i \t hadn = %i \n", hadp, hadn);

  //     // Draw the histograms and get bin contents
  //     TCanvas* c1 = new TCanvas("c1",Form("Sum of weights in target U, B +, configuration %i", config),500,0,  900, 638 );
  //     c1->Divide(3,4);

  //     for (Int_t i = 0; i < nbinxbj; i++){
  //       c1->cd(i+1);
  //       histoUP[i]->Draw();
  //       gPad->SetLogx();
  //       c1->Update(); 
  //     }
  //     c1->Print(Form("/comp01/na50/helena/compass/%i/A1/%s/2ndorder/consecutive/semi-inclusive/WUP_%i_%s.ps", year, week, config));
  //     TCanvas* c3 = new TCanvas("c3",Form("Sum of weights in target D, B +, configuration %i", config),500,0,  900, 638 );
  //     c3->Divide(3,4);
  //     for (Int_t i = 0; i < nbinxbj; i++){
  //       c3->cd(i+1);
  //       histoDP[i]->Draw();
  //       gPad->SetLogx();
  //       c3->Update(); 
  //     }
  
  //     c3->Print(Form("/comp01/na50/helena/compass/%i/A1/%s/2ndorder/consecutive/semi-inclusive/WDP_%i_%s.ps", year, week, config));
  //     TCanvas* c5 = new TCanvas("c5",Form("Sum of weights*Beta in target U, B +, configuration %i", config),500,0,  900, 638 );
  //     c5->Divide(3,4);
  //     for (Int_t i = 0; i < nbinxbj; i++){
  //       c5->cd(i+1);
  //       histoBetaUP[i]->Draw();
  //       gPad->SetLogx();
  //       c5->Update(); 
  //     }
  //     c5->Print(Form("/comp01/na50/helena/compass/%i/A1/%s/2ndorder/consecutive/semi-inclusive/WUPBeta_%i_%s.ps", year, week, config));
  //     TCanvas* c6 = new TCanvas("c6",Form("Sum of weights*Beta in target D, B +, configuration %i", config),500,0,  900, 638 );
  //     c6->Divide(3,4);
  //     for (Int_t i = 0; i < nbinxbj; i++){
  //       c6->cd(i+1);
  //       histoBetaDP[i]->Draw();
  //       gPad->SetLogx();
  //       c6->Update(); 
  //     }
  //     c6->Print(Form("/comp01/na50/helena/compass/%i/A1/%s/2ndorder/consecutive/semi-inclusive/WDPBeta_%i_%s.ps", year, week, config));

  Double_t sumUPP[nbinxbj], sumDPP[nbinxbj]; 
  Double_t sumUPN[nbinxbj], sumDPN[nbinxbj];
  Double_t sumBetaUPP[nbinxbj], sumBetaDPP[nbinxbj];
  Double_t sumBetaUPN[nbinxbj], sumBetaDPN[nbinxbj];

  for(Int_t i=0; i < nbinxbj; i++ ) {
    sumUPP[i] = 0.0; 
    sumDPP[i] = 0.0;
    sumUPN[i] = 0.0; 
    sumDPN[i] = 0.0;
    sumBetaUPP[i] = 0.0; 
    sumBetaDPP[i] = 0.0;
    sumBetaUPN[i] = 0.0; 
    sumBetaDPN[i] = 0.0;
  }

  for(Int_t i = 0; i < nbinxbj; i++){
    for(Int_t j=0; j<101; j++ ) {
      sumUPP[i] += histoUPP[i]->GetBinContent(j); 
      sumDPP[i] += histoDPP[i]->GetBinContent(j); 
      sumUPN[i] += histoUPN[i]->GetBinContent(j); 
      sumDPN[i] += histoDPN[i]->GetBinContent(j); 
      sumBetaUPP[i] += histoBetaUPP[i]->GetBinContent(j); 
      sumBetaDPP[i] += histoBetaDPP[i]->GetBinContent(j); 
      sumBetaUPN[i] += histoBetaUPN[i]->GetBinContent(j); 
      sumBetaDPN[i] += histoBetaDPN[i]->GetBinContent(j); 
    }
  }
  //////////////////////////////////  End of field positive ///////////////////////////////////////////////////

  //////////////////////////////////////////////////////////////////////////////////////////////////////
  //////////////////////////////////  Begining of field negative //////////////////////////////////////////////


  for (Int_t i=0; i < nbinxbj; i++){
    EvUNP[i]=0; EvDNP[i]=0; EvUNN[i]=0; EvDNN[i]=0; CovUN[i]=0; CovDN[i]=0; CorrUN[i]=0; CorrDN[i]=0;
    UNEP[i]=0; UNEN[i]=0; DNEP[i]=0; DNEN[i]=0;
  }

  Int_t hadp1=0;
  Int_t hadn1=0;

  PrEvN = -1;

	     //  n = 0;
  TChain * t_n= new TChain  ("USRSI");
//   gSystem->Exec(Form("awk '/%i / {print substr($0,5,300)}' /comp02/comp/helena/DB/Grouping_%i_%s.list > mytree_%i_%s_%i.list",nRun[iconfig],year,week[iweek].c_str(),year,week[iweek].c_str(),nRun[iconfig]));
//   nFile = fopen (Form("mytree_%i_%s_%i.list",year,week[iweek].c_str(),nRun[iconfig]), "r"); 
//   while (!feof(nFile)) {
//     fgetc (nFile);
//     n++;
//   }
//   rewind (nFile);
//   stnpor=fgets (nLine , 10000 , nFile);
//   posn = stnpor.find(" ");
//   for (int i=0; i<n/6;i++){
//     sscanf (stnpor.substr (posn,6).c_str(), "%d", &runn[i]);
//     posn = posn+6;
    //    t_n->Add(Form("/hometmp/comp/helena/%i/%s/histo/his_%i*.root*",year,week[iweek].c_str(),runn[i]));
    t_n->Add("his_35259-2-7.root.02");
	     //  }
    //  fclose (nFile);

 
  t_n->SetBranchAddress("inTargetD",&inTargetD);  
  t_n->SetBranchAddress("Run",&Run); 
  t_n->SetBranchAddress("EvN",&EvN); 
  t_n->SetBranchAddress("Nout",&Nout);
  t_n->SetBranchAddress("LastHadron",&LastHadron);
  t_n->SetBranchAddress("Muon",&Muon);
  t_n->SetBranchAddress("Hcharge",&Hcharge);
  t_n->SetBranchAddress("z",&z);
  t_n->SetBranchAddress("like",&like);
  t_n->SetBranchAddress("weight",&weight); 
  t_n->SetBranchAddress("Xbj",&Xbj); 
  t_n->SetBranchAddress("Q2",&Q2); 
  t_n->SetBranchAddress("P_UP",&P_UP);
  t_n->SetBranchAddress("P_DOWN",&P_DOWN);
  t_n->SetBranchAddress("P_CENTER",&P_CENTER);
  t_n->SetBranchAddress("inTargetU",&inTargetU);
  t_n->SetBranchAddress("inTargetC",&inTargetC);
  t_n->SetBranchAddress("TM",&TM);  
  //read all entries and fill the histograms

  nentries = (Int_t)t_n->GetEntries(); 
  for (Int_t j=0; j<nentries; j++) { 
    t_n->GetEntry(j); 
    if( (TM&8) || (TM&256) || (TM&1024) )    TI = true; // O, M, J/psi
    else 							  TI = false;
    if(Hcharge > 0) {
      hadp1++;
    }else if(Hcharge > 0) {
      hadn1++;
    }
    //    if (Nout > 1){
    if  (z > 0.2 && z < 1 && !Muon){
      if(Hcharge > 0) {
	hadp_z2++;
	for (Int_t i=0; i < nbinxbj; i++){
	  if (Xbj >= wxbj[i] && Xbj < wxbj[i+1]){
	    if (Run < 41399) {
	      if(inTargetU){
		histoUNP[i] ->Fill(Xbj,weight); 
		histoBetaUNP[i]->Fill(Xbj,(P_UP)*weight*weight); 
	      }else if (inTargetD){
		histoDNP[i] ->Fill(Xbj,weight);
		histoBetaDNP[i]->Fill(Xbj,(P_DOWN)*weight*weight);
	      }
	    } else {
	      if(inTargetU){
		histoUNP[i] ->Fill(Xbj,weight); 
		histoBetaUNP[i]->Fill(Xbj,(P_UP)*weight*weight);
	      }else if(inTargetD){
		histoUNP[i] ->Fill(Xbj,weight); 
		histoBetaUNP[i]->Fill(Xbj,(P_DOWN)*weight*weight);
	      }else if (inTargetC){
		histoDNP[i] ->Fill(Xbj,weight);
		histoBetaDNP[i]->Fill(Xbj,(P_CENTER)*weight*weight);
	      }
	    }
	  }
	}
      }
      if(Hcharge > 0) {
	hadn_z2++;
	for (Int_t i=0; i < nbinxbj; i++){
	  if (Xbj >= wxbj[i] && Xbj < wxbj[i+1]){
	    if (Run < 41399) {
	      if(inTargetU){
		histoUNN[i] ->Fill(Xbj,weight); 
		histoBetaUNN[i]->Fill(Xbj,(P_UP)*weight*weight); 
	      }else if (inTargetD){
		histoDNN[i] ->Fill(Xbj,weight);
		histoBetaDNN[i]->Fill(Xbj,(P_DOWN)*weight*weight);
	      }
	    } else {
	      if(inTargetU){
		histoUNN[i] ->Fill(Xbj,weight); 
		histoBetaUNN[i]->Fill(Xbj,(P_UP)*weight*weight);
	      }else if(inTargetD){
		histoUNN[i] ->Fill(Xbj,weight); 
		histoBetaUNN[i]->Fill(Xbj,(P_DOWN)*weight*weight);
	      }else if (inTargetC){
		histoDNN[i] ->Fill(Xbj,weight);
		histoBetaDNN[i]->Fill(Xbj,(P_CENTER)*weight*weight);
	      }
	    }
	  } 
	}
      }
    }//   if (z > 0.2 && z < 1){
    //    }//   if (Nout > 1){
    if (Run < 41399) {
      if(inTargetU){//for the correlation matrices
	for (Int_t i=0; i < nbinxbj; i++){
	  if(EvN != PrEvN) {
	    UNEP[i]=0; UNEN[i]=0;
	  }
	  if (Xbj >= wxbj[i] && Xbj < wxbj[i+1]){
	    if(Hcharge < 0 && z > 0.2 && z < 1 && !Muon) {
	      hadp_z2++;
	      UNEP[i]++;
	    }else if((Muon && (TI || (!TI && Nout > 1)))){//else if is muon
	      //	    }else if(Muon){
	      hadn_z2++;
	      UNEN[i]=1;
	    }
	    if(LastHadron){
	      EvUNP[i] += UNEP[i]*weight*UNEP[i]*weight;
	      EvUNN[i] += UNEN[i]*weight*UNEN[i]*weight;
	      CovUN[i] += UNEP[i]*weight*UNEN[i]*weight;
	    }
	  }
	}
      }else if(inTargetD){
	for (Int_t i=0; i < nbinxbj; i++){
	  if(EvN != PrEvN) {
	    DNEP[i]=0; DNEN[i]=0;
	  }
	  if (Xbj >= wxbj[i] && Xbj < wxbj[i+1]){
	    if(Hcharge < 0 && z > 0.2 && z < 1 && !Muon) {
	      hadp_z2++;
	      DNEP[i]++;
	    }else if((Muon && (TI || (!TI && Nout > 1)))){//else if is muon
	      //	    }else if(Muon){
	      hadn_z2++;
	      DNEN[i]=1;
	    }
	    if(LastHadron){
	      EvDNP[i] += DNEP[i]*weight*DNEP[i]*weight;
	      EvDNN[i] += DNEN[i]*weight*DNEN[i]*weight;
	      CovDN[i] += DNEP[i]*weight*DNEN[i]*weight;
	    }
	  }	  
	}
      }
    }else {
      if(inTargetU){//for the correlation matrices
	for (Int_t i=0; i < nbinxbj; i++){
	  if(EvN != PrEvN) {
	    UNEP[i]=0; UNEN[i]=0;
	  }
	  if (Xbj >= wxbj[i] && Xbj < wxbj[i+1]){
	    if(Hcharge < 0 && z > 0.2 && z < 1 && !Muon) {
	      hadp_z2++;
	      UNEP[i]++;
	    }else if((Muon && (TI || (!TI && Nout > 1)))){//else if is muon
	      //	    }else if(Muon){
	      hadn_z2++;
	      UNEN[i]=1;
	    }
	    if(LastHadron){
	      EvUNP[i] += UNEP[i]*weight*UNEP[i]*weight;
	      EvUNN[i] += UNEN[i]*weight*UNEN[i]*weight;
	      CovUN[i] += UNEP[i]*weight*UNEN[i]*weight;
	    }
	  }
	}
      }else if(inTargetD){/////
	for (Int_t i=0; i < nbinxbj; i++){
	  if(EvN != PrEvN) {
	    UNEP[i]=0; UNEN[i]=0;
	  }
	  if (Xbj >= wxbj[i] && Xbj < wxbj[i+1]){
	    if(Hcharge < 0 && z > 0.2 && z < 1 && !Muon) {
	      hadp_z2++;
	      UNEP[i]++;
	    }else if((Muon && (TI || (!TI && Nout > 1)))){//else if is muon
	      //	    }else if(Muon){
	      hadn_z2++;
	      UNEN[i]=1;
	    }
	    if(LastHadron){
	      EvUNP[i] += UNEP[i]*weight*UNEP[i]*weight;
	      EvUNN[i] += UNEN[i]*weight*UNEN[i]*weight;
	      CovUN[i] += UNEP[i]*weight*UNEN[i]*weight;
	    }
	  }	  
	}
      }else if(inTargetC){
	for (Int_t i=0; i < nbinxbj; i++){
	  if(EvN != PrEvN) {
	    DNEP[i]=0; DNEN[i]=0;
	  }
	  if (Xbj >= wxbj[i] && Xbj < wxbj[i+1]){
	    if(Hcharge < 0 && z > 0.2 && z < 1 && !Muon) {
	      hadp_z2++;
	      DNEP[i]++;
	    }else if((Muon && (TI || (!TI && Nout > 1)))){//else if is muon
	      //	    }else if(Muon){
	      hadn_z2++;
	      DNEN[i]=1;
	    }
	    if(LastHadron){
	      EvDNP[i] += DNEP[i]*weight*DNEP[i]*weight;
	      EvDNN[i] += DNEN[i]*weight*DNEN[i]*weight;
	      CovDN[i] += DNEP[i]*weight*DNEN[i]*weight;
	    }
	  }	  
	}
      }
    }
  
    // for next hadron
    PrEvN = EvN;

    if(j == 1) { // search min. run number
      minRn = Run;
    }else if(j==nentries-1){ // search max. run number
      maxRn = Run;
    }
  }
 
  fprintf  (A1out, "hadp1 = %i \t hadn1 = %i \n", hadp1, hadn1);
  //     // Draw the histograms and get bin contents
  //     TCanvas* c7 = new TCanvas("c7",Form("Sum of weights in target U, B -, configuration %i", config),500,0,  900, 638 );
  //     c7->Divide(3,4);
  //     for (Int_t i = 0; i < nbinxbj; i++){
  //       c7->cd(i+1);
  //       histoUN[i]->Draw();
  //       gPad->SetLogx();
  //       c7->Update(); 
  //     }

  //     c7->Print(Form("/comp01/na50/helena/compass/%i/A1/%s/2ndorder/consecutive/semi-inclusive/WUN_%i_%s.ps", year, week, config));
  //     TCanvas* c9 = new TCanvas(Form("c9_%i",config),Form("Sum of weights in target D, B -, configuration %i", config),500,0,  900, 638 );
  //     c9->Divide(3,4);
  //     for (Int_t i = 0; i < nbinxbj; i++){
  //       c9->cd(i+1);
  //       histoDN[i]->Draw();
  //       gPad->SetLogx();
  //       c9->Update(); 
  //     }
  //     c9->Print(Form("/comp01/na50/helena/compass/%i/A1/%s/2ndorder/consecutive/semi-inclusive/WDN_%i_%s.ps", year, week, config));
  //     TCanvas* c11 = new TCanvas("c11",Form("Sum of weights*Beta in target U, B -, configuration %i", config),500,0,  900, 638 );
  //     c11->Divide(3,4);
  //     for (Int_t i = 0; i < nbinxbj; i++){
  //       c11->cd(i+1);
  //       histoBetaUN[i]->Draw();
  //       gPad->SetLogx();
  //       c11->Update(); 
  //     }
  //     c11->Print(Form("/comp01/na50/helena/compass/%i/A1/%s/2ndorder/consecutive/semi-inclusive/WUNBeta_%i_%s.ps", year, week, config));
  //     TCanvas* c12 = new TCanvas("c12",Form("Sum of weights*Beta in target D, B -, configuration %i", config),500,0,  900, 638 );
  //     c12->Divide(3,4);
  //     for (Int_t i = 0; i < nbinxbj; i++){
  //       c12->cd(i+1);
  //       histoBetaDN[i]->Draw();
  //       gPad->SetLogx();
  //       c12->Update(); 
  //     }
  //     c12->Print(Form("/comp01/na50/helena/compass/%i/A1/%s/2ndorder/consecutive/semi-inclusive/WDNBeta_%i_%s.ps", year, week, config));


  Double_t sumUNP[nbinxbj],  sumDNP[nbinxbj];
  Double_t sumUNN[nbinxbj],  sumDNN[nbinxbj];
  Double_t sumBetaUNP[nbinxbj], sumBetaDNP[nbinxbj];
  Double_t sumBetaUNN[nbinxbj], sumBetaDNN[nbinxbj];
 
  for(Int_t i=0; i < nbinxbj; i++ ) {
    sumUNP[i]  = 0.0; 
    sumDNP[i]  = 0.0;
    sumUNN[i]  = 0.0; 
    sumDNN[i]  = 0.0;
    sumBetaUNP[i] = 0.0; 
    sumBetaDNP[i] = 0.0;
    sumBetaUNN[i] = 0.0; 
    sumBetaDNN[i] = 0.0;
  }

  for (Int_t i = 0; i < nbinxbj; i++){
    for(int j=0; j<101; j++ ) {
      sumUNP[i] += histoUNP[i]->GetBinContent(j); 
      sumDNP[i] += histoDNP[i]->GetBinContent(j);
      sumUNN[i] += histoUNN[i]->GetBinContent(j); 
      sumDNN[i] += histoDNN[i]->GetBinContent(j);
      sumBetaUNP[i] += histoBetaUNP[i]->GetBinContent(j); 
      sumBetaDNP[i] += histoBetaDNP[i]->GetBinContent(j); 
      sumBetaUNN[i] += histoBetaUNN[i]->GetBinContent(j); 
      sumBetaDNN[i] += histoBetaDNN[i]->GetBinContent(j); 
    }
  } 

  //////////////////////////////////  End of field negative ///////////////////////////////////////////////////
 
  if (minRp < minRn) {
    RunMin = minRp;
  }    else {
    RunMin = minRn;
  }      

  if (maxRp > maxRn) {
    RunMax = maxRp;
  }    else {
    RunMax = maxRn;
  }
  Double_t UPP[nbinxbj], DPP[nbinxbj], UNP[nbinxbj], DNP[nbinxbj];
  Double_t UPN[nbinxbj], DPN[nbinxbj], UNN[nbinxbj], DNN[nbinxbj];

  for (Int_t i = 0; i < nbinxbj; i++){
    UPP[i]=histoUPP[i]->GetEntries();
    DPP[i]=histoDPP[i]->GetEntries();
    UNP[i]=histoUNP[i]->GetEntries();
    DNP[i]=histoDNP[i]->GetEntries();
    UPN[i]=histoUPN[i]->GetEntries();
    DPN[i]=histoDPN[i]->GetEntries();
    UNN[i]=histoUNN[i]->GetEntries();
    DNN[i]=histoDNN[i]->GetEntries();
    fprintf  (A1out, "UPP[%i]\t= %i \t\t DPP[%i]\t= %i \t\t UNP[%i]\t= %i \t\t DNP[%i]\t= %i \n", i, int(UPP[i]), i, int(DPP[i]), i, int(UNP[i]), i, int(DNP[i]));
  }

  fprintf  (A1out,"\n");
  fprintf  (A1out,"\n");
  for (Int_t i = 0; i < nbinxbj; i++){
    fprintf  (A1out,"UPN[%i]\t= %i \t\t DPN[%i]\t= %i \t\t UNN[%i]\t= %i \t\t DNN[%i]\t= %i \n", i, int(UPN[i]), i, int(DPN[i]), i, int(UNN[i]), i, int(DNN[i]));
  }
  fprintf  (A1out,"\n");
  fprintf  (A1out,"\n");

  for (Int_t i = 0; i < nbinxbj; i++){
    fprintf  (A1out, "sumUPP[%i]\t= %f \t sumDPP[%i]\t= %f \t sumUNP[%i]\t= %f \t sumDNP[%i]\t= %f \n", i, sumUPP[i], i, sumDPP[i], i, sumUNP[i], i, sumDNP[i]);
  }
  fprintf  (A1out,"\n");
  fprintf  (A1out,"\n");

  for (Int_t i = 0; i < nbinxbj; i++){
    fprintf  (A1out, "sumUPN[%i]\t= %f \t sumDPN[%i]\t= %f \t sumUNN[%i]\t= %f \t sumDNN[%i]\t= %f \n", i, sumUPN[i], i, sumDPN[i], i, sumUNN[i], i, sumDNN[i]);
  }
  fprintf  (A1out,"\n");
  fprintf  (A1out,"\n");
 
  for (Int_t i = 0; i < nbinxbj; i++){
    fprintf  (A1out, "sumBetaUPP[%i] = %f \t sumBetaDPP[%i] = %f \t\t sumBetaUNP[%i] = %f  \t  sumBetaDNP[%i] = %f \n", i, sumBetaUPP[i], i, sumBetaDPP[i], i, sumBetaUNP[i], i, sumBetaDNP[i]);
  }
  fprintf  (A1out,"\n");
  fprintf  (A1out,"\n");
 
  for (Int_t i = 0; i < nbinxbj; i++){
    fprintf  (A1out, "sumBetaUPN[%i] = %f \t sumBetaDPN[%i] = %f \t\t sumBetaUNN[%i] = %f  \t  sumBetaDNN[%i] = %f \n", i, sumBetaUPN[i], i, sumBetaDPN[i], i, sumBetaUNN[i], i, sumBetaDNN[i]);
  }
  fprintf  (A1out,"\n");
  fprintf  (A1out,"\n");
  fprintf (hadout,"\n");
 
  for(int i=0; i<nbinxbj; i++ ) {
    HADP[i] = int(UPP[i] + DPP[i] + UNP[i] + DNP[i]);
    HADN[i] = int(UPN[i] + DPN[i] + UNN[i] + DNN[i]);
    fprintf (hadout, "HADP[%i] = %f \t HADN[%i] = %f  \t\t year = %i \t period = %s \t config = %i \n", i,HADP[i],i, HADN[i], year, week[iweek].c_str(), iconfig );
  }
 
  //////////////////////////////////  Begining of Asymmetry calculations //////////////////////////////////////////////
  Double_t deltap[nbinxbj], deltan[nbinxbj], err_deltap[nbinxbj], err_deltan[nbinxbj], CovTerm[nbinxbj];
  Double_t a1p[nbinxbj], a2p[nbinxbj], ap[nbinxbj], b1p[nbinxbj], b2p[nbinxbj], bp[nbinxbj], cp[nbinxbj]; //A1+
  Double_t a1n[nbinxbj], a2n[nbinxbj], an[nbinxbj], b1n[nbinxbj], b2n[nbinxbj], bn[nbinxbj], cn[nbinxbj]; //A1-

  Double_t betaUPP[nbinxbj], betaDPP[nbinxbj], betaUNP[nbinxbj], betaDNP[nbinxbj];
  Double_t betaUPN[nbinxbj], betaDPN[nbinxbj], betaUNN[nbinxbj], betaDNN[nbinxbj];

  Double_t d1p[nbinxbj], d2p[nbinxbj], ep[nbinxbj]; //A1+
  Double_t d1n[nbinxbj], d2n[nbinxbj], en[nbinxbj]; //A1-
  Double_t err_e1p[nbinxbj], err_e2p[nbinxbj], err_e3p[nbinxbj], err_ep[nbinxbj];
  Double_t err_e1n[nbinxbj], err_e2n[nbinxbj], err_e3n[nbinxbj], err_en[nbinxbj];
  Double_t A1[nbinxbj], err_A1[nbinxbj];
  for(Int_t i=0; i < nbinxbj; i++ ) {
    A1p[i] = 0.0; 
    err_A1p[i] = 0.0;
    CorrAV[i] = 0.0; 
    err_CorrAV[i] = 0.0;
    A1[i] = 0.0; 
    err_A1[i] = 0.0;
    deltap[i] = 0.0;
    deltan[i] = 0.0;
    err_deltap[i] = 0.0;
    err_deltan[i] = 0.0;
    CovTerm[i]= 0.0;
    betaUPP[i] = 0.0;
    betaUNP[i] = 0.0;
    betaDPP[i] = 0.0;
    betaDNP[i] = 0.0;
    betaUPN[i] = 0.0;
    betaUNN[i] = 0.0;
    betaDPN[i] = 0.0;
    betaDNN[i] = 0.0;
    a1p[i] = 0.0;
    a2p[i] = 0.0;
    ap[i] = 0.0;
    b1p[i] = 0.0;
    b2p[i] = 0.0;
    bp[i] = 0.0;
    cp[i] = 0.0;
    a1n[i] = 0.0;
    a2n[i] = 0.0;
    an[i] = 0.0;
    b1n[i] = 0.0;
    b2n[i] = 0.0;
    bn[i] = 0.0;
    cn[i] = 0.0;
    d1p[i] = 0.0;
    d2p[i] = 0.0;
    ep[i] = 0.0;
    d1n[i] = 0.0;
    d2n[i] = 0.0;
    en[i] = 0.0;

    err_e1p[i] = 0.0;
    err_e2p[i] = 0.0;
    err_e3p[i] = 0.0;
    err_ep[i] = 0.0;
    err_A1p[i] = 0.0;
    err_e1n[i] = 0.0;
    err_e2n[i] = 0.0;
    err_e3n[i] = 0.0;
    err_en[i] = 0.0;
  }


  for(Int_t i=0; i < nbinxbj; i++ ) {
    x[i] = 0.0; 
    err_x[i] = 0.0;
  }
  for (Int_t i = 0; i < nbinxbj; i++){
    x[i] = histoUPP[i]->GetMean(1);
    err_x[i] = histoUPP[i]->GetMeanError(1);
    hout <<"x["<<i<<"] = "<<x[i]<<endl; 
    hout <<"err_x["<<i<<"] = "<<err_x[i] <<endl;
    hout <<" "<<endl;
    GetAcceptance[i] = aA + bA*x[i];
  }

  for (Int_t i = 0; i < nbinxbj; i++){
    fprintf  (A1out, "CovUP[%i]\t= %f \t CovDP[%i]\t= %f \t CovUN[%i]\t= %f \t CovDN[%i]\t= %f \n", i, CovUP[i], i, CovDP[i], i, CovUN[i], i, CovDN[i]);
  }
  fprintf  (A1out,"\n");
  fprintf  (A1out,"\n");

  for (Int_t i = 0; i < nbinxbj; i++){
    fprintf  (A1out, "SigmaUPP[%i]\t= %f \t SigmaDPP[%i]\t= %f \t SigmaUNP[%i]\t= %f \t SigmaDNP[%i]\t= %f \n", i, sqrt(EvUPP[i]), i, sqrt(EvDPP[i]), i, sqrt(EvUNP[i]), i, sqrt(EvDNP[i]));
    fprintf  (A1out, "SigmaUPN[%i]\t= %f \t SigmaDPN[%i]\t= %f \t SigmaUNN[%i]\t= %f \t SigmaDNN[%i]\t= %f \n", i, sqrt(EvUPN[i]), i, sqrt(EvDPN[i]), i, sqrt(EvUNN[i]), i, sqrt(EvDNN[i]));
  }
  fprintf  (A1out,"\n");
  fprintf  (A1out,"\n");


  for (Int_t i = 0; i < nbinxbj; i++){//2nd order weighting method 
    if (sumUNP[i]==0 || sumDPP[i]==0||sumUPP[i]==0 || sumDNP[i]==0||sumUNN[i]==0 || sumDPN[i]==0||sumUPN[i]==0 || sumDNN[i]==0) continue;

    deltap[i] = (sumUPP[i] * sumDNP[i]) / (sumUNP[i] * sumDPP[i]); //A1+
    deltan[i] = (sumUPN[i] * sumDNN[i]) / (sumUNN[i] * sumDPN[i]); //A1-

    err_deltap[i] = sqrt((sumDNP[i] / (sumUNP[i] * sumDPP[i]))*(sumDNP[i] / (sumUNP[i] * sumDPP[i]))*EvUPP[i] + (sumUPP[i] / (sumUNP[i] * sumDPP[i]))*(sumUPP[i] / (sumUNP[i] * sumDPP[i]))*EvDNP[i] + (-(sumUPP[i] * sumDNP[i]) / (sumUNP[i]*sumUNP[i] * sumDPP[i]))*(-(sumUPP[i] * sumDNP[i]) / (sumUNP[i]*sumUNP[i] * sumDPP[i]))*EvUNP[i] +(-(sumUPP[i] * sumDNP[i]) / (sumUNP[i] * sumDPP[i]*sumDPP[i]))*(-(sumUPP[i] * sumDNP[i]) / (sumUNP[i] * sumDPP[i]*sumDPP[i]))*EvDPP[i]) ; //

    err_deltan[i] = sqrt((sumDNN[i] / (sumUNN[i] * sumDPN[i]))*(sumDNN[i] / (sumUNN[i] * sumDPN[i]))*EvUPN[i] + (sumUPN[i] / (sumUNN[i] * sumDPN[i]))*(sumUPN[i] / (sumUNN[i] * sumDPN[i]))*EvDNN[i] + (-(sumUPN[i] * sumDNN[i]) / (sumUNN[i]*sumUNN[i] * sumDPN[i]))*(-(sumUPN[i] * sumDNN[i]) / (sumUNN[i]*sumUNN[i] * sumDPN[i]))*EvUNN[i] +(-(sumUPN[i] * sumDNN[i]) / (sumUNN[i] * sumDPN[i]*sumDPN[i]))*(-(sumUPN[i] * sumDNN[i]) / (sumUNN[i] * sumDPN[i]*sumDPN[i]))*EvDPN[i]) ; //

    hout <<"deltap["<<i<<"] = "<<deltap[i]<<endl;
    hout <<"err_deltap["<<i<<"] = "<<err_deltap[i]<<endl;
    hout <<"deltan["<<i<<"] = "<<deltan[i]<<endl;
    hout <<"err_deltan["<<i<<"] = "<<err_deltan[i]<<endl;

    betaUPP[i] = sumBetaUPP[i] / sumUPP[i]; 
    betaDPP[i] = sumBetaDPP[i] / sumDPP[i]; 
    betaUNP[i] = sumBetaUNP[i] / sumUNP[i]; 
    betaDNP[i] = sumBetaDNP[i] / sumDNP[i];
    betaUPN[i] = sumBetaUPN[i] / sumUPN[i]; 
    betaDPN[i] = sumBetaDPN[i] / sumDPN[i]; 
    betaUNN[i] = sumBetaUNN[i] / sumUNN[i]; 
    betaDNN[i] = sumBetaDNN[i] / sumDNN[i];

    hout <<"betaUPP["<<i<<"] = "<<betaUPP[i]<<endl;
    hout <<"betaDPP["<<i<<"] = "<<betaDPP[i]<<endl;
    hout <<"betaUNP["<<i<<"] = "<<betaUNP[i]<<endl;
    hout <<"betaDNP["<<i<<"] = "<<betaDNP[i]<<endl;
    hout <<"betaUPN["<<i<<"] = "<<betaUPN[i]<<endl;
    hout <<"betaDPN["<<i<<"] = "<<betaDPN[i]<<endl;
    hout <<"betaUNN["<<i<<"] = "<<betaUNN[i]<<endl;
    hout <<"betaDNN["<<i<<"] = "<<betaDNN[i]<<endl;


    a1p[i] = betaUNP[i] * betaDPP[i];
    a2p[i] = betaUPP[i] * betaDNP[i];
    ap[i]  = deltap[i] * a1p[i] - a2p[i];
    a1n[i] = betaUNN[i] * betaDPN[i];
    a2n[i] = betaUPN[i] * betaDNN[i];
    an[i]  = deltan[i] * a1n[i] - a2n[i];

    hout <<"a1p["<<i<<"] = "<<a1p[i]<<endl;
    hout <<"a2p["<<i<<"] = "<<a2p[i]<<endl;
    hout <<"ap["<<i<<"] = "<<ap[i]<<endl;
    hout <<"a1n["<<i<<"] = "<<a1n[i]<<endl;
    hout <<"a2n["<<i<<"] = "<<a2n[i]<<endl;
    hout <<"an["<<i<<"] = "<<an[i]<<endl;

    b1p[i] = betaUNP[i] + betaDPP[i];
    b2p[i] = betaUPP[i] + betaDNP[i];
    bp[i]  = deltap[i] * b1p[i] - b2p[i];
    b1n[i] = betaUNN[i] + betaDPN[i];
    b2n[i] = betaUPN[i] + betaDNN[i];
    bn[i]  = deltan[i] * b1n[i] - b2n[i];

    hout <<"b1p["<<i<<"] = "<<b1p[i]<<endl;
    hout <<"b2p["<<i<<"] = "<<b2p[i]<<endl;
    hout <<"bp["<<i<<"] = "<<bp[i]<<endl;
    hout <<"b1n["<<i<<"] = "<<b1n[i]<<endl;
    hout <<"b2n["<<i<<"] = "<<b2n[i]<<endl;
    hout <<"bn["<<i<<"] = "<<bn[i]<<endl;

    cp[i] = deltap[i] - 1;
    cn[i] = deltan[i] - 1;
    hout <<"cp["<<i<<"] = "<<cp[i]<<endl;
    hout <<"cn["<<i<<"] = "<<cn[i]<<endl;  

    A1p[i] = (bp[i] + sqrt(bp[i]*bp[i] - 4*ap[i]*cp[i])) / (2*ap[i]);
    // A1n[i] = (bn[i] + sqrt(bn[i]*bn[i] - 4*an[i]*cn[i])) / (2*an[i]);

    hout <<"bp[i]*bp[i] - 4*ap[i]*cp[i])["<<i<<"] = "<<bp[i]*bp[i] - 4*ap[i]*cp[i]<<endl;
    hout <<"A1p["<<i<<"] = "<<A1p[i]<<endl;
    hout <<"bn[i]*bn[i] - 4*an[i]*cn[i])["<<i<<"] = "<<bn[i]*bn[i] - 4*an[i]*cn[i]<<endl;
    //    hout <<"A1n["<<i<<"] = "<<A1n[i]<<endl;

   
    //    rfactor[i] = ((UPN[i]+DPN[i]+UNN[i]+DNN[i])/(UPP[i]+DPP[i]+UNP[i]+DNP[i]))/GetAcceptance[i];

    /// err_rfactor[i] =sqrt( ((1/(UPP[i]+DPP[i]+UNP[i]+DNP[i]))*(1/(UPP[i]+DPP[i]+UNP[i]+DNP[i]))*(UPN[i]+DPN[i]+UNN[i]+DNN[i]))/(GetAcceptance[i]*GetAcceptance[i]) + (((UPN[i]+DPN[i]+UNN[i]+DNN[i])/((UPP[i]+DPP[i]+UNP[i]+DNP[i])*(UPP[i]+DPP[i]+UNP[i]+DNP[i])))*((UPN[i]+DPN[i]+UNN[i]+DNN[i])/((UPP[i]+DPP[i]+UNP[i]+DNP[i])*(UPP[i]+DPP[i]+UNP[i]+DNP[i])))*(UPP[i]+DPP[i]+UNP[i]+DNP[i]))/(GetAcceptance[i]*GetAcceptance[i]) - 2*(1/UPP[i])*(UPN[i]/(UPP[i]*UPP[i])*CovUP[i]) - 2*(1/DPP[i])*(DPN[i]/(DPP[i]*DPP[i])*CovDP[i]) - 2*(1/UNP[i])*(UNN[i]/(UNP[i]*UNP[i])*CovUN[i]) - 2*(1/DNP[i])*(DNN[i]/(DNP[i]*DNP[i])*CovDN[i]) );


    //    A1[i] = (1/(1-rfactor[i]))*(A1p[i] - rfactor[i]*A1n[i]);
    //hout <<"A1["<<i<<"] = "<<A1[i]<<endl;
    ///////////// for statistical error calculations

    d1p[i] = b1p[i]*b1p[i] - 4*a1p[i];
    d2p[i] = 2*b1p[i]*b2p[i] - 4*a1p[i] - 4*a2p[i];
    ep[i]  = b2p[i]*b2p[i] - 4*a2p[i];
    d1n[i] = b1n[i]*b1n[i] - 4*a1n[i];
    d2n[i] = 2*b1n[i]*b2n[i] - 4*a1n[i] - 4*a2n[i];
    en[i]  = b2n[i]*b2n[i] - 4*a2n[i];

    err_e1p[i] = (0.5* (1/ sqrt(deltap[i]*deltap[i]*d1p[i] - deltap[i]*d2p[i] + ep[i])) * (2*deltap[i]*d1p[i] - d2p[i]) + b1p[i]) * (2*deltap[i]*a1p[i] - 2*a2p[i]);
 
    err_e2p[i] = 2*a1p[i]*( sqrt(deltap[i]*deltap[i]*d1p[i] - deltap[i]*d2p[i] + ep[i]) + deltap[i]*b1p[i] - b2p[i]);

    err_e3p[i] = (2*deltap[i]*a1p[i] - 2*a2p[i])*(2*deltap[i]*a1p[i] - 2*a2p[i]);

    err_ep[i] = (err_e1p[i] - err_e2p[i]) / err_e3p[i];
    
    err_A1p[i] = sqrt(err_ep[i] * err_ep[i] * err_deltap[i]*err_deltap[i]); //A1+

    err_e1n[i] = (0.5* (1/ sqrt(deltan[i]*deltan[i]*d1n[i] - deltan[i]*d2n[i] + en[i])) * (2*deltan[i]*d1n[i] - d2n[i]) + b1n[i]) * (2*deltan[i]*a1n[i] - 2*a2n[i]);
 
    err_e2n[i] = 2*a1n[i]*( sqrt(deltan[i]*deltan[i]*d1n[i] - deltan[i]*d2n[i] + en[i]) + deltan[i]*b1n[i] - b2n[i]);

    err_e3n[i] = (2*deltan[i]*a1n[i] - 2*a2n[i])*(2*deltan[i]*a1n[i] - 2*a2n[i]);

    err_en[i] = (err_e1n[i] - err_e2n[i]) / err_e3n[i];
    
    //   err_A1n[i] = sqrt(err_en[i] * err_en[i] * err_deltan[i]*err_deltan[i]); //A1-

    if (!(EvUPP[i]==0 || EvUPN[i]==0 || EvDPP[i]==0 || EvDPN[i]==0 ||EvUNP[i]==0 || EvUNN[i]==0 || EvDNP[i]==0 || EvDNN[i]==0)) {
      CorrUP[i] = CovUP[i]/(sqrt(EvUPP[i])*sqrt(EvUPN[i]));
      CorrDP[i] = CovDP[i]/(sqrt(EvDPP[i])*sqrt(EvDPN[i]));
      CorrUN[i] = CovUN[i]/(sqrt(EvUNP[i])*sqrt(EvUNN[i]));
      CorrDN[i] = CovDN[i]/(sqrt(EvDNP[i])*sqrt(EvDNN[i]));   
      CorrAV[i] = 0.25*(CorrUP[i]+CorrDP[i]+CorrUN[i]+CorrDN[i]);
      err_CorrAV[i] = 0.01;
    }

    hout <<"errA1p["<<i<<"] = "<<err_A1p[i]<<endl;
    //    hout <<"errA1n["<<i<<"] = "<<err_A1n[i]<<endl;


    //    err_A1[i] = sqrt((1/(1-rfactor[i]))*(1/(1-rfactor[i]))*err_A1p[i]*err_A1p[i]+rfactor[i]*rfactor[i]*(1/(1-rfactor[i]))*(1/(1-rfactor[i]))*err_A1n[i]*err_A1n[i] + ((1/(1-rfactor[i]))*(A1p[i] - (2*rfactor[i]-1)*A1n[i]))*((1/(1-rfactor[i]))*(A1p[i] - (2*rfactor[i]-1)*A1n[i]))*err_rfactor[i]*err_rfactor[i] -2*(1/(1-rfactor[i]))*(1/(1-rfactor[i]))*rfactor[i]*(CovUP[i]/(sumUPP[i]*sumUPN[i])+CovDP[i]/(sumDPP[i]*sumDPN[i])+CovUN[i]/(sumUNP[i]*sumUNN[i])+CovDN[i]/(sumDNP[i]*sumDNN[i])));
    
    //    CovTerm[i] = -2*(1/(1-rfactor[i]))*(1/(1-rfactor[i]))*rfactor[i]*(CovUP[i]/(sumUPP[i]*sumUPN[i])+CovDP[i]/(sumDPP[i]*sumDPN[i])+CovUN[i]/(sumUNP[i]*sumUNN[i])+CovDN[i]/(sumDNP[i]*sumDNN[i]));

    hout <<"errA1["<<i<<"] = "<<err_A1[i]<<endl;
    hout <<" "<<endl;
    /// plots
	
  }//2nd order weighting method 
  for (Int_t i = 0; i < nbinxbj; i++){
    fprintf  (A1out, "CorrUP[%i]\t= %f \t CorrDP[%i]\t= %f \t CorrUN[%i]\t= %f \t CorrDN[%i]\t= %f \n", i, CorrUP[i], i, CorrDP[i], i, CorrUN[i], i, CorrDN[i]);
  }
  fprintf  (A1out,"\n");
  fprintf  (A1out,"\n");
  for (Int_t i = 0; i < nbinxbj; i++){
    fprintf  (A1out, "CorrAV[%i]\t= %f \t= err_CorrAV[%i]\t= %f\n", i, CorrAV[i], i, err_CorrAV[i]);
  }

  fprintf  (A1out,"\n");
  fprintf  (A1out,"\n");
  for (Int_t i = 0; i < nbinxbj; i++){
    fprintf  (A1out, "rfactor[%i]\t= %f \n", i, rfactor[i] );
  }

  fprintf  (A1out,"\n");
  fprintf  (A1out,"\n");
  for (Int_t i = 0; i < nbinxbj; i++){
    fprintf  (A1out, "CovTerm[%i]\t= %f \n", i, CovTerm[i]);
  }
  fprintf  (A1out,"\n");
  fprintf  (A1out,"\n");

  //   for (Int_t i = 0; i < nbinxbj; i++){
  //     fprintf  (A1out, "A1[%i]\t= %f errA1[%i]\t= %f \n", i, A1[i], i,err_A1[i] );
  //   }
  fprintf  (A1out,"\n");
  fprintf  (A1out,"\n");
  //   for (Int_t i = 0; i < nbinxbj; i++){
  //     if (err_A1[i]!=0){
  //       fprintf  (A1out, "CovTerm[%i]/errA1[%i]\t= %f \n", i, i, CovTerm[i]/err_A1[i]);
  //     }
  //   }
  fprintf  (A1out,"\n");
  fprintf  (A1out,"\n");
  fclose(A1out);
  fclose(hadout);

  //////////////////////////////////////  End of Asymmetry calculations //////////////////////////////////////////////

  for(int i=0; i<nbinxbj; i++ ) {
    histoUPP[i]->Delete();
    histoDPP[i]->Delete();
    histoUNN[i]->Delete();
    histoDNN[i]->Delete();
    histoUNP[i]->Delete();
    histoDNP[i]->Delete();
    histoUPN[i]->Delete();
    histoDPN[i]->Delete();
    histoBetaUPP[i]->Delete(); 
    histoBetaDPP[i]->Delete(); 
    histoBetaUNP[i]->Delete(); 
    histoBetaDNP[i]->Delete(); 
    histoBetaUPN[i]->Delete(); 
    histoBetaDPN[i]->Delete(); 
    histoBetaUNN[i]->Delete(); 
    histoBetaDNN[i]->Delete(); 
  } 


}