{
c1 = new TCanvas("rapidity","rapidity",200,100,700,500);
//
double event, event_previous = 0, id, px, py,pz,e, eta;
h1=new TH1F("h1", "multiplicity", 100,0,500.);
h2=new TH1F("h2", "pseudo-rapidity", 200,-10.,10.);
h3=new TH1F("h3", "pt", 200,0.,10.);
//
int n=0;
ifstream fin("output00.txt");
while (fin >> event >> id >> px>>py>>pz>>e ) {
double pt = sqrt(px*px + py*py);
if(pt<0.1) continue;  // cut out those with pt<100 MeV/C
double p=sqrt(px*px+py*py+pz*pz);
h3->Fill(pt);
eta=999.;
double theta = atan2(pt,pz);
eta = -log(tan(theta/2.));
//
//if(fabs(e-e*pz/p)>1.e-15) {
//  if( (e+e*pz/p)/(e-e*pz/p)>1.e-15)
//} 
if(n<20) cout<<"ievent,n, (px,py,pz) = "<<event<<" "<<n<<" ("<<
px<<","<<py<<","<<pz<<"), eta="<<eta<<endl;
if(event!=event_previous) {
h1->Fill(n);
n = 0;
}
event_previous = event;
h2->Fill(eta);
n++;
}
fin.close();
cout<<"number of evets="<<event<<" , number of particles="<<n<<endl;
//
h1->GetXaxis()->SetTitle("charged particle multiplicity");
h1->Draw();
c1->Print("multiplicity_900.png");
//
h2->Scale(1/event*10.);  // factor 10. for bin width correction
h2->GetXaxis()->SetTitle("#eta (pseudo rapidity)");
h2->GetYaxis()->SetTitle("charged particle per #eta");
h2->Draw();
c1->Print("eta_900.png");
//
h3->GetXaxis()->SetTitle("P_{t} (GeV/c))");
gPad->SetLogy();
h3->Draw();
c1->Print("pt_900.png");
gPad->SetLogy(0);
}