DAQ/macros/TPCPedestalCalibration.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
0018 /*************************************************************/
0019 /*                TPC Pedestal Calibration                   */
0020 /*               Thomas Marshall,Aditya Dash                 */
0021 /*        rosstom@ucla.edu,aditya55@physics.ucla.edu         */
0022 /*************************************************************/
0023
0024 void TPCPedestalCalibration(vector<string> inputFiles = {"/sphenix/user/rosstom/test/testFiles/TPC_ebdc00_pedestal-00010131-0000.prdf_TPCRawDataTree.root"}){
0025   for (int fileNum = 0; fileNum < inputFiles.size(); fileNum++){
0026     string sectorNum = inputFiles[fileNum];
0027     size_t pos = sectorNum.find("ebdc");
0028     sectorNum.erase(pos+6);
0029     sectorNum.erase(sectorNum.begin(),sectorNum.begin()+pos+4);
0030     Int_t sector;
0031     if (sectorNum.substr(0,1) == "0")
0032     {
0033       sectorNum.erase(sectorNum.begin(),sectorNum.end()-1);
0034       sector = std::stoi(sectorNum);
0035     }
0036     else sector = std::stoi(sectorNum);
0037
0038     Int_t 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};
0039     Int_t slot_arr[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};
0040
0041     TFile *pedestalInput = TFile::Open(inputFiles[fileNum].c_str());
0042     TTreeReader myReader("SampleTree", pedestalInput);
0043
0044     TTreeReaderValue<Int_t> nSamples_in(myReader, "nSamples");
0045     TTreeReaderValue<Int_t> fee_in(myReader, "fee");
0046     TTreeReaderArray<UShort_t> adcSamples_in(myReader, "adcSamples");
0047     TTreeReaderValue<Int_t> Channel_in(myReader, "Channel");
0048
0049 //Adding the output tree
0050     TString* outputfilename=new TString("pedestalCalibration_TPC_ebdc"+sectorNum+".root");
0051     TFile* outputfile=new TFile(outputfilename->Data(),"recreate");
0052     TTree* outputTree=new TTree("outputTree","outputTree");
0053     Int_t isAlive=1; // 1 if the channel is working properly, 0 if no signal(number of samples is 0, adc value is 0 or nan), pedestal above 200 or below 10
0054     Float_t pedMean; // average pedestal value over all samples
0055     Float_t pedStd; // pedestal standard deviation over all samples
0056     Int_t chan; // channel number
0057     Int_t fee; // fee number
0058     Int_t module; // module number (e.g. R1, R2, R3)
0059     Int_t slot; // slot number
0060     outputTree->Branch("isAlive",&isAlive,"isAlive/I");
0061     outputTree->Branch("pedMean",&pedMean,"pedMean/F");
0062     outputTree->Branch("pedStd",&pedStd,"pedStd/F");
0063     outputTree->Branch("sector",&sector,"sector/I");
0064     outputTree->Branch("fee",&fee,"fee/I");
0065     outputTree->Branch("chan",&chan,"chan/I");
0066     outputTree->Branch("module",&module,"module/I");
0067     outputTree->Branch("slot",&slot,"slot/I");
0068
0069     Float_t ave_adc_fee_channel[26][256];
0070     Float_t std_adc_fee_channel[26][256];
0071     Float_t counts_adc_fee_channel[26][256];
0072     Int_t alive_array_fee_channel[26][256];
0073     for(int fee_no=0;fee_no<26;fee_no++){
0074         for(int channel_no=0;channel_no<256;channel_no++)
0075     {
0076             ave_adc_fee_channel[fee_no][channel_no]=0.0;
0077             std_adc_fee_channel[fee_no][channel_no]=0.0;
0078             counts_adc_fee_channel[fee_no][channel_no]=0.0;
0079         alive_array_fee_channel[fee_no][channel_no]=1;
0080         }
0081     }
0082
0083     while(myReader.Next())
0084         {
0085             if(*nSamples_in==0)
0086         {
0087         alive_array_fee_channel[*fee_in][*Channel_in]=0;
0088         continue;
0089         }
0090
0091             bool dead = false;
0092             for(int adc_sam_no=0;adc_sam_no<*nSamples_in;adc_sam_no++)
0093             {
0094                 if (adcSamples_in[adc_sam_no] == 0 || TMath::IsNaN(float(adcSamples_in[adc_sam_no])))
0095                 {
0096             alive_array_fee_channel[*fee_in][*Channel_in]=0;
0097         }
0098             }
0099             if (dead) {continue;}
0100
0101             for(int adc_sam_no=0;adc_sam_no<*nSamples_in;adc_sam_no++){
0102                 ave_adc_fee_channel[*fee_in][*Channel_in]+=adcSamples_in[adc_sam_no];
0103                 std_adc_fee_channel[*fee_in][*Channel_in]+=pow(adcSamples_in[adc_sam_no],2);
0104                 counts_adc_fee_channel[*fee_in][*Channel_in]+=1.0;
0105             }
0106
0107         }
0108
0109     for(int fee_no=0;fee_no<26;fee_no++)
0110     {
0111           for(int channel_no=0;channel_no<256;channel_no++)
0112           {
0113             if(counts_adc_fee_channel[fee_no][channel_no] != 0.0)
0114             {
0115                 float temp1 = ave_adc_fee_channel[fee_no][channel_no]/counts_adc_fee_channel[fee_no][channel_no];
0116                 float temp2 = std_adc_fee_channel[fee_no][channel_no]/counts_adc_fee_channel[fee_no][channel_no];
0117                 ave_adc_fee_channel[fee_no][channel_no] = temp1;
0118                 std_adc_fee_channel[fee_no][channel_no] = temp2;
0119             }
0120             else
0121         {
0122                 ave_adc_fee_channel[fee_no][channel_no] = 0.0;
0123                 std_adc_fee_channel[fee_no][channel_no] = 0.0;
0124         alive_array_fee_channel[fee_no][channel_no]=0;
0125             }
0126
0127             if(ave_adc_fee_channel[fee_no][channel_no] > 200 || ave_adc_fee_channel[fee_no][channel_no] < 10)
0128             {
0129         alive_array_fee_channel[fee_no][channel_no]=0;
0130             }
0131
0132         pedMean=ave_adc_fee_channel[fee_no][channel_no];
0133         pedStd=sqrt(std_adc_fee_channel[fee_no][channel_no] - pow(ave_adc_fee_channel[fee_no][channel_no],2));
0134         isAlive=alive_array_fee_channel[fee_no][channel_no];
0135         chan=channel_no;
0136             fee=fee_no;
0137         module=mod_arr[fee_no];
0138         slot=slot_arr[fee_no];
0139         outputTree->Fill();
0140           }
0141       }
0142       outputfile->cd();
0143       outputTree->Write();
0144       outputfile->Close();
0145       pedestalInput->Close();
0146   }
0147 }