LaserClusterQA/src/laserClusterQA.cc

0001 #include "laserClusterQA.h"
0002
0003 #include <trackbase/LaserClusterContainerv1.h>
0004 #include <trackbase/LaserClusterv1.h>
0005 #include <trackbase/TpcDefs.h>
0006
0007 #include <tpc/TpcDistortionCorrection.h>
0008 #include <tpc/TpcDistortionCorrectionContainer.h>
0009
0010 #include <ffaobjects/EventHeader.h>
0011
0012 #include <fun4all/Fun4AllReturnCodes.h>
0013
0014 #include <phool/PHCompositeNode.h>
0015 #include <phool/getClass.h>
0016 #include <phool/phool.h>
0017
0018 #include <TFile.h>
0019 #include <TH1.h>
0020 #include <TH2.h>
0021 #include <TVector3.h>
0022
0023 #include <string>
0024
0025
0026
0027 laserClusterQA::laserClusterQA(const std::string& name)
0028 : SubsysReco(name)
0029 {
0030 }
0031
0032 int laserClusterQA::InitRun(PHCompositeNode* topNode)
0033 {
0034
0035     h = new TH2D("h",Form("Pedestal %d;#phi_{reco};R_{reco} [cm]",m_pedestal),1000,0.036,0.044,1000,43,45);
0036
0037
0038     adcSpectrum = new TH1D("adcSpectrum",Form("Pedestal %d;Hit ADC",m_pedestal),1001,-0.5,1000.5);
0039
0040     adcSpecClus = new TH2D("adcSpecClus",Form("Pedestal %d;Hit ADC;Cluster R_{reco} [cm]",m_pedestal),1001,-0.5,1000.5,1000,43,45);
0041     adcSpecHit = new TH2D("adcSpecHit",Form("Pedestal %d;Hit ADC;Hit R_{reco} [cm]",m_pedestal),1001,-0.5,1000.5,1000,43,45);
0042
0043     for(int i=0; i<5; i++){
0044         hSpot[i] = new TH2D(Form("h_%d",i),Form("Spot %d Pedestal %d;#phi_{reco};R_{reco} [cm]",i,m_pedestal),1000,0.036,0.044,1000,43,45);
0045         adcSpecSpot[i] = new TH1D(Form("adcSpecSpot_%d",i),Form("Spot %d Pedestal %d;Hit ADC",i,m_pedestal),1001,-0.5,1000.5);
0046     }
0047
0048     return Fun4AllReturnCodes::EVENT_OK;
0049
0050 }
0051
0052
0053 int laserClusterQA::process_event(PHCompositeNode* topNode)
0054 {
0055
0056     EventHeader *eventHeader = findNode::getClass<EventHeader>(topNode, "EventHeader");
0057     if(!eventHeader)
0058     {
0059         std::cout << PHWHERE << " EventHeader Node missing, abort" << std::endl;
0060         return Fun4AllReturnCodes::ABORTRUN;
0061     }
0062
0063     LaserClusterContainer *lcc = findNode::getClass<LaserClusterContainer>(topNode, "LASER_CLUSTER");
0064     if (!lcc)
0065     {
0066         std::cout << PHWHERE << "LASER_CLUSTER Node missing, abort." << std::endl;
0067         return Fun4AllReturnCodes::ABORTRUN;
0068     }
0069
0070     m_dcc_in_module_edge = findNode::getClass<TpcDistortionCorrectionContainer>(topNode, "TpcDistortionCorrectionContainerModuleEdge");
0071     if (m_dcc_in_module_edge)
0072     {
0073         std::cout << "TpcCentralMembraneMatching:   found TPC distortion correction container module edge" << std::endl;
0074     }
0075
0076     m_dcc_in_static = findNode::getClass<TpcDistortionCorrectionContainer>(topNode, "TpcDistortionCorrectionContainerStatic");
0077     if (m_dcc_in_static)
0078     {
0079         std::cout << "TpcCentralMembraneMatching:   found TPC distortion correction container static" << std::endl;
0080     }
0081
0082     std::cout << "working on event " << eventHeader->get_EvtSequence() << std::endl;
0083
0084     int n1700Event = 0;
0085
0086     auto clusrange = lcc->getClusters();
0087     for (auto cmitr = clusrange.first;
0088         cmitr != clusrange.second;
0089         ++cmitr)
0090     {
0091         const auto& [cmkey, cmclus_orig] = *cmitr;
0092         LaserCluster* cmclus = cmclus_orig;
0093         bool side = (bool) TpcDefs::getSide(cmkey);
0094         const unsigned int nhits = cmclus->getNhits();
0095
0096         bool clus1700 = false;
0097
0098         if(cmclus->getNLayers()<=1 || cmclus->getNIPhi()<=2)
0099         {
0100             continue;
0101         }
0102
0103         Acts::Vector3 tmp_pos(cmclus->getX(), cmclus->getY(), cmclus->getZ());
0104         if (m_dcc_in_module_edge)
0105         {
0106             tmp_pos = m_distortionCorrection.get_corrected_position(tmp_pos, m_dcc_in_module_edge);
0107         }
0108         if (m_dcc_in_static)
0109         {
0110             tmp_pos = m_distortionCorrection.get_corrected_position(tmp_pos, m_dcc_in_static);
0111         }
0112
0113         TVector3 pos(tmp_pos[0], tmp_pos[1], tmp_pos[2]);
0114
0115         double R = pos.Perp();
0116         double phi = pos.Phi();
0117         int spot = -1;
0118
0119         if(side && R > 43.0 && R < 45.0 && phi > 0.036 && phi < 0.044)
0120         {
0121             n1700Event++;
0122             n1700Total++;
0123             clus1700 = true;
0124             h->Fill(phi,R);
0125
0126             if(R > 200*phi + 36.14)
0127             {
0128                 spot = 0;
0129             }
0130             else if(R > -200*phi + 51.96 && phi > 0.0398)
0131             {
0132                 spot = 1;
0133             }
0134             else if(R < -175.44*phi + 50.701)
0135             {
0136                 spot = 3;
0137             }
0138             else if(R < 152.6*phi + 37.661)
0139             {
0140                 spot = 4;
0141             }
0142             else
0143             {
0144                 spot = 2;
0145             }
0146
0147             if(spot != -1)
0148             {
0149                 hSpot[spot]->Fill(phi,R);
0150             }
0151         }
0152
0153         for(int h=0; h<(int)nhits; h++)
0154         {
0155             adcSpectrum->Fill(cmclus->getHitAdc(h));
0156
0157             if(!clus1700)
0158             {
0159                 continue;
0160             }
0161
0162             adcSpecClus->Fill(cmclus->getHitAdc(h),R);
0163             adcSpecHit->Fill(cmclus->getHitAdc(h), sqrt(pow(cmclus->getHitX(h),2) + pow(cmclus->getHitY(h),2)) );
0164
0165             if(spot != -1)
0166             {
0167                 adcSpecSpot[spot]->Fill(cmclus->getHitAdc(h));
0168             }
0169         }
0170
0171     }
0172
0173     std::cout << "Found " << n1700Event << " in this event and " << n1700Total << " total" << std::endl;
0174
0175     return Fun4AllReturnCodes::EVENT_OK;
0176 }
0177
0178 int laserClusterQA::End(PHCompositeNode* /*topNode*/)
0179 {
0180     TFile *outfile = new TFile(m_output.c_str(),"RECREATE");
0181     outfile->cd();
0182
0183     h->Write();
0184     for(int i=0; i<5; i++)
0185     {
0186         hSpot[i]->Write();
0187     }
0188
0189     adcSpectrum->Write();
0190     adcSpecClus->Write();
0191     adcSpecHit->Write();
0192
0193     for(int i=0; i<5; i++)
0194     {
0195         adcSpecSpot[i]->Write();
0196     }
0197
0198     outfile->Close();
0199
0200     return Fun4AllReturnCodes::EVENT_OK;
0201 }