TPCGainDebug/src/TPCGainDebug.cc

0001 #include "TPCGainDebug.h"
0002
0003 #include <trackbase/ActsGeometry.h>
0004 // #include <trackbase/ClusterErrorPara.h>
0005 // #include <trackbase/TrkrCluster.h>
0006 #include <trackbase/TrkrDefs.h>
0007 #include <trackbase/TrkrHit.h>
0008 #include <trackbase/TrkrHitSetContainer.h>
0009 #include <trackbase/TrackFitUtils.h>
0010
0011 #include <trackbase_historic/TrackSeedContainer_v1.h>
0012
0013 #include <trackbase/TpcDefs.h>
0014
0015 // #include <trackbase/TrkrClusterContainer.h>
0016 // #include <trackbase/TrkrClusterHitAssoc.h>
0017 // #include <trackbase/TrkrClusterIterationMapv1.h>
0018 #include <trackbase/TrkrHitSet.h>
0019
0020 #include <globalvertex/GlobalVertex.h>
0021 #include <globalvertex/GlobalVertexMap.h>
0022 #include <globalvertex/SvtxVertex.h>
0023 #include <globalvertex/SvtxVertexMap.h>
0024
0025 #include <trackbase_historic/TrackAnalysisUtils.h>
0026 #include <trackbase_historic/ActsTransformations.h>
0027 #include <trackbase_historic/SvtxTrack.h>
0028 #include <trackbase_historic/SvtxTrackMap.h>
0029 #include <trackbase_historic/TrackSeed.h>
0030
0031 #include <g4detectors/PHG4TpcCylinderGeom.h>
0032 #include <g4detectors/PHG4TpcCylinderGeomContainer.h>
0033 #include <trackermillepedealignment/HelicalFitter.h>
0034
0035 #include <fun4all/Fun4AllReturnCodes.h>
0036 #include <fun4all/SubsysReco.h>
0037
0038 #include <phool/PHTimer.h>
0039 #include <phool/getClass.h>
0040 #include <phool/phool.h>
0041 #include <phool/recoConsts.h>
0042
0043 #include <TFile.h>
0044 #include <TH1F.h>
0045 #include <TH1I.h>
0046 #include <TNtuple.h>
0047 #include <TVector3.h>
0048 #include <TEnv.h>
0049
0050 #include <cmath>
0051 #include <iomanip>
0052 #include <iostream>
0053 #include <iterator>
0054 #include <map>
0055 #include <memory>  // for shared_ptr
0056 #include <set>     // for _Rb_tree_cons...
0057 #include <utility>
0058 #include <vector>
0059
0060 using namespace std;
0061
0062 // TF1* cleverGaus(TH1* h, const char* title="h", Float_t c = 2.5, bool quietMode=true);
0063 TPCGainDebug::TPCGainDebug(const string& /*name*/, const string& filename, const int runnumber)
0064   : SubsysReco("TPCGainDebug")
0065   , _ievent(0)
0066   , _iseed(0)
0067   , m_fSeed(NAN)
0068   , _save_ntuple(true)
0069   , _do_gain(true)
0070   , _save_debugHist(true)
0071   , _do_pedSigmaCut(true)
0072   // , _do_gain(true)
0073   , _save_sharkFin(false)
0074   , _isSharkFin(false)
0075   // , _isSharkFinForRun(false)
0076   // , _ntp_hit(nullptr)
0077   , _filename(filename)
0078   , _runnumber(runnumber)
0079   , _tfile(nullptr)
0080   , _timer(nullptr)
0081   , AllHitADCThreshold(-99999)
0082 {
0083 }
0084
0085 TPCGainDebug::~TPCGainDebug()
0086 {
0087   if (Verbosity() >= 1)
0088     cout << "TPCGainDebug:: DONE "  << endl;
0089   delete _timer;
0090   delete gaintree;
0091   delete _ntp_hit;
0092   // delete pedhist;
0093 }
0094
0095 int TPCGainDebug::Init(PHCompositeNode* /*topNode*/)
0096 {
0097   _ievent = 0;
0098   m_run = _runnumber;
0099
0100   // pedhist = new TH1F("pedhist", "pedhist", 251, -0.5, 1000.5);
0101
0102   _tfile = new TFile(_filename.c_str(), "RECREATE");
0103   _tfile->SetCompressionLevel(7);
0104
0105   if (_save_ntuple || _save_sharkFin){
0106     _ntp_hit = new TTree("ntp_allhit",Form("all hits > %.2f", AllHitADCThreshold));
0107     _ntp_hit->Branch("run",&m_run);
0108     _ntp_hit->Branch("event",&m_event);
0109     _ntp_hit->Branch("side",&a_side);
0110     _ntp_hit->Branch("sector",&a_sector);
0111     _ntp_hit->Branch("layer",&a_layer);
0112     _ntp_hit->Branch("phibin",&a_phibin);
0113     _ntp_hit->Branch("phi",&a_phi);
0114     _ntp_hit->Branch("pedmean",&a_pedmean);
0115     _ntp_hit->Branch("pedwidth",&a_pedwidth);
0116     _ntp_hit->Branch("tbin",&a_tbin);
0117     _ntp_hit->Branch("adc",&a_adc);
0118     _ntp_hit->Branch("x",&a_x);
0119     _ntp_hit->Branch("y",&a_y);
0120     _ntp_hit->Branch("z",&a_z);
0121     if (_save_sharkFin){
0122       _ntp_hit->Branch("isSharkFin",&a_isSharkFin);
0123     }
0124     // _ntp_hit = new TNtuple("ntp_allhit", "svtxhit => max truth",
0125     //     "run:event:seed:hitID:e:adc:layer:sector:side:"
0126     //     "phibin:tbin:phi:x:y:z"
0127     //     );
0128     a_tbin.reserve(360);
0129     a_adc.reserve(360);
0130     a_x.reserve(360);
0131     a_y.reserve(360);
0132     a_z.reserve(360);
0133   }
0134
0135   if(_save_sharkFin){
0136     h1I_event_shark = new TH1I("h1I_event_shark", ";event number;", 25000, 0, 25000);
0137
0138   }
0139   if(_do_gain){
0140     gaintree = new TTree("tgain","TPC Gain Tree");
0141     gaintree->Branch("run",&m_run);
0142     gaintree->Branch("event",&m_event);
0143     gaintree->Branch("side",&m_side);
0144     gaintree->Branch("sector",&m_sector);
0145     gaintree->Branch("layer",&m_layer);
0146     gaintree->Branch("phibin",&m_phibin);
0147     gaintree->Branch("tbin",&m_tbin);
0148     gaintree->Branch("adcsum",&m_adcsum);
0149     gaintree->Branch("pathlength",&m_pathlength);
0150     gaintree->Branch("phisize",&m_phisize);
0151     // gaintree->Branch("errorcode",&m_errorcode);
0152   }
0153
0154   // float tbin_i = 100;
0155   // float tbin_f = 250;
0156   // h1F_tbin = new TH1F("h1F_tbin", ";tbin;", tbin_f-tbin_i,tbin_i,tbin_f);
0157   // h1F_tbin_prev = (TH1F*) h1F_tbin->Clone("h1F_tbin_prev");
0158   // h1F_tbin_next = (TH1F*) h1F_tbin->Clone("h1F_tbin_next");
0159   _timer = new PHTimer("_eval_timer");
0160   _timer->stop();
0161
0162   return Fun4AllReturnCodes::EVENT_OK;
0163 }
0164
0165 int TPCGainDebug::InitRun(PHCompositeNode*)
0166 {
0167   cout << "Fun4AllServer::BeginRun: InitRun for TPCGainDebug" << endl;
0168   return Fun4AllReturnCodes::EVENT_OK;
0169 }
0170
0171 int TPCGainDebug::process_event(PHCompositeNode* topNode)
0172 {
0173   if ((Verbosity() > 1))
0174   {
0175     cout << "TPCGainDebug::process_event test" << endl;
0176   }
0177
0178   if(_ievent<startevt || _ievent>endevt){
0179     ++_ievent;
0180     return Fun4AllReturnCodes::DISCARDEVENT;
0181   }
0182   _isSharkFin = false;
0183   sf_side = 0;
0184   sf_layer = 0;
0185   sf_phibin = 0;
0186
0187   if ((Verbosity() > 1) && (_ievent % 100 == 0))
0188   {
0189     cout << "TPCGainDebug::process_event - Event = " << _ievent << endl;
0190   }
0191
0192   recoConsts* rc = recoConsts::instance();
0193   if (rc->FlagExist("RANDOMSEED"))
0194   {
0195     _iseed = rc->get_IntFlag("RANDOMSEED");
0196     m_fSeed = _iseed;
0197   }
0198   else
0199   {
0200     _iseed = 0;
0201     m_fSeed = NAN;
0202   }
0203   if (Verbosity() > 1)
0204   {
0205     cout << "TPCGainDebug::process_event - Seed = " << _iseed << endl;
0206   }
0207
0208   ResetEvent(topNode);
0209   //-----------------------------------
0210   // calculate gain!
0211   //-----------------------------------
0212   CalculateGain(topNode);
0213
0214   //---------------------------
0215   // fill the NTuples
0216   //---------------------------
0217   if(_save_sharkFin && _isSharkFin){
0218     if (Verbosity() > 1)
0219     {
0220       cout << "TPCGainDebug::save shark fin "  << endl;
0221     }
0222
0223     fillOutputNtuplesAllhits(topNode);
0224     h1I_event_shark->Fill(_ievent+_eventOffset-1);
0225   }
0226
0227   //--------------------------------------------------
0228   // Print out the ancestry information for this event
0229   //--------------------------------------------------
0230
0231   // printOutputInfo(topNode);
0232
0233   _ievent++;
0234   return Fun4AllReturnCodes::EVENT_OK;
0235 }
0236
0237 int TPCGainDebug::End(PHCompositeNode* /*topNode*/)
0238 {
0239   _tfile->cd();
0240   if (_do_gain){ gaintree->Write();}
0241   if (_save_ntuple || _save_sharkFin){ _ntp_hit->Write();}
0242   if (_save_sharkFin){
0243     h1I_event_shark->Write();
0244   }
0245
0246   TEnv outConfig;
0247   outConfig.SetValue("do_pedSigmaCut", _do_pedSigmaCut);
0248   outConfig.SetValue("adcSigmaThreshold", adcSigmaThreshold);
0249   outConfig.SetValue("adcThreshold", adcThreshold);
0250   outConfig.SetValue("ntp_allhit_adcThreshold", AllHitADCThreshold);
0251   outConfig.SetValue("clusterPhibinArea", clusterPhibinArea);
0252   outConfig.SetValue("clusterTbinArea", clusterTbinArea);
0253   outConfig.SetValue("seedPhibinDistance", seedPhibinDistance);
0254   outConfig.SetValue("seedTbinDistance", seedTbinDistance);
0255   outConfig.SetValue("AdcClockPeriod", AdcClockPeriod);
0256   outConfig.SetValue("binWidthForPedEst", binWidthForPedEst);
0257   if(!_save_sharkFin)
0258     outConfig.Write("config", TObject::kOverwrite);
0259   _tfile->Close();
0260
0261   delete _tfile;
0262
0263   if (Verbosity() > 1)
0264   {
0265     cout << "========================= TPCGainDebug::End() ============================" << endl;
0266     cout << " " << _ievent << " events of output written to: " << _filename << endl;
0267     cout << "===========================================================================" << endl;
0268   }
0269
0270   return Fun4AllReturnCodes::EVENT_OK;
0271 }
0272
0273 //////////////////////////////////////////////////////////
0274 void TPCGainDebug::CalculateGain(PHCompositeNode* topNode)
0275 {
0276   if (Verbosity() > 1)
0277   {
0278     cout << "TPCGainDebug::CalculateGain() entered" << endl;
0279   }
0280
0281   // auto m_tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0282
0283   if (Verbosity() > 1)
0284   {
0285     cout << "Calculating TPC gain, time reset" << endl;
0286     _timer->restart();
0287   }
0288   // need things off of the DST...
0289   TrkrHitSetContainer* hitmap = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0290   if (!hitmap)
0291   {
0292     std::cout << PHWHERE << "ERROR: Can't find hitmap TRKR_HITSET" << std::endl;
0293     return;
0294   }
0295
0296   PHG4TpcCylinderGeomContainer* geom_container =
0297     findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
0298   if (!geom_container)
0299   {
0300     std::cout << PHWHERE << "ERROR: Can't find node CYLINDERCELLGEOM_SVTX" << std::endl;
0301     return;
0302   }
0303   auto m_tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0304
0305   TH1F pedhist("pedhist", "pedhist", 325, -0.5, 1299.5);
0306   // TH1F pedhist("pedhist", "pedhist", 251, -0.5, 1000.5);
0307   // pedhist->Reset();
0308   float hpedestal = 0;
0309   float hpedwidth = 0;
0310   m_event = _ievent+_eventOffset;
0311   for(int iside=0; iside< maxside; ++iside){
0312     for(int isector=0; isector< maxsector; ++isector){
0313       resizeHitsPerLayer();
0314       clearHitsPerLayer();
0315       std::vector<unsigned int> noisyChannel;
0316       for (int ilayer=0; ilayer < maxlayer; ++ilayer){
0317         const TrkrDefs::hitsetkey hitset_key = TpcDefs::genHitSetKey(ilayer, isector, iside);
0318
0319         if (TrkrDefs::getTrkrId(hitset_key) != TrkrDefs::TrkrId::tpcId){
0320           if (Verbosity() > 1){
0321             cout << "The hitsetkey is not associated with TPC" << endl;
0322           }
0323           continue;
0324         }
0325
0326         TrkrHitSet* hitset = hitmap->findHitSet(hitset_key);
0327         if(!hitset){
0328           if (Verbosity() > 3){
0329             cout << "No hits found in TPC (side " << iside << ", sector " << isector << ", layer " << ilayer << ")" << endl;
0330           }
0331           continue;
0332         }
0333
0334         PHG4TpcCylinderGeom* GeoLayer_local = geom_container->GetLayerCellGeom(ilayer);
0335         double radius = GeoLayer_local->get_radius();
0336         unsigned short NTBins = (unsigned short) GeoLayer_local->get_zbins();
0337         if (Verbosity() > 2){
0338           if(NTBins!=360) cout << "NTBins = " << NTBins << endl;
0339           cout << "YJTEST (side " << iside << ", sector " << isector << ", layer " << ilayer << ")"  << ", phibins: "
0340             << GeoLayer_local->get_phibins() << ", "
0341             << GeoLayer_local->get_binning() << ", "
0342             << endl;
0343         }
0344
0345         // for (uint16_t sampleNum = 0; sampleNum < sampleMax; sampleNum++)
0346         // {
0347         //   auto hit_key_temp = TpcDefs::genHitKey(currentPhiBin, (unsigned int) sampleNum);
0348         //   auto hit_temp = hitset->getHit(hit_key_temp);
0349
0350         //   uint16_t adc = hit_temp->getAdc();
0351         //   pedhist.Fill(adc);
0352         // }
0353
0354         if (Verbosity() > 2){
0355           cout << "TrkrHitSet Start!! "<< endl;
0356         }
0357
0358         int currentPhiBin = -1;
0359         float currentPhi = -999;
0360         bool isFirst = true;
0361         int testn = 0;
0362         TrkrHitSet::ConstRange hitrangei = hitset->getHits();
0363         for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
0364             hitr != hitrangei.second;
0365             ++hitr)
0366         {
0367           TrkrDefs::hitkey hit_key = hitr->first;
0368           TrkrHit* hit = hitr->second;
0369           float _phibin = (float) TpcDefs::getPad(hit_key);
0370           float _phi = GeoLayer_local->get_phicenter(_phibin);
0371
0372           // calculate pedestals
0373           if(currentPhiBin != _phibin){
0374
0375             //fill all hits tree when
0376             if(a_tbin.size()!=0 && !isFirst && !_save_sharkFin){
0377               a_side = iside;
0378               a_sector = isector;
0379               a_layer = ilayer;
0380               a_phibin = currentPhiBin;
0381               a_phi = currentPhi;
0382               a_pedmean = hpedestal;
0383               a_pedwidth = hpedwidth;
0384
0385               _ntp_hit->Fill();
0386
0387               a_tbin.clear();
0388               a_adc.clear();
0389               a_x.clear();
0390               a_y.clear();
0391               a_z.clear();
0392               if (Verbosity() > 2)
0393               {
0394                 cout << "TPCGainDebug:: STORE TEST all hits in (side, sector, layer = " << iside << ", " << isector << ", " << ilayer << ") phibin = " << currentPhiBin << ", testn = " << testn<< endl;
0395               }
0396             }
0397             testn = 0;
0398
0399             //this is the start of a new phibin!
0400             currentPhiBin = _phibin;
0401             currentPhi = _phi;
0402             hpedestal = 0;
0403             hpedwidth = 0;
0404             pedhist.Reset();
0405
0406             //fill ADC distribution
0407             uint16_t finaladc = 0;
0408             for (uint16_t sampleNum = 0; sampleNum < sampleMax; sampleNum++)
0409             {
0410               auto hit_key_temp = TpcDefs::genHitKey(currentPhiBin, (unsigned int) sampleNum);
0411               auto hit_temp = hitset->getHit(hit_key_temp);
0412
0413               uint16_t adc = hit_temp->getAdc();
0414               pedhist.Fill(adc);
0415               if(sampleNum==sampleMax-1) finaladc=adc;
0416             }
0417             int hmaxbin = pedhist.GetMaximumBin();
0418
0419             // calc peak position
0420             double adc_sum = 0.0;
0421             double ibin_sum = 0.0;
0422             double ibin2_sum = 0.0;
0423
0424             // calculate pedestal mean and sigma
0425             if(pedhist.GetMaximum()==360 || pedhist.GetEntries()==0){
0426               hpedestal = finaladc;
0427               hpedwidth = 0;
0428             } else{
0429               for (int isum = -1*binWidthForPedEst; isum <= binWidthForPedEst; isum++)
0430               {
0431                 float val = pedhist.GetBinContent(hmaxbin + isum);
0432                 float center = pedhist.GetBinCenter(hmaxbin + isum);
0433                 ibin_sum += center * val;
0434                 ibin2_sum += center * center * val;
0435                 adc_sum += val;
0436               }
0437               hpedestal = ibin_sum / adc_sum;
0438               hpedwidth = sqrt(ibin2_sum / adc_sum - (hpedestal * hpedestal));
0439             }
0440
0441             //noisy channel masking
0442             if(hpedwidth<0.5 || hpedestal<10){
0443               unsigned int temp_noisyKey = _phibin*100 + ilayer*1;
0444               noisyChannel.push_back(temp_noisyKey);
0445               if (Verbosity() > 2)
0446               {
0447                 cout << "TPCGainDebug:: noisy channel of " << std::fixed << temp_noisyKey << " phibin, layer, sector, side = " << _phibin << ", " << ilayer << ", " << isector << ", " << iside << endl;
0448               }
0449             }
0450
0451             //////////////////////////////////////////////
0452             // check if this event is shark fin
0453
0454             if(_save_sharkFin){
0455               int sharkFinBinPrevTbin = 0;
0456               // bool signChange = false;
0457               int sharkFinBinCount = 0;
0458               for (uint16_t sampleNum = 0; sampleNum < sampleMax; sampleNum++)
0459               {
0460                 auto hit_key_temp = TpcDefs::genHitKey(currentPhiBin, (unsigned int) sampleNum);
0461                 auto hit_temp = hitset->getHit(hit_key_temp);
0462
0463                 uint16_t adc = hit_temp->getAdc();
0464                 if((adc < hpedestal-hpedwidth*5) && hpedwidth!=0 && hpedwidth > 1){
0465                   if(sharkFinBinPrevTbin+2 == sampleNum){
0466                     sharkFinBinCount++;
0467                   }
0468                   sharkFinBinPrevTbin = sampleNum;
0469
0470                   if (Verbosity() > 1)
0471                   {
0472                     cout << "TPCGainDebug1:: sharkFin  phibin, layer, sector, side = " << _phibin << ", " << ilayer << ", " << isector << ", " << iside << " :: pedestal, width, adc = " << hpedestal << ", " << hpedwidth << ", " << adc << " :: sharkFinCount, sharkFinBinPrevTbin, tbin = " << sharkFinBinCount << ", " << sharkFinBinPrevTbin << ", " << sampleNum << endl;
0473                   }
0474                 }
0475
0476                 if(sharkFinBinPrevTbin+1 == sampleNum){
0477                   if (Verbosity() > 1)
0478                   {
0479                     cout << "TPCGainDebug2:: sharkFin  phibin, layer, sector, side = " << _phibin << ", " << ilayer << ", " << isector << ", " << iside << " :: pedestal, width, adc = " << hpedestal << ", " << hpedwidth << ", " << adc << " :: sharkFinCount, sharkFinBinPrevTbin, tbin = " << sharkFinBinCount << ", " << sharkFinBinPrevTbin << ", " << sampleNum << endl;
0480                   }
0481                   // if( (adc > hpedestal-hpedwidth*5) && hpedwidth!=0 && hpedwidth > 1) {
0482                   //   signChange= true;
0483                   // } else{
0484                   //   signChange= false;
0485                   //   sharkFinBinCount = 0;
0486                   // }
0487                   if( !((adc > hpedestal+hpedwidth*5) && hpedwidth!=0 && hpedwidth > 1)) {
0488                     sharkFinBinCount = 0;
0489                   }
0490                 }
0491
0492                 if(sharkFinBinCount >= 3){
0493                   if (Verbosity() > 1)
0494                   {
0495                     cout << "TPCGainDebug3:: sharkFin  phibin, layer, sector, side = " << _phibin << ", " << ilayer << ", " << isector << ", " << iside << " :: pedestal, width, adc = " << hpedestal << ", " << hpedwidth << ", " << adc << " :: sharkFinCount, sharkFinBinPrevTbin, tbin = " << sharkFinBinCount << ", " << sharkFinBinPrevTbin << ", " << sampleNum << endl;
0496                   }
0497                 // if(signChange == true && sharkFinBinCount >= 3){
0498                   _isSharkFin=true;
0499                   sf_side = iside;
0500                   sf_layer = ilayer;
0501                   sf_phibin = _phibin;
0502
0503                   if (Verbosity() > 1)
0504                   {
0505                     cout << "TPCGainDebug:: sharkFin found " << endl;
0506                   }
0507                   return;
0508                 }
0509               }
0510             }
0511             // if(isSharkFin) break;
0512           }//if(currentPhiBin != _phibin) // the calculation of pedestal should be done once per each phibin
0513
0514           /// pedestal subtraction
0515           // float adc = hit->getAdc();
0516           float adc = hit->getAdc() - hpedestal;
0517           ihit _hit;
0518           _hit.adc = adc;
0519           _hit.phibin = _phibin;
0520           _hit.tbin = (float) TpcDefs::getTBin(hit_key);
0521
0522           if (Verbosity() >= 3)
0523           {
0524             cout << "TPCGainDebug:: TEST all hits in (side, sector, layer = " << iside << ", " << isector << ", " << ilayer << ") = adc " << adc << ", phibin = " << _phibin << endl;
0525           }
0526
0527           if((_save_ntuple && adc>AllHitADCThreshold) || (_save_sharkFin && (adc < -1*hpedwidth*adcSigmaThreshold))){
0528             float x = NAN;
0529             float y = NAN;
0530             float z = NAN;
0531
0532             double zdriftlength = _hit.tbin * m_tGeometry->get_drift_velocity() * AdcClockPeriod;
0533             // convert z drift length to z position in the TPC
0534             //    cout << " tbin: " << tbin << " vdrift " <<m_tGeometry->get_drift_velocity() << " l drift: " << zdriftlength  <<endl;
0535             double m_tdriftmax = AdcClockPeriod * NTBins / 2.0;
0536             double clusz = (m_tdriftmax * m_tGeometry->get_drift_velocity()) - zdriftlength;
0537             if (iside == 0)
0538             {
0539               clusz = -clusz;
0540             }
0541             z = clusz;
0542             x = radius * cos(_phi);
0543             y = radius * sin(_phi);
0544
0545             if (Verbosity() > 3)
0546             {
0547               cout << "TPCGainDebug:: TEST all hits in (side, sector, layer = " << iside << ", " << isector << ", " << ilayer << ") phibin = " << currentPhiBin
0548                 << ", tbin = " << _hit.tbin
0549                 << ", adc = " << _hit.adc
0550                 << endl;
0551             }
0552             a_tbin.push_back(_hit.tbin);
0553             a_adc.push_back(hit->getAdc());
0554             a_x.push_back(x);
0555             a_y.push_back(y);
0556             a_z.push_back(z);
0557             testn++;
0558           }//if(currentPhiBin != _phibin)
0559           isFirst=false;
0560
0561
0562           if(!_do_gain) continue;
0563
0564           // remove some noise channels for the gain seed finding!
0565           if(hpedwidth<=0.5) continue; // to remove stuck channels (all same ADC values throughout tbin)
0566           if(_do_pedSigmaCut){
0567             if(adc <= (hpedwidth*adcSigmaThreshold)) continue;
0568             // int temp_noisyKey = phibin*100000 + ilayer*1000 + isector*10 + iside;
0569             unsigned int temp_noisyKey = _phibin*100 + ilayer*1;
0570             if(std::find(noisyChannel.begin(), noisyChannel.end(), temp_noisyKey) != noisyChannel.end()){
0571               if (Verbosity() > 1)
0572               {
0573                 cout << "TPCGainDebug:: skip this seed for this noisy channel of " << temp_noisyKey << endl;
0574               }
0575               continue;
0576             }
0577           } else{
0578             if(adc<=adcThreshold) continue;
0579           }
0580           if (Verbosity() > 2)
0581           {
0582             cout << "TPCGainDebug:: TEST all hits in (side, sector, layer = " << iside << ", " << isector << ", " << ilayer << ") phibin = " << currentPhiBin << ", tbin = " << _hit.tbin <<  ", adc = " << adc;
0583             cout << ", hpedestal, hpedwidth = " << hpedestal << ", " << hpedwidth << endl;
0584           }
0585
0586
0587           // if the seed is in other position, add it
0588           // (taking into account of multiple tracks in one sector in a given event)
0589           // also find the local maximum
0590           bool isLocal = true;
0591           int temp_seed_size = hpl_seed[ilayer].size();
0592           if(temp_seed_size==0){
0593             hpl_seed[ilayer].push_back(_hit);
0594           } else{
0595             for(int iseed=0; iseed<temp_seed_size; ++iseed){
0596               if(  (abs(hpl_seed[ilayer][iseed].phibin - _hit.phibin) < seedPhibinDistance)
0597                   && (abs(hpl_seed[ilayer][iseed].tbin - _hit.tbin) < seedTbinDistance) )
0598               {
0599                 isLocal = true;
0600                 if(hpl_seed[ilayer][iseed].adc < _hit.adc){
0601                   hpl_seed[ilayer][iseed].adc = _hit.adc;
0602                   hpl_seed[ilayer][iseed].phibin = _hit.phibin;
0603                   hpl_seed[ilayer][iseed].tbin = _hit.tbin;
0604                 }
0605               } else{
0606                 isLocal = false;
0607               }
0608             } // seed loop
0609             if(!isLocal)
0610               hpl_seed[ilayer].push_back(_hit);
0611           }
0612
0613         }// hits
0614
0615         //fill all hits tree for the last phibin in (layer, sector, side)
0616         if(a_tbin.size()!=0 && !_save_sharkFin){
0617           a_side = iside;
0618           a_sector = isector;
0619           a_layer = ilayer;
0620           a_phibin = currentPhiBin;
0621           a_phi = currentPhi;
0622           a_pedmean = hpedestal;
0623           a_pedwidth = hpedwidth;
0624
0625           _ntp_hit->Fill();
0626
0627           a_tbin.clear();
0628           a_adc.clear();
0629           a_x.clear();
0630           a_y.clear();
0631           a_z.clear();
0632           if (Verbosity() > 3)
0633           {
0634             cout << "TPCGainDebug:: STORE TEST all hits in (side, sector, layer = " << iside << ", " << isector << ", " << ilayer << ") phibin = " << currentPhiBin << ", testn = " << testn<< endl;
0635           }
0636         }
0637         testn = 0;
0638
0639       }//layer
0640
0641       if(!_do_gain) continue;
0642
0643       // print out all gain seeds!
0644       if (Verbosity() > 1){
0645         cout << " >>> DONE Finding seeds "
0646           << " for side " << iside << ", sector " << isector << endl;
0647         for (int ilayer=0; ilayer < maxlayer; ++ilayer){
0648           if(hpl_seed[ilayer].size() ==0) continue;
0649           for(unsigned int it=0; it< hpl_seed[ilayer].size(); it++){
0650             cout << " >>> Seed adc, phibin, tbin for layer " << ilayer << " = "
0651               << hpl_seed[ilayer][it].adc << ", "
0652               << hpl_seed[ilayer][it].phibin << ", "
0653               << hpl_seed[ilayer][it].tbin << ", "
0654               << endl;
0655           }
0656         }
0657       }
0658
0659       //-----------------------------------
0660       // calculate gain
0661       //-----------------------------------
0662       for (int ilayer=1; ilayer < maxlayer-1; ++ilayer){
0663         unsigned int seedsize_prev = hpl_seed[ilayer-1].size();
0664         unsigned int seedsize_next = hpl_seed[ilayer+1].size();
0665         if(seedsize_prev==0 || seedsize_next==0) continue;
0666
0667         PHG4TpcCylinderGeom* GeoLayer_local = geom_container->GetLayerCellGeom(ilayer);
0668         float radius = GeoLayer_local->get_radius();
0669
0670         std::vector<igain> gain;
0671         for (unsigned int j=0; j < seedsize_prev; ++j){//n-1
0672           for (unsigned int k=0; k < seedsize_next; ++k){//n+1
0673             //////////////// track should be perpendicular to the sectors, no tangential tracks
0674             if(abs(hpl_seed[ilayer+1][k].phibin - hpl_seed[ilayer-1][j].phibin)>seedPhibinDistance) continue;
0675             if(abs(hpl_seed[ilayer+1][k].tbin - hpl_seed[ilayer-1][j].tbin)>seedTbinDistance) continue;
0676
0677             // h1F_tbin->Reset();
0678             // h1F_tbin_prev->Reset();
0679             // h1F_tbin_next->Reset();
0680
0681             // if these seeds in layer n-1 and n+1 seem like from different tracks,
0682             float phibin_n = (hpl_seed[ilayer+1][k].phibin + hpl_seed[ilayer-1][j].phibin)/2.;
0683             float tbin_n = (hpl_seed[ilayer+1][k].tbin + hpl_seed[ilayer-1][j].tbin)/2.;
0684
0685             // look for if we already calculated
0686             bool gainfound=false;
0687             int temp_gain_size = gain.size();
0688             if(temp_gain_size >= 1){
0689               for(int ig = 0; ig < temp_gain_size; ig++){
0690                 if((abs(gain[ig].phibin-phibin_n) < seedPhibinDistance) && (abs(gain[ig].tbin-tbin_n) < seedTbinDistance)){
0691                   gainfound = true;
0692                   break;
0693                 }
0694               }
0695             }
0696             if(gainfound) continue;
0697
0698
0699             //////////////////////////////
0700             // find maxADC per each pad around the seeds in n-th layer
0701             const TrkrDefs::hitsetkey hitset_key = TpcDefs::genHitSetKey(ilayer, isector, iside);
0702             if (TrkrDefs::getTrkrId(hitset_key) != TrkrDefs::TrkrId::tpcId){
0703               continue;
0704             }
0705             TrkrHitSet* hitset = hitmap->findHitSet(hitset_key);
0706             if(!hitset){
0707               continue;
0708             }
0709
0710             if (Verbosity() > 2){
0711               cout << "TrkrHitSet Start!! "<< endl;
0712             }
0713
0714             std::map<float, ADCInfo> maxADCperPad; // pad, pair of (maxADC, tbin), highest ADC value over samples (time bins)
0715             // std::map<float, std::pair<float, float>> maxADCperPad; // pad, pair of (maxADC, tbin), highest ADC value over samples (time bins)
0716             TrkrHitSet::ConstRange hitrangei = hitset->getHits();
0717
0718             int currentPhiBin = -1;
0719             for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
0720                 hitr != hitrangei.second;
0721                 ++hitr)
0722             {
0723               TrkrDefs::hitkey hit_key = hitr->first;
0724               TrkrHit* hit = hitr->second;
0725               float phibin = (float) TpcDefs::getPad(hit_key);
0726               if((phibin < (phibin_n - clusterPhibinArea)) || (phibin > (phibin_n + clusterPhibinArea))) continue;
0727               unsigned int temp_noisyKey = phibin*100 + ilayer*1;
0728               if(std::find(noisyChannel.begin(), noisyChannel.end(), temp_noisyKey) != noisyChannel.end()){
0729                 continue;
0730                 if (Verbosity() > 1)
0731                 {
0732                   cout << "TPCGainDebug:: skip this cluster for this noisy channel of " << temp_noisyKey << endl;
0733                 }
0734               }
0735
0736               // calculate pedestals
0737               if(currentPhiBin != phibin){
0738
0739                 //this is the start of a new phibin!
0740                 currentPhiBin = phibin;
0741                 hpedestal = 0;
0742                 hpedwidth = 0;
0743                 pedhist.Reset();
0744
0745                 //fill ADC distribution
0746                 uint16_t finaladc = 0;
0747                 for (uint16_t sampleNum = 0; sampleNum < sampleMax; sampleNum++)
0748                 {
0749                   auto hit_key_temp = TpcDefs::genHitKey(currentPhiBin, (unsigned int) sampleNum);
0750                   auto hit_temp = hitset->getHit(hit_key_temp);
0751
0752                   uint16_t adc = hit_temp->getAdc();
0753                   pedhist.Fill(adc);
0754                   if(sampleNum==sampleMax-1) finaladc=adc;
0755                 }
0756                 int hmaxbin = pedhist.GetMaximumBin();
0757
0758                 // calc peak position
0759                 double adc_sum = 0.0;
0760                 double ibin_sum = 0.0;
0761                 double ibin2_sum = 0.0;
0762
0763                 // calculate pedestal mean and sigma
0764                 if(pedhist.GetMaximum()==360 || pedhist.GetEntries()==0){
0765                   hpedestal = finaladc;
0766                   hpedwidth = 0;
0767                 } else{
0768                   for (int isum = -1*binWidthForPedEst; isum <= binWidthForPedEst; isum++)
0769                   {
0770                     float val = pedhist.GetBinContent(hmaxbin + isum);
0771                     float center = pedhist.GetBinCenter(hmaxbin + isum);
0772                     ibin_sum += center * val;
0773                     ibin2_sum += center * center * val;
0774                     adc_sum += val;
0775                   }
0776                   hpedestal = ibin_sum / adc_sum;
0777                   hpedwidth = sqrt(ibin2_sum / adc_sum - (hpedestal * hpedestal));
0778                 }
0779               }
0780               if(hpedwidth==0) continue;
0781
0782               float tbin = (float) TpcDefs::getTBin(hit_key);
0783               if((tbin < (tbin_n - clusterTbinArea)) || (tbin > (tbin_n + clusterTbinArea))) continue;
0784
0785               ////////////////////////////////////////////
0786               // subtract pedestal from ADC in nth layer
0787               float adc = hit->getAdc()-hpedestal;
0788               ////////////////////////////////////////////
0789
0790               if (maxADCperPad.find(phibin) == maxADCperPad.end() || adc > maxADCperPad[phibin].maxADC) {
0791                 ADCInfo tempADCInfo;
0792                 tempADCInfo.maxADC = adc;
0793                 tempADCInfo.tbin = tbin;
0794                 tempADCInfo.isSignal = (hpedwidth != 0 && adc>=hpedwidth*adcSigmaThreshold_phisize) ? 1 : 0;
0795                 maxADCperPad[phibin] = tempADCInfo;
0796                 // maxADCperPad[phibin] = std::make_pair(adc, tbin);
0797               }
0798
0799               // adcsum += adc;
0800
0801               // if(_save_debugHist){
0802               //   h1F_tbin->Fill(tbin,adc);
0803               // }
0804               if (Verbosity() > 2)
0805               {
0806                 cout << "TPCGainDebug:: all hits around the seeds in (side, sector, layer = " << iside << ", " << isector << ", " << ilayer << ") = adc " << adc << ", phibin = " << phibin << ", tbin = " << tbin << endl;
0807               }
0808             }// hits
0809
0810             if(maxADCperPad.empty()){
0811               continue; // it can't be... but let's be safe
0812             }
0813
0814             //////////////////////////////
0815             // find position of max of maxADC among all pads around the seeds in n-th layer (find the highest phibin, tbin positionin a localized phi and time window)
0816             // also add up ADC values over all pads (phi bins) around the seeds
0817             float adcsum = 0;
0818             float maxADCofAll = 0;
0819             int padcenter = -1;
0820             int tcenter = -1;
0821             int phisize_ = 0;
0822             for (auto it = maxADCperPad.begin(); it != maxADCperPad.end(); ++it) {
0823               float adc_ = it->second.maxADC;
0824               // float adc_ = it->second.first;
0825               adcsum += adc_;
0826               if (Verbosity() > 2)
0827               {
0828                 cout << "TPCGainDebug:: max ADC hit around the seeds in (side, sector, layer = " << iside << ", " << isector << ", " << ilayer << ") for pad " << it->first << " =  adc " << adc_ << ", tbin = " << it->second.tbin << endl;
0829               }
0830               if(maxADCofAll < adc_){
0831                 maxADCofAll = adc_;
0832                 padcenter = it->first;
0833                 tcenter = it->second.tbin;
0834               }
0835               if(it->second.isSignal==1) phisize_++;
0836             }
0837             // if(adcsum==0) continue;
0838
0839             //////////////////////////////
0840             //calculate pathlength
0841             float dphibin = abs(hpl_seed[ilayer+1][k].phibin - hpl_seed[ilayer-1][j].phibin);
0842             float dz = abs(hpl_seed[ilayer+1][k].tbin - hpl_seed[ilayer-1][j].tbin) * m_tGeometry->get_drift_velocity() * AdcClockPeriod;
0843             float drphi = radius*abs(GeoLayer_local->get_phicenter(padcenter)-GeoLayer_local->get_phicenter((padcenter-dphibin<0) ? padcenter+dphibin : padcenter-dphibin));
0844             float dR = 0;// (ilayer <=22 ? 0.6 : 1.2);
0845             if(ilayer==7) dR=0.687;
0846             else if(ilayer<=22 && ilayer%2==0) dR=0.595;
0847             else if(ilayer<=22 && ilayer%2!=0) dR=0.534;
0848             else if(ilayer>=23 && ilayer<=38) dR=1.012;
0849             else if(ilayer>=39) dR=1.088;
0850             else dR=999999;
0851             float _pathlength = dz*dz + drphi*drphi + dR*dR;
0852
0853             if (Verbosity() > 1)
0854             {
0855               cout << "TPCGainDebug:: max of max ADC hit around the seeds in (side, sector, layer = " << iside << ", " << isector << ", " << ilayer << ") ::::: pad = " << padcenter << ", adc = " << maxADCofAll << ", tbin = " << tcenter << ", pathlength = " << _pathlength << endl;
0856             }
0857
0858             // if(_save_debugHist && adcsum!=0){
0859             //   _tfile->cd();
0860             //   h1F_tbin->Write(Form("h1F_tbin_side%d_sector%d_layer%d_event%d", iside, isector, ilayer,_ievent));
0861             //   h1F_tbin_prev->Write(Form("h1F_tbin_side%d_sector%d_layer%d_event%d_prev", iside, isector, ilayer,_ievent));
0862             //   h1F_tbin_next->Write(Form("h1F_tbin_side%d_sector%d_layer%d_event%d_next", iside, isector, ilayer,_ievent));
0863             //   // _tfile->Close();
0864             // }
0865
0866
0867             igain temp_gain;
0868             temp_gain.adc=adcsum;
0869             temp_gain.phibin=padcenter;
0870             temp_gain.tbin=tcenter;
0871             temp_gain.pathlength=_pathlength;
0872             temp_gain.phisize=phisize_;
0873             // temp_gain.errorcode=_errorcode;
0874
0875             gain.push_back(temp_gain);
0876
0877             // delete hitset;
0878
0879           }//n+1 layer
0880         }//n-1 layer
0881
0882         int gain_size= gain.size();
0883         if(gain_size==0) continue;
0884
0885         for(int ig = 0; ig < gain_size; ig++){
0886           m_side.push_back(iside);
0887           m_sector.push_back(isector);
0888           m_layer.push_back(ilayer);
0889           m_phibin.push_back(gain[ig].phibin);
0890           m_tbin.push_back(gain[ig].tbin);
0891           m_adcsum.push_back(gain[ig].adc);
0892           m_pathlength.push_back(gain[ig].pathlength);
0893           m_phisize.push_back(gain[ig].phisize);
0894           // m_errorcode.push_back(_errorcode);
0895
0896           if (Verbosity() > 2)
0897           {
0898             cout << "TPCGainDebug::FINAL ADC sum in layer " << ilayer << " (side, sector, phibin, tbin = " <<  iside << ", " << isector << ", " << gain[ig].phibin << ", " << gain[ig].tbin << ") = " << gain[ig].adc << endl;
0899           }
0900         }//gain size
0901
0902       }//nth layer
0903     }//sector
0904   }//side
0905
0906   if(_do_gain && m_side.size()!=0)
0907     gaintree->Fill();
0908
0909   // if(_save_ntuple && a_side.size()!=0)
0910   //   _ntp_hit->Fill();
0911   // delete hitmap;
0912   // delete geom_container;
0913
0914   return;
0915 }
0916
0917 //////////////////////////////////////////////////////////
0918 void TPCGainDebug::fillOutputNtuplesAllhits(PHCompositeNode* topNode)
0919 {
0920   if (Verbosity() > 1)
0921   {
0922     cout << "TPCGainDebug::fillOutputNtuplesAllhits() entered" << endl;
0923   }
0924
0925   // need things off of the DST...
0926   TrkrHitSetContainer* hitmap = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0927   if (!hitmap)
0928   {
0929     std::cout << PHWHERE << "ERROR: Can't find hitmap TRKR_HITSET" << std::endl;
0930     return;
0931   }
0932
0933   PHG4TpcCylinderGeomContainer* geom_container =
0934     findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
0935   if (!geom_container)
0936   {
0937     std::cout << PHWHERE << "ERROR: Can't find node CYLINDERCELLGEOM_SVTX" << std::endl;
0938     return;
0939   }
0940   auto m_tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0941
0942   TH1F pedhist("pedhist", "pedhist", 325, -0.5, 1299.5);
0943   // TH1F pedhist("pedhist", "pedhist", 251, -0.5, 1000.5);
0944   // pedhist->Reset();
0945   float hpedestal = 0;
0946   float hpedwidth = 0;
0947   m_event = _ievent+_eventOffset;
0948   for(int iside=0; iside< maxside; ++iside){
0949     for(int isector=0; isector< maxsector; ++isector){
0950       for (int ilayer=0; ilayer < maxlayer; ++ilayer){
0951         const TrkrDefs::hitsetkey hitset_key = TpcDefs::genHitSetKey(ilayer, isector, iside);
0952
0953         if (TrkrDefs::getTrkrId(hitset_key) != TrkrDefs::TrkrId::tpcId){
0954           if (Verbosity() > 1){
0955             cout << "The hitsetkey is not associated with TPC" << endl;
0956           }
0957           continue;
0958         }
0959
0960         TrkrHitSet* hitset = hitmap->findHitSet(hitset_key);
0961         if(!hitset){
0962           if (Verbosity() > 3){
0963             cout << "No hits found in TPC (side " << iside << ", sector " << isector << ", layer " << ilayer << ")" << endl;
0964           }
0965           continue;
0966         }
0967
0968         PHG4TpcCylinderGeom* GeoLayer_local = geom_container->GetLayerCellGeom(ilayer);
0969         double radius = GeoLayer_local->get_radius();
0970         unsigned short NTBins = (unsigned short) GeoLayer_local->get_zbins();
0971         if (Verbosity() > 2){
0972           if(NTBins!=360) cout << "NTBins = " << NTBins << endl;
0973           cout << "YJTEST (side " << iside << ", sector " << isector << ", layer " << ilayer << ")"  << ", phibins: "
0974             << GeoLayer_local->get_phibins() << ", "
0975             << GeoLayer_local->get_binning() << ", "
0976             << endl;
0977         }
0978
0979         if (Verbosity() > 2){
0980           cout << "TrkrHitSet Start!! "<< endl;
0981         }
0982
0983         int currentPhiBin = -1;
0984         float currentPhi = -999;
0985         bool isFirst = true;
0986         int testn = 0;
0987         TrkrHitSet::ConstRange hitrangei = hitset->getHits();
0988         for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
0989             hitr != hitrangei.second;
0990             ++hitr)
0991         {
0992           TrkrDefs::hitkey hit_key = hitr->first;
0993           TrkrHit* hit = hitr->second;
0994           float _phibin = (float) TpcDefs::getPad(hit_key);
0995           float _phi = GeoLayer_local->get_phicenter(_phibin);
0996
0997           // calculate pedestals
0998           if(currentPhiBin != _phibin){
0999
1000             //fill all hits tree when
1001             if(a_tbin.size()!=0 && !isFirst){
1002               a_side = iside;
1003               a_sector = isector;
1004               a_layer = ilayer;
1005               a_phibin = currentPhiBin;
1006               a_phi = currentPhi;
1007               a_pedmean = hpedestal;
1008               a_pedwidth = hpedwidth;
1009               if(_save_sharkFin){
1010                 if(sf_side == iside && sf_layer == ilayer && sf_phibin == currentPhiBin) a_isSharkFin = 1;
1011                 else a_isSharkFin = 0;
1012               }
1013
1014               _ntp_hit->Fill();
1015
1016               a_tbin.clear();
1017               a_adc.clear();
1018               a_x.clear();
1019               a_y.clear();
1020               a_z.clear();
1021               if (Verbosity() > 2)
1022               {
1023                 cout << "TPCGainDebug:: STORE TEST all hits in (side, sector, layer = " << iside << ", " << isector << ", " << ilayer << ") phibin = " << currentPhiBin << ", testn = " << testn<< endl;
1024               }
1025             }
1026             testn = 0;
1027
1028             //this is the start of a new phibin!
1029             currentPhiBin = _phibin;
1030             currentPhi = _phi;
1031             hpedestal = 0;
1032             hpedwidth = 0;
1033             pedhist.Reset();
1034
1035             //fill ADC distribution
1036             uint16_t finaladc = 0;
1037             for (uint16_t sampleNum = 0; sampleNum < sampleMax; sampleNum++)
1038             {
1039               auto hit_key_temp = TpcDefs::genHitKey(currentPhiBin, (unsigned int) sampleNum);
1040               auto hit_temp = hitset->getHit(hit_key_temp);
1041
1042               uint16_t adc = hit_temp->getAdc();
1043               pedhist.Fill(adc);
1044               if(sampleNum==sampleMax-1) finaladc=adc;
1045             }
1046             int hmaxbin = pedhist.GetMaximumBin();
1047
1048             // calc peak position
1049             double adc_sum = 0.0;
1050             double ibin_sum = 0.0;
1051             double ibin2_sum = 0.0;
1052
1053             // calculate pedestal mean and sigma
1054             if(pedhist.GetMaximum()==360 || pedhist.GetEntries()==0){
1055               hpedestal = finaladc;
1056               hpedwidth = 0;
1057             } else{
1058               for (int isum = -1*binWidthForPedEst; isum <= binWidthForPedEst; isum++)
1059               {
1060                 float val = pedhist.GetBinContent(hmaxbin + isum);
1061                 float center = pedhist.GetBinCenter(hmaxbin + isum);
1062                 ibin_sum += center * val;
1063                 ibin2_sum += center * center * val;
1064                 adc_sum += val;
1065               }
1066               hpedestal = ibin_sum / adc_sum;
1067               hpedwidth = sqrt(ibin2_sum / adc_sum - (hpedestal * hpedestal));
1068             }
1069           }//if(currentPhiBin != _phibin) // the calculation of pedestal should be done once per each phibin
1070
1071           /// pedestal subtraction
1072           float adc = hit->getAdc() - hpedestal;
1073           ihit _hit;
1074           _hit.adc = adc;
1075           _hit.phibin = _phibin;
1076           _hit.tbin = (float) TpcDefs::getTBin(hit_key);
1077
1078           if (Verbosity() >= 3)
1079           {
1080             cout << "TPCGainDebug:: TEST all hits in (side, sector, layer = " << iside << ", " << isector << ", " << ilayer << ") = adc " << adc << ", phibin = " << _phibin << endl;
1081           }
1082
1083           float x = NAN;
1084           float y = NAN;
1085           float z = NAN;
1086
1087           double zdriftlength = _hit.tbin * m_tGeometry->get_drift_velocity() * AdcClockPeriod;
1088           double m_tdriftmax = AdcClockPeriod * NTBins / 2.0;
1089           double clusz = (m_tdriftmax * m_tGeometry->get_drift_velocity()) - zdriftlength;
1090           if (iside == 0)
1091           {
1092             clusz = -clusz;
1093           }
1094           z = clusz;
1095           x = radius * cos(_phi);
1096           y = radius * sin(_phi);
1097
1098           if (Verbosity() > 3)
1099           {
1100             cout << "TPCGainDebug:: TEST all hits in (side, sector, layer = " << iside << ", " << isector << ", " << ilayer << ") phibin = " << currentPhiBin
1101               << ", tbin = " << _hit.tbin
1102               << ", adc = " << _hit.adc
1103               << endl;
1104           }
1105           a_tbin.push_back(_hit.tbin);
1106           a_adc.push_back(hit->getAdc());
1107           a_x.push_back(x);
1108           a_y.push_back(y);
1109           a_z.push_back(z);
1110           testn++;
1111
1112           isFirst=false;
1113
1114           if (Verbosity() > 2)
1115           {
1116             cout << "TPCGainDebug:: TEST all hits in (side, sector, layer = " << iside << ", " << isector << ", " << ilayer << ") phibin = " << currentPhiBin << ", tbin = " << _hit.tbin <<  ", adc = " << adc;
1117             cout << ", hpedestal, hpedwidth = " << hpedestal << ", " << hpedwidth << endl;
1118           }
1119         }//hits
1120
1121         //fill all hits tree for the last phibin in (layer, sector, side)
1122         if(a_tbin.size()!=0){
1123           a_side = iside;
1124           a_sector = isector;
1125           a_layer = ilayer;
1126           a_phibin = currentPhiBin;
1127           a_phi = currentPhi;
1128           a_pedmean = hpedestal;
1129           a_pedwidth = hpedwidth;
1130           if(_save_sharkFin){
1131             if(sf_side == iside && sf_layer == ilayer && sf_phibin == currentPhiBin) a_isSharkFin = 1;
1132             else a_isSharkFin = 0;
1133           }
1134
1135           _ntp_hit->Fill();
1136
1137           a_tbin.clear();
1138           a_adc.clear();
1139           a_x.clear();
1140           a_y.clear();
1141           a_z.clear();
1142         }
1143       }//layer
1144     }//sector
1145   }//side
1146 }
1147
1148 // void TPCGainDebug::fillOutputNtuplesAllhits(PHCompositeNode* topNode)
1149 // {
1150 //   if (Verbosity() > 1)
1151 //   {
1152 //     cout << "TPCGainDebug::fillOutputNtuplesAllhits() entered" << endl;
1153 //   }
1154
1155 //   if (_ntp_hit)
1156 //   {
1157 //     auto m_tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
1158
1159 //     if (Verbosity() > 1)
1160 //     {
1161 //       cout << "Filling ntp_hit " << endl;
1162 //       _timer->restart();
1163 //     }
1164 //     // need things off of the DST...
1165 //     TrkrHitSetContainer* hitmap = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
1166
1167 //     PHG4TpcCylinderGeomContainer* geom_container =
1168 //       findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
1169 //     if (!geom_container)
1170 //     {
1171 //       std::cout << PHWHERE << "ERROR: Can't find node CYLINDERCELLGEOM_SVTX" << std::endl;
1172 //       return;
1173 //     }
1174
1175 //     // for(float i=7; i<55; i++){
1176 //     //   PHG4TpcCylinderGeom* geom = geom_container->GetLayerCellGeom(i);
1177 //     //   std::cout << "PHG4TpcCylinderGeom - layer: " << geom->get_layer() << std::endl;
1178 //     //   std::cout
1179 //     //     << "  binnig: " << geom->get_binning()
1180 //     //     << ", radius: " << geom->get_radius()
1181 //     //     << ", nzbins: " << geom->get_zbins()
1182 //     //     << ", zmin: " << geom->get_zmin()
1183 //     //     << ", zstep: " << geom->get_zstep()
1184 //     //     << ", nphibins: " << geom->get_phibins()
1185 //     //     << ", phimin: " << geom->get_phimin()
1186 //     //     << ", phistep: " << geom->get_phistep()
1187 //     //     << ", thickness: " << geom->get_thickness()
1188 //     //     << std::endl;
1189
1190 //     //   // std::cout << "  sector_R_bias: " << geom.sector_R_bias << std::endl;
1191 //     //   // std::cout << "  sector_Phi_bias: " << geom.sector_Phi_bias << std::endl;
1192 //     //   // std::cout << "  sector_min_Phi: " << geom.sector_min_Phi << std::endl;
1193 //     //   // std::cout << "  sector_max_Phi: " << geom.sector_max_Phi << std::endl;
1194 //     //   // cout << GeoLayer_local;
1195 //     // }
1196
1197 //     if (!hitmap)
1198 //     {
1199 //       std::cout << PHWHERE << "ERROR: Can't find hitmap TRKR_HITSET" << std::endl;
1200 //       return;
1201 //     }
1202
1203 //     TrkrHitSetContainer::ConstRange all_hitsets = hitmap->getHitSets();
1204 //     for (TrkrHitSetContainer::ConstIterator iter = all_hitsets.first;
1205 //         iter != all_hitsets.second;
1206 //         ++iter)
1207 //     {
1208 //       const TrkrDefs::hitsetkey hitset_key = iter->first;
1209 //       TrkrHitSet* hitset = iter->second;
1210 //       // get all hits for this hitset
1211 //       TrkrHitSet::ConstRange hitrangei = hitset->getHits();
1212 //       for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
1213 //           hitr != hitrangei.second;
1214 //           ++hitr)
1215 //       {
1216 //         TrkrDefs::hitkey hit_key = hitr->first;
1217 //         TrkrHit* hit = hitr->second;
1218 //         float run = _runnumber;
1219 //         float event = _ievent;
1220 //         float hitID = hit_key;
1221 //         float e = hit->getEnergy();
1222 //         float adc = hit->getAdc();
1223 //         float layer_local = TrkrDefs::getLayer(hitset_key);
1224 //         float sector = TpcDefs::getSectorId(hitset_key);
1225 //         float side = TpcDefs::getSide(hitset_key);
1226
1227 //         if(adc<AllHitADCThreshold) continue;
1228 //         float phibin = NAN;
1229 //         float tbin = NAN;
1230 //         float phi = NAN;
1231 //         float phi_center = NAN;
1232 //         float x = NAN;
1233 //         float y = NAN;
1234 //         float z = NAN;
1235
1236 //         if (TrkrDefs::getTrkrId(hitset_key) == TrkrDefs::TrkrId::tpcId)
1237 //         {
1238 //           PHG4TpcCylinderGeom* GeoLayer_local = geom_container->GetLayerCellGeom(layer_local);
1239 //           double radius = GeoLayer_local->get_radius();
1240 //           phibin = (float) TpcDefs::getPad(hit_key);
1241 //           tbin = (float) TpcDefs::getTBin(hit_key);
1242 //           phi = GeoLayer_local->get_phicenter(phibin);
1243
1244 //           double zdriftlength = tbin * m_tGeometry->get_drift_velocity() * AdcClockPeriod;
1245 //           // convert z drift length to z position in the TPC
1246 //           //    cout << " tbin: " << tbin << " vdrift " <<m_tGeometry->get_drift_velocity() << " l drift: " << zdriftlength  <<endl;
1247 //           unsigned short NTBins = (unsigned short) GeoLayer_local->get_zbins();
1248 //           double m_tdriftmax = AdcClockPeriod * NTBins / 2.0;
1249 //           double clusz = (m_tdriftmax * m_tGeometry->get_drift_velocity()) - zdriftlength;
1250 //           if (side == 0)
1251 //           {
1252 //             clusz = -clusz;
1253 //           }
1254 //           z = clusz;
1255 //           phi_center = GeoLayer_local->get_phicenter(phibin);
1256 //           x = radius * cos(phi_center);
1257 //           y = radius * sin(phi_center);
1258 //         }
1259
1260 //         float hit_data[] = {
1261 //           run,
1262 //           event,
1263 //           (float) _iseed,
1264 //           hitID,
1265 //           e,
1266 //           adc,
1267 //           layer_local,
1268 //           sector,
1269 //           side,
1270 //           (float) phibin,
1271 //           (float) tbin,
1272 //           phi,
1273 //           x,
1274 //           y,
1275 //           z,
1276 //         };
1277
1278 //         _ntp_hit->Fill(hit_data);
1279 //       }
1280 //     }
1281 //     if (Verbosity() > 1)
1282 //     {
1283 //       _timer->stop();
1284 //       cout << "hit time:                " << _timer->get_accumulated_time() / 1000. << " sec" << endl;
1285 //     }
1286 //   }
1287 //   return;
1288 // }
1289
1290 int TPCGainDebug::ResetEvent(PHCompositeNode *)
1291 {
1292   m_side.clear();
1293   m_sector.clear();
1294   m_layer.clear();
1295   m_phibin.clear();
1296   m_tbin.clear();
1297   m_adcsum.clear();
1298   m_pathlength.clear();
1299   m_phisize.clear();
1300
1301   a_tbin.clear();
1302   a_adc.clear();
1303   a_x.clear();
1304   a_y.clear();
1305   a_z.clear();
1306
1307   return Fun4AllReturnCodes::EVENT_OK;
1308 }
1309
1310 //bool SampleFit_PowerLawDoubleExp(        //
1311 //    const std::vector<double> &samples,  //
1312 //    double &peak,                        //
1313 //    double &peak_sample,                 //
1314 //    double &pedestal,                    //
1315 //    std::map<int, double> &parameters_io,
1316 //    const int verbosity)
1317 //{
1318 //  static const int n_parameter = 7;
1319
1320 //  // inital guesses
1321 //  int peakPos = 0.;
1322
1323 //  //  assert(samples.size() == n_samples);
1324 //  const int n_samples = samples.size();
1325
1326 //  TGraph gpulse(n_samples);
1327 //  for (int i = 0; i < n_samples; i++)
1328 //  {
1329 //    (gpulse.GetX())[i] = i;
1330
1331 //    (gpulse.GetY())[i] = samples[i];
1332 //  }
1333
1334 //  //Saturation correction - Abhisek
1335 //  //  for (int ipoint = 0; ipoint < gpulse.GetN(); ipoint++)
1336 //  //    if ((gpulse.GetY())[ipoint] >= ((1 << 10) - 10)  // drop point if touching max or low limit on ADCs
1337 //  //        or (not isnormal((gpulse.GetY())[ipoint])))
1338 //  //    {
1339 //  //      gpulse.RemovePoint(ipoint);
1340 //  //      ipoint--;
1341 //  //    }
1342
1343 //  pedestal = gpulse.GetY()[0];  //(double) PEDESTAL;
1344 //  double peakval = pedestal;
1345 //  const double risetime = 1.5;
1346
1347 //  for (int iSample = 0; iSample < n_samples - risetime * 3; iSample++)
1348 //  {
1349 //    if (abs(gpulse.GetY()[iSample] - pedestal) > abs(peakval - pedestal))
1350 //    {
1351 //      peakval = gpulse.GetY()[iSample];
1352 //      peakPos = iSample;
1353 //    }
1354 //  }
1355 //  peakval -= pedestal;
1356
1357 //  if (verbosity)
1358 //  {
1359 //    cout << "SampleFit_PowerLawDoubleExp - "
1360 //         << "pedestal = " << pedestal << ", "
1361 //         << "peakval = " << peakval << ", "
1362 //         << "peakPos = " << peakPos << endl;
1363 //  }
1364
1365
1366
1367
1368 //double SignalShape_PowerLawExp(double *x, double *par)
1369 //{
1370 //  double pedestal = par[4];
1371 //  //                        + ((x[0] - 1.5 * par[1]) > 0) * par[5];  // quick fix on exting tails on the signal function
1372 //  if (x[0] < par[1])
1373 //    return pedestal;
1374 //  //double  signal = (-1)*par[0]*pow((x[0]-par[1]),par[2])*exp(-(x[0]-par[1])*par[3]);
1375 //  double signal = par[0] * pow((x[0] - par[1]), par[2]) * exp(-(x[0] - par[1]) * par[3]);
1376 //  return pedestal + signal;
1377 //}
1378
1379 //double SignalShape_PowerLawDoubleExp(double *x, double *par)
1380 //{
1381 //  double pedestal = par[4];
1382 //  //                        + ((x[0] - 1.5 * par[1]) > 0) * par[5];  // quick fix on exting tails on the signal function
1383 //  if (x[0] < par[1])
1384 //    return pedestal;
1385 //  //double  signal = (-1)*par[0]*pow((x[0]-par[1]),par[2])*exp(-(x[0]-par[1])*par[3]);
1386 //  //  peak / pow(fits.GetParameter(2) / fits.GetParameter(3), fits.GetParameter(2)) * exp(fits.GetParameter(2)) = fits.GetParameter(0);  // exact peak height is (p0*Power(p2/p3,p2))/Power(E,p2)
1387 //  //  fits.GetParameter(2) / peak_shift =  fits.GetParameter(3);  // signal peak time
1388
1389 //  double signal =                                                                                         //
1390 //    par[0]                                                                                              //
1391 //    * pow((x[0] - par[1]), par[2])                                                                      //
1392 //    * (((1. - par[5]) / pow(par[3], par[2]) * exp(par[2])) * exp(-(x[0] - par[1]) * (par[2] / par[3]))  //
1393 //        + (par[5] / pow(par[6], par[2]) * exp(par[2])) * exp(-(x[0] - par[1]) * (par[2] / par[6]))       //
1394 //      );
1395 //  return pedestal + signal;
1396 //}
1397 //TF1* cleverGaus(TH1* h, const char* title, Float_t c, bool quietMode)
1398 //{
1399 //  if ( h->GetEntries() == 0 )
1400 //  {
1401 //    TF1 *fit0 = new TF1(title,"gaus",-1,1);
1402 //    fit0->SetParameters(0,0,0);
1403 //    return fit0;
1404 //  }
1405
1406 //  Int_t peakBin  = h->GetMaximumBin();
1407 //  Double_t peak =  h->GetBinCenter(peakBin);
1408 //  Double_t sigma = h->GetRMS();
1409
1410 //  TF1 *fit1 = new TF1(title,"gaus",peak-c*sigma,peak+c*sigma);
1411 //  fit1->SetParameter(1, peak);
1412 //  //fit1->SetParameter(1, 0.0005);
1413 //  //if (quietMode) h->Fit(fit1,"Q R");
1414 //  if (quietMode) h->Fit(fit1,"LL M O Q R");
1415 //  else    h->Fit(fit1,"LL M O R");
1416 //  return fit1;
1417 //}