macros/prelimEvtDisplay/TPCEventDisplay.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 <boost/format.hpp>
0020 #include <boost/math/special_functions/sign.hpp>
0021
0022 /*************************************************************/
0023 /*               TPC Raw Data Event Display                  */
0024 /*               Thomas Marshall,Aditya Dash                 */
0025 /*        rosstom@ucla.edu,aditya55@physics.ucla.edu         */
0026 /*************************************************************/
0027
0028 TVector2 getBinPosition(int sector, int inFee, int channel, TTree* R1, TTree* R2, TTree* R3)
0029 {
0030   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};
0031   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};
0032   TTreeReader* myReader;
0033   if (mod_arr[inFee] == 1)
0034   {
0035     myReader = new TTreeReader(R1);
0036   }
0037   else if (mod_arr[inFee] == 2)
0038   {
0039     myReader = new TTreeReader(R2);
0040   }
0041   else if (mod_arr[inFee] == 3)
0042   {
0043     myReader = new TTreeReader(R3);
0044   }
0045   TTreeReaderValue<Float_t> FEE(*myReader, "FEE");
0046   TTreeReaderValue<Int_t> FEE_Chan(*myReader, "FEE_Chan");
0047   //TTreeReaderValue<Double_t> x(myReader, "x");
0048   //TTreeReaderValue<Double_t> y(myReader, "y");
0049   TTreeReaderValue<Double_t> phi(*myReader, "phi");
0050   //TTreeReaderValue<Double_t> PadX(myReader, "PadX");
0051   //TTreeReaderValue<Double_t> PadY(myReader, "PadY");
0052   TTreeReaderValue<Double_t> PadR(*myReader, "PadR");
0053   //TTreeReaderValue<Double_t> PadPhi(myReader, "PadPhi");
0054
0055   double phiOffset;
0056   if (sector < 12) phiOffset = 2*M_PI*((double)sector/12.) - 2*M_PI*(1./12.)/2;
0057   else phiOffset = 2*M_PI*(((double)sector-12.)/12.) - 2*M_PI*(1./12.)/2;
0058
0059   while(myReader->Next())
0060   {
0061     if (slot[inFee]-1 == *FEE && channel == *FEE_Chan)
0062     {
0063       double x = (*PadR/10)*std::cos(*phi+phiOffset);
0064       double y = (*PadR/10)*std::sin(*phi+phiOffset);
0065       TVector2 pos;
0066       pos.SetX(x); pos.SetY(y);
0067       return pos;
0068     }
0069   }
0070   TVector2 pos;
0071   pos.SetX(0); pos.SetY(0);
0072   return pos;
0073 }
0074
0075
0076 void TPCEventDisplay(const float adcCut = 100.,
0077                      const string &outfile = "/sphenix/user/rosstom/test/TPCEventDisplay_10684",
0078                      const string &indir = "/sphenix/user/rosstom/test/testFiles/",
0079                      const string &runName = "beam-00010684",
0080              const bool writeAllHits = true){
0081
0082   //output nTuple
0083   TString* outputfilename=new TString(outfile+".root");
0084   TFile* outputfile=new TFile(outputfilename->Data(),"recreate");
0085   TTree* hitTree=new TTree("hitTree","x,y,z position and max ADC value for all samples passing cut");
0086   double x_Pos, y_Pos, z_Pos;
0087   int max_adc;
0088   hitTree->Branch("x_Pos",&x_Pos,"x_Pos/D");
0089   hitTree->Branch("y_Pos",&y_Pos,"y_Pos/D");
0090   hitTree->Branch("z_Pos",&z_Pos,"z_Pos/D");
0091   hitTree->Branch("max_adc",&max_adc,"max_adc/I");
0092
0093   std::ofstream outdata;
0094   outdata.open((outfile+".json").c_str());
0095   if( !outdata ) { // file couldn't be opened
0096       cerr << "ERROR: file could not be opened" << endl;
0097       exit(1);
0098   }
0099
0100   string runNumberString = runName;
0101   size_t pos = runNumberString.find("000");
0102   runNumberString.erase(runNumberString.begin(),runNumberString.begin()+pos+3);
0103
0104   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;
0105
0106   outdata << "    \"META\": {\n       \"HITS\": {\n          \"INNERTRACKER\": {\n              \"type\": \"3D\",\n              \"options\": {\n              \"size\": 5,\n              \"color\": 16777215\n              } \n          },\n" << endl;
0107   outdata << "          \"TRACKHITS\": {\n              \"type\": \"3D\",\n              \"options\": {\n              \"size\": 5,\n              \"transparent\": 0.5,\n              \"color\": 16777215\n              } \n          },\n" << endl;
0108   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;
0109   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;
0110   outdata << "    \"TRACKHITS\": [\n\n ";
0111
0112   TTree *R1 = new TTree("R1","TPC Sector Mapping for R1");
0113   TTree *R2 = new TTree("R2","TPC Sector Mapping for R1");
0114   TTree *R3 = new TTree("R3","TPC Sector Mapping for R1");
0115
0116   R1->ReadFile("/sphenix/u/rosstom/calibrations/TPC/Mapping/PadPlane/AutoPad-R1-RevA.sch.ChannelMapping.csv",
0117                "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",',');
0118   R2->ReadFile("/sphenix/u/rosstom/calibrations/TPC/Mapping/PadPlane/AutoPad-R2-RevA-Pads.sch.ChannelMapping.csv",
0119                "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",',');
0120   R3->ReadFile("/sphenix/u/rosstom/calibrations/TPC/Mapping/PadPlane/AutoPad-R3-RevA.sch.ChannelMapping.csv",
0121                "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",',');
0122
0123   bool firstClus = true;
0124
0125   //int framevals[12] = {77,67,69,71,70,79,77,70,69,71,64,77};
0126
0127   for (int q = 0; q < 24; q++)
0128   {
0129     cout << "Processing sector " << q << endl;
0130
0131     string sectorNumber;
0132     if (q < 10) sectorNumber = "0"+std::to_string(q);
0133     else sectorNumber = std::to_string(q);
0134
0135     string fileName = indir+"TPC_ebdc"+sectorNumber+"_"+runName+"-0005.prdf_TPCRawDataTree.root";
0136
0137     TFile *TPCFile = TFile::Open(fileName.c_str());
0138     TTreeReader myReader("SampleTree", TPCFile);
0139
0140     TTreeReaderValue<Int_t> nSamples(myReader, "nSamples");
0141     TTreeReaderValue<Int_t> fee(myReader, "fee");
0142     TTreeReaderArray<UShort_t> adcSamples(myReader, "adcSamples");
0143     TTreeReaderValue<Int_t> Channel(myReader, "Channel");
0144     TTreeReaderValue<Int_t> BCO(myReader, "BCO");
0145     TTreeReaderValue<Int_t> frame(myReader, "frame");
0146     TTreeReaderValue<Int_t> checksumError(myReader, "checksumError");
0147
0148     while(myReader.Next())
0149     {
0150       //if (*frame != framevals[q] && *frame != framevals[q]+1) continue;
0151
0152       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};
0153
0154       if (*nSamples <= 1) continue;
0155
0156       float adcMax = 0;
0157       float adcMaxPos = 0;
0158       if (writeAllHits && *checksumError == 0)
0159       {
0160         for(int adc_sam_no=0;adc_sam_no<*nSamples;adc_sam_no++)
0161         {
0162           if (adcSamples[adc_sam_no] - adcSamples[0] > adcCut)
0163           {
0164             max_adc = adcSamples[adc_sam_no];
0165             double zPos;
0166             if (q < 12)
0167             { //0.4 cm per sample, 50 ns/sample, 8 cm/us drift speed
0168               zPos = 105. - 0.4*(double)adc_sam_no; // - (*BCO - globalEventBCO)*(8./9.4);
0169               //if (zPos < 0. || zPos > 105.) continue;
0170         }
0171             else
0172             {
0173               zPos = -105. + 0.4*(double)adc_sam_no; // + (*BCO - globalEventBCO)*(8./9.4);
0174               //if (zPos > 0. || zPos < 105.) continue;
0175             }
0176
0177             TVector2 binXY = getBinPosition(q,*fee,*Channel,R1,R2,R3);
0178
0179             if (binXY.X() == 0 || binXY.Y() == 0) continue;
0180
0181             x_Pos=binXY.X();
0182             y_Pos=binXY.Y();
0183             z_Pos=zPos;
0184             max_adc=adcMax;
0185             hitTree->Fill();
0186
0187             stringstream spts;
0188
0189             if (firstClus) firstClus = false;
0190             else spts << ",";
0191
0192         spts << "{ \"x\": ";
0193             spts << binXY.X();
0194             spts << ", \"y\": ";
0195             spts << binXY.Y();
0196             spts << ", \"z\": ";
0197             spts << zPos;
0198             spts << ", \"e\": ";
0199             spts << adcMax;
0200             spts << "}";
0201
0202             outdata << (boost::format("%1%") % spts.str());
0203             spts.clear();
0204             spts.str("");
0205           }
0206         }
0207       }
0208       else
0209       {
0210         for(int adc_sam_no=0;adc_sam_no<*nSamples;adc_sam_no++)
0211         {
0212           if (adc_sam_no == 0)
0213       {
0214         adcMax = adcSamples[adc_sam_no];
0215         adcMaxPos = adc_sam_no;
0216       }
0217           else
0218           {
0219         if (adcSamples[adc_sam_no] > adcMax)
0220         {
0221           adcMax = adcSamples[adc_sam_no];
0222               adcMaxPos = adc_sam_no;
0223         }
0224           }
0225         }
0226         if (adcMax > adcCut)
0227         {
0228       double zPos;
0229           //adcMaxPos = 0;
0230           if (q < 12)
0231           { //0.4 cm per sample, 50 ns/sample, 8 cm/us drift speed
0232             zPos = 105. - 0.4*(double)adcMaxPos; // - (*BCO - globalEventBCO)*(8./9.4);
0233             if (zPos < 0. || zPos > 105.) continue;
0234       }
0235           else
0236           {
0237             zPos = -105. + 0.4*(double)adcMaxPos; // + (*BCO - globalEventBCO)*(8./9.4);
0238             if (zPos > 0. || zPos < -105.) continue;
0239           }
0240
0241           TVector2 binXY = getBinPosition(q,*fee,*Channel,R1,R2,R3);
0242
0243           if (binXY.X() == 0 || binXY.Y() == 0) continue;
0244
0245           x_Pos=binXY.X();
0246           y_Pos=binXY.Y();
0247           z_Pos=zPos;
0248           max_adc=adcMax;
0249           hitTree->Fill();
0250
0251           stringstream spts;
0252
0253           if (firstClus) firstClus = false;
0254           else spts << ",";
0255
0256           spts << "{ \"x\": ";
0257           spts << binXY.X();
0258           spts << ", \"y\": ";
0259           spts << binXY.Y();
0260           spts << ", \"z\": ";
0261           spts << zPos;
0262           spts << ", \"e\": ";
0263           spts << adcMax;
0264           spts << "}";
0265
0266           outdata << (boost::format("%1%") % spts.str());
0267           spts.clear();
0268           spts.str("");
0269         }
0270       }
0271     }
0272   }
0273   outdata << "],\n    \"JETS\": [\n         ]\n    }," << endl;
0274   outdata << "\"TRACKS\": {" << endl;
0275   outdata <<"\""<<"INNERTRACKER"<<"\": [";
0276   outdata << "]" << endl;
0277   outdata << "}" << endl;
0278   outdata << "}" << endl;
0279
0280   outdata.close();
0281
0282   outputfile->cd();
0283   hitTree->Write();
0284   outputfile->Close();
0285 }