File indexing completed on 2024-04-06 12:30:17
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
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
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
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
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
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) {
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
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
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
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
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) {
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
0651 return tmp;
0652 }
