gROOT->Reset();
#include "Riostream.h"

void readfile(){
dptdst= new TH1F("ptdst","pt of K0s MC gen",100,0.0,0.2);

 TH1F dptdstr ("dptdstr", "pt of K0s MC reco",100,0.0,0.25);
 TH1F dptrho ("ptrho"," t ",100,0.,0.1);

   const Int_t n = 101;





  Float_t x,y,z;
   char s[164],h[28];
  ifstream in;
  Int_t nlines = 0;
   in.open("/Users/bhaag/desktop/root/data/MCrecoNewR/theoryComparison/starlight/xnxn01int.out");


float index=0.;
float pxyz1[3];
float pxyz2[3];
float x=0;

while (1) {
  in >> s;
  // cout<<s<<endl;
  // if(s=="TRACK:") cout<<" track found "<<endl;
  x=-20000.;

  sscanf(s,"%f",&x);
  //  cout<<" from scanf "<<h<<" "<<x<<endl;
 




  if(x==8.){

    index+=x;

  int iii=0; 

    while (1) {
      in >> s;
      sscanf(s,"%f",&x);
      //cout<<" inside loop "<<s<<" number "<<x;
      if(iii==0.)
	pxyz1[0]=x;
	if(iii==1.)
	pxyz1[1]=x;
      if(iii==2.)
	pxyz1[2]=x;
      iii++;
      if(iii==3) break;
    }





  }

 if(x==9.){

   index+=x;

  int iii=0; 
  
    while (1) {
      in >> s;
      sscanf(s,"%f",&x);
      //cout<<" inside loop "<<s<<" number "<<x;
      if(iii==0.)
	pxyz2[0]=x;
	if(iii==1.)
	pxyz2[1]=x;
      if(iii==2.)
	pxyz2[2]=x;
      iii++;
      if(iii==3) break;
    }





  }

 if(index>17) index=0;

 // cout<<index<<endl;
 if(index==17){
   //cout<<"event"<<endl;
    //dptdst->Fill(sqrt(pxyz1[0]*pxyz1[0]+pxyz1[1]*pxyz1[1]));
    //dptdst->Fill(sqrt(pxyz2[0]*pxyz2[0]+pxyz2[1]*pxyz2[1]));
    //dptdst->Fill(sqrt(pxyz2[0]*pxyz2[0]+pxyz2[1]*pxyz2[1]));
    dptdst->Fill(  sqrt((pxyz1[0]+pxyz2[0])*(pxyz1[0]+pxyz2[0])+(pxyz1[1]+pxyz2[1])*(pxyz1[1]+pxyz2[1])) );
    dptrho->Fill(  (pxyz1[0]+pxyz2[0])*(pxyz1[0]+pxyz2[0])+(pxyz1[1]+pxyz2[1])*(pxyz1[1]+pxyz2[1]) );


    /*cout<<" 8 pt "<<sqrt(pxyz1[0]*pxyz1[0]+pxyz1[1]*pxyz1[1])<<"    "<<pxyz1[0]<<"   "<<pxyz1[1]<<endl;    
    cout<<" 9 pt "<<sqrt(pxyz2[0]*pxyz2[0]+pxyz2[1]*pxyz2[1])<<"    "<<pxyz2[0]<<"   "<<pxyz2[1]<<endl;    
    cout<<" pt^2 of rho "<<
(pxyz1[0]+pxyz2[0])*(pxyz1[0]+pxyz2[0])+(pxyz1[1]+pxyz2[1])*(pxyz1[1]+pxyz2[1])<<endl;*/
   index=0.;
 }
 

      if (!in.good()) break;
      nlines++;
      /* if(nlines>45) break;
       */
}
//mytree->Draw("pPairs.pPt>>dptdstr");

TCanvas *c2 = new TCanvas("c2","tSpec");
//dptrho->GetXaxis()->SetLimits(0.0,0.01);
	
   const Int_t n = 101;

   Double_t pt[n] = {0.0000000, 2.0000000, 4.0000000, 6.0000000, 8.0000000, 10.0000000, 12.0000000, 14.0000000, 16.0000000,
18.0000000, 20.0000000, 22.0000000, 24.0000000, 26.0000000, 28.0000000, 30.0000000, 32.0000000, 34.0000000, 36.0000000, 38.0000000, 40.0000000, 42.0000000, 44.0000000, 46.0000000,
48.0000000, 50.0000000, 52.0000000, 54.0000000, 56.0000000, 58.0000000, 60.0000000, 62.0000000,
64.0000000, 66.0000000, 68.0000000, 70.0000000, 72.0000000, 74.0000000, 76.0000000, 78.0000000,
80.0000000, 82.0000000, 84.0000000, 86.0000000, 88.0000000, 90.0000000, 92.0000000, 94.0000000,
96.0000000, 98.0000000, 100.0000000, 102.0000000, 104.0000000, 106.0000000, 108.0000000,
110.0000000, 112.0000000, 114.0000000, 116.0000000, 118.0000000, 120.0000000, 122.0000000,
124.0000000, 126.0000000, 128.0000000, 130.0000000, 132.0000000, 134.0000000, 136.0000000, 138.0000000,
140.0000000, 142.0000000, 144.0000000, 146.0000000, 148.0000000, 150.0000000, 152.0000000,
154.0000000, 156.0000000, 158.0000000, 160.0000000, 162.0000000, 164.0000000, 166.0000000,
168.0000000, 170.0000000, 172.0000000, 174.0000000, 176.0000000, 178.0000000, 180.0000000,
182.0000000, 184.0000000, 186.0000000, 188.0000000, 190.0000000, 192.0000000, 194.0000000,
196.0000000, 198.0000000, 200.0000000 };

   Double_t x1[n] = {0.0000000, 0.0133314, 0.0477901, 0.0918247, 0.1361367, 0.1761786, 0.2102764, 0.2381471, 0.2602134, 0.2770273, 0.2891507, 0.2971844, 0.3016474, 0.3029967, 0.3016852, 0.2981032, 0.2925932, 0.2854916, 0.2770904, 0.2676442, 0.2573995, 0.2465688, 0.2353330, 0.2238635, 0.2123050, 0.2007753, 0.1893829, 0.1782154, 0.1673393, 0.1568130, 0.1466799, 0.1369677, 0.1276986, 0.1188850, 0.1105289, 0.1026292, 0.0951795, 0.0881672, 0.0815788, 0.0753988, 0.0696084, 0.0641896, 0.0591242, 0.0543928,0.0499772, 0.0458602, 0.0420243, 0.0384537, 0.0351334, 0.0320487, 0.0291862, 0.0265334, 0.0240782, 0.0218092, 0.0197160, 0.0177881, 0.0160160, 0.0143903, 0.0129021, 0.0115427, 0.0103039, 0.0091778, 0.0081564, 0.0072325, 0.0063989, 0.0056486, 0.0049752, 0.0043723, 0.0038339, 0.0033545, 0.0029287, 0.0025514, 0.0022180, 0.0019242, 0.0016658, 0.0014392, 0.0012410, 0.0010680, 0.0009173, 0.0007863, 0.0006728, 0.0005747, 0.0004899, 0.0004169, 0.0003542, 0.0003004, 0.0002543, 0.0002149, 0.0001813, 0.0001527, 0.0001283, 0.0001077, 0.0000902, 0.0000755, 0.0000630, 0.0000525, 0.0000437, 0.0000363, 0.0000301, 0.0000249, 0.0000206 };



for(Int_t j =0; j < n; j++){
	x1[j]=x1[j]*31.341;
}

for(Int_t m =0; m < n; m++){
	pt[m]=pt[m]/1000;
}

TH1F *hencken = new TH1F("hencken","noExcitation",n-1,0.0,0.1);
   for (Int_t k=0;k<=100;k++) {
      hencken->SetBinContent(k+1,x1[k]);
   }
   hencken->Draw();

Double_t hintegral = hencken.Integral(0,100);
Double_t dptintegral = dptrho.Integral(0,100);
Double_t dptdstint = dptdst.Integral(0,50);

hencken->Scale((1/hintegral));
dptrho->Scale((1/dptintegral));
dptdst->Scale((1/dptdstint));

    Float_t a[n];
	for(Int_t g=0;g<n;g++){
		a[g]=dptdst->GetBinContent(g+1);
		Float_t num = 0;
		num = dptdst->GetBinCenter(g+1);
		//printf("%f %f \n", num, a[g]);
		printf("%f \n", a[g]);
		}
new TCanvas();
hencken->Draw();
dptrho->Draw("same");

 TFile myfile("dptrho.root","RECREATE");
   dptdst->Write();
   dptrho->Write();
   hencken->Write();
   myfile.Close();

}