Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TPCEventDisplay_Updated.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TPCEventDisplay_Updated.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 <tpc/TpcMap.h>
20 
21 #include <boost/format.hpp>
22 #include <boost/math/special_functions/sign.hpp>
23 
24 R__LOAD_LIBRARY(libtpc.so)
25 /*************************************************************/
26 /* TPC Raw Data Event Display */
27 /* Thomas Marshall,Aditya Dash */
28 /* rosstom@ucla.edu,aditya55@physics.ucla.edu */
29 /*************************************************************/
30 
31 void TPCEventDisplay_Updated(const float adcCut = 10.,
32  const string &outfile = "TPCEventDisplay_25926_100frame_test",
33  const string &indir = "./",
34  const string &runName = "cosmics-00025926",
35  const bool writeAllHits = true){
36 
37  //output nTuple
38  TString* outputfilename=new TString(outfile+".root");
39  TFile* outputfile=new TFile(outputfilename->Data(),"recreate");
40  TTree* hitTree=new TTree("hitTree","x,y,z position and max ADC value for all samples passing cut");
41  double x_Pos, y_Pos, z_Pos, frameVal;
42  int adcMinusPedestal, sectorNum, layer, side, phiElement, timeBin, pad;
43  hitTree->Branch("x_Pos",&x_Pos,"x_Pos/D");
44  hitTree->Branch("y_Pos",&y_Pos,"y_Pos/D");
45  hitTree->Branch("z_Pos",&z_Pos,"z_Pos/D");
46  hitTree->Branch("frameVal",&frameVal,"frameVal/D");
47  hitTree->Branch("adcMinusPedestal",&adcMinusPedestal,"adcMinusPedestal/I");
48  hitTree->Branch("sectorNum",&sectorNum,"sectorNum/I");
49  hitTree->Branch("layer",&layer,"layer/I");
50  hitTree->Branch("side",&side,"side/I");
51  hitTree->Branch("phiElement",&phiElement,"phiElement/I");
52  hitTree->Branch("timeBin",&timeBin,"timeBin/I");
53  hitTree->Branch("pad",&pad,"pad/I");
54 
55  std::ofstream outdata;
56  outdata.open((outfile+".json").c_str());
57  if( !outdata ) { // file couldn't be opened
58  cerr << "ERROR: file could not be opened" << endl;
59  exit(1);
60  }
61 
62  string runNumberString = runName;
63  size_t pos = runNumberString.find("000");
64  runNumberString.erase(runNumberString.begin(),runNumberString.begin()+pos+3);
65 
66  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 \"runstats\": [ \n \"sPHENIX Time Projection Chamber\", \"2023-08-23, Run 25926 - All EBDCs, 100 frames\", \"Cosmics Data\"] \n },\n" << endl;
67 
68  outdata << " \"META\": {\n \"HITS\": {\n \"INNERTRACKER\": {\n \"type\": \"3D\",\n \"options\": {\n \"size\": 2,\n \"color\": 16777215\n } \n },\n" << endl;
69  outdata << " \"TRACKHITS\": {\n \"type\": \"3D\",\n \"options\": {\n \"size\": 2,\n \"transparent\": 0.5,\n \"color\": 16777215\n } \n },\n" << endl;
70  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;
71  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;
72  outdata << " \"TRACKHITS\": [\n\n ";
73 
74  bool firstClus = true;
75 
76  //int framevals[12] = {77,67,69,71,70,79,77,70,69,71,64,77};
77  //int framevals[24] = {922,843,839,875,863,890,834,849,853,815,837,855,0,0,0,0,0,0,0,0,0,0,0,0}; //sector 3
78  //int framevals[24] = {0,0,0,875,0,0,0,0,0,0,0,0,1152,1071,1076,1083,1211,1137,1054,1158,1095,67,1532,1169}; //sector 12
79  //int framevals[24] = {23,21,23,21,21,20,20,23,22,19,24,24,23,22,23,24,24,23,25,21,25,22,-1,-1};
80  //int framevals[24] = {88,87,91,92,86,91,91,91,90,88,90,88,98,96,-1,87,112,92,89,86,105,95,280,251};
81  //int framevals[24] = {9,16,15,13,9,13,16,14,13,13,11,14,14,12,10,14,13,14,13,14,14,14,13,13};
82  //int framevals[24] = {20,30,27,25,19,31,36,27,26,25,22,24,24,24,22,25,23,26,23,26,24,24,25,25};
83  //int framevals[24] = {38,48,-1,43,37,52,64,50,44,41,39,40,47,40,47,43,41,46,39,48,42,42,50,44};
84  //int framevals[24] = {26,36,33,31,25,38,47,34,33,31,28,30,34,30,-1,31,29,34,29,36,32,32,34,32};
85  //int framevals[24] = {4,4,6,4,4,4,6,6,4,4,4,4,4,4,4,5,5,5,4,4,6,4,4,4};
86  //int framevals[24] = {11,18,17,16,11,17,19,16,16,16,14,16,16,14,12,16,15,16,15,16,16,16,15,15};
87  //int framevals[24] = {28,38,35,34,27,40,51,38,36,33,30,32,37,32,36,33,31,36,31,38,34,34,38,34};
88  //int framevals[24] = {30,40,37,37,29,42,54,41,38,35,32,34,39,34,38,35,33,40,33,40,36,36,42,36};
89  //int framevals[24] = {38,48,37,43,37,52,64,50,44,41,39,40,47,40,47,43,41,46,39,48,42,42,50,44};
90  //int framevals[24] = {44,55,50,52,43,62,76,58,51,48,45,46,58,46,54,53,47,52,45,54,50,50,57,53};
91  //int framevals[24] = {80,87,83,85,75,93,76,91,83,81,75,74,92,74,91,84,77,80,85,85,82,82,93,86};
92  //int framevals[24] = {84,91,89,89,79,97,76,97,87,85,79,74,96,80,97,90,83,86,90,90,86,86,97,90};
93 
94  TpcMap M;
95  M.setMapNames("AutoPad-R1-RevA.sch.ChannelMapping.csv", "AutoPad-R2-RevA-Pads.sch.ChannelMapping.csv", "AutoPad-R3-RevA.sch.ChannelMapping.csv");
96 
97  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};
98  // int FEE_map[26] = {3, 2, 5, 3, 4, 0, 2, 1, 3, 4, 5, 7, 6, 2, 0, 1, 0, 1, 4, 5, 11, 9, 10, 8, 6, 7};
99 
100  // updated FEE mapping from Tom's description
101  int FEE_map[26] = {
102  4,
103  5,
104  0,
105  2,
106  1,
107  11,
108  9,
109  10,
110  8,
111  7,
112  6,
113  0,
114  1,
115  3,
116  7,
117  6,
118  5,
119  4,
120  3,
121  2,
122  0,
123  2,
124  1,
125  3,
126  5,
127  4
128  };
129 
130  float pads_per_sector[3] = {96, 128, 192};
131 
132  //ULong64_t BCOVal = 128330850912;
133 
134  for (int q = 0; q < 24; q++)
135  {
136  //if (q == 17) continue;
137  cout << "Processing sector " << q << endl;
138 
139  sectorNum = q;
140 
141  string sectorNumber;
142  if (q < 10) sectorNumber = "0"+std::to_string(q);
143  else sectorNumber = std::to_string(q);
144 
145  string fileName = indir+"TPC_ebdc"+sectorNumber+"_"+runName+"-0000.prdf_test.root";
146 
147  TFile *TPCFile = TFile::Open(fileName.c_str());
148  TTree *sampleTree = (TTree*)TPCFile->Get("SampleTree");
149  TTree *taggerTree = (TTree*)TPCFile->Get("TaggerTree");
150 
151  Int_t nSamples, fee, Channel, frame, checksumError, taggerFrame;
152  ULong64_t BCO;
153  Double_t xPos, yPos;
154  UShort_t adcSamples[5000];
155  sampleTree->SetBranchAddress("nSamples",&nSamples);
156  sampleTree->SetBranchAddress("fee",&fee);
157  sampleTree->SetBranchAddress("adcSamples",&adcSamples);
158  sampleTree->SetBranchAddress("Channel",&Channel);
159  //sampleTree->SetBranchAddress("BCO",&BCO);
160  sampleTree->SetBranchAddress("frame",&frame);
161  sampleTree->SetBranchAddress("checksumError",&checksumError);
162  sampleTree->SetBranchAddress("xPos",&xPos);
163  sampleTree->SetBranchAddress("yPos",&yPos);
164 
165  taggerTree->SetBranchAddress("bco",&BCO);
166  taggerTree->SetBranchAddress("frame",&taggerFrame);
167 
168  double channelMean[26][256];
169  double channelSigma[26][256];
170  double channelSamples[26][256];
171  int nHitsChannel[26][256];
172 
173  for(int fee_no=0;fee_no<26;fee_no++)
174  {
175  for(int channel_no=0;channel_no<256;channel_no++)
176  {
177  channelMean[fee_no][channel_no] = 0.;
178  channelSigma[fee_no][channel_no] = 0.;
179  channelSamples[fee_no][channel_no] = 0.;
180  nHitsChannel[fee_no][channel_no] = 0;
181  }
182  }
183 
184  //TTreeReader myReader("SampleTree", TPCFile);
185 
186  /*
187  TTreeReaderValue<Int_t> nSamples(myReader, "nSamples");
188  TTreeReaderValue<Int_t> fee(myReader, "fee");
189  TTreeReaderArray<UShort_t> adcSamples(myReader, "adcSamples");
190  TTreeReaderValue<Int_t> Channel(myReader, "Channel");
191  TTreeReaderValue<Int_t> BCO(myReader, "BCO");
192  TTreeReaderValue<Int_t> frame(myReader, "frame");
193  TTreeReaderValue<Int_t> checksumError(myReader, "checksumError");
194  TTreeReaderValue<double> xPosition(myReader, "xPos");
195  TTreeReaderValue<double> yPosition(myReader, "yPos");
196  */
197 
198  std::vector<int> frameValsToCheck;
199 
200  for(int i = 0; i < taggerTree->GetEntries(); i++)
201  {
202  taggerTree->GetEntry(i);
203  //if ((BCO < BCOVal + 5) && (BCO > BCOVal - 5))
204  //if (BCO == BCOVal)
205  //{
206  // frameValsToCheck.push_back(taggerFrame);
207  //}
208  }
209 
210  for(int i = 0; i < sampleTree->GetEntries(); i++)
211  {
212  sampleTree->GetEntry(i);
213  if (checksumError == 0)
214  {
215  for(int adc_sam_no=0;adc_sam_no<10;adc_sam_no++)
216  {
217  channelSamples[fee][Channel] += 1.;
218  channelMean[fee][Channel] += adcSamples[adc_sam_no];
219  channelSigma[fee][Channel] += pow(adcSamples[adc_sam_no],2);
220  }
221  }
222  }
223 
224  for(int fee_no=0;fee_no<26;fee_no++)
225  {
226  for(int channel_no=0;channel_no<256;channel_no++)
227  {
228  double temp1 = channelMean[fee_no][channel_no]/channelSamples[fee_no][channel_no];
229  double temp2 = channelSigma[fee_no][channel_no]/channelSamples[fee_no][channel_no];
230  channelMean[fee_no][channel_no] = temp1;
231  channelSigma[fee_no][channel_no] = sqrt(temp2 - pow(temp1,2));
232  }
233  }
234 
235  //int nHitsInEntry[1000];
236  for(int i = 0; i < sampleTree->GetEntries(); i++)
237  {
238  sampleTree->GetEntry(i);
239  //if (frame != framevals[q] && frame != framevals[q]+1) continue;
240  bool recordHits = false;
241  for (auto & vals : frameValsToCheck)
242  {
243  if (frame < vals + 2 && frame > vals - 1)
244  {
245  recordHits = true;
246  break;
247  }
248  }
249  //if (!recordHits) continue;
250  if (checksumError == 0 && channelMean[fee][Channel] > 20. && channelMean[fee][Channel] < 200.)
251  {
252  for(int adc_sam_no=0;adc_sam_no<nSamples;adc_sam_no++)
253  {
254  if ((adcSamples[adc_sam_no] - channelMean[fee][Channel] > channelSigma[fee][Channel]*4) && (adcSamples[adc_sam_no] - channelMean[fee][Channel] > adcCut))
255  {
256  nHitsChannel[fee][Channel] += 1;
257  //nHitsInEntry[frame] += 1;
258  }
259  }
260  }
261  }
262 
263  for(int i = 0; i < sampleTree->GetEntries(); i++)
264  {
265  sampleTree->GetEntry(i);
266  //if (*frame > 100) continue;
267  //if (q == 2 && (*frame < framevals[q] - 5 || *frame > framevals[q]+5)) continue;
268  //if (frame != framevals[q] && frame != framevals[q]+1) continue;
269  if (nHitsChannel[fee][Channel] > 50) continue;
270  //if (nHitsInEntry[frame] > 100) continue;
271  if (nSamples <= 1) continue;
272  bool recordHits = false;
273  for (auto & vals : frameValsToCheck)
274  {
275  if (frame < vals + 2 && frame > vals - 1)
276  {
277  recordHits = true;
278  break;
279  }
280  }
281  //if (!recordHits) continue;
282 
283  float adcMax = 0;
284  float adcMaxPos = 0;
285  if (checksumError == 0 && channelMean[fee][Channel] > 20. && channelMean[fee][Channel] < 200.)
286  {
287  for(int adc_sam_no=0;adc_sam_no<nSamples;adc_sam_no++)
288  {
289  if ((adcSamples[adc_sam_no] - channelMean[fee][Channel] > channelSigma[fee][Channel]*5) && (adcSamples[adc_sam_no] - channelMean[fee][Channel] > adcCut))
290  {
291  adcMinusPedestal = adcSamples[adc_sam_no] - channelMean[fee][Channel];
292  double zPos;
293  if (q < 12)
294  { //0.4 cm per sample, 50 ns/sample, 8 cm/us drift speed
295  zPos = 105. - (105./230.)*((double)adc_sam_no - 10.); // - (*BCO - globalEventBCO)*(8./9.4);
296  //if (zPos < 0. || zPos > 105.) continue;
297  side = 0;
298  phiElement = q;
299  }
300  else
301  {
302  zPos = -105. + (105./230.)*((double)adc_sam_no - 10.); // + (*BCO - globalEventBCO)*(8./9.4);
303  //if (zPos > 0. || zPos < 105.) continue;
304  side = 1;
305  phiElement = q - 12;
306  }
307 
308  int feeM = FEE_map[fee];
309  if(mod_arr[fee]==2) feeM += 6;
310  else if(mod_arr[fee]==3) feeM += 14;
311  layer = M.getLayer(feeM, Channel);
312  pad = M.getPad(feeM, Channel) + (phiElement)*pads_per_sector[mod_arr[fee]-1];
313 
314 
315  if(layer!=0)
316  {
317  double R = M.getR(feeM, Channel);
318  double phi = (sectorNum<12? -M.getPhi(feeM, Channel) : M.getPhi(feeM, Channel) ) + (sectorNum - side*12. )* M_PI / 6. ;
319  R /= 10.; //convert mm to cm
320 
321  xPos = R*cos(phi);
322  yPos = R*sin(phi);
323  }
324 
325  x_Pos = xPos;
326  y_Pos = yPos;
327  z_Pos = zPos;
328  frameVal = frame;
329  timeBin = adc_sam_no;
330  hitTree->Fill();
331 
332  stringstream spts;
333 
334  if (firstClus) firstClus = false;
335  else spts << ",";
336 
337  spts << "{ \"x\": ";
338  spts << x_Pos;
339  spts << ", \"y\": ";
340  spts << y_Pos;
341  spts << ", \"z\": ";
342  spts << z_Pos;
343  spts << ", \"e\": ";
344  spts << adcMinusPedestal;
345  spts << "}";
346 
347  outdata << (boost::format("%1%") % spts.str());
348  spts.clear();
349  spts.str("");
350  }
351  }
352  }
353  }
354  TPCFile->Close();
355  }
356  outdata << "],\n \"JETS\": [\n ]\n }," << endl;
357  outdata << "\"TRACKS\": {" << endl;
358  outdata <<"\""<<"INNERTRACKER"<<"\": [";
359  outdata << "]" << endl;
360  outdata << "}" << endl;
361  outdata << "}" << endl;
362 
363  outdata.close();
364 
365  outputfile->cd();
366  hitTree->Write();
367  outputfile->Close();
368  std::cout << "Done" << std::endl;
369 }