Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:15:08 

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_BCO(const float adcCut = 10., 
0032                      const string &outfile = "TPCEventDisplay_25926_BCO", 
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\", \"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 
0099   //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};
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 
0131   float pads_per_sector[3] = {96, 128, 192};  
0132 
0133   //ULong64_t BCOVal = 128330850912;
0134   //  ULong64_t BCOVal = 128330850911;
0135   ULong64_t BCOVal =128434038131;
0136 
0137   for (int q = 0; q < 24; q++)
0138   {
0139     //if (q == 17) continue;
0140     cout << "Processing sector " << q << endl;
0141   
0142     sectorNum = q;
0143 
0144     string sectorNumber;
0145     if (q < 10) sectorNumber = "0"+std::to_string(q);
0146     else sectorNumber = std::to_string(q);   
0147 
0148     string fileName = indir+"TPC_ebdc"+sectorNumber+"_"+runName+"-0001.prdf_test.root";
0149 
0150     TFile *TPCFile = TFile::Open(fileName.c_str());
0151     TTree *sampleTree = (TTree*)TPCFile->Get("SampleTree");
0152     TTree *taggerTree = (TTree*)TPCFile->Get("TaggerTree");
0153 
0154     Int_t nSamples, fee, Channel, frame, packet, checksumError, taggerFrame;
0155     ULong64_t BCO;
0156     Double_t xPos, yPos;
0157     UShort_t adcSamples[5000];
0158     sampleTree->SetBranchAddress("nSamples",&nSamples);
0159     sampleTree->SetBranchAddress("fee",&fee);
0160     sampleTree->SetBranchAddress("adcSamples",&adcSamples);
0161     sampleTree->SetBranchAddress("Channel",&Channel);
0162     //sampleTree->SetBranchAddress("BCO",&BCO);
0163     sampleTree->SetBranchAddress("frame",&frame);
0164     sampleTree->SetBranchAddress("checksumError",&checksumError);
0165     sampleTree->SetBranchAddress("xPos",&xPos);
0166     sampleTree->SetBranchAddress("yPos",&yPos);
0167 
0168     taggerTree->SetBranchAddress("bco",&BCO);
0169     taggerTree->SetBranchAddress("packet",&packet);
0170     taggerTree->SetBranchAddress("frame",&taggerFrame);
0171 
0172     double channelMean[26][256];
0173     double channelSigma[26][256];
0174     double channelSamples[26][256];
0175     int nHitsChannel[26][256];
0176 
0177     for(int fee_no=0;fee_no<26;fee_no++)
0178     {
0179       for(int channel_no=0;channel_no<256;channel_no++)
0180       {
0181         channelMean[fee_no][channel_no] = 0.;
0182         channelSigma[fee_no][channel_no] = 0.;
0183         channelSamples[fee_no][channel_no] = 0.;
0184         nHitsChannel[fee_no][channel_no] = 0;
0185       }
0186     }
0187 
0188     //TTreeReader myReader("SampleTree", TPCFile);
0189 
0190     /*
0191     TTreeReaderValue<Int_t> nSamples(myReader, "nSamples");
0192     TTreeReaderValue<Int_t> fee(myReader, "fee");
0193     TTreeReaderArray<UShort_t> adcSamples(myReader, "adcSamples");
0194     TTreeReaderValue<Int_t> Channel(myReader, "Channel");
0195     TTreeReaderValue<Int_t> BCO(myReader, "BCO");
0196     TTreeReaderValue<Int_t> frame(myReader, "frame");
0197     TTreeReaderValue<Int_t> checksumError(myReader, "checksumError");
0198     TTreeReaderValue<double> xPosition(myReader, "xPos");
0199     TTreeReaderValue<double> yPosition(myReader, "yPos");
0200     */ 
0201     
0202     std::vector<int> frameValsToCheck;
0203 
0204     for(int i = 0; i < taggerTree->GetEntries(); i++)
0205     {
0206       taggerTree->GetEntry(i);
0207       if ((BCO < BCOVal + 5) && (BCO > BCOVal - 5))
0208       //if (BCO == BCOVal)
0209       {
0210     cout <<"Sector "<<q<<" found packet = "<<packet<<" taggerFrame = "<<taggerFrame<<" with BCO = "<<BCO<<endl;
0211         frameValsToCheck.push_back(taggerFrame);
0212       } 
0213     } 
0214 
0215     for(int i = 0; i < sampleTree->GetEntries(); i++)
0216     {
0217       sampleTree->GetEntry(i);
0218       if (checksumError == 0)
0219       {
0220         for(int adc_sam_no=0;adc_sam_no<10;adc_sam_no++)
0221         {
0222           channelSamples[fee][Channel] += 1.;
0223           channelMean[fee][Channel] += adcSamples[adc_sam_no]; 
0224           channelSigma[fee][Channel] += pow(adcSamples[adc_sam_no],2); 
0225         } 
0226       } 
0227     }
0228    
0229     for(int fee_no=0;fee_no<26;fee_no++)
0230     {
0231       for(int channel_no=0;channel_no<256;channel_no++)
0232       {
0233         double temp1 = channelMean[fee_no][channel_no]/channelSamples[fee_no][channel_no];
0234         double temp2 = channelSigma[fee_no][channel_no]/channelSamples[fee_no][channel_no];
0235         channelMean[fee_no][channel_no] = temp1;
0236         channelSigma[fee_no][channel_no] = sqrt(temp2 - pow(temp1,2));
0237       }
0238     }
0239 
0240     //int nHitsInEntry[1000];
0241     for(int i = 0; i < sampleTree->GetEntries(); i++)
0242     {
0243       sampleTree->GetEntry(i);
0244       //if (frame != framevals[q] && frame != framevals[q]+1) continue;
0245       bool recordHits = false;
0246       for (auto & vals : frameValsToCheck)
0247       {
0248         if (frame < vals + 2 && frame > vals - 1) 
0249     {
0250           recordHits = true;
0251       break;
0252     }
0253       }
0254       //if (!recordHits) continue;
0255       if (checksumError == 0 && channelMean[fee][Channel] > 20. && channelMean[fee][Channel] < 200.)
0256       {
0257         for(int adc_sam_no=0;adc_sam_no<nSamples;adc_sam_no++)
0258         {
0259           if ((adcSamples[adc_sam_no] - channelMean[fee][Channel] > channelSigma[fee][Channel]*4) && (adcSamples[adc_sam_no] - channelMean[fee][Channel] > adcCut))
0260           {
0261             nHitsChannel[fee][Channel] += 1;
0262             //nHitsInEntry[frame] += 1;
0263           }
0264         }
0265       }
0266     }
0267 
0268     int sample_count = 0;
0269     for(int i = 0; i < sampleTree->GetEntries(); i++)
0270     {
0271       sampleTree->GetEntry(i);
0272       //if (*frame > 100) continue;
0273       //if (q == 2 && (*frame < framevals[q] - 5 || *frame > framevals[q]+5)) continue;
0274       //if (frame != framevals[q] && frame != framevals[q]+1) continue;
0275       if (nHitsChannel[fee][Channel] > 50) continue;
0276       //if (nHitsInEntry[frame] > 100) continue;
0277       if (nSamples <= 1) continue;
0278       bool recordHits = false;
0279       for (auto & vals : frameValsToCheck)
0280       {
0281         if (frame < vals + 2 && frame > vals - 1) 
0282     {
0283           recordHits = true;
0284       break;
0285     }
0286       }
0287       if (!recordHits) continue;
0288       ++sample_count;
0289 
0290       float adcMax = 0;
0291       float adcMaxPos = 0;
0292       if (checksumError == 0 && channelMean[fee][Channel] > 20. && channelMean[fee][Channel] < 200.)
0293       {
0294         for(int adc_sam_no=0;adc_sam_no<nSamples;adc_sam_no++)
0295         {
0296           if ((adcSamples[adc_sam_no] - channelMean[fee][Channel] > channelSigma[fee][Channel]*5) && (adcSamples[adc_sam_no] - channelMean[fee][Channel] > adcCut))
0297           {
0298             adcMinusPedestal = adcSamples[adc_sam_no] - channelMean[fee][Channel];
0299             double zPos;
0300             if (q < 12)
0301             { //0.4 cm per sample, 50 ns/sample, 8 cm/us drift speed
0302               zPos = 105. - (105./230.)*((double)adc_sam_no - 10.); // - (*BCO - globalEventBCO)*(8./9.4);
0303               //if (zPos < 0. || zPos > 105.) continue;
0304               side = 0;
0305               phiElement = q;
0306         }
0307             else
0308             {
0309               zPos = -105. + (105./230.)*((double)adc_sam_no - 10.); // + (*BCO - globalEventBCO)*(8./9.4);
0310               //if (zPos > 0. || zPos < 105.) continue;
0311               side = 1;
0312               phiElement = q - 12;
0313             }
0314 
0315             int feeM = FEE_map[fee];
0316             if(mod_arr[fee]==2) feeM += 6;
0317             else if(mod_arr[fee]==3) feeM += 14;
0318             layer = M.getLayer(feeM, Channel);
0319             pad = M.getPad(feeM, Channel) + (phiElement)*pads_per_sector[mod_arr[fee]-1]; 
0320 
0321         
0322         if(layer!=0)
0323           {
0324         double R = M.getR(feeM, Channel);
0325         double phi = (sectorNum<12? -M.getPhi(feeM, Channel) : M.getPhi(feeM, Channel) ) + (sectorNum - side*12. )* M_PI / 6. ;
0326         R /= 10.; //convert mm to cm  
0327  
0328         xPos = R*cos(phi);
0329         yPos = R*sin(phi);
0330           }
0331 
0332             x_Pos = xPos;
0333             y_Pos = yPos;            
0334             z_Pos = zPos;
0335             frameVal = frame;
0336             timeBin = adc_sam_no;
0337             hitTree->Fill();
0338 
0339             stringstream spts;
0340 
0341             if (firstClus) firstClus = false;
0342             else spts << ","; 
0343             
0344         spts << "{ \"x\": ";
0345             spts << x_Pos;
0346             spts << ", \"y\": ";
0347             spts << y_Pos;
0348             spts << ", \"z\": ";
0349             spts << z_Pos;
0350             spts << ", \"e\": ";
0351             spts << adcMinusPedestal;
0352             spts << "}";
0353 
0354             outdata << (boost::format("%1%") % spts.str());
0355             spts.clear();
0356             spts.str("");
0357           }
0358         }
0359       }
0360 
0361     }
0362 
0363     cout <<"sample_count = "<<sample_count<<endl;
0364 
0365     TPCFile->Close();
0366   }
0367   outdata << "],\n    \"JETS\": [\n         ]\n    }," << endl;
0368   outdata << "\"TRACKS\": {" << endl;
0369   outdata <<"\""<<"INNERTRACKER"<<"\": [";
0370   outdata << "]" << endl;
0371   outdata << "}" << endl;
0372   outdata << "}" << endl; 
0373  
0374   outdata.close();    
0375 
0376   outputfile->cd();
0377   hitTree->Write();
0378   outputfile->Close();
0379   std::cout << "Done" << std::endl; 
0380 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!