DAQ/macros/TPCEventDisplay_Updated.C

0001 #include <cstdlib>
0002 #include <iostream>
0003 #include <map>
0004 #include <string>
0005 #include <vector>
0006
0007 #include "TChain.h"
0008 #include "TFile.h"
0009 #include "TTree.h"
0010 #include "TString.h"
0011 #include "TObjString.h"
0012 #include "TSystem.h"
0013 #include "TROOT.h"
0014 #include "TTreeReader.h"
0015 #include "TTreeReaderValue.h"
0016 #include "TTreeReaderArray.h"
0017 #include "TVector3.h"
0018
0019 #include <tpc/TpcMap.h>
0020
0021 #include <boost/format.hpp>
0022 #include <boost/math/special_functions/sign.hpp>
0023
0024 R__LOAD_LIBRARY(libtpc.so)
0025 /*************************************************************/
0026 /*               TPC Raw Data Event Display                  */
0027 /*               Thomas Marshall,Aditya Dash                 */
0028 /*        rosstom@ucla.edu,aditya55@physics.ucla.edu         */
0029 /*************************************************************/
0030
0031 void TPCEventDisplay_Updated(const float adcCut = 10.,
0032                      const string &outfile = "TPCEventDisplay_25926_100frame_test",
0033                      const string &indir = "./",
0034                      const string &runName = "cosmics-00025926",
0035              const bool writeAllHits = true){
0036
0037   //output nTuple
0038   TString* outputfilename=new TString(outfile+".root");
0039   TFile* outputfile=new TFile(outputfilename->Data(),"recreate");
0040   TTree* hitTree=new TTree("hitTree","x,y,z position and max ADC value for all samples passing cut");
0041   double x_Pos, y_Pos, z_Pos, frameVal;
0042   int adcMinusPedestal, sectorNum, layer, side, phiElement, timeBin, pad;
0043   hitTree->Branch("x_Pos",&x_Pos,"x_Pos/D");
0044   hitTree->Branch("y_Pos",&y_Pos,"y_Pos/D");
0045   hitTree->Branch("z_Pos",&z_Pos,"z_Pos/D");
0046   hitTree->Branch("frameVal",&frameVal,"frameVal/D");
0047   hitTree->Branch("adcMinusPedestal",&adcMinusPedestal,"adcMinusPedestal/I");
0048   hitTree->Branch("sectorNum",&sectorNum,"sectorNum/I");
0049   hitTree->Branch("layer",&layer,"layer/I");
0050   hitTree->Branch("side",&side,"side/I");
0051   hitTree->Branch("phiElement",&phiElement,"phiElement/I");
0052   hitTree->Branch("timeBin",&timeBin,"timeBin/I");
0053   hitTree->Branch("pad",&pad,"pad/I");
0054
0055   std::ofstream outdata;
0056   outdata.open((outfile+".json").c_str());
0057   if( !outdata ) { // file couldn't be opened
0058       cerr << "ERROR: file could not be opened" << endl;
0059       exit(1);
0060   }
0061
0062   string runNumberString = runName;
0063   size_t pos = runNumberString.find("000");
0064   runNumberString.erase(runNumberString.begin(),runNumberString.begin()+pos+3);
0065
0066   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;
0067
0068   outdata << "    \"META\": {\n       \"HITS\": {\n          \"INNERTRACKER\": {\n              \"type\": \"3D\",\n              \"options\": {\n              \"size\": 2,\n              \"color\": 16777215\n              } \n          },\n" << endl;
0069   outdata << "          \"TRACKHITS\": {\n              \"type\": \"3D\",\n              \"options\": {\n              \"size\": 2,\n              \"transparent\": 0.5,\n              \"color\": 16777215\n              } \n          },\n" << endl;
0070   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;
0071   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;
0072   outdata << "    \"TRACKHITS\": [\n\n ";
0073
0074   bool firstClus = true;
0075
0076   //int framevals[12] = {77,67,69,71,70,79,77,70,69,71,64,77};
0077   //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
0078   //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
0079   //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};
0080   //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};
0081   //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};
0082   //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};
0083   //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};
0084   //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};
0085   //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};
0086   //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};
0087   //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};
0088   //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};
0089   //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};
0090   //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};
0091   //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};
0092   //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};
0093
0094   TpcMap M;
0095   M.setMapNames("AutoPad-R1-RevA.sch.ChannelMapping.csv", "AutoPad-R2-RevA-Pads.sch.ChannelMapping.csv", "AutoPad-R3-RevA.sch.ChannelMapping.csv");
0096
0097   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};
0098   // 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};
0099
0100   // updated FEE mapping from Tom's description
0101   int FEE_map[26] = {
0102     4,
0103     5,
0104     0,
0105     2,
0106     1,
0107     11,
0108     9,
0109     10,
0110     8,
0111     7,
0112     6,
0113     0,
0114     1,
0115     3,
0116     7,
0117     6,
0118     5,
0119     4,
0120     3,
0121     2,
0122     0,
0123     2,
0124     1,
0125     3,
0126     5,
0127     4
0128   };
0129
0130   float pads_per_sector[3] = {96, 128, 192};
0131
0132   //ULong64_t BCOVal = 128330850912;
0133
0134   for (int q = 0; q < 24; q++)
0135   {
0136     //if (q == 17) continue;
0137     cout << "Processing sector " << q << endl;
0138
0139     sectorNum = q;
0140
0141     string sectorNumber;
0142     if (q < 10) sectorNumber = "0"+std::to_string(q);
0143     else sectorNumber = std::to_string(q);
0144
0145     string fileName = indir+"TPC_ebdc"+sectorNumber+"_"+runName+"-0000.prdf_test.root";
0146
0147     TFile *TPCFile = TFile::Open(fileName.c_str());
0148     TTree *sampleTree = (TTree*)TPCFile->Get("SampleTree");
0149     TTree *taggerTree = (TTree*)TPCFile->Get("TaggerTree");
0150
0151     Int_t nSamples, fee, Channel, frame, checksumError, taggerFrame;
0152     ULong64_t BCO;
0153     Double_t xPos, yPos;
0154     UShort_t adcSamples[5000];
0155     sampleTree->SetBranchAddress("nSamples",&nSamples);
0156     sampleTree->SetBranchAddress("fee",&fee);
0157     sampleTree->SetBranchAddress("adcSamples",&adcSamples);
0158     sampleTree->SetBranchAddress("Channel",&Channel);
0159     //sampleTree->SetBranchAddress("BCO",&BCO);
0160     sampleTree->SetBranchAddress("frame",&frame);
0161     sampleTree->SetBranchAddress("checksumError",&checksumError);
0162     sampleTree->SetBranchAddress("xPos",&xPos);
0163     sampleTree->SetBranchAddress("yPos",&yPos);
0164
0165     taggerTree->SetBranchAddress("bco",&BCO);
0166     taggerTree->SetBranchAddress("frame",&taggerFrame);
0167
0168     double channelMean[26][256];
0169     double channelSigma[26][256];
0170     double channelSamples[26][256];
0171     int nHitsChannel[26][256];
0172
0173     for(int fee_no=0;fee_no<26;fee_no++)
0174     {
0175       for(int channel_no=0;channel_no<256;channel_no++)
0176       {
0177         channelMean[fee_no][channel_no] = 0.;
0178         channelSigma[fee_no][channel_no] = 0.;
0179         channelSamples[fee_no][channel_no] = 0.;
0180         nHitsChannel[fee_no][channel_no] = 0;
0181       }
0182     }
0183
0184     //TTreeReader myReader("SampleTree", TPCFile);
0185
0186     /*
0187     TTreeReaderValue<Int_t> nSamples(myReader, "nSamples");
0188     TTreeReaderValue<Int_t> fee(myReader, "fee");
0189     TTreeReaderArray<UShort_t> adcSamples(myReader, "adcSamples");
0190     TTreeReaderValue<Int_t> Channel(myReader, "Channel");
0191     TTreeReaderValue<Int_t> BCO(myReader, "BCO");
0192     TTreeReaderValue<Int_t> frame(myReader, "frame");
0193     TTreeReaderValue<Int_t> checksumError(myReader, "checksumError");
0194     TTreeReaderValue<double> xPosition(myReader, "xPos");
0195     TTreeReaderValue<double> yPosition(myReader, "yPos");
0196     */
0197
0198     std::vector<int> frameValsToCheck;
0199
0200     for(int i = 0; i < taggerTree->GetEntries(); i++)
0201     {
0202       taggerTree->GetEntry(i);
0203       //if ((BCO < BCOVal + 5) && (BCO > BCOVal - 5))
0204       //if (BCO == BCOVal)
0205       //{
0206       //  frameValsToCheck.push_back(taggerFrame);
0207       //}
0208     }
0209
0210     for(int i = 0; i < sampleTree->GetEntries(); i++)
0211     {
0212       sampleTree->GetEntry(i);
0213       if (checksumError == 0)
0214       {
0215         for(int adc_sam_no=0;adc_sam_no<10;adc_sam_no++)
0216         {
0217           channelSamples[fee][Channel] += 1.;
0218           channelMean[fee][Channel] += adcSamples[adc_sam_no];
0219           channelSigma[fee][Channel] += pow(adcSamples[adc_sam_no],2);
0220         }
0221       }
0222     }
0223
0224     for(int fee_no=0;fee_no<26;fee_no++)
0225     {
0226       for(int channel_no=0;channel_no<256;channel_no++)
0227       {
0228         double temp1 = channelMean[fee_no][channel_no]/channelSamples[fee_no][channel_no];
0229         double temp2 = channelSigma[fee_no][channel_no]/channelSamples[fee_no][channel_no];
0230         channelMean[fee_no][channel_no] = temp1;
0231         channelSigma[fee_no][channel_no] = sqrt(temp2 - pow(temp1,2));
0232       }
0233     }
0234
0235     //int nHitsInEntry[1000];
0236     for(int i = 0; i < sampleTree->GetEntries(); i++)
0237     {
0238       sampleTree->GetEntry(i);
0239       //if (frame != framevals[q] && frame != framevals[q]+1) continue;
0240       bool recordHits = false;
0241       for (auto & vals : frameValsToCheck)
0242       {
0243         if (frame < vals + 2 && frame > vals - 1)
0244     {
0245           recordHits = true;
0246       break;
0247     }
0248       }
0249       //if (!recordHits) continue;
0250       if (checksumError == 0 && channelMean[fee][Channel] > 20. && channelMean[fee][Channel] < 200.)
0251       {
0252         for(int adc_sam_no=0;adc_sam_no<nSamples;adc_sam_no++)
0253         {
0254           if ((adcSamples[adc_sam_no] - channelMean[fee][Channel] > channelSigma[fee][Channel]*4) && (adcSamples[adc_sam_no] - channelMean[fee][Channel] > adcCut))
0255           {
0256             nHitsChannel[fee][Channel] += 1;
0257             //nHitsInEntry[frame] += 1;
0258           }
0259         }
0260       }
0261     }
0262
0263     for(int i = 0; i < sampleTree->GetEntries(); i++)
0264     {
0265       sampleTree->GetEntry(i);
0266       //if (*frame > 100) continue;
0267       //if (q == 2 && (*frame < framevals[q] - 5 || *frame > framevals[q]+5)) continue;
0268       //if (frame != framevals[q] && frame != framevals[q]+1) continue;
0269       if (nHitsChannel[fee][Channel] > 50) continue;
0270       //if (nHitsInEntry[frame] > 100) continue;
0271       if (nSamples <= 1) continue;
0272       bool recordHits = false;
0273       for (auto & vals : frameValsToCheck)
0274       {
0275         if (frame < vals + 2 && frame > vals - 1)
0276     {
0277           recordHits = true;
0278       break;
0279     }
0280       }
0281       //if (!recordHits) continue;
0282
0283       float adcMax = 0;
0284       float adcMaxPos = 0;
0285       if (checksumError == 0 && channelMean[fee][Channel] > 20. && channelMean[fee][Channel] < 200.)
0286       {
0287         for(int adc_sam_no=0;adc_sam_no<nSamples;adc_sam_no++)
0288         {
0289           if ((adcSamples[adc_sam_no] - channelMean[fee][Channel] > channelSigma[fee][Channel]*5) && (adcSamples[adc_sam_no] - channelMean[fee][Channel] > adcCut))
0290           {
0291             adcMinusPedestal = adcSamples[adc_sam_no] - channelMean[fee][Channel];
0292             double zPos;
0293             if (q < 12)
0294             { //0.4 cm per sample, 50 ns/sample, 8 cm/us drift speed
0295               zPos = 105. - (105./230.)*((double)adc_sam_no - 10.); // - (*BCO - globalEventBCO)*(8./9.4);
0296               //if (zPos < 0. || zPos > 105.) continue;
0297               side = 0;
0298               phiElement = q;
0299         }
0300             else
0301             {
0302               zPos = -105. + (105./230.)*((double)adc_sam_no - 10.); // + (*BCO - globalEventBCO)*(8./9.4);
0303               //if (zPos > 0. || zPos < 105.) continue;
0304               side = 1;
0305               phiElement = q - 12;
0306             }
0307
0308             int feeM = FEE_map[fee];
0309             if(mod_arr[fee]==2) feeM += 6;
0310             else if(mod_arr[fee]==3) feeM += 14;
0311             layer = M.getLayer(feeM, Channel);
0312             pad = M.getPad(feeM, Channel) + (phiElement)*pads_per_sector[mod_arr[fee]-1];
0313
0314
0315         if(layer!=0)
0316           {
0317         double R = M.getR(feeM, Channel);
0318         double phi = (sectorNum<12? -M.getPhi(feeM, Channel) : M.getPhi(feeM, Channel) ) + (sectorNum - side*12. )* M_PI / 6. ;
0319         R /= 10.; //convert mm to cm
0320
0321         xPos = R*cos(phi);
0322         yPos = R*sin(phi);
0323           }
0324
0325             x_Pos = xPos;
0326             y_Pos = yPos;
0327             z_Pos = zPos;
0328             frameVal = frame;
0329             timeBin = adc_sam_no;
0330             hitTree->Fill();
0331
0332             stringstream spts;
0333
0334             if (firstClus) firstClus = false;
0335             else spts << ",";
0336
0337         spts << "{ \"x\": ";
0338             spts << x_Pos;
0339             spts << ", \"y\": ";
0340             spts << y_Pos;
0341             spts << ", \"z\": ";
0342             spts << z_Pos;
0343             spts << ", \"e\": ";
0344             spts << adcMinusPedestal;
0345             spts << "}";
0346
0347             outdata << (boost::format("%1%") % spts.str());
0348             spts.clear();
0349             spts.str("");
0350           }
0351         }
0352       }
0353     }
0354     TPCFile->Close();
0355   }
0356   outdata << "],\n    \"JETS\": [\n         ]\n    }," << endl;
0357   outdata << "\"TRACKS\": {" << endl;
0358   outdata <<"\""<<"INNERTRACKER"<<"\": [";
0359   outdata << "]" << endl;
0360   outdata << "}" << endl;
0361   outdata << "}" << endl;
0362
0363   outdata.close();
0364
0365   outputfile->cd();
0366   hitTree->Write();
0367   outputfile->Close();
0368   std::cout << "Done" << std::endl;
0369 }