Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TPCEventDisplay.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TPCEventDisplay.C
1 #include <cstdlib>
2 #include <iostream>
3 #include <map>
4 #include <string>
5 #include <vector>
6 
7 #include "TChain.h"
8 #include "TFile.h"
9 #include "TTree.h"
10 #include "TString.h"
11 #include "TObjString.h"
12 #include "TSystem.h"
13 #include "TROOT.h"
14 #include "TTreeReader.h"
15 #include "TTreeReaderValue.h"
16 #include "TTreeReaderArray.h"
17 #include "TVector3.h"
18 
19 #include <boost/format.hpp>
20 #include <boost/math/special_functions/sign.hpp>
21 
22 /*************************************************************/
23 /* TPC Raw Data Event Display */
24 /* Thomas Marshall,Aditya Dash */
25 /* rosstom@ucla.edu,aditya55@physics.ucla.edu */
26 /*************************************************************/
27 
28 TVector2 getBinPosition(int sector, int inFee, int channel, TTree* R1, TTree* R2, TTree* R3)
29 {
30  int mod_arr[26] = {2,2,1,1,1,3,3,3,3,3,3,2,2,1,2,2,1,1,2,2,3,3,3,3,3,3};
31  float slot[26] = {5,6,1,3,2,12,10,11,9,8,7,1,2,4,8,7,6,5,4,3,1,3,2,4,6,5};
32  TTreeReader* myReader;
33  if (mod_arr[inFee] == 1)
34  {
35  myReader = new TTreeReader(R1);
36  }
37  else if (mod_arr[inFee] == 2)
38  {
39  myReader = new TTreeReader(R2);
40  }
41  else if (mod_arr[inFee] == 3)
42  {
43  myReader = new TTreeReader(R3);
44  }
45  TTreeReaderValue<Float_t> FEE(*myReader, "FEE");
46  TTreeReaderValue<Int_t> FEE_Chan(*myReader, "FEE_Chan");
47  //TTreeReaderValue<Double_t> x(myReader, "x");
48  //TTreeReaderValue<Double_t> y(myReader, "y");
49  TTreeReaderValue<Double_t> phi(*myReader, "phi");
50  //TTreeReaderValue<Double_t> PadX(myReader, "PadX");
51  //TTreeReaderValue<Double_t> PadY(myReader, "PadY");
52  TTreeReaderValue<Double_t> PadR(*myReader, "PadR");
53  //TTreeReaderValue<Double_t> PadPhi(myReader, "PadPhi");
54 
55  double phiOffset;
56  if (sector < 12) phiOffset = 2*M_PI*((double)sector/12.) - 2*M_PI*(1./12.)/2;
57  else phiOffset = 2*M_PI*(((double)sector-12.)/12.) - 2*M_PI*(1./12.)/2;
58 
59  while(myReader->Next())
60  {
61  if (slot[inFee]-1 == *FEE && channel == *FEE_Chan)
62  {
63  double x = (*PadR/10)*std::cos(*phi+phiOffset);
64  double y = (*PadR/10)*std::sin(*phi+phiOffset);
65  TVector2 pos;
66  pos.SetX(x); pos.SetY(y);
67  return pos;
68  }
69  }
70  TVector2 pos;
71  pos.SetX(0); pos.SetY(0);
72  return pos;
73 }
74 
75 
76 void TPCEventDisplay(const float adcCut = 100.,
77  const string &outfile = "/sphenix/user/rosstom/test/TPCEventDisplay_10684",
78  const string &indir = "/sphenix/user/rosstom/test/testFiles/",
79  const string &runName = "beam-00010684",
80  const bool writeAllHits = true){
81 
82  //output nTuple
83  TString* outputfilename=new TString(outfile+".root");
84  TFile* outputfile=new TFile(outputfilename->Data(),"recreate");
85  TTree* hitTree=new TTree("hitTree","x,y,z position and max ADC value for all samples passing cut");
86  double x_Pos, y_Pos, z_Pos;
87  int max_adc;
88  hitTree->Branch("x_Pos",&x_Pos,"x_Pos/D");
89  hitTree->Branch("y_Pos",&y_Pos,"y_Pos/D");
90  hitTree->Branch("z_Pos",&z_Pos,"z_Pos/D");
91  hitTree->Branch("max_adc",&max_adc,"max_adc/I");
92 
93  std::ofstream outdata;
94  outdata.open((outfile+".json").c_str());
95  if( !outdata ) { // file couldn't be opened
96  cerr << "ERROR: file could not be opened" << endl;
97  exit(1);
98  }
99 
100  string runNumberString = runName;
101  size_t pos = runNumberString.find("000");
102  runNumberString.erase(runNumberString.begin(),runNumberString.begin()+pos+3);
103 
104  outdata << "{\n \"EVENT\": {\n \"runid\":" << runNumberString << ", \n \"evtid\": 1, \n \"time\": 0, \n \"timeStr\": \"2023-05-30, 15:23:30 EST\", \n \"type\": \"Cosmics\", \n \"s_nn\": 0, \n \"B\": 0.0,\n \"pv\": [0,0,0] \n },\n" << endl;
105 
106  outdata << " \"META\": {\n \"HITS\": {\n \"INNERTRACKER\": {\n \"type\": \"3D\",\n \"options\": {\n \"size\": 5,\n \"color\": 16777215\n } \n },\n" << endl;
107  outdata << " \"TRACKHITS\": {\n \"type\": \"3D\",\n \"options\": {\n \"size\": 5,\n \"transparent\": 0.5,\n \"color\": 16777215\n } \n },\n" << endl;
108  outdata << " \"JETS\": {\n \"type\": \"JET\",\n \"options\": {\n \"rmin\": 0,\n \"rmax\": 78,\n \"emin\": 0,\n \"emax\": 30,\n \"color\": 16777215,\n \"transparent\": 0.5 \n }\n }\n }\n }\n," << endl;
109  outdata << " \"HITS\": {\n \"CEMC\":[{\"eta\": 0, \"phi\": 0, \"e\": 0}\n ],\n \"HCALIN\": [{\"eta\": 0, \"phi\": 0, \"e\": 0}\n ],\n \"HCALOUT\": [{\"eta\": 0, \"phi\": 0, \"e\": 0}\n \n ],\n\n" << endl;
110  outdata << " \"TRACKHITS\": [\n\n ";
111 
112  TTree *R1 = new TTree("R1","TPC Sector Mapping for R1");
113  TTree *R2 = new TTree("R2","TPC Sector Mapping for R1");
114  TTree *R3 = new TTree("R3","TPC Sector Mapping for R1");
115 
116  R1->ReadFile("/sphenix/u/rosstom/calibrations/TPC/Mapping/PadPlane/AutoPad-R1-RevA.sch.ChannelMapping.csv",
117  "Entry/I:Radius/I:Pad/I:U/I:G/I:Pin/C:PinColID/I:PinRowID/I:PadName/C:FEE/F:FEE_Connector/C:FEE_Chan/I:phi/D:x/D:y/D:PadX/D:PadY/D:PadR/D:PadPhi/D",',');
118  R2->ReadFile("/sphenix/u/rosstom/calibrations/TPC/Mapping/PadPlane/AutoPad-R2-RevA-Pads.sch.ChannelMapping.csv",
119  "Entry/I:Radius/I:Pad/I:U/I:G/I:Pin/C:PinColID/I:PinRowID/I:PadName/C:FEE/F:FEE_Connector/C:FEE_Chan/I:phi/D:x/D:y/D:PadX/D:PadY/D:PadR/D:PadPhi/D",',');
120  R3->ReadFile("/sphenix/u/rosstom/calibrations/TPC/Mapping/PadPlane/AutoPad-R3-RevA.sch.ChannelMapping.csv",
121  "Entry/I:Radius/I:Pad/I:U/I:G/I:Pin/C:PinColID/I:PinRowID/I:PadName/C:FEE/F:FEE_Connector/C:FEE_Chan/I:phi/D:x/D:y/D:PadX/D:PadY/D:PadR/D:PadPhi/D",',');
122 
123  bool firstClus = true;
124 
125  //int framevals[12] = {77,67,69,71,70,79,77,70,69,71,64,77};
126 
127  for (int q = 0; q < 24; q++)
128  {
129  cout << "Processing sector " << q << endl;
130 
131  string sectorNumber;
132  if (q < 10) sectorNumber = "0"+std::to_string(q);
133  else sectorNumber = std::to_string(q);
134 
135  string fileName = indir+"TPC_ebdc"+sectorNumber+"_"+runName+"-0005.prdf_TPCRawDataTree.root";
136 
137  TFile *TPCFile = TFile::Open(fileName.c_str());
138  TTreeReader myReader("SampleTree", TPCFile);
139 
140  TTreeReaderValue<Int_t> nSamples(myReader, "nSamples");
141  TTreeReaderValue<Int_t> fee(myReader, "fee");
142  TTreeReaderArray<UShort_t> adcSamples(myReader, "adcSamples");
143  TTreeReaderValue<Int_t> Channel(myReader, "Channel");
144  TTreeReaderValue<Int_t> BCO(myReader, "BCO");
145  TTreeReaderValue<Int_t> frame(myReader, "frame");
146  TTreeReaderValue<Int_t> checksumError(myReader, "checksumError");
147 
148  while(myReader.Next())
149  {
150  //if (*frame != framevals[q] && *frame != framevals[q]+1) continue;
151 
152  int mod_arr[26]={2,2,1,1,1,3,3,3,3,3,3,2,2,1,2,2,1,1,2,2,3,3,3,3,3,3};
153 
154  if (*nSamples <= 1) continue;
155 
156  float adcMax = 0;
157  float adcMaxPos = 0;
158  if (writeAllHits && *checksumError == 0)
159  {
160  for(int adc_sam_no=0;adc_sam_no<*nSamples;adc_sam_no++)
161  {
162  if (adcSamples[adc_sam_no] - adcSamples[0] > adcCut)
163  {
164  max_adc = adcSamples[adc_sam_no];
165  double zPos;
166  if (q < 12)
167  { //0.4 cm per sample, 50 ns/sample, 8 cm/us drift speed
168  zPos = 105. - 0.4*(double)adc_sam_no; // - (*BCO - globalEventBCO)*(8./9.4);
169  //if (zPos < 0. || zPos > 105.) continue;
170  }
171  else
172  {
173  zPos = -105. + 0.4*(double)adc_sam_no; // + (*BCO - globalEventBCO)*(8./9.4);
174  //if (zPos > 0. || zPos < 105.) continue;
175  }
176 
177  TVector2 binXY = getBinPosition(q,*fee,*Channel,R1,R2,R3);
178 
179  if (binXY.X() == 0 || binXY.Y() == 0) continue;
180 
181  x_Pos=binXY.X();
182  y_Pos=binXY.Y();
183  z_Pos=zPos;
184  max_adc=adcMax;
185  hitTree->Fill();
186 
187  stringstream spts;
188 
189  if (firstClus) firstClus = false;
190  else spts << ",";
191 
192  spts << "{ \"x\": ";
193  spts << binXY.X();
194  spts << ", \"y\": ";
195  spts << binXY.Y();
196  spts << ", \"z\": ";
197  spts << zPos;
198  spts << ", \"e\": ";
199  spts << adcMax;
200  spts << "}";
201 
202  outdata << (boost::format("%1%") % spts.str());
203  spts.clear();
204  spts.str("");
205  }
206  }
207  }
208  else
209  {
210  for(int adc_sam_no=0;adc_sam_no<*nSamples;adc_sam_no++)
211  {
212  if (adc_sam_no == 0)
213  {
214  adcMax = adcSamples[adc_sam_no];
215  adcMaxPos = adc_sam_no;
216  }
217  else
218  {
219  if (adcSamples[adc_sam_no] > adcMax)
220  {
221  adcMax = adcSamples[adc_sam_no];
222  adcMaxPos = adc_sam_no;
223  }
224  }
225  }
226  if (adcMax > adcCut)
227  {
228  double zPos;
229  //adcMaxPos = 0;
230  if (q < 12)
231  { //0.4 cm per sample, 50 ns/sample, 8 cm/us drift speed
232  zPos = 105. - 0.4*(double)adcMaxPos; // - (*BCO - globalEventBCO)*(8./9.4);
233  if (zPos < 0. || zPos > 105.) continue;
234  }
235  else
236  {
237  zPos = -105. + 0.4*(double)adcMaxPos; // + (*BCO - globalEventBCO)*(8./9.4);
238  if (zPos > 0. || zPos < -105.) continue;
239  }
240 
241  TVector2 binXY = getBinPosition(q,*fee,*Channel,R1,R2,R3);
242 
243  if (binXY.X() == 0 || binXY.Y() == 0) continue;
244 
245  x_Pos=binXY.X();
246  y_Pos=binXY.Y();
247  z_Pos=zPos;
248  max_adc=adcMax;
249  hitTree->Fill();
250 
251  stringstream spts;
252 
253  if (firstClus) firstClus = false;
254  else spts << ",";
255 
256  spts << "{ \"x\": ";
257  spts << binXY.X();
258  spts << ", \"y\": ";
259  spts << binXY.Y();
260  spts << ", \"z\": ";
261  spts << zPos;
262  spts << ", \"e\": ";
263  spts << adcMax;
264  spts << "}";
265 
266  outdata << (boost::format("%1%") % spts.str());
267  spts.clear();
268  spts.str("");
269  }
270  }
271  }
272  }
273  outdata << "],\n \"JETS\": [\n ]\n }," << endl;
274  outdata << "\"TRACKS\": {" << endl;
275  outdata <<"\""<<"INNERTRACKER"<<"\": [";
276  outdata << "]" << endl;
277  outdata << "}" << endl;
278  outdata << "}" << endl;
279 
280  outdata.close();
281 
282  outputfile->cd();
283  hitTree->Write();
284  outputfile->Close();
285 }