Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TPCEventDisplay_Updated_BCO.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TPCEventDisplay_Updated_BCO.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_BCO(const float adcCut = 10.,
32  const string &outfile = "TPCEventDisplay_25926_BCO",
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\", \"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 
99  //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};
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 
131  float pads_per_sector[3] = {96, 128, 192};
132 
133  //ULong64_t BCOVal = 128330850912;
134  // ULong64_t BCOVal = 128330850911;
135  ULong64_t BCOVal =128434038131;
136 
137  for (int q = 0; q < 24; q++)
138  {
139  //if (q == 17) continue;
140  cout << "Processing sector " << q << endl;
141 
142  sectorNum = q;
143 
144  string sectorNumber;
145  if (q < 10) sectorNumber = "0"+std::to_string(q);
146  else sectorNumber = std::to_string(q);
147 
148  string fileName = indir+"TPC_ebdc"+sectorNumber+"_"+runName+"-0001.prdf_test.root";
149 
150  TFile *TPCFile = TFile::Open(fileName.c_str());
151  TTree *sampleTree = (TTree*)TPCFile->Get("SampleTree");
152  TTree *taggerTree = (TTree*)TPCFile->Get("TaggerTree");
153 
154  Int_t nSamples, fee, Channel, frame, packet, checksumError, taggerFrame;
155  ULong64_t BCO;
156  Double_t xPos, yPos;
157  UShort_t adcSamples[5000];
158  sampleTree->SetBranchAddress("nSamples",&nSamples);
159  sampleTree->SetBranchAddress("fee",&fee);
160  sampleTree->SetBranchAddress("adcSamples",&adcSamples);
161  sampleTree->SetBranchAddress("Channel",&Channel);
162  //sampleTree->SetBranchAddress("BCO",&BCO);
163  sampleTree->SetBranchAddress("frame",&frame);
164  sampleTree->SetBranchAddress("checksumError",&checksumError);
165  sampleTree->SetBranchAddress("xPos",&xPos);
166  sampleTree->SetBranchAddress("yPos",&yPos);
167 
168  taggerTree->SetBranchAddress("bco",&BCO);
169  taggerTree->SetBranchAddress("packet",&packet);
170  taggerTree->SetBranchAddress("frame",&taggerFrame);
171 
172  double channelMean[26][256];
173  double channelSigma[26][256];
174  double channelSamples[26][256];
175  int nHitsChannel[26][256];
176 
177  for(int fee_no=0;fee_no<26;fee_no++)
178  {
179  for(int channel_no=0;channel_no<256;channel_no++)
180  {
181  channelMean[fee_no][channel_no] = 0.;
182  channelSigma[fee_no][channel_no] = 0.;
183  channelSamples[fee_no][channel_no] = 0.;
184  nHitsChannel[fee_no][channel_no] = 0;
185  }
186  }
187 
188  //TTreeReader myReader("SampleTree", TPCFile);
189 
190  /*
191  TTreeReaderValue<Int_t> nSamples(myReader, "nSamples");
192  TTreeReaderValue<Int_t> fee(myReader, "fee");
193  TTreeReaderArray<UShort_t> adcSamples(myReader, "adcSamples");
194  TTreeReaderValue<Int_t> Channel(myReader, "Channel");
195  TTreeReaderValue<Int_t> BCO(myReader, "BCO");
196  TTreeReaderValue<Int_t> frame(myReader, "frame");
197  TTreeReaderValue<Int_t> checksumError(myReader, "checksumError");
198  TTreeReaderValue<double> xPosition(myReader, "xPos");
199  TTreeReaderValue<double> yPosition(myReader, "yPos");
200  */
201 
202  std::vector<int> frameValsToCheck;
203 
204  for(int i = 0; i < taggerTree->GetEntries(); i++)
205  {
206  taggerTree->GetEntry(i);
207  if ((BCO < BCOVal + 5) && (BCO > BCOVal - 5))
208  //if (BCO == BCOVal)
209  {
210  cout <<"Sector "<<q<<" found packet = "<<packet<<" taggerFrame = "<<taggerFrame<<" with BCO = "<<BCO<<endl;
211  frameValsToCheck.push_back(taggerFrame);
212  }
213  }
214 
215  for(int i = 0; i < sampleTree->GetEntries(); i++)
216  {
217  sampleTree->GetEntry(i);
218  if (checksumError == 0)
219  {
220  for(int adc_sam_no=0;adc_sam_no<10;adc_sam_no++)
221  {
222  channelSamples[fee][Channel] += 1.;
223  channelMean[fee][Channel] += adcSamples[adc_sam_no];
224  channelSigma[fee][Channel] += pow(adcSamples[adc_sam_no],2);
225  }
226  }
227  }
228 
229  for(int fee_no=0;fee_no<26;fee_no++)
230  {
231  for(int channel_no=0;channel_no<256;channel_no++)
232  {
233  double temp1 = channelMean[fee_no][channel_no]/channelSamples[fee_no][channel_no];
234  double temp2 = channelSigma[fee_no][channel_no]/channelSamples[fee_no][channel_no];
235  channelMean[fee_no][channel_no] = temp1;
236  channelSigma[fee_no][channel_no] = sqrt(temp2 - pow(temp1,2));
237  }
238  }
239 
240  //int nHitsInEntry[1000];
241  for(int i = 0; i < sampleTree->GetEntries(); i++)
242  {
243  sampleTree->GetEntry(i);
244  //if (frame != framevals[q] && frame != framevals[q]+1) continue;
245  bool recordHits = false;
246  for (auto & vals : frameValsToCheck)
247  {
248  if (frame < vals + 2 && frame > vals - 1)
249  {
250  recordHits = true;
251  break;
252  }
253  }
254  //if (!recordHits) continue;
255  if (checksumError == 0 && channelMean[fee][Channel] > 20. && channelMean[fee][Channel] < 200.)
256  {
257  for(int adc_sam_no=0;adc_sam_no<nSamples;adc_sam_no++)
258  {
259  if ((adcSamples[adc_sam_no] - channelMean[fee][Channel] > channelSigma[fee][Channel]*4) && (adcSamples[adc_sam_no] - channelMean[fee][Channel] > adcCut))
260  {
261  nHitsChannel[fee][Channel] += 1;
262  //nHitsInEntry[frame] += 1;
263  }
264  }
265  }
266  }
267 
268  int sample_count = 0;
269  for(int i = 0; i < sampleTree->GetEntries(); i++)
270  {
271  sampleTree->GetEntry(i);
272  //if (*frame > 100) continue;
273  //if (q == 2 && (*frame < framevals[q] - 5 || *frame > framevals[q]+5)) continue;
274  //if (frame != framevals[q] && frame != framevals[q]+1) continue;
275  if (nHitsChannel[fee][Channel] > 50) continue;
276  //if (nHitsInEntry[frame] > 100) continue;
277  if (nSamples <= 1) continue;
278  bool recordHits = false;
279  for (auto & vals : frameValsToCheck)
280  {
281  if (frame < vals + 2 && frame > vals - 1)
282  {
283  recordHits = true;
284  break;
285  }
286  }
287  if (!recordHits) continue;
288  ++sample_count;
289 
290  float adcMax = 0;
291  float adcMaxPos = 0;
292  if (checksumError == 0 && channelMean[fee][Channel] > 20. && channelMean[fee][Channel] < 200.)
293  {
294  for(int adc_sam_no=0;adc_sam_no<nSamples;adc_sam_no++)
295  {
296  if ((adcSamples[adc_sam_no] - channelMean[fee][Channel] > channelSigma[fee][Channel]*5) && (adcSamples[adc_sam_no] - channelMean[fee][Channel] > adcCut))
297  {
298  adcMinusPedestal = adcSamples[adc_sam_no] - channelMean[fee][Channel];
299  double zPos;
300  if (q < 12)
301  { //0.4 cm per sample, 50 ns/sample, 8 cm/us drift speed
302  zPos = 105. - (105./230.)*((double)adc_sam_no - 10.); // - (*BCO - globalEventBCO)*(8./9.4);
303  //if (zPos < 0. || zPos > 105.) continue;
304  side = 0;
305  phiElement = q;
306  }
307  else
308  {
309  zPos = -105. + (105./230.)*((double)adc_sam_no - 10.); // + (*BCO - globalEventBCO)*(8./9.4);
310  //if (zPos > 0. || zPos < 105.) continue;
311  side = 1;
312  phiElement = q - 12;
313  }
314 
315  int feeM = FEE_map[fee];
316  if(mod_arr[fee]==2) feeM += 6;
317  else if(mod_arr[fee]==3) feeM += 14;
318  layer = M.getLayer(feeM, Channel);
319  pad = M.getPad(feeM, Channel) + (phiElement)*pads_per_sector[mod_arr[fee]-1];
320 
321 
322  if(layer!=0)
323  {
324  double R = M.getR(feeM, Channel);
325  double phi = (sectorNum<12? -M.getPhi(feeM, Channel) : M.getPhi(feeM, Channel) ) + (sectorNum - side*12. )* M_PI / 6. ;
326  R /= 10.; //convert mm to cm
327 
328  xPos = R*cos(phi);
329  yPos = R*sin(phi);
330  }
331 
332  x_Pos = xPos;
333  y_Pos = yPos;
334  z_Pos = zPos;
335  frameVal = frame;
336  timeBin = adc_sam_no;
337  hitTree->Fill();
338 
339  stringstream spts;
340 
341  if (firstClus) firstClus = false;
342  else spts << ",";
343 
344  spts << "{ \"x\": ";
345  spts << x_Pos;
346  spts << ", \"y\": ";
347  spts << y_Pos;
348  spts << ", \"z\": ";
349  spts << z_Pos;
350  spts << ", \"e\": ";
351  spts << adcMinusPedestal;
352  spts << "}";
353 
354  outdata << (boost::format("%1%") % spts.str());
355  spts.clear();
356  spts.str("");
357  }
358  }
359  }
360 
361  }
362 
363  cout <<"sample_count = "<<sample_count<<endl;
364 
365  TPCFile->Close();
366  }
367  outdata << "],\n \"JETS\": [\n ]\n }," << endl;
368  outdata << "\"TRACKS\": {" << endl;
369  outdata <<"\""<<"INNERTRACKER"<<"\": [";
370  outdata << "]" << endl;
371  outdata << "}" << endl;
372  outdata << "}" << endl;
373 
374  outdata.close();
375 
376  outputfile->cd();
377  hitTree->Write();
378  outputfile->Close();
379  std::cout << "Done" << std::endl;
380 }