LaserClusterQA/src/laserQA_multiple.cc

0001 #include "laserQA_multiple.h"
0002
0003 #include <trackbase/LaserClusterContainer.h>
0004 #include <trackbase/LaserCluster.h>
0005 #include <trackbase/TpcDefs.h>
0006
0007
0008 #include <ffaobjects/EventHeader.h>
0009
0010 #include <fun4all/Fun4AllReturnCodes.h>
0011
0012 #include <phool/PHCompositeNode.h>
0013 #include <phool/getClass.h>
0014 #include <phool/phool.h>
0015
0016 //#include <TFile.h>
0017 //#include <TTree.h>
0018 #include <TVector3.h>
0019
0020 #include <string>
0021
0022
0023
0024 laserQA_multiple::laserQA_multiple(const std::string& name)
0025 : SubsysReco(name)
0026 {
0027 }
0028
0029 int laserQA_multiple::InitRun(PHCompositeNode* topNode)
0030 {
0031
0032   outputFile = new TFile(m_output.c_str(),"RECREATE");
0033
0034   //hitHist = new TH3D("hitHist",";layer;iphi;it;adc",48,6.5,57.5,1501,-0.5,1500.5,426,-0.5,425.5);
0035
0036   T = new TTree("T","T");
0037   T->Branch("event",&event);
0038   T->Branch("hits",&hitHist);
0039   T->Branch("clusters",&clustHist);
0040   T->Branch("adc",&adc);
0041   T->Branch("layer",&layer);
0042   T->Branch("iphi",&iphi);
0043   T->Branch("it",&it);
0044   T->Branch("phi",&phi);
0045   T->Branch("R",&R);
0046   T->Branch("fitMode",&fitMode);
0047
0048   //T->Branch("laserCluster",&lc);
0049
0050   return Fun4AllReturnCodes::EVENT_OK;
0051
0052 }
0053
0054
0055 int laserQA_multiple::process_event(PHCompositeNode* topNode)
0056 {
0057
0058   LaserClusterContainer *lcc = findNode::getClass<LaserClusterContainer>(topNode, "LASER_CLUSTER");
0059   if (!lcc)
0060   {
0061     std::cout << PHWHERE << "LASER_CLUSTER Node missing, abort." << std::endl;
0062     return Fun4AllReturnCodes::ABORTRUN;
0063   }
0064
0065   EventHeader *eventHeader = findNode::getClass<EventHeader>(topNode, "EventHeader");
0066   if (!eventHeader)
0067   {
0068     std::cout << PHWHERE << " EventHeader Node missing, abort" << std::endl;
0069     return Fun4AllReturnCodes::ABORTRUN;
0070   }
0071
0072   event = eventHeader->get_EvtSequence();
0073
0074
0075   if(lcc->size() < 1000) return Fun4AllReturnCodes::EVENT_OK;
0076
0077   //hitHist = nullptr;
0078
0079   //bool fillTree = false;
0080
0081   //std::cout << "good event " << std::endl;
0082
0083
0084   std::cout << "there are " << lcc->size() << " clusters in this event" << std::endl;
0085
0086   std::vector<std::vector<TrkrDefs::cluskey>> clusters;
0087
0088   auto clusrange = lcc->getClusters();
0089   for (auto cmitr = clusrange.first;
0090        cmitr != clusrange.second;
0091        ++cmitr)
0092   {
0093
0094     const auto& [cmkey, cmclus_orig] = *cmitr;
0095     //LaserCluster* cmclus = cmclus_orig;
0096     //      bool side = (bool) TpcDefs::getSide(cmkey);
0097     //      const unsigned int nhits = cmclus->getNhits();
0098
0099     //layer = cmclus->getLayer();
0100     //if(cmclus->getNIPhi() <= 2 || (cmclus->getNLayers() == 1 && (layer == 22 || layer == 38 || layer == 54))) continue;
0101
0102     clusters.push_back(std::vector<TrkrDefs::cluskey>{cmkey});
0103
0104     if(clusters.size() == 131)
0105     {
0106       std::cout << "Cluster 131 hits:" << std::endl;
0107       LaserCluster* cmclus = cmclus_orig;
0108       for(unsigned int i=0; i<cmclus->getNhits(); i++){
0109     std::cout << "   hit " << i << ": (" << cmclus->getHitLayer(i) << ", " << cmclus->getHitIPhi(i) << ", " << cmclus->getHitIT(i) << ", " << cmclus->getHitAdc(i) << ")" << std::endl;
0110       }
0111     }
0112
0113   }
0114
0115   for(unsigned int i=0; i<clusters.size(); i++)
0116   {
0117
0118     LaserCluster* cmclus = lcc->findCluster(clusters[i][0]);
0119     float li = cmclus->getLayer();
0120     float pi = cmclus->getIPhi();
0121     float ti = cmclus->getIT();
0122
0123     for(unsigned int j=0; j<clusters.size(); j++)
0124     {
0125       if(i == j) continue;
0126       if(TpcDefs::getSide(clusters[i][0]) != TpcDefs::getSide(clusters[j][0])) continue;
0127       if(TpcDefs::getSectorId(clusters[i][0]) != TpcDefs::getSectorId(clusters[j][0])) continue;
0128       LaserCluster* cmclus2 = lcc->findCluster(clusters[j][0]);
0129       float lj = cmclus2->getLayer();
0130       float pj = cmclus2->getIPhi();
0131       float tj = cmclus2->getIT();
0132
0133       if(fabs(li - lj) < 5 && fabs(pi - pj) < 5 && fabs(ti - tj) < 8)
0134       {
0135     clusters[i].push_back(clusters[j][0]);
0136       }
0137
0138     }
0139   }
0140
0141
0142   for(unsigned int i=0; i<clusters.size(); i++)
0143   {
0144     LaserCluster* cmclusI = lcc->findCluster(clusters[i][0]);
0145     layer = cmclusI->getLayer();
0146     iphi = cmclusI->getIPhi();
0147     it = cmclusI->getIT();
0148     adc = cmclusI->getAdc();
0149
0150     fitMode = cmclusI->getFitMode();
0151
0152     hitHist = new TH3D("hitHist",";layer;iphi;it;adc",9,round(layer)-4.5,round(layer)+4.5,9,round(iphi)-4.5,round(iphi)+4.5,15,round(it)-7.5,round(it)+7.5);
0153     clustHist = new TH3D("clustHist",";layer;iphi;it;adc",9,round(layer)-4.5,round(layer)+4.5,9,round(iphi)-4.5,round(iphi)+4.5,15,round(it)-7.5,round(it)+7.5);
0154
0155     for(unsigned int j=0; j<clusters[i].size(); j++)
0156     {
0157
0158       LaserCluster* cmclus = lcc->findCluster(clusters[i][j]);
0159
0160       unsigned int nHits = cmclus->getNhits();
0161
0162       for(unsigned int h=0; h<nHits; h++)
0163       {
0164     //std::cout << "   hit num: " << i << std::endl;
0165     hitHist->Fill(cmclus->getHitLayer(h), cmclus->getHitIPhi(h), cmclus->getHitIT(h), cmclus->getHitAdc(h));
0166     clustHist->Fill(cmclus->getHitLayer(h), cmclus->getHitIPhi(h), cmclus->getHitIT(h), j+1);
0167       }
0168       //hitHist = (TH3D*)tmpHitHist->Clone();
0169
0170     }
0171
0172     T->Fill();
0173
0174     if(hitHist)
0175     {
0176       delete hitHist;
0177     }
0178
0179     if(clustHist)
0180     {
0181       delete clustHist;
0182     }
0183
0184   }
0185
0186   return Fun4AllReturnCodes::EVENT_OK;
0187 }
0188
0189 int laserQA_multiple::End(PHCompositeNode* /*topNode*/)
0190 {
0191     outputFile->cd();
0192
0193     T->Write();
0194
0195     outputFile->Close();
0196
0197     return Fun4AllReturnCodes::EVENT_OK;
0198 }