#include <stdio.h>

double f1(int id,double t,double *x[],double *v[]);
double f2(int id,double t,double *x[],double *v[],double *B[]);

int imax=600E3;double endKickPar=0.5;
int typ=12;double facB=1;double TK=260;double start_t0=0E-9;//weak_0.48 mod-3 redo 20190624 CONFIRM kick+17
double vin[3]={0,0,0};//pitch_-0.8mrad 
double x0=33.3*1E-2;//
double z0=0.*1E-2;//m
double tau=1E20;double vy0=16.299E-3;//hiromi20170214
double pi=TMath::Pi();

double f1(int id,double t,double x[3],double v[3])
{
        return v[id];
}

double f2(int id,double t,double x[3],double v[3],double B[3])
{
  double ans[3];
 ans[0]=v[1]*B[2]-v[2]*B[1];
 ans[1]=v[2]*B[0]-v[0]*B[2];
 ans[2]=v[0]*B[1]-v[1]*B[0];
        return ans[id];
}
//////////////KICK-TOOL
int kickMag(double t,double x[3],double BRkick0,double BYkick0,double &BRkick,double &BYkick){
 BRkick0=BRkick0*1E-4;
 BYkick0=BYkick0*1E-4;
    double theta=t/(double)(TK*1E-9)*2*pi;
    BRkick=BRkick0*sin(theta);
    BYkick=BYkick0*sin(theta);
    return 1;
}

void UniBkick()//Feb2018_*.txt, Gauss
{
//double y0=0.6;//
double y0=0.2758;//20200427
  int home=0;
  int w_frag=1;//write all trajectory 
  FILE *ofile,*infile;
  char oname[256];
  char iname[256];
  char sdmy[256];
  double ldmy;
  int idmy;
      sprintf(oname,"kickStop_UniBKick_hs.txt"); 
      ofile=fopen(oname,"w"); 
  double xtmp,ytmp,ztmp,vxtmp,vytmp,vztmp;
  double vxtmp0,vytmp0,vztmp0,vtmp0;
  int ttmp;
  double timeStop;
        double t=0;
        double x[3],v[3],dt;
        double sp[3],dsp[3],mom[3],betaSpin[3],betaB;
        double k1[6][2],k2[6][2],k3[6][2],k4[6][2];
        double B[3]={0,3,0};
        double mom0=300;//MeV/c
        double beta_mu;
        double mu_mass=105.66;//MeV
        double Etot=sqrt(pow(mu_mass,2)+pow(mom0,2));
        double gamma_mu=Etot/(double)mu_mass;
        double light_c=TMath::C();
        double q=1.60217733*1E-19;//[c]
        double GeVkG=1E-26/(double)5.61;
        beta_mu=sqrt(pow(gamma_mu,2)-1)/(double)gamma_mu;
        double v0=light_c*beta_mu;//[m/sec]
        printf("beta_mu=%e gamma_mu=%e\n",beta_mu,gamma_mu);//getchar();

        double par1;
        par1=mu_mass*1E-3*GeVkG*gamma_mu;
    //par1=e_mass*1E-3*GeVkG*gamma_e;
        double amu=0.0011659208;
        double par2=q/(double)par1;
        double par3=q/(double)(mu_mass*1E-3*GeVkG)*amu*gamma_mu/(double)(gamma_mu+1);

        double R0=v0/(double)(B[1]*par2);
        printf("par2=%lf R0=%e\n",par2,R0);//getchar();
        double pitch=-15E-3;
        double tmp=1/(double)(2*par2)*pitch*2*pi/(double)TK*1E9;
        printf("A=1/2*(par1)=%e\n",1/(double)(2*par1));
        printf("kick omega=%e\n",2*pi/(double)TK*1E9);
        printf("tmp=%e\n",tmp);getchar();

        x[0]=x0;                  //
        x[1]=y0;                  //
        x[2]=z0;                  //
        v[0]=0;                  //
        v[1]=v0*pitch;                  //
        v[2]=v0*sqrt(1-pow(pitch,2));                  //

        //dt=1E-9*1E-1;                //
        //dt=1E-9*1E-2;                //
        //dt=5E-9*1E-3;                //
        dt=1E-9*1E-3;                //

/////////////set kickB-field/////////
   double start_kick;
   start_kick=start_t0;//sec
   double end_kick=endKickPar*TK*1E-9;//sec
   double BRkick,BYkick;
   double BRkick0,BYkick0;
   double kickTime;
   double int_time;
   int kick_status;

        double xx[128*100],zz[128*100],yy[128*100],px[128*100],py[128*100],pz[128*100],p0;
        double vy[128*100];
        double vv;
        double tg[128*100],brg_kick[128*100],brg[128*100],brg_weak[128*100];
        double tg2[128*100];
        double dsx[128*100],dsy[128*100],dsz[128*100];
        double sx[128*100],sy[128*100],sz[128*100],tt[128*100];
        double spinmom[128*100];
        int ig2=0;
        double sp1[3],sp2[3],dsp1[3],dsp2[3],sp3[3],dsp3[3];
        int ig=0;
        double x1[3],x2[3],v1[3],v2[3],x3[3],v3[3];
for(int ii=0;ii<imax;ii++) {
           double RM=sqrt(x[0]*x[0]+x[2]*x[2]);
           double YM=x[1];
           if(ii==0){
                 double V=sqrt(pow(v[0],2)+pow(v[1],2)+pow(v[2],2));
                 sp[0]=-1*v[0]/(double)V;
                 sp[1]=-1*v[1]/(double)V;
                 sp[2]=-1*v[2]/(double)V;
                   for(int iip=0;iip<3;++iip)betaSpin[iip]=v[iip]/(double)light_c;
                   betaB=betaSpin[0]*B[0]+betaSpin[1]*B[1]+betaSpin[2]*B[2];
                 dsp[0]=dt*f2(0,t,sp,sp,B)*par2*(1+amu*gamma_mu)-dt*betaB*par3*f2(0,t,sp,sp,betaSpin);
                 dsp[1]=dt*f2(1,t,sp,sp,B)*par2*(1+amu*gamma_mu)-dt*betaB*par3*f2(1,t,sp,sp,betaSpin);
                 dsp[2]=dt*f2(2,t,sp,sp,B)*par2*(1+amu*gamma_mu)-dt*betaB*par3*f2(2,t,sp,sp,betaSpin);
                 printf("initial x=%e %e %e\n",x[0],x[1],x[2]);//getchar();
            }
/////////////Get kick-B-spatial/////////
            BRkick=0;BYkick=0;
            BRkick0=0;;
            BYkick0=0;;
            int kick_spatial=0;
            //if(sqrt(pow(x[1],2))<0.6){
            if(sqrt(pow(x[1],2))<y0){
              //BRkick0=6;
              BRkick0=sqrt(pow(tmp,2))*1E4;//Gauss
              BYkick0=0;
            }
///////////Get kickB-field/////////
           kick_status=0;
           if(t>=start_kick && t<(start_kick+end_kick)){
              int_time=t-start_kick;
              kick_status=kickMag(int_time,x,BRkick0,BYkick0,BRkick,BYkick);
           }
           if( t>=(start_kick+end_kick))kick_status=2;
/////////////Get B-field/////////
           double BR=0;double BY=3;
           B[0]=(BR+BRkick)*x[0]/(double)RM;
           B[2]=(BR+BRkick)*x[2]/(double)RM;
           B[1]=BY+BYkick;
////////////Get Bofield END////////
                k1[0][0]=dt*f1(0,t,x,v);
                k1[0][1]=dt*f2(0,t,x,v,B)*par2;
                k1[1][0]=dt*f1(1,t,x,v);
                k1[1][1]=dt*f2(1,t,x,v,B)*par2;
                k1[2][0]=dt*f1(2,t,x,v);
                k1[2][1]=dt*f2(2,t,x,v,B)*par2;

                for(int i=0;i<3;++i)x1[i]=x[i]+k1[i][0]/2.0;
                for(int i=0;i<3;++i)v1[i]=v[i]+k1[i][1]/2.0;
                for(int i=0;i<3;++i)k2[i][0]=dt*f1(i,t+dt/2.0,x1,v1);
                for(int i=0;i<3;++i)k2[i][1]=dt*f2(i,t+dt/2.0,x1,v1,B)*par2;

                for(int i=0;i<3;++i)x2[i]=x[i]+k2[i][0]/2.0;
                for(int i=0;i<3;++i)v2[i]=v[i]+k2[i][1]/2.0;
                for(int i=0;i<3;++i)k3[i][0]=dt*f1(i,t+dt/2.0,x2,v2);
                for(int i=0;i<3;++i)k3[i][1]=dt*f2(i,t+dt/2.0,x2,v2,B)*par2;

                for(int i=0;i<3;++i)x3[i]=x[i]+k3[i][0];
                for(int i=0;i<3;++i)v3[i]=v[i]+k3[i][1];
                for(int i=0;i<3;++i)k4[i][0]=dt*f1(i,t+dt,x3,v3);
                for(int i=0;i<3;++i)k4[i][1]=dt*f2(i,t+dt,x3,v3,B)*par2;

                for(int i=0;i<3;++i)x[i]=x[i]+(k1[i][0]+2.0*k2[i][0]+2.0*k3[i][0]+k4[i][0])/6.0;
                for(int i=0;i<3;++i)v[i]=v[i]+(k1[i][1]+2.0*k2[i][1]+2.0*k3[i][1]+k4[i][1])/6.0;

/////spin
                for(int iip=0;iip<3;++iip)betaSpin[iip]=v[iip]/(double)light_c;
                betaB=betaSpin[0]*B[0]+betaSpin[1]*B[1]+betaSpin[2]*B[2];
                k1[3][1]=dt*f2(0,t,sp,sp,B)*par2*(1+amu*gamma_mu)-dt*betaB*par3*f2(0,t,sp,sp,betaSpin);
                k1[4][1]=dt*f2(1,t,sp,sp,B)*par2*(1+amu*gamma_mu)-dt*betaB*par3*f2(1,t,sp,sp,betaSpin);
                k1[5][1]=dt*f2(2,t,sp,sp,B)*par2*(1+amu*gamma_mu)-dt*betaB*par3*f2(2,t,sp,sp,betaSpin);
                  //printf("t=%lf y=%lf sp[0]=%lf sp[1]=%lf sp[2]=%lf\n",t,x[1],sp[0],sp[1],sp[2]); 
                  //printf("B[0]=%lf B[1]=%lf B[2]=%lf\n",B[0],B[1],B[2]); 
                  //printf("par2*(1+amu*gamma_mu)=%lf\n",par2*(1+amu*gamma_mu)); 
                //  printf("t=%e k1[3][1]=%lf k1[4][1]=%lf k1[5][1]=%lf\n",t,k1[3][1],k1[4][1],k1[5][1]); 
                //  getchar();
                for(int iip=0;iip<3;++iip)sp1[iip]=sp[iip]+k1[iip+3][1]/2.0;
                for(int iip=0;iip<3;++iip){
                  k2[iip+3][1]=dt*f2(iip,t+dt/2.0,sp1,sp1,B)*par2*(1+amu*gamma_mu)-dt*betaB*par3*f2(iip,t,sp1,sp1,betaSpin);
                  //printf(" sp1[%d]=%e k2[%d][1]=%e \n",iip,sp1[iip],iip+3,k2[iip+3][1]); 
                }
              //  getchar();

                for(int iip=0;iip<3;++iip)sp2[iip]=sp[iip]+k2[iip+3][1]/2.0;
                for(int iip=0;iip<3;++iip){
                   k3[iip+3][1]=dt*f2(iip,t+dt/2.0,sp2,sp2,B)*par2*(1+amu*gamma_mu)-dt*betaB*par3*f2(iip,t,sp2,sp2,betaSpin);
         //         printf("sp2[%d]=%e k3[%d][1]=%e \n",iip,sp2[iip],iip+3,k3[iip+3][1]); 
                }
              //  getchar();

                for(int iip=0;iip<3;++iip)sp3[iip]=sp[iip]+k3[iip+3][1];
                for(int iip=0;iip<3;++iip){
                   k4[iip+3][1]=dt*f2(iip,t+dt,sp3,sp3,B)*par2*(1+amu*gamma_mu)-dt*betaB*par3*f2(iip,t,sp3,sp3,betaSpin);
        //          printf("sp3[%d]=%e k4[%d][1]=%e \n",iip,sp3[iip],iip+3,k4[iip+3][1]); 
                }
              //  getchar();

                for(int iip=0;iip<3;++iip){
                  sp[iip]=sp[iip]+(k1[iip+3][1]+2.0*k2[iip+3][1]+2.0*k3[iip+3][1]+k4[iip+3][1])/6.0;
                  //printf("sp[%d]=%e\n",iip,sp[iip]); 
                }
               // printf("sp=%e\n",sqrt(sp[0]*sp[0]+sp[1]*sp[1]+sp[2]*sp[2]));
              //  getchar(); 

//----------------spin end -----------------
               if( ii%100==0 && ig<4*2000){  //for half sin
                  xx[ig]=x[0];
                  zz[ig]=x[2];
                  yy[ig]=x[1];
                  vv=sqrt(pow(v[0],2)+pow(v[1],2)+pow(v[2],2));
                  vy[ig]=v[1]/(double)vv;
                  tg[ig]=t;
                  px[ig]=(v[0]/vv);
                  py[ig]=(v[1]/vv);
                  pz[ig]=(v[2]/vv);
                  sx[ig]=sp[0];
                  sy[ig]=sp[1];
                  sz[ig]=sp[2];
                  spinmom[ig]=(sx[ig]*px[ig]+sy[ig]*py[ig]+sz[ig]*pz[ig])*1;
                  brg[ig]=(BR+BRkick)*2E1;
                  brg_kick[ig]=(BRkick)*2E1;
                  brg_weak[ig]=(BR)*2E1;
                  //printf("now is ii=%d ig=%d x=%lf %lf %lf B=%e %e %e::vy=%e kick=%d\n",
                        //ii,ig,x[0],x[1],x[2],B[0],B[1],B[2],vy[ig],kick_status);
                  if( x[1]>0.95){
                        printf("pos y too big %lf\n",x[1]);
                         break;
                  }
                  ig++;
                  if(ii>=0)fprintf(ofile,"%e,%e,%e,%e,%e,%e,%e\n",t,x[0],x[1],x[2],v[0],v[1],v[2]);
               }
               t=t+dt;
}//END(for imax)


if( w_frag==1 )fclose(ofile);
///////Graph
gROOT->SetStyle("Plain");
gROOT->ForceStyle();
gStyle->SetOptStat(0);
gStyle->SetCanvasBorderMode(0);
gStyle->SetPadColor(0);
gStyle->SetTitleFillColor(0);
gStyle->SetStatColor(0);
gStyle->SetCanvasColor(0);
gStyle->SetFrameLineColor(1);

TGraph *gr_xz,*gr_xy,*gr_yy,*gr_vy,*gr_BR,*gr_BRkick,*gr_BRall;
TGraph *gr_tpx,*gr_tpy,*gr_tpz,*gr_tsx,*gr_tsy,*gr_tsz,*gr_tsp;
 gr_xz = new TGraph(ig,xx,zz);
 gr_xy = new TGraph(ig,xx,yy);//??
 gr_yy = new TGraph(ig,yy,vy);
 gr_vy = new TGraph(ig,tg,vy);
 gr_BR = new TGraph(ig,tg,brg_weak);
 gr_BRkick = new TGraph(ig,tg,brg_kick);
 gr_BRall = new TGraph(ig,tg,brg);

 gr_tpx = new TGraph(ig,tg,px);
 gr_tpy = new TGraph(ig,tg,py);
 gr_tpz = new TGraph(ig,tg,pz);
 gr_tsx = new TGraph(ig,tg,sx);
 gr_tsy = new TGraph(ig,tg,sy);
 gr_tsz = new TGraph(ig,tg,sz);
 gr_tsp = new TGraph(ig,tg,spinmom);

 TCanvas *c = new TCanvas("c","check",200,100,500,500);
 c->SetLeftMargin(0.13);
 TH2F* frame = new TH2F("frame"," ",10, -0.4, 0.4, 100, -0.4, 0.4);
 frame->SetTitle("");
 frame->GetXaxis()->SetLabelSize(0.04);
 frame->GetXaxis()->SetTitleSize(0.05);
 frame->GetYaxis()->SetLabelSize(0.04);
 frame->GetYaxis()->SetTitleSize(0.05);
 frame->GetXaxis()->SetTitle("posx[m]");
 frame->GetYaxis()->SetTitle("posz[m]");
 frame->GetYaxis()->SetTitleOffset(1.2);
 frame->SetStats(0);
 frame->Draw();
   gr_xz->SetMarkerStyle(9);
   gr_xz->SetMarkerColor(1);
   gr_xz->SetLineWidth(2);
   gr_xz->Draw("pl");
 c->SetGrid();

 TCanvas *cxy = new TCanvas("cxy","check x-y",200,100,500,500);
 cxy->SetLeftMargin(0.13);
 TH2F* framexy = new TH2F("framexy"," ",10, -0.4, 0.4, 100, -0.1, 0.4);
 framexy->SetTitle("");
 framexy->GetXaxis()->SetLabelSize(0.04);
 framexy->GetXaxis()->SetTitleSize(0.05);
 framexy->GetYaxis()->SetLabelSize(0.04);
 framexy->GetYaxis()->SetTitleSize(0.05);
 framexy->GetXaxis()->SetTitle("posx[m]");
 framexy->GetYaxis()->SetTitle("posy[m]");
 framexy->GetYaxis()->SetTitleOffset(1.2);
 framexy->SetStats(0);
 framexy->Draw();
   gr_xy->SetMarkerStyle(9);
   gr_xy->SetMarkerColor(1);
   gr_xy->SetLineWidth(3);
   gr_xy->Draw("pl");
 cxy->SetGrid();

 TCanvas *cyy = new TCanvas("cyy","check y-vy",200,100,500,500);
 cyy->SetLeftMargin(0.15);
// TH2F* frameyy = new TH2F("frameyy"," ",10, -1.4, 1.4, 100, -0.1, 0.1);
 //TH2F* frameyy = new TH2F("frameyy"," ",10, -0.5, 1.0, 100, -0.3, 0.01);
 TH2F* frameyy = new TH2F("frameyy"," ",10, -0.5, 1.0, 100, -0.3, 0.3);
// TH2F* frameyy = new TH2F("frameyy"," ",10, -2E-2, 2E-2, 100, -1E-3, 1E-3);
// TH2F* frameyy = new TH2F("frameyy"," ",10, -5E-2, 5E-2, 100, -5E-3, 5E-3);
 frameyy->SetTitle("");
 frameyy->GetXaxis()->SetLabelSize(0.04);
 frameyy->GetXaxis()->SetTitleSize(0.05);
 frameyy->GetYaxis()->SetLabelSize(0.04);
 frameyy->GetYaxis()->SetTitleSize(0.05);
 frameyy->GetXaxis()->SetTitle("posy[m]");
 frameyy->GetYaxis()->SetTitle("pitch=vy/v[rad]");
 frameyy->GetYaxis()->SetTitleOffset(1.3);
 frameyy->SetStats(0);
 frameyy->Draw();
   gr_yy->SetMarkerStyle(9);
   gr_yy->SetMarkerColor(1);
   gr_yy->SetLineWidth(3);
   gr_yy->Draw("pl");
 cyy->SetGrid();

 TCanvas *cte = new TCanvas("ct","check time",200,100,500,500);
 cte->SetLeftMargin(0.2);
 cte->SetRightMargin(0.2);
 //TH2F* framet = new TH2F("framet"," ",10, 0., 1.1E-6, 100, -0.01, 0.01);
 TH2F* framet = new TH2F("framet"," ",10, 0., 1.1E-6, 100, -0.1, 0.1);
 framet->SetTitle("");
 framet->GetXaxis()->SetLabelSize(0.04);
 framet->GetXaxis()->SetTitleSize(0.05);
 framet->GetYaxis()->SetLabelSize(0.04);
 framet->GetYaxis()->SetTitleSize(0.05);
 framet->GetXaxis()->SetTitle("[nsec]");
 framet->GetYaxis()->SetTitle("vy/v[rad], b");
 framet->GetYaxis()->SetTitleOffset(1.3);
 framet->SetStats(0);
 framet->Draw();

 TF1 *f1=new TF1("f1","x",-0.01*0.05*1E4,0.01*0.05*1E4);
   TGaxis *A1 = new TGaxis(1.1E-6,-0.01,1.1E-6,0.01,"f1",4,"-");
   A1->SetTitle("B gauss");
   A1->SetTitleOffset(-0.26);
   A1->SetLabelOffset(-0.12);
   A1->SetLabelColor(2);
   A1->SetTitleColor(2);
   A1->Draw();
   gr_vy->SetMarkerStyle(9);
   gr_vy->SetMarkerColor(1);
   gr_vy->SetLineWidth(1);
   gr_vy->SetLineColor(1);
   gr_vy->Draw("pl");

   gr_BR->SetMarkerStyle(9);
   gr_BR->SetMarkerColor(2);
   gr_BR->SetLineWidth(1);
   gr_BR->SetLineColor(2);
   gr_BR->Draw("pl");

   gr_BRkick->SetMarkerStyle(9);
   gr_BRkick->SetMarkerColor(4);
   gr_BRkick->SetLineWidth(1);
   gr_BRkick->SetLineColor(4);
   gr_BRkick->Draw("pl");

   gr_BRall->SetMarkerStyle(9);
   gr_BRall->SetMarkerColor(3);
   gr_BRall->SetLineWidth(1);
   gr_BRall->SetLineColor(3);
   gr_BRall->Draw("pl");
 cte->SetGrid();

 TCanvas *cSpin = new TCanvas("cSpin","check spin time",200,100,900,500);
 cSpin->SetLeftMargin(0.1);
 cSpin->SetRightMargin(0.1);
 //TH2F* frameS = new TH2F("frameS"," ",10, 0., 1.1E-7, 100, -1, 1);
 //TH2F* frameS = new TH2F("frameS"," ",10, 0., 3E-7, 100, -1, 1);
 TH2F* frameS = new TH2F("frameS"," ",10, 0., 1.2E-6, 100, -1, 1);
 frameS->SetTitle("");
 frameS->GetXaxis()->SetLabelSize(0.04);
 frameS->GetXaxis()->SetTitleSize(0.05);
 frameS->GetYaxis()->SetLabelSize(0.04);
 frameS->GetYaxis()->SetTitleSize(0.05);
 frameS->GetXaxis()->SetTitle("[sec]");
 frameS->GetYaxis()->SetTitle("spin, momentum");
 frameS->GetYaxis()->SetTitleOffset(0.9);
 frameS->SetStats(0);
 frameS->Draw();

   gr_tpx->SetMarkerStyle(9);
   gr_tpx->SetMarkerColor(7);
   gr_tpx->SetLineWidth(1);
   gr_tpx->SetLineColor(7);
//   gr_tpx->Draw("l");

   gr_tpy->SetMarkerStyle(9);
   gr_tpy->SetMarkerColor(8);
   gr_tpy->SetLineWidth(2);
   gr_tpy->SetLineColor(8);
  gr_tpy->Draw("pl");

   gr_tpz->SetMarkerStyle(9);
   gr_tpz->SetMarkerColor(kMagenta+2);
   gr_tpz->SetLineWidth(1);
   gr_tpz->SetLineColor(kMagenta+2);
//   gr_tpz->Draw("pl");

   gr_tsx->SetMarkerStyle(7);
   gr_tsx->SetMarkerColor(1);
   gr_tsx->SetLineWidth(1);
   gr_tsx->SetLineColor(6);
   gr_tsx->SetLineWidth(2);
//   gr_tsx->Draw("l");
   gr_tsy->SetMarkerStyle(7);
   gr_tsy->SetMarkerColor(7);
   gr_tsy->SetLineWidth(2);
   gr_tsy->SetLineColor(2);
   gr_tsy->Draw("l");
   gr_tsz->SetMarkerStyle(7);
   gr_tsz->SetMarkerColor(4);
   gr_tsz->SetLineWidth(2);
   gr_tsz->SetLineColor(4);
//   gr_tsz->Draw("l");
 //
   gr_tsp->SetMarkerStyle(7);
   gr_tsp->SetMarkerColor(1);
   gr_tsp->SetLineWidth(1);
   gr_tsp->SetLineColor(1);
   gr_tsp->Draw("pl");
 TLegend *l1;
 l1=new TLegend(0.75,0.2,0.99,0.3);
 l1->SetTextSize(0.05);
 l1->SetFillColor(0);
// l1->AddEntry(gr_tsx,"sx","l");
// l1->AddEntry(gr_tsz,"sy","l");
 l1->AddEntry(gr_tsy,"sz","l");
// l1->AddEntry(gr_tpx,"px","lp");
// l1->AddEntry(gr_tpz,"py","lp");
 l1->AddEntry(gr_tpy,"pz","lp");
 l1->AddEntry(gr_tsp,"spin*mom","lp");
 l1->Draw("same");
 cSpin->SetGrid();
}