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 "TH1I.h"
0018 #include "TCanvas.h"
0019 #include "TApplication.h"
0020 #include "TRefArray.h"
0021 #include "TStyle.h"
0022 #include "TGraph.h"
0023 
0024 static unsigned int antiproton=12, proton=13, neutron=14, heavy=15, ions=16;
0025 
0026 void AnalyseH2TB(char element[6], char list[20], char ene[6], char part[4], int sav=0, int nMax=-1, bool debug=false) {
0027 
0028   char *g4ver = "9.2.ref01P";
0029   bool detail = true;
0030 
0031   int  energy = atoi(ene);
0032   char fname[120];
0033   sprintf (fname, "%s%s_%s_%sGeV.root", element, list, part, ene);
0034   char ofile[130];
0035   sprintf (ofile, "histo/histo_%s", fname);
0036 
0037   double rhol = rhoL(element);
0038   double atwt = atomicWt(element);
0039   std::vector<std::string> types = types();
0040   std::cout << fname << " rhoL " << rhol << " atomic weight " << atwt << "\n";
0041   
0042   TFile *fout = new TFile(ofile, "recreate");
0043   TH1F *hiKE0[20], *hiKE1[20], *hiKE2[20], *hiCT0[20], *hiCT1[20], *hiCT2[20];
0044   TH1I *hiMulti[20];
0045   TH1F *hiParticle[5][20], *hiTotalKE[20], *hiMomInclusive[20];
0046   TH1F *hiSumP[20], *hiPT2[20], *hiEP2[4], *baryon1, *baryon2;
0047   TH1F *hProton[2], *hNeutron[2], *hHeavy[2], *hIon[2], *hBaryon[2];;
0048   char name[80], title[180], ctype[20], ytitle[20];
0049   double xbin;
0050 
0051   for (unsigned int ii=0; ii<=(types.size()); ii++) {
0052     if      (ii == 0) sprintf (ctype, "All Particles");
0053     else              sprintf (ctype, "%s", types[ii-1].c_str());
0054     
0055     sprintf (title, "%s GeV %s on %s: Inclusive Mom Dist of %s (%s %s)", ene, part, element, ctype, g4ver, list);
0056     sprintf (name, "hiMomInclusive_%s%s%sGeV(%s)", element, list, ene, ctype);
0057     hiMomInclusive[ii] = new TH1F (name, title, 6000, 0., 300.);
0058 
0059     for (unsigned int jj=0; jj<5; jj++) {
0060       sprintf (title, "Particle %i : %s in %s at %s GeV (%s)",jj, ctype, element, ene, list);
0061       sprintf (name, "Particle%i_KE%s%s%sGeV(%s)",jj,element, list, ene, ctype);
0062       if (ii==ions) hiParticle[jj][ii] = new TH1F (name, title, 50000, 0., 10.);
0063       else          hiParticle[jj][ii] = new TH1F (name, title, 15500, 0., 310.);
0064       if (debug) std::cout << "hiParticle[" << jj << "][" << ii << "] = " << hiParticle[jj][ii] << " " <<  name << " Particle KE Energy  " << title << "\n"; 
0065     }
0066     if (debug) std::cout << "Ctype " << ctype << "\n";
0067 
0068     sprintf (title, "%s in %s at %s GeV (%s)", ctype, element, ene, list);
0069     sprintf (name, "KE0%s%s%sGeV(%s)", element, list, ene, ctype);
0070     hiKE0[ii] = new TH1F (name, title, 15000, 0., 300.);
0071     hiKE0[ii]->GetXaxis()->SetTitle("Kinetic Energy (GeV)");
0072     xbin = hiKE0[ii]->GetBinWidth(1);
0073     sprintf (ytitle, "Events/%6.3f GeV", xbin);
0074     hiKE0[ii]->GetYaxis()->SetTitle(ytitle);
0075     if (debug) std::cout << "hiKE0[" << ii << "] = " << hiKE0[ii] << " " <<  name << " KE Energy  " << title << "\n";
0076 
0077     sprintf (name, "CT0%s%s%sGeV(%s)", element, list, ene, ctype);
0078     hiCT0[ii] = new TH1F (name, title, 100, -1.0, 1.0.);
0079     hiCT0[ii]->GetXaxis()->SetTitle("cos (#theta)");
0080     xbin = hiCT0[ii]->GetBinWidth(1);
0081     sprintf (ytitle, "Events/%6.3f", xbin);
0082     hiCT0[ii]->GetYaxis()->SetTitle(ytitle);
0083     if (debug) std::cout << "hiCT0[" << ii << "] = " << hiCT0[ii] << " " <<  name << " cos(T#eta) " << title << "\n";
0084 
0085     sprintf (title, "%s (Elastic) in %s at %s GeV (%s)", ctype, element, ene, list);
0086     sprintf (name, "KE1%s%s%sGeV(%s)", element, list, ene, ctype);
0087     hiKE1[ii] = new TH1F (name, title, 15000, 0., 300.);
0088     hiKE1[ii]->GetXaxis()->SetTitle("Kinetic Energy (GeV)");
0089     xbin = hiKE1[ii]->GetBinWidth(1);
0090     sprintf (ytitle, "Events/%6.3f GeV", xbin);
0091     hiKE1[ii]->GetYaxis()->SetTitle(ytitle);
0092     if (debug) std::cout << "hiKE1[" << ii << "] = " << hiKE1[ii] << " " <<  name << " KE Energy  " << title << "\n";
0093 
0094     sprintf (name, "CT1%s%s%sGeV(%s)", element, list, ene, ctype);
0095     hiCT1[ii] = new TH1F (name, title, 100, -1.0, 1.0.);
0096     hiCT1[ii]->GetXaxis()->SetTitle("cos (#theta)");
0097     xbin = hiCT1[ii]->GetBinWidth(1);
0098     sprintf (ytitle, "Events/%6.3f", xbin);
0099     hiCT1[ii]->GetYaxis()->SetTitle(ytitle);
0100     if (debug) std::cout << "hiCT1[" << ii << "] = " << hiCT1[ii] << " " <<  name << " cos(T#eta) " << title << "\n";
0101 
0102     sprintf (title, "%s (InElastic) in %s at %s GeV (%s)", ctype, element, ene, list);
0103     sprintf (name, "KE2%s%s%sGeV(%s)", element, list, ene, ctype);
0104     hiKE2[ii] = new TH1F (name, title, 15000, 0., 300.);
0105     hiKE2[ii]->GetXaxis()->SetTitle("Kinetic Energy (GeV)");
0106     xbin = hiKE2[ii]->GetBinWidth(1);
0107     sprintf (ytitle, "Events/%6.3f GeV", xbin);
0108     hiKE2[ii]->GetYaxis()->SetTitle(ytitle);
0109     if (debug) std::cout << "hiKE2[" << ii << "] = " << hiKE2[ii] << " " <<  name << " KE Energy  " << title << "\n";
0110 
0111     sprintf (name, "CT2%s%s%sGeV(%s)", element, list, ene, ctype);
0112     hiCT2[ii] = new TH1F (name, title, 100, -1.0, 1.0.);
0113     hiCT2[ii]->GetXaxis()->SetTitle("cos (#theta)");
0114     xbin = hiCT2[ii]->GetBinWidth(1);
0115     sprintf (ytitle, "Events/%6.3f", xbin);
0116     hiCT2[ii]->GetYaxis()->SetTitle(ytitle);
0117     if (debug) std::cout << "hiCT2[" << ii << "] = " << hiCT2[ii] << " " <<  name << " cos(T#eta) " << title << "\n";
0118 
0119     sprintf (name, "PT2%s%s%sGeV(%s)", element, list, ene, ctype);
0120     hiPT2[ii] = new TH1F (name, title, 15000, 0.0, 3000.0.);
0121     hiPT2[ii]->GetXaxis()->SetTitle("p_T (GeV)");
0122     xbin = hiCT2[ii]->GetBinWidth(1);
0123     sprintf (ytitle, "Events/(%6.3f GeV)", xbin);
0124     hiPT2[ii]->GetYaxis()->SetTitle(ytitle);
0125     if (debug) std::cout << "hiPT2[" << ii << "] = " << hiPT2[ii] << " " <<  name << " pT " << title << "\n";
0126 
0127     sprintf (name, "Multi%s%s%sGeV(%s)", element, list, ene, ctype);
0128     sprintf (title,"%s multiplicity in %s at %s GeV (%s)", ctype, element, ene, list);
0129     hiMulti[ii] = new TH1I (name, title, 101, -1, 100);
0130     hiMulti[ii]->GetXaxis()->SetTitle("Multiplicity");
0131     if (debug) std::cout << "hiMulti[" << ii << "] = " << hiMulti[ii] << " " <<  name << " Multiplicity\n";
0132 
0133     sprintf (name, "TotalKE%s%s%sGeV(%s)", element, list, ene, ctype);
0134     sprintf (title,"%s (inelastic) in %s at %s GeV (%s)", ctype, element, ene, list);
0135     hiTotalKE[ii] = new TH1F (name, title, 15500, 0, 310);
0136     sprintf (title, "Total KE carried by %s", ctype);
0137     hiTotalKE[ii]->GetXaxis()->SetTitle(title);
0138     if (debug) std::cout << "hiTotalKE[" << ii << "] = " << hiTotalKE[ii] << " " <<  name << " " << title << "\n";
0139 
0140     sprintf(name, "hiSumMomentum%s%s%sGeV(%s)", element, list, ene, ctype);
0141     sprintf (title, "%s GeV %s on %s: Total Mom carried by %s (%s %s)", ene, part, element, ctype, g4ver, list);
0142     hiSumP[ii] = new TH1F (name, title, 6000, 0., 310.);
0143     if (debug) std::cout << "hiSumP[" << ii << "] = " << hiSumP[ii] << " " <<  name << " " << title << "\n";
0144   }
0145   hiEP2[0] = new TH1F ("sumPX", "Sum px", 2000., -100., 100.);
0146   hiEP2[0]->GetXaxis()->SetTitle("Momentum balance (x)");
0147   hiEP2[0]->GetYaxis()->SetTitle("Events");
0148   hiEP2[1] = new TH1F ("sumPY", "Sum py", 2000., -100., 100.);
0149   hiEP2[1]->GetXaxis()->SetTitle("Momentum balance (y)");
0150   hiEP2[1]->GetYaxis()->SetTitle("Events");
0151   hiEP2[2] = new TH1F ("sumPZ", "Sum pz", 2000., -100., 100.);
0152   hiEP2[2]->GetXaxis()->SetTitle("Momentum balance (z)");
0153   hiEP2[2]->GetYaxis()->SetTitle("Events");
0154   hiEP2[3] = new TH1F ("sumE",  "Sum E",  2000., -100., 100.);
0155   hiEP2[3]->GetXaxis()->SetTitle("Energy balance");
0156   hiEP2[3]->GetYaxis()->SetTitle("Events");
0157 
0158   if (detail) {
0159     for (int i=0; i<2; i++) {
0160       sprintf(title, "proton%i_%s%s%sGeV(%s)", i, element, list, ene, ctype);
0161       hProton[i] = new TH1F(title, title, 15500, 0., 310.);
0162       hProton[i]->GetXaxis()->SetTitle("Kinetic Energy (GeV)");
0163 
0164       sprintf(title, "neutron%i_%s%s%sGeV(%s)", i, element, list, ene, ctype);
0165       hNeutron[i] = new TH1F(title, title, 15500, 0., 310.);
0166       hNeutron[i]->GetXaxis()->SetTitle("Kinetic Energy (GeV)");
0167 
0168       sprintf(title, "heavy%i_%s%s%sGeV(%s)", i, element, list, ene, ctype);
0169       hHeavy[i] = new TH1F(title, title, 15500, 0., 310.);
0170       hHeavy[i]->GetXaxis()->SetTitle("Kinetic Energy (GeV)");
0171 
0172       sprintf(title, "ion%i_%s%s%sGeV(%s)", i, element, list, ene, ctype);
0173       hIon[i] = new TH1F(title, title, 50000, 0., 10.);
0174       hIon[i]->GetXaxis()->SetTitle("Kinetic Energy (GeV)");
0175 
0176       sprintf(title, "baryon%i_%s%s%sGeV(%s)", i, element, list, ene, ctype);
0177       hBaryon[i] = new TH1F(title, title, 15500, 0., 310.);
0178       hBaryon[i]->GetXaxis()->SetTitle("Kinetic Energy (GeV)");
0179     }
0180 
0181     sprintf(title, "baryonX_%s%s%sGeV", element, list, ene);
0182     baryon1 = new TH1F("baryon1", title, 15500, 0., 310.);
0183     sprintf(title, "baryonY_%s%s%sGeV", element, list, ene);
0184     baryon2 = new TH1F("baryon2", title, 15500, 0., 310.);
0185   }
0186 
0187   //  sprintf(fname, "out/%s", fname);
0188   std::cout << "Reading from " << fname << std::endl;
0189 
0190   TFile *file = new TFile(fname);
0191   TTree *tree = (TTree *) file->Get("T1");
0192   
0193   if (!tree) {
0194     std::cout << "Cannot find Tree T1 in file " << fname << "\n";
0195   } else {
0196     std::cout << "Tree T1 found with " << tree->GetEntries() << " entries\n";
0197     int nentry = tree->GetEntries();
0198     int ninter=0, elastic=0, inelastic=0;
0199     if (nMax > 0 && nMax < nentry) nentry = nMax;
0200 
0201     for (int i=0; i<nentry; i++) {
0202       if (i%1000 == 0) std::cout << "Start processing event " << i << "\n";
0203       std::vector<int>                     *nsec, *procids;
0204       std::vector<double>                  *px, *py, *pz, *mass;
0205       std::vector<std::string>             *procs;
0206       tree->SetBranchAddress("NumberSecondaries", &nsec);
0207       tree->SetBranchAddress("ProcessID",         &procids);
0208       tree->SetBranchAddress("SecondaryPx",       &px);
0209       tree->SetBranchAddress("SecondaryPy",       &py);
0210       tree->SetBranchAddress("SecondaryPz",       &pz);
0211       tree->SetBranchAddress("SecondaryMass",     &mass);
0212       tree->GetEntry(i);
0213 
0214       
0215       if ((*nsec).size() > 0) {
0216 
0217     if (debug) std::cout << "Entry " << i << " no. of secondaries " << (*nsec)[0] << "\n";
0218     ninter++;
0219     bool isItElastic = false;
0220     if ((*procids)[0] == 17) {elastic++; isItElastic = true;}
0221     else                     inelastic++;
0222     
0223     if (ninter <3) {
0224       std::cout << "Interaction " << ninter << "/" << i+1 << " Type "
0225             << (*procids)[0]  << " with " << (*nsec)[0] << " secondaries\n";
0226       for (int k=0; k<(*nsec)[0]; k++)
0227         std::cout << " Secondary " << k << " Px " << (*px)[k] << " Py " << (*py)[k] << " Pz " << (*pz)[k] << " Mass " << (*mass)[k] << "\n";
0228     }
0229     
0230     std::list<double> pProton, pNeutron, pIon, pHeavy;
0231     int counter[20];
0232     double sumKE[20], sumPx[20], sumPy[20], sumPz[20];
0233     for(unsigned int nct=0; nct<=types.size(); nct++) {
0234       counter[nct] = 0;   sumKE[nct] = 0;
0235       sumPx[nct]   = 0.0; sumPy[nct] = 0.0; sumPz[nct] = 0.0;
0236     }
0237 
0238     double sumpx=0, sumpy=-energy, sumpz=0, sume=-energy;
0239     int num = (*nsec)[0];
0240     if (debug) std::cout << "Secondaries: " << num << "\n";
0241     std::vector<double> partEne(num,0.0); 
0242     std::vector<int>    partType(num,999);
0243     double              kBaryon=0;
0244 
0245     for (int k=0; k<(*nsec)[0]; k++) {
0246       int type = type((*mass)[k]);
0247       double m  = (((*mass)[k]) >=0 ? ((*mass)[k]): -((*mass)[k]));
0248       m        /= 1000.;
0249       double pl = ((*py)[k])/1000.;
0250       double p1 = ((*px)[k])/1000.;
0251       double p2 = ((*pz)[k])/1000.;
0252       double pt = (p1*p1 + p2*p2);
0253       double pp = (pt+pl*pl);
0254       double ke = (sqrt (pp + m*m) - m);
0255       pp        = sqrt (pp);
0256       double cth= (pp == 0. ? -2. : (pl/pp));
0257       pt        = sqrt(pt);
0258 
0259       if      (type == proton)  pProton.push_back(ke);
0260       else if (type == neutron) pNeutron.push_back(ke);
0261       else if (type == heavy)   pHeavy.push_back(ke);
0262       else if (type == ions)    pIon.push_back(ke);
0263       if (type == proton || type == neutron || type == heavy) kBaryon += ke;
0264 
0265       sumpx += p1;
0266       sumpy += pl;
0267       sumpz += p2;
0268       sume  += (ke+m);
0269       partEne[k] = ke; partType[k] = type;
0270       hiMomInclusive[0]->Fill(pp);
0271       hiMomInclusive[type]->Fill(pp);
0272       if (debug) std::cout << "Entry " << i << " Secondary " << k << " Mass " << (*mass)[k] << " Type " << type << " Cth " << cth << " KE " << ke << "\n";
0273       hiKE0[0]->Fill(ke);
0274       hiCT0[0]->Fill(cth);
0275       hiKE0[type]->Fill(ke);
0276       hiCT0[type]->Fill(cth);
0277       if (isItElastic) {
0278         hiKE1[0]->Fill(ke);
0279         hiCT1[0]->Fill(cth);
0280         hiKE1[type]->Fill(ke);
0281         hiCT1[type]->Fill(cth);
0282       } else {
0283         hiKE2[0]->Fill(ke);
0284         hiCT2[0]->Fill(cth);
0285         hiPT2[0]->Fill(pt);
0286         hiKE2[type]->Fill(ke);
0287         hiCT2[type]->Fill(cth);
0288         hiPT2[type]->Fill(pt);
0289         counter[0]    += 1;
0290         counter[type] += 1;
0291         sumKE[0]      += ke;
0292         sumKE[type]   += ke;
0293         sumPx[0]      += (*px)[k];
0294         sumPy[0]      += (*py)[k];
0295         sumPz[0]      += (*pz)[k];
0296         sumPx[type]   += (*px)[k];
0297         sumPy[type]   += (*py)[k];
0298         sumPz[type]   += (*pz)[k];
0299       }
0300     } // loop over particles
0301     
0302     if (debug) std::cout << "Entry " << i << "  :: Elastic (?) " << isItElastic << "\n";
0303 
0304     if( !isItElastic ) {
0305       for (unsigned int nct=0; nct<=(types.size()); nct++) {
0306         double sumP = std::sqrt(sumPx[nct]*sumPx[nct] + sumPy[nct]*sumPy[nct] + sumPz[nct]*sumPz[nct]);
0307         if(debug) {
0308           if(nct<1) sdt::cout <<  sumP << "   "; 
0309           else      sdt::cout << types[nct-1] << " " <<  sumP << "   "; 
0310         }
0311         hiSumP[nct]->Fill( sumP );
0312 
0313         hiMulti[nct]->Fill(counter[nct]);
0314         hiTotalKE[nct]->Fill(sumKE[nct]);
0315         if (debug) {
0316           for (unsigned int nct=0; nct<=types.size(); nct++) 
0317         std::cout << "  [" << nct <<"]:" << counter[nct] << " KE " << sumKE[nct] << "\n";
0318         }
0319         hiEP2[0]->Fill(sumpx);
0320         hiEP2[1]->Fill(sumpy);
0321         hiEP2[2]->Fill(sumpz);
0322         hiEP2[3]->Fill(sume);
0323       }
0324 
0325       if (detail) {
0326         list<double>::iterator iter;
0327         double kMaxB=0;
0328         if (pProton.size() > 0) {
0329           pProton.sort();
0330           iter = pProton.end(); iter--;
0331           double pMax = *iter;
0332           kMaxB = pMax;
0333           hProton[0]->Fill(pMax);
0334           double pNext= 0;
0335           if (pProton.size() > 1) {
0336         iter--; pNext = *iter;
0337         hProton[1]->Fill(pNext);
0338           }
0339           if (debug) std::cout << "Proton " << pProton.size() << " " << pMax << " " << pNext << "\n";
0340         }
0341 
0342         if (pNeutron.size() > 0) {
0343           pNeutron.sort();
0344           iter = pNeutron.end(); iter--;
0345           double pMax = *iter;
0346           if (pMax > kMaxB) kMaxB = pMax;
0347           hNeutron[0]->Fill(pMax);
0348           double pNext= 0;
0349           if (pNeutron.size() > 1) {
0350         iter--; pNext = *iter;
0351         hNeutron[1]->Fill(pNext);
0352           }
0353           if (debug) std::cout << "Neutron " << pNeutron.size() << " " << pMax << " " << pNext << "\n";
0354         }
0355 
0356         if (pHeavy.size() > 0) {
0357           pHeavy.sort();
0358           iter = pHeavy.end(); iter--;
0359           double pMax = *iter;
0360           if (pMax > kMaxB) kMaxB = pMax;
0361           hHeavy[0]->Fill(pMax);
0362           double pNext= 0;
0363           if (pHeavy.size() > 1) {
0364         iter--; pNext = *iter;
0365         hHeavy[1]->Fill(pNext);
0366           }
0367         }
0368 
0369         if (kMaxB > 0) {
0370           hBaryon[0]->Fill(kMaxB);
0371           hBaryon[1]->Fill(kBaryon-kMaxB);
0372           if (debug) std::cout << "Baryon " << kMaxB << " " << kBaryon-kMaxB << "\n";
0373         }
0374 
0375         if (pIon.size() > 0) {
0376           pIon.sort();
0377           iter = pIon.end(); iter--;
0378           double pMax = *iter;
0379           hIon[0]->Fill(pMax);
0380           double pNext= 0;
0381           if (pIon.size() > 1) {
0382         iter--; pNext = *iter;
0383         hIon[1]->Fill(pNext);
0384           }
0385         }
0386       
0387         for(int ipart=0; ipart<partEne.size(); ipart++){
0388           for(int jpart=ipart+1; jpart<partEne.size(); jpart++){
0389         if(partEne[ipart] < partEne[jpart]){
0390           double tempE = partEne[ipart];    int tempI       = partType[ipart];
0391           partEne[ipart] = partEne[jpart];  partType[ipart] = partType[jpart];
0392           partEne[jpart] = tempE;           partType[jpart] = tempI;
0393         }
0394           }
0395         }
0396 
0397         int nPart[20]; for(unsigned int ii=0; ii<20; ii++) nPart[ii]=0;
0398         int nbaryon=0; 
0399         bool firstBaryon  = false;
0400         bool secondBaryon = false;
0401         bool first = true;
0402         for(int unsigned ii=0; ii<partEne.size(); ii++){
0403           nPart[partType[ii]] += 1;
0404           if(partType[ii]==proton || partType[ii]==neutron){
0405         if (first)     baryon1->Fill(partEne[ii]);
0406         else           baryon2->Fill(partEne[ii]);
0407         first = false;
0408           }
0409           
0410           if(partType[ii]==antiproton || partType[ii]==proton || partType[ii]==neutron || partType[ii]==heavy || partType[ii]==ions) {
0411         nbaryon++;
0412         if(nbaryon == 1) firstBaryon  = true;
0413         if(nbaryon == 2) secondBaryon = true;
0414           }
0415 
0416           unsigned int ip = partType[ii];
0417           if(firstBaryon) {
0418         if (debug && partEne[ii]>200) std::cout << "nbaryon " << nbaryon << " ii " << ii << "  partType " << partType[ii] << " " << partEne[ii] << "\n";
0419         if(nPart[ip]==1) hiParticle[0][ip]->Fill(partEne[ii]);
0420           }
0421 
0422           firstBaryon = false;
0423           if(secondBaryon){
0424         if (debug && partEne[ii]>200) std::cout << "nbaryon " << nbaryon << " ii " << ii << "  partType " << partType[ii] << " " << partEne[ii] << "\n";
0425         hiParticle[1][ip]->Fill(partEne[ii]);
0426           }
0427           secondBaryon = false;
0428         } 
0429       }
0430 
0431     } // ifInElastic
0432       }
0433     } // loop over entries
0434 
0435     gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite);
0436     gStyle->SetPadColor(kWhite);    gStyle->SetFrameBorderMode(0);
0437     gStyle->SetFrameBorderSize(1);  gStyle->SetFrameFillColor(0);
0438     gStyle->SetFrameFillStyle(0);   gStyle->SetFrameLineColor(1);
0439     gStyle->SetFrameLineStyle(1);   gStyle->SetFrameLineWidth(1);
0440     gStyle->SetOptLogy(1);          gStyle->SetTitleOffset(1.2,"Y");
0441 
0442     if (sav < 0) {
0443       TCanvas *cc1[20], *cc2[20], *cc3; 
0444       cc3 = new TCanvas("c_multiplicity", "c_multiplicity", 800, 800);
0445       TLegend *leg = new TLegend(0.5, 0.5, 0.8, 0.8);
0446       for (unsigned int iia=0; iia<=(types.size()); iia++) {
0447     if      (iia == 0) sprintf (ctype, "All Particles");
0448     else               sprintf (ctype, "%s", types[iia-1].c_str());
0449 
0450     sprintf (title, "%s in %s at %s GeV (%s)", ctype, element, ene, list);
0451     sprintf(name, "C-KE%i", iia);
0452     cc1[iia] = new TCanvas(name,title,800,600); cc1[iia]->Divide(2,2);
0453     cc1[iia]->cd(1); if (hiKE0[iia]->GetEntries() > 0) hiKE0[iia]->Draw(); 
0454     cc1[iia]->cd(3); if (hiKE1[iia]->GetEntries() > 0) hiKE1[iia]->Draw();
0455     cc1[iia]->cd(4); if (hiKE2[iia]->GetEntries() > 0) hiKE2[iia]->Draw(); 
0456     sprintf(name, "C-CT%i", iia);
0457     cc2[iia] = new TCanvas(name,title,800,600); cc2[iia]->Divide(2,2);
0458     cc2[iia]->cd(1); if (hiCT0[iia]->GetEntries() > 0) hiCT0[iia]->Draw(); 
0459     cc2[iia]->cd(3); if (hiCT1[iia]->GetEntries() > 0) hiCT1[iia]->Draw();
0460     cc2[iia]->cd(4); if (hiCT2[iia]->GetEntries() > 0) hiCT2[iia]->Draw(); 
0461     cc3->cd();
0462     hiMulti[iia]->SetLineColor(iia+1);
0463     if(iia>=9) hiMulti[iia]->SetLineColor(iia+2);
0464     if(iia==0) hiMulti[iia]->Draw();
0465     else       hiMulti[iia]->Draw("sames");
0466     leg->AddEntry(hiMulti[iia], title, "l");
0467       }
0468       cc3->cd();
0469       leg->Draw("same");
0470 
0471       if (detail) {
0472     TCanvas *cc4 = new TCanvas("c_Nucleon", "c_Nucleon", 800, 800);
0473     cc4->Divide(2,2);
0474     cc4->cd(1); hProton[0]->Draw();
0475     cc4->cd(2); hProton[1]->Draw();
0476     cc4->cd(3); hNeutron[0]->Draw();
0477     cc4->cd(4); hNeutron[1]->Draw();
0478 
0479     TCanvas *cc5 = new TCanvas("c_Ion", "c_Ion", 800, 800);
0480     cc5->Divide(2,2);
0481     cc5->cd(1); hHeavy[0]->Draw();
0482     cc5->cd(2); hHeavy[1]->Draw();
0483     cc5->cd(3); hIon[0]->Draw();
0484     cc5->cd(4); hIon[1]->Draw();
0485     
0486     TCanvas *cc6 = new TCanvas("c_Baryon", "c_Baryon", 800, 400);
0487     cc6->Divide(2,1);
0488     cc6->cd(1); hBaryon[0]->Draw();
0489     cc6->cd(2); hBaryon[1]->Draw();
0490       }
0491 
0492     } else {
0493 
0494       std::cout << "Writing histograms to " << ofile << "\n";
0495       fout->cd();
0496       fout->Write();
0497     }
0498 
0499     std::cout << ninter << " interactions seen in " << nentry << " trials\n"
0500           << "Elastic/Inelastic " << elastic << "/" << inelastic << "\n";
0501     if( nentry-ninter != 0 ) {
0502       double sigma = atwt*10000.*log((double)(nentry)/(double)(nentry-ninter))/(rhol*6.023);
0503       double dsigma    = sigma/sqrt(double(ninter));
0504       double sigmaEl   = sigma*((double)(elastic))/((double)(ninter));
0505       double dsigmaEl  = sigmaEl/sqrt(double(max(1,elastic)));
0506       double sigmaInel = sigma*((double)(inelastic))/((double)(ninter));
0507       double dsigmaInel= sigmaInel/sqrt(double(max(1,inelastic)));
0508       std::cout << "Total     " << sigma << " +- " << dsigma 
0509         << " mb (" << ninter << " events)\n"
0510         << "Elastic   " << sigmaEl<< " +- " << dsigmaEl
0511         << " mb (" << elastic << " events)\n"
0512         << "Inelasric " << sigmaInel << " +- " << dsigmaInel
0513         << " mb (" << inelastic << " events)\n";
0514     }
0515 
0516   } // if tree is found
0517 } //end analysis()
0518 
0519 double rhoL(char element[6]) {
0520   
0521   double tmp=0;
0522   if      (element == "Brass") tmp = 8.50 * 0.40;
0523   else if (element == "PbWO4") tmp = 8.28 * 0.30;
0524   else if (element == "Fe")    tmp = 7.87 * 0.30;
0525   else if (element == "H")     tmp = 0.0708 * 12.0;
0526   else if (element == "D")     tmp = 0.162  * 6.0;
0527   return tmp;
0528 }
0529 
0530 double atomicWt(char element[6]) {
0531   
0532   double tmp=0;
0533   if      (element == "Brass") tmp = 64.228;
0534   else if (element == "PbWO4") tmp = 455.036;
0535   else if (element == "Fe")    tmp = 55.85;
0536   else if (element == "H")     tmp = 1.0079;
0537   else if (element == "D")     tmp = 2.01;
0538   return tmp;
0539 }
0540 
0541 int type(double mass) {
0542 
0543   double m  = (mass >=0 ? mass: -mass);
0544   int    tmp=0;
0545   if      (m < 0.01)   {tmp = 1;}
0546   else if (m < 1.00)   {if (mass < 0) tmp = 2; else tmp = 3;}
0547   else if (m < 115.00) {if (mass < 0) tmp = 4; else tmp = 5;}
0548   else if (m < 135.00) tmp = 6;
0549   else if (m < 140.00) {if (mass < 0) tmp = 7; else tmp = 8;}
0550   else if (m < 495.00) {if (mass < 0) tmp = 9; else tmp = 10;}
0551   else if (m < 500.00) tmp = 11;
0552   else if (m < 938.50) {if (mass < 0) tmp = 12; else tmp = 13;}
0553   else if (m < 940.00) tmp = 14;
0554   else if (m < 1850.0) {tmp = 15;}
0555   else                 {tmp = 16;}
0556   //  std::cout << "Mass " << mass << " type " << tmp << "\n";
0557   return tmp;
0558 }
0559 
0560 std::vector<std::string> types() {
0561 
0562   std::vector<string> tmp;
0563   tmp.push_back("Photon/Neutrino");
0564   tmp.push_back("Electron");
0565   tmp.push_back("Positron");
0566   tmp.push_back("MuMinus");
0567   tmp.push_back("MuPlus");
0568   tmp.push_back("Pizero");
0569   tmp.push_back("Piminus");
0570   tmp.push_back("Piplus");
0571   tmp.push_back("Kminus");
0572   tmp.push_back("Kiplus");
0573   tmp.push_back("Kzero");
0574   tmp.push_back("AntiProton");
0575   tmp.push_back("Proton");
0576   tmp.push_back("Neutron/AntiNeutron");
0577   tmp.push_back("Heavy Hadrons");
0578   tmp.push_back("Ions");
0579 
0580   return tmp;
0581 }
0582