Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:29:37

0001 #include <iostream>
0002 #include <fstream>
0003 #include <iomanip>
0004 #include <string>
0005 #include <list>
0006 
0007 #include <math.h>
0008 #include <vector>
0009 
0010 #include "Rtypes.h"
0011 #include "TROOT.h"
0012 #include "TRint.h"
0013 #include "TObject.h"
0014 #include "TFile.h"
0015 #include "TTree.h"
0016 #include "TH1F.h"
0017 #include "TCanvas.h"
0018 #include "TApplication.h"
0019 #include "TRefArray.h"
0020 #include "TStyle.h"
0021 #include "TGraph.h"
0022 
0023 void AnalyseITEP(char element[2], char list[20], char ene[6], char part[2]="p", int scan=1, char plot='Y') {
0024 
0025   static double massP = 938.272;
0026   static double massN = 939.565;
0027   static double pi  = 3.1415926;
0028   static double deg = pi/180.;
0029   static double dcth= 0.05, de=0.02;
0030   double massp = massP;
0031   std::string fpart("Proton");
0032   if (part == "n" || part == "N") { 
0033     massp = massN;
0034     fpart = "Neutron";
0035   }
0036   char fname[40];
0037   sprintf (fname, "%s%s%sGeV.root", element, list, ene);
0038 
0039   double rhol = rhoL(element);
0040   double atwt = atomicWt(element);
0041   cout << fname << " rhoL " << rhol << " atomic weight " << atwt << "\n";
0042   std::vector<double> angles = angleScan(scan);
0043   std::vector<double> energies = energyScan(part);
0044 
0045   char name[60], title[160];
0046   sprintf (title, "All %s", fpart.c_str());
0047   TH1F *hiK0 = new TH1F ("hiK0", title,  800, 0.,8.);
0048   TH1F *hiC0 = new TH1F ("hiC0", title,  100,-1.,1.);
0049   sprintf (title, "Elastc Scattered %s", fpart.c_str());
0050   TH1F *hiK1 = new TH1F ("hiK1", title,  800, 0.,8.);
0051   TH1F *hiC1 = new TH1F ("hiC1", title,  100,-1.,1.);
0052   sprintf (title, "Inelastc Scattered %s", fpart.c_str());
0053   TH1F *hiK2 = new TH1F ("hiK2", title,  800, 0.,8.);
0054   TH1F *hiC2 = new TH1F ("hiC2", title,  100,-1.,1.);
0055   std::vector<double> cthmin, cthmax;
0056   TH1F *hiKE1[30], *hiKE2[30];
0057   unsigned int ii=0;
0058   for (ii=0; ii<angles.size(); ii++) {
0059     double cth = cos(angles[ii]);
0060     cthmin.push_back(cth-0.5*dcth);
0061     cthmax.push_back(cth+0.5*dcth);
0062     sprintf (name, "KE1%s%s%sGeV%5.1f", element, list, ene, angles[ii]/deg);
0063     sprintf (title, "p+%s at %s GeV (%s) (#theta = %8.2f)", element, ene, list, angles[ii]/deg);
0064     hiKE1[ii] = new TH1F (name, title, 800, 0., 8.);
0065     //std::cout << "hiKE1[" << ii << "] = " << hiKE1[ii] << " " <<  name << "   " << title << "\n";
0066     sprintf (name, "KE2%s%s%sGeV%5.1f", element, list, ene, angles[ii]/deg);
0067     sprintf (title, "p+%s at %s GeV (%s) (#theta = %8.2f)", element, ene, list, angles[ii]/deg);
0068     hiKE2[ii] = new TH1F (name, title, 800, 0., 8.);
0069     //std::cout << "hiKE2[" << ii << "] = " << hiKE2[ii] << " " <<  name << "   " << title << "\n";
0070   }
0071 
0072   std::vector<double> emin, emax;
0073   TH1F *hiCT1[30], *hiCT2[30];
0074   for (ii=0; ii<energies.size(); ii++) {
0075     double en = energies[ii];
0076     emin.push_back(en-0.5*de);
0077     emax.push_back(en+0.5*de);
0078     sprintf (name, "CT1%s%s%sGeV%4.2f", element, list, ene, energies[ii]);
0079     sprintf (title, "p+%s at %s GeV (%s) (KE = %6.2f GeV)", element, ene, list, energies[ii]);
0080     hiCT1[ii] = new TH1F (name, title, 80, -1., 1.);
0081     //std::cout << "hiCT1[" << ii << "] = " << hiCT1[ii] << " " <<  name << "   " << title << "\n";
0082     sprintf (name, "CT2%s%s%sGeV%4.2f", element, list, ene, energies[ii]);
0083     sprintf (title, "p+%s at %s GeV (%s) (KE = %6.2f GeV)", element, ene, list, energies[ii]);
0084     hiCT2[ii] = new TH1F (name, title, 80, -1., 1.);
0085     //std::cout << "hiCT2[" << ii << "] = " << hiCT2[ii] << " " <<  name << "   " << title << "\n";
0086   }
0087 
0088   TFile *file = new TFile(fname);
0089   TTree *tree = (TTree *) file->Get("T1");
0090   int interval = 100000;
0091   if (plot == 'N' || plot == 'n') interval = 100000;
0092 
0093   if (!tree) {
0094     std::cout << "Cannot find Tree T1 in file " << fname << "\n";
0095   } else {
0096     std::cout << "Tree T1 found with " << tree->GetEntries() << " entries\n";
0097     int nentry = tree->GetEntries();
0098     int ninter=0, elastic=0, inelastic=0;
0099     for (int i=0; i<nentry; i++) {
0100       if (i%interval == 0) std:cout << "Started with event # " << i << "\n";
0101       std::vector<int>                     *nsec, *procids;
0102       std::vector<double>                  *px, *py, *pz, *mass;
0103       std::vector<std::string>             *procs;
0104       tree->SetBranchAddress("NumberSecondaries", &nsec);
0105       tree->SetBranchAddress("ProcessID",         &procids);
0106       //      tree->SetBranchAddress("ProcessNames",      &procs);
0107       tree->SetBranchAddress("SecondaryPx",       &px);
0108       tree->SetBranchAddress("SecondaryPy",       &py);
0109       tree->SetBranchAddress("SecondaryPz",       &pz);
0110       tree->SetBranchAddress("SecondaryMass",     &mass);
0111       tree->GetEntry(i);
0112       if ((*nsec).size() > 0) {
0113     ninter++;
0114     bool isItElastic = false;
0115     if ((*procids)[0] == 17) {elastic++; isItElastic = true;}
0116     else                     inelastic++;
0117 
0118     if (ninter <3) {
0119       std::cout << "Interaction " << ninter << "/" << i+1 << " Type "
0120             << (*procids)[0]  << " with " << (*nsec)[0] << " secondaries\n";
0121       for (int k=0; k<(*nsec)[0]; k++)
0122         std::cout << " Secondary " << k << " Px " << (*px)[k] << " Py " << (*py)[k] << " Pz " << (*pz)[k] << " Mass " << (*mass)[k] << "\n";
0123     }
0124 
0125     for (int k=0; k<(*nsec)[0]; k++) {
0126       if (abs((*mass)[k]-massp) < 0.01) { // This is the required particle
0127         double pl = (*py)[k];
0128         double pt = ((*px)[k])*((*px)[k])+((*pz)[k])*((*pz)[k]);
0129         double pp = (pt+pl*pl);
0130         double ke = (sqrt (pp + massp*massp) - massp)/1000.;
0131         pp        = sqrt (pp);
0132         double cth= (pp == 0. ? -2. : (pl/pp));
0133         double wt = (pp == 0. ?  0. : (1000./pp));
0134         // std::cout << "Entry " << i << " Secondary " << k << " Cth " << cth << " KE " << ke << " WT " << wt << "\n";
0135         hiK0->Fill(ke);
0136         hiC0->Fill(cth);
0137         if (isItElastic) {
0138           hiK1->Fill(ke);
0139           hiC1->Fill(cth);
0140         } else {
0141           hiK2->Fill(ke);
0142           hiC2->Fill(cth);
0143           for (ii=0; ii<angles.size(); ii++) {
0144         if (cth > cthmin[ii] && cth <= cthmax[ii]) {
0145           // std::cout << " Loop " << ii << " Limit " << cthmin[ii] << " " << cthmax[ii] << " " << hiKE1[ii] << " " << hiKE2[ii] << "\n";
0146           hiKE1[ii]->Fill(ke);
0147           hiKE2[ii]->Fill(ke,wt);
0148         }
0149           }
0150           for (ii=0; ii<energies.size(); ii++) {
0151         if (ke > emin[ii] && ke <= emax[ii]) {
0152           // std::cout << " Loop " << ii << " Limit " << emin[ii] << " " << emax[ii] << " " << hiCT1[ii] << " " << hiCT2[ii] << "\n";
0153           hiCT1[ii]->Fill(cth);
0154           hiCT2[ii]->Fill(cth,wt);
0155         }
0156           }
0157         }
0158       }
0159     }
0160       }
0161     }
0162 
0163     std::cout << ninter << " interactions seen in " << nentry << " trials\n";
0164     double sigma = atwt*10000.*log((double)(nentry)/(double)(nentry-ninter))/(rhol*6.023);
0165     double dsigma    = sigma/sqrt(double(max(ninter,1)));
0166     double sigmaEl   = sigma*((double)(elastic))/((double)(max(ninter,1)));
0167     double dsigmaEl  = sigmaEl/sqrt(double(max(elastic,1)));
0168     double sigmaInel = sigma*((double)(inelastic))/((double)(max(ninter,1)));
0169     double dsigmaInel= sigmaInel/sqrt(double(max(inelastic,1)));
0170     std::cout << "Total     " << sigma << " +- " << dsigma 
0171           << " mb (" << ninter << " events)\n"
0172           << "Elastic   " << sigmaEl<< " +- " << dsigmaEl
0173           << " mb (" << ninter << " events)\n"
0174           << "Inelastic " << sigmaInel << " +- " << dsigmaInel
0175           << " mb (" << ninter << " events)\n";
0176   }
0177 
0178   gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite);
0179   gStyle->SetPadColor(kWhite);    gStyle->SetFrameBorderMode(0);
0180   gStyle->SetFrameBorderSize(1);  gStyle->SetFrameFillColor(0);
0181   gStyle->SetFrameFillStyle(0);   gStyle->SetFrameLineColor(1);
0182   gStyle->SetFrameLineStyle(1);   gStyle->SetFrameLineWidth(1);
0183   gStyle->SetOptLogy(1);          gStyle->SetTitleOffset(1.2,"Y");
0184     
0185   sprintf (title, "Kinetic Energy of %s (GeV)", fpart.c_str());
0186   hiK0->GetXaxis()->SetTitle(title);
0187   hiK1->GetXaxis()->SetTitle(title);
0188   hiK2->GetXaxis()->SetTitle(title);
0189   if (plot != 'N' && plot != 'n') {
0190     TCanvas *c1 = new TCanvas("c1","K.E.",800,600); c1->Divide(2,2);
0191     c1->cd(1); hiK1->Draw(); c1->cd(2); hiK2->Draw(); c1->cd(3); hiK0->Draw();
0192   }
0193   sprintf (title, "cos (#theta) of scattered %s", fpart.c_str());
0194   hiC0->GetXaxis()->SetTitle(title);
0195   hiC1->GetXaxis()->SetTitle(title);
0196   hiC2->GetXaxis()->SetTitle(title);
0197   if (plot != 'N' && plot != 'n') {
0198     TCanvas *c2 = new TCanvas("c2","cos#theta",800,600); c2->Divide(2,2);
0199     c2->cd(1); hiC1->Draw(); c2->cd(2); hiC2->Draw(); c2->cd(3); hiC0->Draw();
0200   }
0201   TCanvas *cc[30];
0202   TH1F    *hiKE0[30];
0203   for (ii=0; ii<angles.size(); ii++) {
0204     double xbin = hiKE1[ii]->GetBinWidth(1);
0205     sprintf (title, "Kinetic Energy of %s (GeV)", fpart.c_str());
0206     hiKE1[ii]->GetXaxis()->SetTitle(title);
0207     sprintf (title, "Events/%6.3f GeV", xbin);
0208     hiKE1[ii]->GetYaxis()->SetTitle(title);
0209     double xbin  = hiKE2[ii]->GetBinWidth(1);
0210     double scale = sigmaInel/(((double)(max(inelastic,1)))*xbin*2.*pi*dcth);
0211     std::cout << "Bin " << ii << " Angle " << angles[ii]/deg << " Bin " << xbin << " Scale " << scale << " " << title << "\n";
0212     sprintf (title, "Kinetic Energy of %s (GeV)", fpart.c_str());
0213     hiKE2[ii]->GetXaxis()->SetTitle(title);
0214     sprintf (title, "Events (scaled by #frac{1}{p})/%6.3f GeV", xbin);
0215     hiKE2[ii]->GetYaxis()->SetTitle(title);
0216     sprintf (name, "KE0%s%s%sGeV%5.1f", element, list, ene, angles[ii]/deg);
0217     hiKE0[ii] = (TH1F*)hiKE2[ii]->Clone();
0218     hiKE0[ii]->SetName(name); hiKE0[ii]->Scale(scale);
0219     hiKE0[ii]->GetYaxis()->SetTitle("E#frac{d^{3}#sigma}{dp^{3}} (mb/GeV^{2})");
0220       
0221     if (plot != 'N' && plot != 'n') {
0222       sprintf(name, "Canvas%i", ii);
0223       sprintf (title, "p+%s at %s GeV (%s) (#theta = %8.2f)", element, ene, list, angles[ii]/deg);
0224       cc[ii] = new TCanvas(name,title,800,600); cc[ii]->Divide(2,2);
0225     
0226       std::cout << "hiKE1: " << hiKE1[ii]->GetName() << " " << hiKE1[ii]->GetEntries() << " " << hiKE1[ii] << "\n";
0227       cc[ii]->cd(1); hiKE1[ii]->Draw();
0228       std::cout << "hiKE0: " << hiKE0[ii]->GetName() << " " << hiKE0[ii]->GetEntries() << " " << hiKE0[ii] << "\n";
0229       cc[ii]->cd(2); hiKE0[ii]->Draw();
0230       std::cout << "hiKE2: " << hiKE2[ii]->GetName() << " " << hiKE2[ii]->GetEntries() << " " << hiKE2[ii] << "\n";
0231       cc[ii]->cd(3); hiKE2[ii]->Draw(); 
0232     }
0233   }
0234     
0235   TCanvas *ct[30];
0236   TH1F    *hiCT0[30];
0237   for (ii=0; ii<energies.size(); ii++) {
0238     double xbin = hiCT1[ii]->GetBinWidth(1);
0239     sprintf (title, "Events/%6.3f", xbin);
0240     hiCT1[ii]->GetXaxis()->SetTitle("cos (#theta)");
0241     hiCT1[ii]->GetYaxis()->SetTitle(title);
0242     double xbin  = hiCT2[ii]->GetBinWidth(1);
0243     double scale = sigmaInel/(((double)(max(inelastic,1)))*xbin*2.*pi*de);
0244     std::cout << "Bin " << ii << " KE " << energies[ii] << " GeV Bin " << xbin << " Scale " << scale << " " << title << "\n";
0245     sprintf (title, "Events (scaled by #frac{1}{p})/%6.3f", xbin);
0246     hiCT2[ii]->GetXaxis()->SetTitle("cos (#theta)");
0247     hiCT2[ii]->GetYaxis()->SetTitle(title);
0248     sprintf (name, "CT0%s%s%sGeV%4.2f", element, list, ene, energies[ii]);
0249     hiCT0[ii] = (TH1F*)hiCT2[ii]->Clone();
0250     hiCT0[ii]->SetName(name); hiCT0[ii]->Scale(scale);
0251     hiCT0[ii]->GetYaxis()->SetTitle("E#frac{d^{3}#sigma}{dp^{3}} (mb/GeV^{2})");
0252       
0253     if (plot != 'N' && plot != 'n') {
0254       sprintf(name, "Canvas0%i", ii);
0255       sprintf (title, "p+%s at %s GeV (%s) (KE = %6.2f GeV)", element, ene, list, energies[ii]);
0256       ct[ii] = new TCanvas(name,title,800,600); ct[ii]->Divide(2,2);
0257       std::cout << "hiCT1: " << hiCT1[ii]->GetName() << " " << hiCT1[ii]->GetEntries() << " " << hiCT1[ii] << "\n";
0258       ct[ii]->cd(1); hiCT1[ii]->Draw();
0259       std::cout << "hiCT0: " << hiCT0[ii]->GetName() << " " << hiCT0[ii]->GetEntries() << " " << hiCT0[ii] << "\n";
0260       ct[ii]->cd(2); hiCT0[ii]->Draw();
0261       std::cout << "hiCT2: " << hiCT2[ii]->GetName() << " " << hiCT2[ii]->GetEntries() << " " << hiCT2[ii] << "\n";
0262       ct[ii]->cd(3); hiCT2[ii]->Draw(); 
0263     }
0264   }
0265 
0266   char ofile[40];
0267   sprintf (ofile, "%s%s%sGeV_%i.root", element, list, ene, scan);
0268   TFile f(ofile, "recreate");
0269   hiK0->Write(); hiK1->Write(); hiK2->Write();
0270   hiC0->Write(); hiC1->Write(); hiC2->Write();
0271   for (ii=0; ii<angles.size(); ii++) {
0272     hiKE1[ii]->Write(); hiKE0[ii]->Write(); hiKE2[ii]->Write();
0273   }
0274   for (ii=0; ii<energies.size(); ii++) {
0275     hiCT1[ii]->Write(); hiCT0[ii]->Write(); hiCT2[ii]->Write();
0276   }
0277   f.Close();
0278   std::cout << "o/p saved in file " << ofile << "\n";
0279 
0280   file->Close();
0281 }
0282 
0283 void AnalyseBNL(char element[2], char list[10], char ene[6], char part[4]="p", char plot='Y') {
0284   
0285   static double pi   = 3.1415926;
0286   static double deg  = pi/180.;
0287   int           indx = 13;
0288   std::string fpart("Proton");
0289   if      (part == "n"    || part == "N")    {fpart = "Neutron";    indx = 14;}
0290   else if (part == "pbar" || part == "PBAR") {fpart = "\bar{p}";    indx = 12;}
0291   else if (part == "kz"   || part == "KZ")   {fpart = "K^{0}";      indx = 11;}
0292   else if (part == "kp"   || part == "KP")   {fpart = "K^{+}";      indx = 10;}
0293   else if (part == "km"   || part == "KM")   {fpart = "K^{-}";      indx = 9;}
0294   else if (part == "pip"  || part == "PIP")  {fpart = "#pi^{+}";    indx = 8;}
0295   else if (part == "pim"  || part == "PIM")  {fpart = "#pi^{-}";    indx = 7;}
0296   char fname[40];
0297   sprintf (fname, "%s%s%sGeV.root", element, list, ene);
0298   
0299   double rhol = rhoL(element);
0300   double atwt = atomicWt(element);
0301   cout << fname << " rhoL " << rhol << " atomic weight " << atwt << "\n";
0302   std::vector<double> rapidities = rapidityScan();
0303   
0304   char name[60], title[160];
0305   sprintf (title, "All %s", fpart.c_str());
0306   TH1F *hiK0 = new TH1F ("hiK0", title,  500, 0.,5.);
0307   TH1F *hiC0 = new TH1F ("hiC0", title,  100,-1.,1.);
0308   sprintf (title, "Elastc Scattered %s", fpart.c_str());
0309   TH1F *hiK1 = new TH1F ("hiK1", title,  500, 0.,5.);
0310   TH1F *hiC1 = new TH1F ("hiC1", title,  100,-1.,1.);
0311   sprintf (title, "Inelastc Scattered %s", fpart.c_str());
0312   TH1F *hiK2 = new TH1F ("hiK2", title,  500, 0.,5.);
0313   TH1F *hiC2 = new TH1F ("hiC2", title,  100,-1.,1.);
0314   std::vector<double> ymin, ymax;
0315   TH1F *hiKE1[30], *hiKE2[30];
0316   unsigned int ii=0;
0317   for (ii=0; ii<rapidities.size()-1; ii++) {
0318     ymin.push_back(rapidities[ii]);
0319     ymax.push_back(rapidities[ii+1]);
0320     double yv = 0.5*(rapidities[ii]+rapidities[ii+1]);
0321     sprintf (name, "KE1%s%s%sGeVy%4.2f", element, list, ene, yv);
0322     sprintf (title, "p+%s at %s GeV (%s) (y = %8.2f)", element, ene, list, yv);
0323     hiKE1[ii] = new TH1F (name, title, 800, 0., 8.);
0324     std::cout << "hiKE1[" << ii << "] = " << hiKE1[ii] << " " <<  name << "   " << title << "\n";
0325     sprintf (name, "KE2%s%s%sGeVy%4.2f", element, list, ene, yv);
0326     sprintf (title, "p+%s at %s GeV (%s) (y = %8.2f)", element, ene, list, yv);
0327     hiKE2[ii] = new TH1F (name, title, 800, 0., 8.);
0328     std::cout << "hiKE2[" << ii << "] = " << hiKE2[ii] << " " <<  name << "   " << title << "\n";
0329   }
0330   
0331   TFile *file = new TFile(fname);
0332   TTree *tree = (TTree *) file->Get("T1");
0333   int interval = 100000;
0334   if (plot == 'N' || plot == 'n') interval = 100000;
0335   
0336   if (!tree) {
0337     std::cout << "Cannot find Tree T1 in file " << fname << "\n";
0338   } else {
0339     std::cout << "Tree T1 found with " << tree->GetEntries() << " entries\n";
0340     int nentry = tree->GetEntries();
0341     int ninter=0, elastic=0, inelastic=0;
0342     for (int i=0; i<nentry; i++) {
0343       if (i%interval == 0) std:cout << "Started with event # " << i << "\n";
0344       std::vector<int>                     *nsec, *procids;
0345       std::vector<double>                  *px, *py, *pz, *mass;
0346       std::vector<std::string>             *procs;
0347       tree->SetBranchAddress("NumberSecondaries", &nsec);
0348       tree->SetBranchAddress("ProcessID",         &procids);
0349       //      tree->SetBranchAddress("ProcessNames",      &procs);
0350       tree->SetBranchAddress("SecondaryPx",       &px);
0351       tree->SetBranchAddress("SecondaryPy",       &py);
0352       tree->SetBranchAddress("SecondaryPz",       &pz);
0353       tree->SetBranchAddress("SecondaryMass",     &mass);
0354       tree->GetEntry(i);
0355       if ((*nsec).size() > 0) {
0356     ninter++;
0357     bool isItElastic = false;
0358     if ((*procids)[0] == 17) {elastic++; isItElastic = true;}
0359     else                     inelastic++;
0360     
0361     if (ninter <10) {
0362       std::cout << "Interaction " << ninter << "/" << i+1 << " Type "
0363             << (*procids)[0]  << " with " << (*nsec)[0] << " secondaries\n";
0364       for (int k=0; k<(*nsec)[0]; k++)
0365         std::cout << " Secondary " << k << " Px " << (*px)[k] << " Py " << (*py)[k] << " Pz " << (*pz)[k] << " Mass " << (*mass)[k] << "\n";
0366     }
0367     
0368     for (int k=0; k<(*nsec)[0]; k++) {
0369       double massp = (*mass)[k];
0370       int    type  = type(massp);
0371       if (type == indx) { // This is the required particle
0372         double pl = (*py)[k];
0373         double pt = ((*px)[k])*((*px)[k])+((*pz)[k])*((*pz)[k]);
0374         double pp = (pt+pl*pl);
0375         double mt = sqrt (pt + massp*massp);
0376         double mtp= (mt - std::abs(massp))/1000.;
0377         double ee = sqrt(pp+massp*massp);
0378         double yv = 0.5*log((ee+pl)/(ee-pl));
0379         pp        = sqrt (pp);
0380         double cth= (pp == 0. ? -2. : (pl/pp));
0381         double wt = (mt == 0. ?  0. : (1000./mt));
0382         if (ninter <10) std::cout << "Entry " << i << " Secondary " << k << " yv " << yv << " mtp " << mtp << " WT " << wt << "\n";
0383         hiK0->Fill(mtp);
0384         hiC0->Fill(cth);
0385         if (isItElastic) {
0386           hiK1->Fill(mtp);
0387           hiC1->Fill(cth);
0388         } else {
0389           hiK2->Fill(mtp);
0390           hiC2->Fill(cth);
0391           for (ii=0; ii<ymin.size(); ii++) {
0392         if (yv > ymin[ii] && yv <= ymax[ii]) {
0393           if (ninter <10) std::cout << " Loop " << ii << " Limit " << ymin[ii] << " " << ymax[ii] << " " << hiKE1[ii] << " " << hiKE2[ii] << "\n";
0394           hiKE1[ii]->Fill(mtp);
0395           hiKE2[ii]->Fill(mtp,wt);
0396         }
0397           }
0398         }
0399       }
0400     }
0401       }
0402     }
0403     
0404     std::cout << ninter << " interactions seen in " << nentry << " trials\n";
0405     double sigma = atwt*10000.*log((double)(nentry)/(double)(nentry-ninter))/(rhol*6.023);
0406     double dsigma    = sigma/sqrt(double(max(ninter,1)));
0407     double sigmaEl   = sigma*((double)(elastic))/((double)(max(ninter,1)));
0408     double dsigmaEl  = sigmaEl/sqrt(double(max(elastic,1)));
0409     double sigmaInel = sigma*((double)(inelastic))/((double)(max(ninter,1)));
0410     double dsigmaInel= sigmaInel/sqrt(double(max(inelastic,1)));
0411     std::cout << "Total     " << sigma << " +- " << dsigma 
0412           << " mb (" << ninter << " events)\n"
0413           << "Elastic   " << sigmaEl<< " +- " << dsigmaEl
0414           << " mb (" << ninter << " events)\n"
0415           << "Inelastic " << sigmaInel << " +- " << dsigmaInel
0416           << " mb (" << ninter << " events)\n";
0417   }
0418 
0419   gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite);
0420   gStyle->SetPadColor(kWhite);    gStyle->SetFrameBorderMode(0);
0421   gStyle->SetFrameBorderSize(1);  gStyle->SetFrameFillColor(0);
0422   gStyle->SetFrameFillStyle(0);   gStyle->SetFrameLineColor(1);
0423   gStyle->SetFrameLineStyle(1);   gStyle->SetFrameLineWidth(1);
0424   gStyle->SetOptLogy(1);          gStyle->SetTitleOffset(1.2,"Y");
0425     
0426   sprintf (title, "Reduced Transverse mass of %s (GeV)", fpart.c_str());
0427   hiK0->GetXaxis()->SetTitle(title);
0428   hiK1->GetXaxis()->SetTitle(title);
0429   hiK2->GetXaxis()->SetTitle(title);
0430   if (plot != 'N' && plot != 'n') {
0431     TCanvas *c1 = new TCanvas("c1","Reduced Transverse Mass",800,600); 
0432     c1->Divide(2,2);
0433     c1->cd(1); hiK1->Draw(); c1->cd(2); hiK2->Draw(); c1->cd(3); hiK0->Draw();
0434   }
0435   sprintf (title, "cos (#theta) of scattered %s", fpart.c_str());
0436   hiC0->GetXaxis()->SetTitle(title);
0437   hiC1->GetXaxis()->SetTitle(title);
0438   hiC2->GetXaxis()->SetTitle(title);
0439   if (plot != 'N' && plot != 'n') {
0440     TCanvas *c2 = new TCanvas("c2","cos#theta",800,600); c2->Divide(2,2);
0441     c2->cd(1); hiC1->Draw(); c2->cd(2); hiC2->Draw(); c2->cd(3); hiC0->Draw();
0442   }
0443   TCanvas *cc[30];
0444   TH1F    *hiKE0[30];
0445   for (ii=0; ii<ymin.size(); ii++) {
0446     double xbin = hiKE1[ii]->GetBinWidth(1);
0447     sprintf (title, "Reduced Transverse Mass of %s (GeV)", fpart.c_str());
0448     hiKE1[ii]->GetXaxis()->SetTitle(title);
0449     sprintf (title, "Events/%6.3f GeV", xbin);
0450     hiKE1[ii]->GetYaxis()->SetTitle(title);
0451     double xbin  = hiKE2[ii]->GetBinWidth(1);
0452     double dy    = (ymax[ii]-ymin[ii]);
0453     double yv    = 0.5*(ymin[ii]+ymax[ii]);
0454     double scale = sigmaInel/(((double)(max(inelastic,1)))*xbin*2.*pi*dy);
0455     std::cout << "Bin " << ii << " yv " << yv << " Bin " << xbin << " Scale " << scale << " " << title << "\n";
0456     sprintf (title, "Reduced Transverse mass of %s (GeV)", fpart.c_str());
0457     hiKE2[ii]->GetXaxis()->SetTitle(title);
0458     sprintf (title, "Events (scaled by #frac{1}{mT})/%6.3f GeV",xbin);
0459     hiKE2[ii]->GetYaxis()->SetTitle(title);
0460     sprintf (name, "KE0%s%s%sGeVy%4.2f", element, list, ene, yv);
0461     hiKE0[ii] = (TH1F*)hiKE2[ii]->Clone();
0462     hiKE0[ii]->SetName(name); hiKE0[ii]->Scale(scale);
0463     hiKE0[ii]->GetYaxis()->SetTitle("E#frac{d^{3}#sigma}{dp^{3}} (mb/GeV^{2})");
0464       
0465     if (plot != 'N' && plot != 'n') {
0466       sprintf(name, "Canvas%i", ii);
0467       sprintf (title, "p+%s at %s GeV (%s) (y = %5.2f)", element,ene,list,yv);
0468       cc[ii] = new TCanvas(name,title,800,600); cc[ii]->Divide(2,2);
0469       std::cout << "hiKE1: " << hiKE1[ii]->GetName() << " " << hiKE1[ii]->GetEntries() << " " << hiKE1[ii] << "\n";
0470       cc[ii]->cd(1); hiKE1[ii]->Draw();
0471       std::cout << "hiKE0: " << hiKE0[ii]->GetName() << " " << hiKE0[ii]->GetEntries() << " " << hiKE0[ii] << "\n";
0472       cc[ii]->cd(2); hiKE0[ii]->Draw();
0473       std::cout << "hiKE2: " << hiKE2[ii]->GetName() << " " << hiKE2[ii]->GetEntries() << " " << hiKE2[ii] << "\n";
0474       cc[ii]->cd(3); hiKE2[ii]->Draw(); 
0475     }
0476   }
0477     
0478   char ofile[40];
0479   sprintf (ofile, "%s%s%sGeV_1.root", element, list, ene);
0480   TFile f(ofile, "recreate");
0481   hiK0->Write(); hiK1->Write(); hiK2->Write();
0482   hiC0->Write(); hiC1->Write(); hiC2->Write();
0483   for (ii=0; ii<ymin.size(); ii++) {
0484     hiKE1[ii]->Write(); hiKE0[ii]->Write(); hiKE2[ii]->Write();
0485   }
0486   f.Close();
0487   std::cout << "o/p saved in file " << ofile << "\n";
0488 
0489   file->Close();
0490 }
0491 
0492 double rhoL(char element[2]) {
0493 
0494   double tmp=0;
0495   if      (element == "H")   tmp = 0.0708 * 800.;
0496   else if (element == "Be")  tmp = 1.848 * 80.;
0497   else if (element == "C")   tmp = 2.265 * 80.;
0498   else if (element == "Al")  tmp = 2.700 * 80.;
0499   else if (element == "Ti")  tmp = 4.530 * 40.;
0500   else if (element == "Fe")  tmp = 7.870 * 30.;
0501   else if (element == "Cu")  tmp = 8.960 * 30.;
0502   else if (element == "Nb")  tmp = 8.550 * 30.;
0503   else if (element == "Cd")  tmp = 8.630 * 30.;
0504   else if (element == "Sn")  tmp = 7.310 * 35.;
0505   else if (element == "Ta")  tmp = 16.65 * 20.;
0506   else if (element == "Au")  tmp = 18.85 * 20.;
0507   else if (element == "Pb")  tmp = 11.35 * 30.;
0508   else if (element == "U")   tmp = 18.95 * 20.;
0509   return tmp;
0510 }
0511 
0512 double atomicWt(char element[2]) {
0513 
0514   double tmp=0;
0515   if      (element == "H")   tmp = 1.00794;
0516   else if (element == "Be")  tmp = 9.0122;
0517   else if (element == "C")   tmp = 12.011;
0518   else if (element == "Al")  tmp = 26.98;
0519   else if (element == "Ti")  tmp = 47.88;
0520   else if (element == "Fe")  tmp = 55.85;
0521   else if (element == "Cu")  tmp = 63.546;
0522   else if (element == "Nb")  tmp = 92.906;
0523   else if (element == "Cd")  tmp = 112.41;
0524   else if (element == "Sn")  tmp = 118.69;
0525   else if (element == "Ta")  tmp = 180.9479;
0526   else if (element == "Au")  tmp = 196.97;
0527   else if (element == "Pb")  tmp = 207.19;
0528   else if (element == "U")   tmp = 238.03;
0529   return tmp;
0530 }
0531 
0532 std::vector<double> angleScan(int scan) {
0533 
0534   static double deg = 3.1415926/180.;
0535   std::vector<double> tmp;
0536   if (scan <= 1) {
0537     tmp.push_back(59.1*deg);
0538     tmp.push_back(89.0*deg);
0539     tmp.push_back(119.0*deg);
0540     tmp.push_back(159.6*deg);
0541   } else {
0542     tmp.push_back(10.1*deg);
0543     tmp.push_back(15.0*deg);
0544     tmp.push_back(19.8*deg);
0545     tmp.push_back(24.8*deg);
0546     tmp.push_back(29.5*deg);
0547     tmp.push_back(34.6*deg);
0548     tmp.push_back(39.6*deg);
0549     tmp.push_back(44.3*deg);
0550     tmp.push_back(49.3*deg);
0551     tmp.push_back(54.2*deg);
0552     tmp.push_back(59.1*deg);
0553     tmp.push_back(64.1*deg);
0554     tmp.push_back(69.1*deg);
0555     tmp.push_back(74.1*deg);
0556     tmp.push_back(79.1*deg);
0557     tmp.push_back(84.1*deg);
0558     tmp.push_back(89.0*deg);
0559     tmp.push_back(98.9*deg);
0560     tmp.push_back(108.9*deg);
0561     tmp.push_back(119.0*deg);
0562     tmp.push_back(129.1*deg);
0563     tmp.push_back(139.1*deg);
0564     tmp.push_back(149.3*deg);
0565     tmp.push_back(159.6*deg);
0566     tmp.push_back(161.4*deg);
0567     tmp.push_back(165.5*deg);
0568     tmp.push_back(169.5*deg);
0569     tmp.push_back(173.5*deg);
0570     tmp.push_back(177.0*deg);
0571   }
0572   std::cout << "Scan " << tmp.size() << " angular regions:\n";
0573   for (unsigned int i=0; i<tmp.size(); i++) {
0574     std::cout << tmp[i]/deg;
0575     if (i == tmp.size()-1) std::cout << " degrees\n";
0576     else                   std::cout << ", ";
0577   }
0578   return tmp;
0579 }
0580 
0581 std::vector<double> energyScan(char part[2]) {
0582 
0583   std::vector<double> tmp;
0584   if (part == "n" || part == "N") {
0585     tmp.push_back(0.01);
0586     tmp.push_back(0.03);
0587     tmp.push_back(0.05);
0588   }
0589   tmp.push_back(0.07);
0590   tmp.push_back(0.09);
0591   tmp.push_back(0.11);
0592   tmp.push_back(0.13);
0593   tmp.push_back(0.15);
0594   tmp.push_back(0.17);
0595   tmp.push_back(0.19);
0596   tmp.push_back(0.21);
0597   tmp.push_back(0.23);
0598   tmp.push_back(0.25);
0599 
0600   std::cout << "Scan " << tmp.size() << " Energy regions:\n";
0601   for (unsigned int i=0; i<tmp.size(); i++) {
0602     std::cout << tmp[i];
0603     if (i == tmp.size()-1) std::cout << " GeV\n";
0604     else                   std::cout << ", ";
0605   }
0606   return tmp;
0607 }
0608 
0609 std::vector<double> rapidityScan() {
0610 
0611   std::vector<double> tmp;
0612   tmp.push_back(0.60);
0613   tmp.push_back(0.80);
0614   tmp.push_back(1.00);
0615   tmp.push_back(1.20);
0616   tmp.push_back(1.40);
0617   tmp.push_back(1.60);
0618   tmp.push_back(1.80);
0619   tmp.push_back(2.00);
0620   tmp.push_back(2.20);
0621   tmp.push_back(2.40);
0622   tmp.push_back(2.60);
0623   tmp.push_back(2.80);
0624   tmp.push_back(3.00);
0625 
0626   std::cout << "Scan " << tmp.size() << " rapidity regions:\n";
0627   for (unsigned int i=0; i<tmp.size(); i++) {
0628     std::cout << tmp[i];
0629     if (i == tmp.size()-1) std::cout << ";\n";
0630     else                   std::cout << ", ";
0631   }
0632   return tmp;
0633 }
0634 
0635 int type(double mass) {
0636 
0637   double m  = (mass >=0 ? mass: -mass);
0638   int    tmp=0;
0639   if      (m < 0.01)   {tmp = 1;}
0640   else if (m < 1.00)   {if (mass < 0) tmp = 2; else tmp = 3;}
0641   else if (m < 115.00) {if (mass < 0) tmp = 4; else tmp = 5;}
0642   else if (m < 135.00) tmp = 6;
0643   else if (m < 140.00) {if (mass < 0) tmp = 7; else tmp = 8;}
0644   else if (m < 495.00) {if (mass < 0) tmp = 9; else tmp = 10;}
0645   else if (m < 500.00) tmp = 11;
0646   else if (m < 938.50) {if (mass < 0) tmp = 12; else tmp = 13;}
0647   else if (m < 940.00) tmp = 14;
0648   else if (m < 1850.0) {tmp = 15;}
0649   else                 {tmp = 16;}
0650   //  std::cout << "Mass " << mass << " type " << tmp << "\n";
0651   return tmp;
0652 }