Double_t fuPtVsY(Double_t *x, Double_t *par)
{
// function y=fu gives p_t vs lab. rapidity
// x[0] - rapidity in the lab.
// par[0]=Pcm, par[1]=mass, par[2]=Ycm
  Double_t Ecm,VL,PL,val;
  Ecm=TMath::Sqrt(par[0]*par[0]+par[1]*par[1]);
  VL=TMath::TanH(x[0]-par[2]);
  PL=VL*Ecm;
  val=0;
  if(TMath::Abs(PL)<par[0]) val=TMath::Sqrt(par[0]*par[0]-PL*PL);

  return val;
}
Double_t fuPvsTheta(Double_t *x, Double_t *par)
{
// par[0]=beta, par[1]=gamma, par[2]=Pcm, par[3]=mass
// y^2 + (x/gamma - beta*Ecm)^2 = Pcm, where y=p_t, x=p_z
// y=p*sin(theta), x=p*cos(theta)
// we want to have p(GeV/c) vs theta(deg.)
  Double_t Ecm=TMath::Sqrt(par[2]*par[2]+par[3]*par[3]);
//  Double_t w=x[0]/par[1]-par[0]*Ecm;
//  Double_t val=TMath::Sqrt(par[2]*par[2]-w*w);
  Double_t kDegRad=3.1415926/180.,a,b,c,det,sinThe,cosThe,val;
  sinThe=TMath::Sin(x[0]*kDegRad);
  cosThe=TMath::Cos(x[0]*kDegRad);
  a=par[1]*par[1]*sinThe*sinThe+cosThe*cosThe;
  b=-2.*par[1]*par[0]*cosThe*Ecm;
  c=par[1]*par[1]*(par[0]*par[0]*Ecm*Ecm-par[2]*par[2]);
  det=b*b-4.*a*c;
  if(det<0) {
    val=0;
  } else {
    val=(-b+TMath::Sqrt(det))/(2.*a);
    if(val<0) val=0;
  }
  return val;
}
void hiau25()
{
  gROOT.Reset();

  TStopwatch timer;
  timer.Start();

  Int_t PdgCode,geantCode,outOfGeantCode=0;
//                 GEANT numbers of particles
//      1. gamma       13. n0          25. n~0         37. D0
//      2. e+          14. p+          26. Lambda~0    38. D~0
//      3. e-          15. p~-         27. Sigma~-     39. D_s+
//      4. neutrino    16. K_S0        28. Sigma~0     40. D_s-
//      5. mu+         17. eta         29. Sigma~+     41. Lambda_c+
//      6. mu-         18. Lambda      30. Xi~0        42. rho+
//      7. pi0         19. Sigma+      31. Xi~+        43. rho-
//      8. pi+         20. Sigma0      32. Omega~+     44. rho0
//      9. pi-         21. Sigma-      33. omega       45. deut
//     10. K_L0        22. Xi0         34. phi         46. t
//     11. K+          23. Xi-         35. D+          47. He4
//     12. K-          24. Omega-      36. D-          48. Geantino
                                                             
  Int_t charge[48]={ 0, 1,-1, 0, 1,-1, 0, 1,-1, 0,
                     1,-1, 0, 1,-1, 0, 0, 0, 1, 0,
                    -1, 0,-1,-1, 0, 0,-1, 0, 1, 0,
                     1, 1, 0, 0, 1,-1, 0, 0, 1,-1,
                     1, 1,-1, 0, 1, 1, 2, 0};
  Int_t baryon[48]={ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                     0, 0, 1, 1,-1, 0, 0, 1, 1, 1,
                     1, 1, 1, 1,-1,-1,-1,-1,-1,-1,
                    -1,-1, 0, 0, 0, 0, 0, 0, 0, 0,
                     0, 0, 0, 0, 0, 0, 0, 0};
//  Float_t massGeant[48]={0.0,0.00051,0.00051,0.0,0.1057,0.1057,0.135,0.1396,0.1396,0.4977,
  Float_t massGeant[48]={0.0,0.00051,0.00051,0.0,0.1057,0.1057,0.138,0.138,0.138,0.498,
//			 0.4936,0.4936,0.9396,0.9383,0.9383,0.4977,0.5488,1.1156,1.1894,1.1926,
			 0.494,0.494,0.938,0.938,0.938,0.498,0.5488,1.1156,1.1894,1.1926,
			 1.1974,1.3149,1.3213,1.6724,0.9396,1.1156,1.1894,1.1926,1.1974,1.3149,
                         1.3213,1.6724,0.782,1.020,1.869,1.869,1.864,1.864,1.969,1.969,
                         2.285,0.769,0.769,0.769,0.0,0.0,0.0,0.0};

  Float_t px,py,pz,p,e,nEvents=50,pt,Charged=0;
  Float_t thedeg=0,yrap,theta,phi;
  nEvents=1000;
  printf(" %i events\n",nEvents);

// Create histograms
  TH1F *hID = new TH1F("hID","HIJING1.381: 10^{3} Au+Au central (b=0) events at E/A=25 GeV",50,0,50);


  Int_t nMom=80, nTheta=90;
  Float_t momMin=0,momMax=80;  
  Float_t thetaMin=0,thetaMax=90;  
  TH1F *hGammaP = new TH1F("hGammaP","HIJING, Au+Au 25 AGeV, primary #gamma",nMom,momMin,momMax);
  TH1F *hGammaTheta = new TH1F("hGammaTheta","HIJING, Au+Au 25 AGeV, primary #gamma",nTheta,thetaMin,thetaMax);
  TH2F *hGammaPvsTheta = new TH2F("hGammaPvsTheta","HIJING, Au+Au 25 AGeV, primary #gamma (1/event/bin)",nTheta,thetaMin,thetaMax,nMom,momMin,momMax);

  Int_t nPt=50, nY=50;
  Float_t ptMin=0,ptMax=10;  
  Float_t yMin=-3,yMax=7;  
  TH1F *hGammaPt = new TH1F("hGammaPt","HIJING, Au+Au 25 AGeV, primary #gamma",nPt,ptMin,ptMax);
  TH1F *hGammaY = new TH1F("hGammaY","HIJING, Au+Au 25 AGeV, primary #gamma",nY,yMin,yMax);
  TH2F *hGammaPtVsY = new TH2F("hGammaPtVsY","HIJING, Au+Au 25 AGeV, primary #gamma (1/event/bin)",nY,yMin,yMax,nPt,ptMin,ptMax);

  Float_t yCM=2.0;
  Float_t deltaYrap=1.0;
  Float_t deltaPt=0.25; // GeV/c
  Float_t pTmax=3.5;
  TH1F *hSigma = new TH1F("hSigma","",pTmax/deltaPt,0,pTmax);

  Float_t nTagSpect=0;

  TFile *file;
// Start loop on files
  Int_t nFiles=1,iEvent=0;
  for(Int_t iFile=0; iFile<nFiles; iFile++) {
    cout << "   FFFFFF     " << iFile+1 << " file     FFFFFF" << endl;

    if(iFile==0) {
      TFile *hfile = new TFile("hiau25-1000.root");
      cout << "   input file : hiau25-1000.root" << endl; 
//     TFile *hfile = new TFile("hipb-10ev.root");
//     cout << "   input file : hipb-10ev.root" << endl; 
    }
    if(iFile==1) {
      TFile *hfile = new TFile("~/RQMD/runB/rqmdB500AuAu25c-2.root");
      cout << "   input file : ~/RQMD/runB/rqmdB500AuAu25c-2.root" << endl; 
    }

    Int_t nparticles = (Int_t) h888->GetEntries();

    printf(" Number of nparticles per event = %i\n",nparticles/nEvents);
    if (nparticles <= 0) return 0;

//Loop over tracks (particles)
    for (Int_t i=0;i<nparticles;i++) {
      h888->GetEvent(i);
      geantCode=h888->GetLeaf("Idpart")->GetValue();
      hID->Fill(geantCode);

      theta=h888->GetLeaf("Theta")->GetValue();
      thedeg=theta;
      phi=h888->GetLeaf("Phi")->GetValue();
      p=h888->GetLeaf("P")->GetValue();
      if(p!=0) { // NOT target spectators 
   
        pt=p*TMath::Sin(theta*TMath::DegToRad());
        pz=p*TMath::Cos(theta*TMath::DegToRad());
        px=pt*TMath::Cos(phi*TMath::DegToRad());
        py=pt*TMath::Sin(phi*TMath::DegToRad());
        e=h888->GetLeaf("E")->GetValue();
        yrap=1.e30;   
        if(TMath::Abs(e != pz)) yrap=0.5*TMath::Log((e+pz)/(e-pz));
        if(geantCode==1) { //gamma
	  hGammaP->Fill(p);
	  hGammaTheta->Fill(theta);
	  hGammaPvsTheta->Fill(theta,p);
	  hGammaPt->Fill(pt);
	  hGammaY->Fill(yrap);
	  hGammaPtVsY->Fill(yrap,pt);
	  if(TMath::Abs(yrap-yCM)<deltaYrap) hSigma->Fill(pt,1.0/pt);
        }
      } else {
      nTagSpect++;
      } //end of if(p!=0) { // NOT target spectators
               
    } //end of for (Int_t i=0;i<nparticles;i++) { 

  } //end of for(Int_t iFile=0; iFile<nFiles; iFile++) {

  cout<<" Mumber of target spectators per event ="<<nTagSpect/nEvents<<endl;

  Int_t xCanvas=40,yCanvas=10,deltaX=5,deltaY=5;
  TCanvas *cID = new TCanvas("cID","HIJING analysis",xCanvas,yCanvas,700,700);
  gPad->SetLogy();
  gPad->SetGridy();
  gPad->SetGridx();
  hID->SetLineStyle(1);
  hID->SetLineWidth(3);
  hID->SetXTitle("GEANT code");
  hID->SetYTitle("multiplicity per event");
  hID->Scale(1./nEvents);
  hID->Draw();

  xCanvas+=deltaX; yCanvas+=deltaY;
  TCanvas *cGammaPY = new TCanvas("cGammaPY","analysis of ECAL hits",xCanvas,yCanvas,700,700);
//  gStyle->SetOptStat(0);
  cGammaPY->Divide(2,2);
  cGammaPY->cd(1);
  hGammaTheta->SetLineStyle(1);
  hGammaTheta->SetLineWidth(3);
  hGammaTheta->SetYTitle("N (1/event/bin)");
  hGammaTheta->SetXTitle("#theta (deg.)");
  hGammaTheta->Scale(1./nEvents);
  hGammaTheta->Draw();
  cGammaPY->cd(2);
  gPad->SetLogy();
  hGammaPt->SetLineStyle(1);
  hGammaPt->SetLineWidth(3);
  hGammaPt->SetYTitle("N (1/event/bin)");
  hGammaPt->SetXTitle("p_{t} (GeV/c)");
  hGammaPt->Scale(1./nEvents);
  hGammaPt->Draw();
  cGammaPY->cd(3);
  gPad->SetLogy();
  hGammaP->SetLineStyle(1);
  hGammaP->SetLineWidth(3);
  hGammaP->SetYTitle("N (1/event/bin)");
  hGammaP->SetXTitle("p (GeV/c)");
  hGammaP->Scale(1./nEvents);
  hGammaP->Draw();
  cGammaPY->cd(4);
  hGammaY->SetLineStyle(1);
  hGammaY->SetLineWidth(3);
  hGammaY->SetYTitle("N (1/event/bin)");
  hGammaY->SetXTitle("y");
  hGammaY->Scale(1./nEvents);
  hGammaY->Draw();
        
  gStyle->SetPalette(1);            
  xCanvas+=deltaX; yCanvas+=deltaY;
  TCanvas *cGammaPvsTheta = new TCanvas("cGammaPvsTheta","analysis of ECAL hits",xCanvas,yCanvas,700,700);
  hGammaPvsTheta->SetXTitle("#theta (deg.)");
  hGammaPvsTheta->SetYTitle("p (GeV/c)");
  hGammaPvsTheta->Scale(1./nEvents);
  hGammaPvsTheta->Draw("COLZ");
  gPad->SetGridx();
  gPad->SetGridy();

  Float_t mass=0.938,Elab,Plab,s,beta,gamma,Pcm,Ycm,Ecm;
  Float_t xmin,xmax;
 
  xmin=thetaMin; xmax=thetaMax;
  Float_t Tlab0=25.,Plab0,Plab1,Plab2,PlabP,PlabT,ElabP,ElabT,PcmP,ThetaLabForThetaCM90;
  Float_t massP,massT,massGamma=0.0,massPi=0.140,massRho0=0.770,massKstar0=0.892,massLambda=1.116,massDelta=1.232,massRes;
  massRes=massGamma;
 // massRes=massPi;
//  massRes=massRho0;
//  massRes=massKstar0;
//  massRes=mass;
//  massRes=massDelta;
  cout<<" mass(GeV) of particle = "<<massRes<<endl;
  cout<<" Tlab0="<<Tlab0<<endl;
  Plab0=TMath::Sqrt(Tlab0*(Tlab0+2*mass));

  Plab1=Plab0; Plab2=0;
  for(Int_t n=1;n<4;n++) { //number of nucleons in proj.
    massP=n*mass;
    PlabP=n*Plab1;
    ElabP=TMath::Sqrt(PlabP*PlabP+massP*massP);
    for(Int_t m=1;m<4;m++) { //number of nucleons in targ.
      cout<<" "<<n<<" "<<m<<" (proj. targ.)"<<endl;
      massT=m*mass;
      PlabT=m*Plab2;
      ElabT=TMath::Sqrt(PlabT*PlabT+massT*massT);
      s=TMath::Sqrt(massP*massP+massT*massT+2*(ElabP*ElabT+PlabP*PlabT));

      beta=(PlabP-PlabT)/(ElabP+ElabT); //=Vcm
      gamma=1/TMath::Sqrt(1-beta*beta);
      PcmP=gamma*(PlabP-beta*ElabP);
      Ycm=TMath::ATanH(beta);
      Pcm=TMath::Sqrt((s*s-(massT+massP+massRes)*(massT+massP+massRes))*(s*s-(massT+massP-massRes)*(massT+massP-massRes)))/(2.*s); // for pi, rho, omega, gamma
//      Pcm=TMath::Sqrt((s*s-(massT+massP-mass+massLambda+massRes)*(massT+massP-mass+massLambda+massRes))*(s*s-(massT+massP-mass+massLambda-massRes)*(massT+massP-mass+massLambda-massRes)))/(2.*s); // for Kstar0
//      Pcm=TMath::Sqrt((s*s-(massT+massP-mass+massRes)*(massT+massP-mass+massRes))*(s*s-(massT+massP-mass-massRes)*(massT+massP-mass-massRes)))/(2.*s); // for proton, Delta
     ThetaLabForThetaCM90=TMath::RadToDeg()*TMath::ATan(Pcm/TMath::Sqrt(Pcm*Pcm+massRes*massRes)/beta/gamma);
      cout<<"     PlabP="<<PlabP<<" PlabT="<<PlabT<<" s="<<s<<" beta="<<beta<<" gamma="<<gamma<<" Ycm="<<Ycm<<" PcmP="<<PcmP<<" ThetaLabForThetaCM90="<<ThetaLabForThetaCM90<<endl;
      cout<<"     momentum of part. in c.m.s="<<Pcm<<endl;

// par[0]=beta, par[1]=gamma, par[2]=Pcm, par[3]=mass
            TF1 *f00 = new TF1("f00F",fuPvsTheta,xmin,xmax,4);
            f00->SetParameters(beta,gamma,Pcm,massRes);
            f00->SetParNames("beta","gamma","Pcm","mass");
	    Int_t iColor;
            iColor=n+m-1;
	    if(iColor==6-1) iColor=6;
            f00->SetLineColor(iColor);
            f00->SetLineStyle(1);
            f00->SetLineWidth(2);
            f00->DrawCopy("same");

    } // end of for(Int_t m=2;m<3;m++) {
  } // end of for(Int_t n=2;n<3;n++) {



        
//  gStyle->SetPalette(1);            
  xCanvas+=deltaX; yCanvas+=deltaY;
  TCanvas *cGammaPtVsY = new TCanvas("cGammaPtVsY","analysis of ECAL hits",xCanvas,yCanvas,700,700);
  hGammaPtVsY->SetXTitle("y");
  hGammaPtVsY->SetYTitle("p_{t} (GeV/c)");
  hGammaPtVsY->Scale(1./nEvents);
  hGammaPtVsY->Draw("COLZ");
  gPad->SetGridx();
  gPad->SetGridy();

  xmin=yMin,xmax=yMax;
  for(Int_t n=1;n<4;n++) { //number of nucleons in proj.
    massP=n*mass;
    PlabP=n*Plab1;
    ElabP=TMath::Sqrt(PlabP*PlabP+massP*massP);
    for(Int_t m=1;m<4;m++) { //number of nucleons in targ.
      cout<<" "<<n<<" "<<m<<" (proj. targ.)"<<endl;
      massT=m*mass;
      PlabT=m*Plab2;
      ElabT=TMath::Sqrt(PlabT*PlabT+massT*massT);
      s=TMath::Sqrt(massP*massP+massT*massT+2*(ElabP*ElabT+PlabP*PlabT));

      beta=(PlabP-PlabT)/(ElabP+ElabT); //=Vcm
      gamma=1/TMath::Sqrt(1-beta*beta);
      PcmP=gamma*(PlabP-beta*ElabP);
      Ycm=TMath::ATanH(beta);
      Pcm=TMath::Sqrt((s*s-(massT+massP+massRes)*(massT+massP+massRes))*(s*s-(massT+massP-massRes)*(massT+massP-massRes)))/(2.*s); // for pi, rho, omega, gamma
//      Pcm=TMath::Sqrt((s*s-(massT+massP-mass+massLambda+massRes)*(massT+massP-mass+massLambda+massRes))*(s*s-(massT+massP-mass+massLambda-massRes)*(massT+massP-mass+massLambda-massRes)))/(2.*s); // for Kstar0
//      Pcm=TMath::Sqrt((s*s-(massT+massP-mass+massRes)*(massT+massP-mass+massRes))*(s*s-(massT+massP-mass-massRes)*(massT+massP-mass-massRes)))/(2.*s); // for proton, Delta
     ThetaLabForThetaCM90=TMath::RadToDeg()*TMath::ATan(Pcm/TMath::Sqrt(Pcm*Pcm+massRes*massRes)/beta/gamma);
      cout<<"     PlabP="<<PlabP<<" PlabT="<<PlabT<<" s="<<s<<" beta="<<beta<<" gamma="<<gamma<<" Ycm="<<Ycm<<" PcmP="<<PcmP<<" ThetaLabForThetaCM90="<<ThetaLabForThetaCM90<<endl;
      cout<<"     momentum of part. in c.m.s="<<Pcm<<endl;

// par[0]=beta, par[1]=gamma, par[2]=Pcm, par[3]=mass
// par[0]=Pcm, par[1]=mass, par[2]=Ycm
            TF1 *f00 = new TF1("f00F",fuPtVsY,xmin,xmax,3);
            f00->SetParameters(Pcm,massRes,Ycm);
            f00->SetParNames("beta","gamma","Pcm","mass");
	    Int_t iColor;
            iColor=n+m-1;
	    if(iColor==6-1) iColor=6;
            f00->SetLineColor(iColor);
            f00->SetLineStyle(1);
            f00->SetLineWidth(2);
            f00->DrawCopy("same");

    } // end of for(Int_t m=2;m<3;m++) {
  } // end of for(Int_t n=2;n<3;n++) {

  xCanvas+=deltaX; yCanvas+=deltaY;
  TCanvas *cGraph = new TCanvas("cGraph","analysis of ECAL hits",xCanvas,yCanvas,700,700);
  gPad->SetLogy();
  TLegend *legD = new TLegend(0.5,0.6,0.85,0.89);
// UrQMD direct photons spectra d^2N/p_t dp_t dy (GeV^{-2}) for Pb+Pb at E_lab=160 AGeV
// PR C57 3271 (1998), hep-ph/9709487
  const Int_t nPoints=9;
  Float_t pT[nPoints]={0.0,0.95,1.25,1.55,1.85,2.15,2.45,2.75,3.05};
  Float_t sig[nPoints]={1e-5,2.4e-1,5.5e-2,1.2e-2,6e-3,3e-3,1e-3,4.3e-4,4e-4};
  TGraph *graph = new TGraph(nPoints,pT,sig);
  graph->SetMarkerColor(2);
  graph->SetMarkerStyle(20);
  graph->SetMaximum(10);
  graph->SetMinimum(1.0e-4);
  graph->GetHistogram()->SetXTitle("p_{t} (GeV/c)");
  graph->GetHistogram()->SetYTitle("d^{2}N/p_{t}dp_{t}dy (GeV^{-2})");
  graph->SetTitle("direct photons");
  legD->AddEntry(graph,"UrQMD, Pb+Pb 160 AGeV","P");
  graph->Draw("AP");
  hSigma->SetLineStyle(1);
  hSigma->SetLineWidth(3);
  hSigma->Scale(1./nEvents/(2.*deltaYrap*deltaPt));
  legD->AddEntry(hSigma,"HIJING, Au+Au 25 AGeV","L");
  hSigma->Draw("same");
  legD->Draw();

  timer.Stop();
  Double_t rtime = timer.RealTime();
  Double_t ctime = timer.CpuTime();

  printf("RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime);

}