#include "TCut.h"
#include "TCanvas.h"
#include "TH1.h"
#include "TF1.h"
#include "TChain.h"
#include "TGraph.h"
#include <math.h>
Double_t  gfpara(TH1F* hnm,Float_t vdn,Float_t vup,Double_t &pk,Double_t &sig)
{
  Double_t gpara[3];
  TF1 *gf=new TF1("gf","gaus");
  gf->SetLineColor(2);
  if(vdn==0.&&vup==0.)
    {
      hnm->Fit("gf");
      gf->GetParameters(&gpara[0]);
      pk=gpara[1];
      sig=gpara[2];
    }
  else
    {
      hnm->Fit("gf","","",vdn,vup);
      gf->GetParameters(&gpara[0]);
      pk=gpara[1];
      sig=gpara[2];
    }
    //YOU MUST RETURN SOMETHING !!!
    return 0;
}
//--main program begin------
void s30crscal(){
TChain chain("ana_data");
chain.Add("S30m020.root");

Double_t pk,sig;
char cutexp[512];
char savefile[512];

TCanvas *myC1= new TCanvas("myC1","S30 ANA HIST",1);
//------ the cut of position-------
	Double_t dntx=-14; Double_t uptx=16;
	sprintf(cutexp,"target1x>=%10.4f && target1x<=%10.4f",dntx,uptx);
	TCut ct1x=cutexp;
	Double_t dnty=-3; Double_t upty=17;
	sprintf(cutexp,"target1y>=%10.4f && target1y<=%10.4f",dnty,upty);
	TCut ct1y=cutexp;
	Double_t dnax=-7; Double_t upax=5;
	sprintf(cutexp,"angle1x>=%10.4f && angle1x<=%10.4f",dnax,upax);
	TCut cax=cutexp;
	Double_t dnay=-5; Double_t upay=5;
	sprintf(cutexp,"angle1y>=%10.4f && angle1y<=%10.4f",dnay,upay);
	TCut cay=cutexp;
	TCut cposi=ct1x && ct1y && cax && cay;
//-------end of the cut of positon-----------
//-----bcpids----
myC1->Divide(3,2);
myC1->cd(1);
	TH1F *S30bcde= new TH1F("S30bcde","BC de S30 no cut",30,160,180);
	chain.Draw("t1si1e>>S30bcde",cposi);
	gfpara(S30bcde,163,174,pk,sig);
		Double_t bcderpk=pk; Double_t bcdersig=sig;
		Double_t bcderdn=bcderpk-2*bcdersig; Double_t bcderup=bcderpk+2*bcdersig;
	sprintf(cutexp,"t1si1e > %10.4f && t1si1e < %10.4f",bcderdn,bcderup);
	TCut Cbcder=cutexp;
myC1->cd(2);
	TH1F *S30bcaoq= new TH1F("S30bcaoq","BC AoQ S30 no cut",50,1.75,2.05);
	chain.Draw("aoq1>>S30bcaoq",cposi);
	gfpara(S30bcaoq,1.83,1.94,pk,sig);
		Double_t bcarpk=pk; Double_t bcarsig=sig;
		Double_t bcardn=bcarpk-2*bcarsig; Double_t bcarup=bcarpk+2*bcarsig;
	sprintf(cutexp,"aoq1 > %10.6f && aoq1 < %10.6f",bcardn,bcarup);
	TCut Cbcaoqr=cutexp; 
	TCut Cbcipidr= Cbcder && Cbcaoqr && cposi;
myC1->cd(3);
	TH1F *S30bcdeC= new TH1F("S30bcdeC","BC Z S30 de with raw cut",30,160,180);
	chain.Draw("t1si1e>>S30bcdeC",Cbcipidr);
	gfpara(S30bcdeC,0.,0.,pk,sig);
		 Double_t bcdepk=pk; Double_t bcdesig=sig;
		 Double_t bcdedn=bcdepk-1.5*bcdesig; Double_t bcdeup=bcdepk+1.5*bcdesig;
	sprintf(cutexp,"t1si1e > %10.6f && t1si1e < %10.6f",bcdedn,bcdeup);
	TCut Cbcdec=cutexp;
	TH1F *S30bcaoqC= new TH1F("S30bcaoqC","BC AoQ S30 with raw cut",50,1.75,2.05);
myC1->cd(4);
	chain.Draw("aoq1>>S30bcaoqC",Cbcipidr);
	gfpara(S30bcaoqC,0.,0.,pk,sig);
		Double_t bcaoqpk=pk; Double_t bcaoqsig=sig;
		Double_t bcadn=bcaoqpk-1.5*bcaoqsig; Double_t bcaup=bcaoqpk+1.5*bcaoqsig;
	sprintf(cutexp,"aoq1 > %10.6f && aoq1 < %10.6f",bcadn,bcaup);
	TCut Cbcaoqc=cutexp;
	TCut Cbcipid=Cbcaoqc && Cbcdec && cposi;	//bcipID: A RAW cut. Name: Cbcipid posicut add

	sprintf(cutexp,"(((t1si1e-%10.6f)/(1.6*%10.6f))**2+((aoq1-%10.6f)/(1.6*%10.6f))**2)<=1.0",bcdepk,bcdesig,bcaoqpk,bcaoqsig);
	TCut Cbcipidt=cutexp;	// the Ellipse cut. Name: Cbcipidt
	TCut bcpidposi=Cbcipidt && cposi;
//--f
	TCut Cbcipidtt=Cbcipidt&&cposi&&("tof2<400");	//tof2 overflow

myC1->Close();
//-------IO correction gate--------@t1Si2-------

TCanvas *myC=new TCanvas("myC","IO CORRECTION GATE",1);
	TH1F *iogatet1si2 = new TH1F("iogatet1si2","IOgate t1si2",300,80,160);
	chain.Draw("t1si2e>>iogatet1si2",bcpidposi);	
//	chain.Draw("t1si2e>>iogatet1si2");	
	gfpara(iogatet1si2,110.,140.,pk,sig);
	Double_t iogpk=pk;	Double_t iogsig=sig;
	Double_t iogdn=iogpk-2.8*iogsig;	Double_t iogup=iogpk+2.8*iogsig;
	sprintf(cutexp,"t1si2e>=%10.4f && t1si2e<=%10.4f",iogdn,iogup);
	TCut iogcut=cutexp;
	TH1F *iogatet1si2c = new TH1F("iogatet1si2c","IOgate t1si2",300,80,160);
	chain.Draw("t1si2e>>iogatet1si2c",bcpidposi&&iogcut,"same");
	iogatet1si2c->SetFillColor(3);
myC->SaveAs("/mnt/LDATA/s30/S30pic/crsection/iogatecut2.osig.eps");
//----The Tof2 cuts-----About the the 2 part of the ac dE and E cuts------
	TCut tof2ac1="tof2<299";
	TCut tof2ac2="tof2>299";
	TCut acipidtl=tof2ac1 && Cbcipidtt;	//the left part of tof2
	TCut acipidtr=tof2ac2 && Cbcipidt &&("tof2<400");	//the right part of tof2

//------AC Pid---C IN Case:iogcut not added--Left Part-----
myC->Close();
myC=new TCanvas("myC","AC dE+E PID",1);
//myC->Divide(2,2);
//myC->cd(1);
gPad->SetLogy();
	TH1F *acz2tl= new TH1F("acz2tl","ac2z L",100,14,18);
	chain.Draw("z2ac>>acz2tl",acipidtl);
	gfpara(acz2tl,15.6,16.3,pk,sig);
	Double_t z2tlpk=pk; Double_t z2tlsig=sig;
	Double_t z2tldn=z2tlpk-2*z2tlsig; Double_t z2tlup=z2tlpk+2*z2tlsig;
	sprintf(cutexp,"z2ac>=%10.6f && z2ac<=%10.6f",z2tldn,z2tlup);
	TCut zactl=cutexp;

	TH1F *acz2tlsc= new TH1F("acz2tlsc","ac2z L-selfc",100,14,18);
	chain.Draw("z2ac>>acz2tlsc",acipidtl&&zactl,"same");
	acz2tlsc->SetFillColor(3);
//	chain.Draw("z2ac:tof2",acipidtl,"cont0");

//myC->cd(2);
//gPad->SetLogy();
	TH1F *eptli= new TH1F("eptli","EPlus ToF L",500,2700,3200);
	chain.Draw("(t2si2e+t2si3e+8.383121*tof2)>>eptli",acipidtl&&zactl);
	gfpara(eptli,3030,3090,pk,sig);
	Double_t eptlpk=pk; Double_t eptlsig=sig;
	Double_t eptldn=eptlpk-3*eptlsig; Double_t eptlup=eptlpk+3*eptlsig;
	sprintf(cutexp,"(t2si2e+t2si3e+8.383121*tof2)>=%10.6f && (t2si2e+t2si3e+8.383121*tof2)<=%10.6f",eptldn,eptlup);
	TCut eptl=cutexp;

	TH1F *ept2lsc= new TH1F("ept2lsc","EPlus ToF L-selfc",500,2700,3200);
	chain.Draw("(t2si2e+t2si3e+8.383121*tof2)>>ept2lsc",acipidtl&&zactl&&eptl,"same");
	ept2lsc->SetFillColor(3);

//---------AC Pid Cuts-----including position cuts------
	TCut aciedetl=acipidtl && eptl && zactl;
	TCut acoedetl=aciedetl && iogcut;
//------AC Pid---C IN Case:iogcut not added--Right Part-----
	TH1F *acz2tr= new TH1F("acz2tr","ac2z R",80,14,18);
	chain.Draw("z2ac>>acz2tr",acipidtr);
	gfpara(acz2tr,15.6,16.5,pk,sig);
	Double_t z2trpk=pk; Double_t z2trsig=sig;
	Double_t z2trdn=z2trpk-2*z2trsig; Double_t z2trup=z2trpk+2*z2trsig;
	sprintf(cutexp,"z2ac>=%10.6f && z2ac<=%10.6f",z2trdn,z2trup);
	TCut zactr=cutexp;



	TH1F *acz2trsc= new TH1F("acz2trsc","ac2z R-selfc",80,14,18);
	chain.Draw("z2ac>>acz2trsc",acipidtr&&zactr,"same");
	acz2trsc->SetFillColor(3);
myC->Close();
myC=new TCanvas("myC","AC dE+E PID",1);

	TH1F *eptr= new TH1F("eptr","EPlus ToF R",1000,2100,3200);
	chain.Draw("(t2si2e+t2si3e+6.89373*tof2)",acipidtr&&zactr);
//break;
	gfpara(eptr,3030,3090,pk,sig);
	Double_t eptrpk=pk; Double_t eptrsig=sig;
	Double_t eptrdn=eptrpk-3*eptrsig; Double_t eptrup=eptrpk+3*eptrsig;
	sprintf(cutexp,"(t2si2e+t2si3e+6.89373*tof2)>=%10.6f && (t2si2e+t2si3e+6.89373*tof2)<=%10.6f",eptrdn,eptrup);
	TCut eptrc=cutexp;	
//break;
	TH1F *ept2rsc= new TH1F("ept2rsc","EPlus ToF R-selfc",500,2700,3200);
	chain.Draw("(t2si2e+t2si3e+6.89373*tof2)>>ept2rsc",acipidtr&&zactr&&eptrc,"same");
	ept2rsc->SetFillColor(3);

//	chain.Draw("z2ac:tof2",acipidtl,"cont0");
//----------the crs dependence of different cuts.----------
	cout<<"begin the crs calculations"<<endl;
	cout<<"to calculate: 1. acpid dependece 2. iogate cut dependence 3. bcpid dependence"<<endl;
	Int_t choice;
//	cin>>choice>>endl;
	choice=2;
	if(choice==1) {		//the acpid dependence

	  Double_t cwidth[50];
	  Double_t inicount[50],inocount[50],outcicount[50],outcocount[50],crs[50];
	  Double_t thick,cs;
	  char hsiname[50],hsititle[50],hsoname[50],hsotitle[50];
//TH1F *acicountx[50], *acocountx[50];
	  for(int i=0;i<50;i++){
	   cwidth[i]=0.;inicount[i]=0.;inocount[i]=0.;outcicount[i]=0.;outcocount[i]=0.;
	   crs[i]=0.;
	  }
	for(int i=0;i<50;i++){
	cs=1.3+(Double_t)i/20;

	cwidth[i]=cs;
	  Double_t z2tldn=z2tlpk-cs*z2tlsig; Double_t z2tlup=z2tlpk+cs*z2tlsig;
	  sprintf(cutexp,"z2ac>=%10.6f && z2ac<=%10.6f",z2tldn,z2tlup);
	TCut stacz=cutexp;
	  Double_t eptldn=eptlpk-cs*eptlsig; Double_t eptlup=eptlpk+cs*eptlsig;
	  sprintf(cutexp,"(t2si2e+t2si3e+8.383121*tof2)>=%10.6f && (t2si2e+t2si3e+8.383121*tof2)<=%10.6f",eptldn,eptlup);
	TCut stacE=cutexp;

	TCut bci=bcpidposi;
	TCut acidep=acipidtl && stacE && stacz;	

	TCut bco=iogcut && bcpidposi;
	TCut acodep=iogcut && acidep;
//	sprintf(hsiname,"iaccw%4.2f",cs);
//	sprintf(hsoname,"oaccw%4.2f",cs);
//bc-count histograms
	sprintf(hsititle,"BC in case count x@target");
	sprintf(hsotitle,"BC out case x@target");
	TH1F *bcin= new TH1F("bcin",hsititle,25,-25,25);
	chain.Draw("target1x>>bcin",bci);
	inicount[i]=bcin->GetEntries();

	TH1F *bcout= new TH1F("bcout",hsotitle,25,-25,25);
	chain.Draw("target1x>>bcout",bco);
	inocount[i]=bcout->GetEntries();
//ac-count histograms
	sprintf(hsititle,"C in acpid cut width%4.2f x@target",cs);
	sprintf(hsotitle,"C out acpid cut width%4.2f x@target",cs);

	TH1F *acicountx= new TH1F("acicountx",hsititle,25,-25,25);
	chain.Draw("target1x>>acicountx",acidep);
	outcicount[i]=acicountx->GetEntries();

	TH1F *acocountx= new TH1F("acocountx",hsotitle,25,-25,25);
	chain.Draw("target1x>>acocountx",acodep);
	outcocount[i]=acocountx->GetEntries();

	crs[i]=-120110000.0*(log((outcicount[i]/inicount[i])/(outcocount[i]/inocount[i])))/(643.0*6.0221);

cout<<"cs="<<cwidth[i]<<"outci"<<outcicount[i]<<"outco"<<outcocount[i]<<"cross section "<<crs[i]<<endl;
		}	//acpid dependence
//---save the results
myC->Close();
TCanvas *myC=new TCanvas("myC","IO_2.0sig ACPID dep Counts",1);
	TGraph *crsacdep=new TGraph(51,cwidth,crs);
	crsacdep->SetMarkerStyle(20);
	crsacdep->SetMarkerColor(2);
	crsacdep->Draw("AP");
//sprintf(savefile,"/mnt/LDATA/s30/S30pic/crsection/acpid_dep_iog2.0.C");
//myC->SaveAS("/mnt/LDATA/s30/S30pic/crsection/acpid_dep_iog2.0.C");
	FILE *fp;
	fp=fopen("countacdep.dat","w");
	for(Int_t i=0;i<50;i++){	
	  fprintf(fp,"%12.4f%12.1f%12.1f%12.1f%12.1f%12.3f\n",cwidth[i],inicount[i],outcicount[i],inocount[i],outcocount[i],crs[i]);
			}
	fclose(fp);


	}	//choice 1 cal. end
	else if(choice==2) {		//the iog cut dependence

myC->Close();
myC=new TCanvas("myC","IO dep AC 3.0sig Counts",1);
	  char hsiname[50],hsititle[50],hsoname[50],hsotitle[50];
	Double_t ioinci[40],iooutci[40],ioinco[40],iooutco[40],iocrs[40];
	Double_t iocs,ciowidth[40];
	for(Int_t i=0;i<40;i++){ioinci[i]=0;iooutci[i]=0;ioinco[i]=0,iooutco[i]=0;iocrs[i]=0;ciowidth[i]=0;}
	for(Int_t i=0;i<40;i++){
	  iocs=1.2+(Double_t)i/20.0;
	  ciowidth[i]=iocs;
	  Double_t iogcsdn=iogpk-iocs*iogsig;	Double_t iogcsup=iogpk+iocs*iogsig;
	  sprintf(cutexp,"t1si2e>=%10.4f && t1si2e<=%10.4f",iogcsdn,iogcsup);
	TCut iogcs=cutexp;
	TCut iobcci=bcpidposi;
	TCut iobcco=bcpidposi && iogcs;
	TCut iocsaci=iobcci && acipidtl && eptl && zactl;
	TCut iocsaco=iobcco && iocsaci;
//-----io dep cin case------------
	sprintf(hsititle,"BC in IO width%4.2f x@target",iocs);
	sprintf(hsotitle,"BC out IO width%4.2f x@target",iocs);
	TH1F *iobci= new TH1F("iobci",hsititle,25,-25,25);
	chain.Draw("target1x>>iobci",iobcci);
	ioinci[i]=iobci->GetEntries();
	TH1F *iobco= new TH1F("iobco",hsotitle,25,-25,25);
	chain.Draw("target1x>>iobco",iobcco);
	ioinco[i]=iobco->GetEntries();
//-----io dep co case----
	sprintf(hsititle,"C in acpid IO width%4.2f x@target",iocs);
	sprintf(hsotitle,"C out acpid IO width%4.2f x@target",iocs);
	TH1F *iogcictx= new TH1F("iogcictx",hsititle,25,-25,25);
	chain.Draw("target1x>>iogcictx",iocsaci);
	iooutci[i]=iogcictx->GetEntries();
	TH1F *iogcoctx= new TH1F("iogcoctx",hsotitle,25,-25,25);
	chain.Draw("target1x>>iogcoctx",iocsaco);
	iooutco[i]=iogcoctx->GetEntries();

	iocrs[i]=-120110000.0*(log((iooutci[i]/ioinci[i])/(iooutco[i]/ioinco[i])))/(643.0*6.0221);

	}
	FILE *fp;
	fp=fopen("iogacz2e3sigdep.dat","w");
	for(Int_t i=0;i<40;i++){	
	fprintf(fp,"%12.4f%12.1f%12.1f%12.1f%12.1f%12.3f\n",ciowidth[i],ioinci[i],iooutci[i],ioinco[i],iooutco[i],iocrs[i]);
	}
	fclose(fp);

myC->Close();
TCanvas *myC=new TCanvas("myC","crs IO dep ACPID sig3.0",1);
	TGraph *ratioio=new TGraph(41,ciowidth,iocrs);
	ratioio->SetMarkerStyle(20);
	ratioio->SetMarkerColor(2);
	ratioio->Draw("AP");


	}	//choice 2 cal. end
	else if(choice==3) {		//the ****dependence


	}	//choice 3 cal. end

	else{}





}//programm main