Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
simRho.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file simRho.cc
1 #include <iostream>
2 #include <fstream>
3 #include "sHelix.h"
4 #include "sChargeMap.h"
5 
6 #include "TStyle.h"
7 #include "TCanvas.h"
8 #include "TMath.h"
9 #include "TRandom3.h"
10 #include "TFile.h"
11 #include "TH1D.h"
12 #include "TH2D.h"
13 #include "TH3D.h"
14 
15 const int npid = 12+1;
16 int apid[npid] = {11, 13, // 2
17  211, // 1
18  321, // 1
19  411, 431, // 2
20  511, // 1
21  2212, // 1
22  3222, 3112, 3312, 3334, // 6
23  0};
24 int acha[npid] = {-1,-1,+1,+1,+1,+1,+1,+1,+1,-1,-1,-1};
25 TString spid[npid] = {"e^{-}",
26  "#mu^{-}",
27  "#pi^{+}",
28  "K^{+}",
29  "D^{+}",
30  "D^{+}_{s}",
31  "B^{+}",
32  "p",
33  "#Sigma^{+}",
34  "#Sigma^{-}",
35  "#Xi^{-}",
36  "#Omega^{-}",
37 
38  "Neutral"};
39 TH1F *hpid = new TH1F("hpid","hpid",npid,-0.5,npid-0.5);
40 TH1F *hrpass = new TH1F("hrpass","hrpass; r [cm]",100,-15,+15);
41 TH1F *hzpass = new TH1F("hzpass","hzpass; z [cm]",100,-15,+15);
42 TH1F *hpidpass = new TH1F("hpidpass","hpidpass",npid,-0.5,npid-0.5);
43 TH1F *hptpass = new TH1F("hptpass","hptpass; pt [GeV]",100,0,10);
44 TH1F *hetapass = new TH1F("hetapass","hetapass; eta",100,-3,+3);
45 TH1F *hlengthpass=new TH1F("hlengthpass","hlengthpass; length [cm]",100,0,150);
46 TH1F *hbxing = new TH1F("hbxing","hbxing",150,-0.5,300-0.5);
47 TH3F *tracking = new TH3F("tracking","Generated tracks: tracking;xi;yi;zi",120,-100,+100,120,-100,+100,120,-100,+100);
48 
49 int main() {
50  for(int n=0; n!=npid; ++n) hpid->GetXaxis()->SetBinLabel(1+n,spid[n].Data());
51  for(int n=0; n!=npid; ++n) hpidpass->GetXaxis()->SetBinLabel(1+n,spid[n].Data());
52 
53  gStyle->SetOptStat(111111111);
54  //================
55  // Gas parameters
56  //float gs_wf_me = 35; // [eV] mean energy for ionisation
57  //float gs_density = 1e-3; // [gr/cm^3]
58  //float gs_np = 14.35; // [cm^-1]
59  float gs_nt = 47.80; // [cm^-1]
60  float gs_ion_mobility = 4; // [cm^2/(V.s)]
61 
62  //================
63  // TPC running conditions
64  float tpc_magnetic_field = 1.5; // [Tesla]
65  float tpc_electric_field = 400; // [V/cm]
66  float tpc_drift_velocity_e = 8; // [cm/us]
67  float tpc_drift_velocity_i = tpc_electric_field*gs_ion_mobility*1e-3; // [cm/ms]
68  int lumi = 189; // 1/(189*106e-9) => 49915.1 Hz
69 
70  // input file
71  std::ifstream inputf("/Users/cperez//ampt/ana/ampt.dat");
72 
73  sChargeMap *map = new sChargeMap(100,30,80,
74  1,0,TMath::TwoPi(),
75  4*int(80/(tpc_drift_velocity_i*(lumi*106e-6)))+1,80,
76  tpc_drift_velocity_e,tpc_drift_velocity_i);
77  // main routine
78  int ncoll;
79  for(ncoll=0; ncoll!=3; ++ncoll) {
80  // =========
81  // COLLISION
82  if(ncoll%100==0) std::cout << ncoll << std::endl;
83  // event wise
84  int eventno, testno, nhadrons, npartA, npartB, npartElaA, npartIneA, npartElaB, npartIneB;
85  float impactPar;
86  // hadron wise
87  int pid;
88  float px, py, pz, mass, x, y, z, t;
89  inputf >> eventno;
90  if(!inputf.good()) {
91  inputf.close();
92  inputf.open("ampt.dat");
93  inputf >> eventno;
94  }
95  inputf >> testno >> nhadrons >> impactPar >> npartA >> npartB >> npartElaA >> npartIneA >> npartElaB >> npartIneB;
96  for(int np=0; np!=nhadrons; ++np) {
97  inputf >> pid >> px >> py >> pz >> mass >> x >> y >> z >> t;
98  if( TMath::AreEqualAbs(t,0,1e-3) ) continue;
99  x *= 1e-13;
100  y *= 1e-13;
101  z *= 1e-13;
102  int fill = npid-1;
103  for(int n=0; n!=npid-1; ++n) if(apid[n]==TMath::Abs(pid)) fill = n;
104  hpid->Fill( fill );
105  if(fill==npid-1) continue;
106  int q;
107  for(int n=0; n!=npid-1; ++n) if(apid[n]==TMath::Abs(pid)) q = acha[n];
108  sHelix a(x,y,z,px,py,pz,q,tpc_magnetic_field);
109  float t1 = a.findFirstInterceptTo(30,80);
110  if(t1>999) continue;
111  float t2 = a.findFirstInterceptTo(80,80);
112  if(t2>999) continue;
113  float length = a.s(t1,t2);
114  if(length<1) continue;
115  hrpass->Fill( sqrt(x*x+y*y) );
116  hzpass->Fill( z );
117  hpidpass->Fill( fill );
118  hptpass->Fill( sqrt(px*px+py*py) );
119  float pmom = sqrt(px*px+py*py+pz*pz);
120  float eta = 1.e30;
121  if (pmom != TMath::Abs(pz)) eta = 0.5*TMath::Log((pmom+pz)/(pmom-pz));
122  hetapass->Fill( eta );
123  hlengthpass->Fill( length );
124  float nprim = gRandom->Poisson(gs_nt*length);
125  float track[100][3];
126  a.breakIntoPieces(t1,t2,track);
127  for(int i=0; i!=100; ++i) {
128  float ri = TMath::Sqrt( track[i][0]*track[i][0] + track[i][1]*track[i][1] );
129  if(ri<30-0.1) break;
130  if(ri>80+0.1) break;
131  float pi = 1;
132  map->Fill(ri, pi, track[i][2], -nprim/100);
133  map->Fill(ri, pi, track[i][2], nprim/100);
134  //if(ncoll==1&&
135  // (
136  // np==622||
137  // np==935
138  // )) {
139  tracking->Fill(track[i][0],track[i][1],track[i][2]);
140  // if(i==0) {
141  // std::cout << "pmom " << pmom << "| eta " << eta;
142  // std::cout << "| x y z " << x << " " << y << " " << z << " " << px << " " << py << " " << pz << std::endl;
143  // std::cout << "t " << t1 << "|" << t2 << std::endl;
144  //}
145  //std::cout << ri << " ";
146  //if(i==99)
147  // std::cout << std::endl;
148  //}
149  }
150  }
151 
152  // ==========
153  // DRIFT TIME
154  int steps = gRandom->Poisson( lumi );
155  hbxing->Fill(steps);
156  float lapse = 106*steps*1e-6; // [ns] -> [ms]
157  map->Propagate(lapse);
158  }
159  map->ScreenShot("map","root",ncoll);
160  inputf.close();
161  map->SaveIonMap("ionmap.root");
162  //map->SaveRho("mark1_0.root",100,1,1600);
163 
164  TFile *qa = new TFile("qa.root","RECREATE");
165  qa->WriteTObject(hpid,hpid->GetName());
166  qa->WriteTObject(hrpass,hrpass->GetName());
167  qa->WriteTObject(hzpass,hzpass->GetName());
168  qa->WriteTObject(hpidpass,hpidpass->GetName());
169  qa->WriteTObject(hptpass,hptpass->GetName());
170  qa->WriteTObject(hetapass,hetapass->GetName());
171  qa->WriteTObject(hlengthpass,hlengthpass->GetName());
172  qa->Close();
173  tracking->SaveAs("tracking.root","root");
174 
175  return 0;
176 }