File indexing completed on 2023-10-25 10:04:23
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
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 }
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 }
0432 }
0433 }
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 }
0517 }
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
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
