EMCal_pi0_Calib_2023/src/CaloAna.cc

0001 #include "CaloAna.h"
0002
0003 // G4Hits includes
0004 #include <TLorentzVector.h>
0005 #include <g4main/PHG4Hit.h>
0006 #include <g4main/PHG4HitContainer.h>
0007 #include <g4main/PHG4Particle.h>
0008
0009 // G4Cells includes
0010 #include <g4detectors/PHG4Cell.h>
0011 #include <g4detectors/PHG4CellContainer.h>
0012
0013 // Tower includes
0014 #include <calobase/RawCluster.h>
0015 #include <calobase/RawClusterContainer.h>
0016 #include <calobase/RawClusterUtility.h>
0017 #include <calobase/RawTower.h>
0018 #include <calobase/RawTowerContainer.h>
0019 #include <calobase/RawTowerGeom.h>
0020 #include <calobase/RawTowerGeomContainer.h>
0021 #include <calobase/TowerInfo.h>
0022 #include <calobase/TowerInfoContainer.h>
0023 #include <calobase/TowerInfoDefs.h>
0024
0025 // Cluster includes
0026 #include <calobase/RawCluster.h>
0027 #include <calobase/RawClusterContainer.h>
0028
0029 #include <fun4all/Fun4AllHistoManager.h>
0030 #include <fun4all/Fun4AllReturnCodes.h>
0031
0032 #include <phool/getClass.h>
0033
0034 #include <globalvertex/GlobalVertex.h>
0035 #include <globalvertex/GlobalVertexMap.h>
0036
0037 // MBD
0038 #include <mbd/BbcGeom.h>
0039 #include <mbd/MbdPmtContainerV1.h>
0040 #include <mbd/MbdPmtHit.h>
0041
0042 #include <TFile.h>
0043 #include <TH1.h>
0044 #include <TH2.h>
0045 #include <TNtuple.h>
0046 #include <TProfile2D.h>
0047 #include <TTree.h>
0048
0049 #include <Event/Event.h>
0050 #include <Event/packet.h>
0051 #include <cassert>
0052 #include <sstream>
0053 #include <string>
0054
0055 #include <iostream>
0056 #include <utility>
0057 #include "TCanvas.h"
0058 #include "TF1.h"
0059 #include "TFile.h"
0060 #include "TH1F.h"
0061 #include"TRandom3.h"
0062
0063 #include <cdbobjects/CDBTTree.h>  // for CDBTTree
0064 #include <ffamodules/CDBInterface.h>
0065 #include <phool/recoConsts.h>
0066
0067 #include <g4main/PHG4TruthInfoContainer.h>
0068
0069 R__LOAD_LIBRARY(libLiteCaloEvalTowSlope.so)
0070
0071 using namespace std;
0072
0073 CaloAna::CaloAna(const std::string& name, const std::string& filename)
0074   : SubsysReco(name)
0075   , detector("HCALIN")
0076   , outfilename(filename)
0077 {
0078   _eventcounter = 0;
0079 }
0080
0081 CaloAna::~CaloAna()
0082 {
0083   delete hm;
0084   delete g4hitntuple;
0085   delete g4cellntuple;
0086   delete towerntuple;
0087   delete clusterntuple;
0088 }
0089
0090 int CaloAna::Init(PHCompositeNode*)
0091 {
0092   hm = new Fun4AllHistoManager(Name());
0093   // create and register your histos (all types) here
0094
0095   outfile = new TFile(outfilename.c_str(), "RECREATE");
0096
0097   // correlation plots
0098   for (int i = 0; i < 96; i++)
0099   {
0100     h_mass_eta_lt[i] = new TH1F(Form("h_mass_eta_lt%d", i), "", 50, 0, 0.5);
0101     h_mass_eta_lt_rw[i] = new TH1F(Form("h_mass_eta_lt_rw%d", i), "", 50, 0, 0.5);
0102     h_pt_eta[i] = new TH1F(Form("h_pt_eta%d", i), "", 100, 0, 10);
0103     h_pt_eta_rw[i] = new TH1F(Form("h_pt_eta_rw%d", i), "", 100, 0, 10);
0104   }
0105
0106   h_cemc_etaphi = new TH2F("h_cemc_etaphi", "", 96, 0, 96, 256, 0, 256);
0107
0108   // 1D distributions
0109   h_InvMass = new TH1F("h_InvMass", "Invariant Mass", 500, 0, 1.0);
0110   h_InvMass_w = new TH1F("h_InvMass_w", "Invariant Mass", 500, 0, 1.0);
0111   h_InvMassMix = new TH1F("h_InvMassMix", "Invariant Mass", 120, 0, 1.2);
0112
0113   // cluster QA
0114   h_etaphi_clus = new TH2F("h_etaphi_clus", "", 140, -1.2, 1.2, 64, -1 * TMath::Pi(), TMath::Pi());
0115   h_clusE = new TH1F("h_clusE", "", 100, 0, 10);
0116
0117   h_emcal_e_eta = new TH1F("h_emcal_e_eta", "", 96, 0, 96);
0118
0119   h_pt1 = new TH1F("h_pt1", "", 100, 0, 5);
0120   h_pt2 = new TH1F("h_pt2", "", 100, 0, 5);
0121
0122   h_nclusters = new TH1F("h_nclusters", "", 1000, 0, 1000);
0123
0124   h_truth_eta = new TH1F("h_truth_eta", "", 100, -1.2, 1.2);
0125   h_truth_e = new TH1F("h_truth_e", "", 100, 0, 10);
0126   h_truth_pt = new TH1F("h_truth_pt", "", 100, 0, 10);
0127   h_delR_recTrth = new TH1F("h_delR_recTrth", "", 500, 0, 5);
0128   h_matched_res = new TH2F("h_matched_res","",100,0,1.5,20,-1,1);
0129
0130
0131   //////////////////////////
0132   // pT rewieghting
0133   frw = new TFile("/sphenix/u/bseidlitz/work/analysis/EMCal_pi0_Calib_2023/macros/rw_pt.root");
0134   for(int i=0; i<96; i++) h_pt_rw[i] = (TH1F*) frw->Get(Form("h_pt_eta%d", i));
0135
0136   rnd = new TRandom3();
0137
0138
0139   return 0;
0140 }
0141
0142 int CaloAna::process_event(PHCompositeNode* topNode)
0143 {
0144   _eventcounter++;
0145
0146   process_towers(topNode);
0147
0148   return Fun4AllReturnCodes::EVENT_OK;
0149 }
0150
0151 int CaloAna::process_towers(PHCompositeNode* topNode)
0152 {
0153   if ((_eventcounter % 1000) == 0) std::cout << _eventcounter << std::endl;
0154
0155   // float emcaldownscale = 1000000 / 800;
0156
0157   float emcal_hit_threshold = 0.5;  // GeV
0158
0159   // cuts
0160   float maxDr = 1.1;
0161   float maxAlpha = 0.6;
0162   float clus_chisq_cut = 4;
0163   float nClus_ptCut = 0.5;
0164   int max_nClusCount = 300;
0165
0166   //-----------------------get vertex----------------------------------------//
0167
0168   float vtx_z = 0;
0169   if (getVtx)
0170   {
0171     GlobalVertexMap* vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0172     if (!vertexmap)
0173     {
0174       // std::cout << PHWHERE << " Fatal Error - GlobalVertexMap node is missing"<< std::endl;
0175       std::cout << "CaloAna GlobalVertexMap node is missing" << std::endl;
0176       // return Fun4AllReturnCodes::ABORTRUN;
0177     }
0178     if (vertexmap && !vertexmap->empty())
0179     {
0180       GlobalVertex* vtx = vertexmap->begin()->second;
0181       if (vtx)
0182       {
0183         vtx_z = vtx->get_z();
0184       }
0185     }
0186   }
0187
0188
0189
0190   //////////////////////////////////////////////
0191   //         towers
0192   vector<float> ht_eta;
0193   vector<float> ht_phi;
0194
0195   // if (!m_vtxCut || abs(vtx_z) > _vz)  return Fun4AllReturnCodes::EVENT_OK;
0196
0197   TowerInfoContainer* towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC");
0198   if (towers)
0199   {
0200     int size = towers->size();  // online towers should be the same!
0201     for (int channel = 0; channel < size; channel++)
0202     {
0203       TowerInfo* tower = towers->get_tower_at_channel(channel);
0204       float offlineenergy = tower->get_energy();
0205       unsigned int towerkey = towers->encode_key(channel);
0206       int ieta = towers->getTowerEtaBin(towerkey);
0207       int iphi = towers->getTowerPhiBin(towerkey);
0208       bool isGood = !(tower->get_isBadChi2());
0209       if (!isGood && offlineenergy > 0.2)
0210       {
0211         ht_eta.push_back(ieta);
0212         ht_phi.push_back(iphi);
0213       }
0214       if (isGood) h_emcal_e_eta->Fill(ieta, offlineenergy);
0215       if (offlineenergy > emcal_hit_threshold)
0216       {
0217         h_cemc_etaphi->Fill(ieta, iphi);
0218       }
0219     }
0220   }
0221
0222   RawClusterContainer* clusterContainer = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_POS_COR_CEMC");
0223   if (!clusterContainer)
0224   {
0225     std::cout << PHWHERE << "funkyCaloStuff::process_event - Fatal Error - CLUSTER_CEMC node is missing. " << std::endl;
0226     return 0;
0227   }
0228
0229   float weight =1;
0230
0231   /////////////////////////////////////////////////
0232   //// Truth info
0233   float wieght = 1;
0234   PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0235   vector <TLorentzVector> truth_photons;
0236   if (truthinfo)
0237   {
0238     PHG4TruthInfoContainer::Range range = truthinfo->GetPrimaryParticleRange();
0239     for (PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter)
0240     {
0241       // Get truth particle
0242       const PHG4Particle* truth = iter->second;
0243       if (!truthinfo->is_primary(truth)) continue;
0244       TLorentzVector myVector;
0245       myVector.SetXYZM(truth->get_px(), truth->get_py(), truth->get_pz(), 0.13497);
0246
0247       float energy = myVector.E();
0248       h_truth_eta->Fill(myVector.Eta());
0249       h_truth_e->Fill(energy, wieght);
0250       weight = myVector.Pt()*TMath::Exp(-3*myVector.Pt());
0251       h_truth_pt->Fill(myVector.Pt(),weight);
0252       // int id =  truth->get_pid();
0253       if (debug) std::cout << "M=" <<  myVector.M()  << "   E=" << energy << "  pt=" << myVector.Pt() << "  eta=" << myVector.Eta() << std::endl;
0254     }
0255
0256     PHG4TruthInfoContainer::Range second_range = truthinfo->GetSecondaryParticleRange();
0257     float m_g4 = 0;
0258
0259     for (PHG4TruthInfoContainer::ConstIterator siter = second_range.first;
0260           siter != second_range.second; ++siter) {
0261       if (m_g4 >= 19999) break;
0262       // Get photons from pi0 decays
0263       const PHG4Particle *truth = siter->second;
0264
0265       if (truth->get_pid() == 22) {
0266         PHG4Particle *parent = truthinfo->GetParticle(truth->get_parent_id());
0267         if (parent->get_pid() == 111 && parent->get_track_id() > 0) {
0268           float phot_pt = sqrt(truth->get_px() * truth->get_px()
0269                             + truth->get_py() * truth->get_py());
0270           //float phot_pz = truth->get_pz();
0271           float phot_e = truth->get_e();
0272           float phot_phi = atan2(truth->get_py(), truth->get_px());
0273           float phot_eta = atanh(truth->get_pz() / sqrt(truth->get_px()*truth->get_px()+truth->get_py()*truth->get_py()+truth->get_pz()*truth->get_pz()));
0274
0275           TLorentzVector photon = TLorentzVector();
0276           photon.SetPtEtaPhiE(phot_pt, phot_eta, phot_phi,phot_e);
0277           truth_photons.push_back(photon);
0278
0279           if (debug) std::cout<< "pt=" <<  phot_pt <<   " e=" << phot_e << " phi=" << phot_phi << " eta=" << phot_eta << endl;
0280        }
0281      }
0282    }
0283   }
0284
0285
0286   //////////////////////////////////////////
0287   // geometry for hot tower/cluster masking
0288   std::string towergeomnodename = "TOWERGEOM_CEMC";
0289   RawTowerGeomContainer* m_geometry = findNode::getClass<RawTowerGeomContainer>(topNode, towergeomnodename);
0290   if (!m_geometry)
0291   {
0292     std::cout << Name() << "::"
0293               << "CreateNodeTree"
0294               << ": Could not find node " << towergeomnodename << std::endl;
0295     throw std::runtime_error("failed to find TOWERGEOM node in RawClusterDeadHotMask::CreateNodeTree");
0296   }
0297
0298   RawClusterContainer::ConstRange clusterEnd = clusterContainer->getClusters();
0299   RawClusterContainer::ConstIterator clusterIter;
0300   RawClusterContainer::ConstIterator clusterIter2;
0301   int nClusCount = 0;
0302   for (clusterIter = clusterEnd.first; clusterIter != clusterEnd.second; clusterIter++)
0303   {
0304     RawCluster* recoCluster = clusterIter->second;
0305
0306     CLHEP::Hep3Vector vertex(0, 0, vtx_z);
0307     CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*recoCluster, vertex);
0308
0309     float clus_pt = E_vec_cluster.perp();
0310     float clus_chisq = recoCluster->get_chi2();
0311
0312     if (clus_pt < nClus_ptCut) continue;
0313     if (clus_chisq > clus_chisq_cut) continue;
0314
0315     nClusCount++;
0316   }
0317
0318   h_nclusters->Fill(nClusCount);
0319
0320   if (nClusCount > max_nClusCount) return Fun4AllReturnCodes::EVENT_OK;
0321
0322   float ptMaxCut = 4;
0323
0324   float pt1ClusCut = 1.3;  // 1.3
0325   float pt2ClusCut = 0.7;  // 0.7
0326
0327   if (nClusCount > 30)
0328   {
0329     pt1ClusCut += 1.4 * (nClusCount - 29) / 200.0;
0330     pt2ClusCut += 1.4 * (nClusCount - 29) / 200.0;
0331   }
0332
0333   float pi0ptcut = 1.22 * (pt1ClusCut + pt2ClusCut);
0334
0335   vector<float> save_pt;
0336   vector<float> save_eta;
0337   vector<float> save_phi;
0338   vector<float> save_e;
0339
0340   float smear = 0.00;
0341
0342   for (clusterIter = clusterEnd.first; clusterIter != clusterEnd.second; clusterIter++)
0343   {
0344     RawCluster* recoCluster = clusterIter->second;
0345
0346     CLHEP::Hep3Vector vertex(0, 0, vtx_z);
0347     CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*recoCluster, vertex);
0348
0349     float clusE = E_vec_cluster.mag();
0350     float clus_eta = E_vec_cluster.pseudoRapidity();
0351     float clus_phi = E_vec_cluster.phi();
0352     float clus_pt = E_vec_cluster.perp();
0353     clus_pt *= rnd->Gaus(1,smear);
0354     float clus_chisq = recoCluster->get_chi2();
0355
0356     h_clusE->Fill(clusE);
0357     if (clusE < 0.2) continue;
0358
0359     if (clus_chisq > clus_chisq_cut) continue;
0360
0361     // loop over the towers in the cluster
0362     RawCluster::TowerConstRange towerCR = recoCluster->get_towers();
0363     RawCluster::TowerConstIterator toweriter;
0364     float lt_e = -1000;
0365     unsigned int lt_eta = -1;
0366     for (toweriter = towerCR.first; toweriter != towerCR.second; ++toweriter)
0367     {
0368       int towereta = m_geometry->get_tower_geometry(toweriter->first)->get_bineta();
0369       int towerphi = m_geometry->get_tower_geometry(toweriter->first)->get_binphi();
0370       unsigned int key = TowerInfoDefs::encode_emcal(towereta, towerphi);
0371       unsigned int channel = towers->decode_key(key);
0372       float energy = towers->get_tower_at_channel(channel)->get_energy();
0373       if (energy > lt_e)
0374       {
0375         lt_e = energy;
0376         lt_eta = towereta;
0377       }
0378
0379     }
0380
0381     if (lt_eta > 95) continue;
0382
0383     h_pt_eta[lt_eta]->Fill(clus_pt);
0384
0385    // float weight = getWeight(lt_eta,clus_pt);
0386     h_pt_eta_rw[lt_eta]->Fill(clus_pt,weight);
0387
0388
0389     h_etaphi_clus->Fill(clus_eta, clus_phi);
0390
0391     TLorentzVector photon1;
0392     photon1.SetPtEtaPhiE(clus_pt, clus_eta, clus_phi, clusE);
0393
0394     for (auto tr_phot : truth_photons){
0395       float delR = photon1.DeltaR(tr_phot);
0396       h_delR_recTrth->Fill(delR);
0397       if (debug) std::cout << delR << " ";
0398       if (delR < 0.015){
0399         float res = photon1.E()/tr_phot.E();
0400         h_matched_res->Fill(res,photon1.Eta());
0401       }
0402     }
0403     if (debug) std::cout << std::endl;
0404
0405     if (clus_pt < pt1ClusCut || clus_pt > ptMaxCut) continue;
0406
0407     for (clusterIter2 = clusterEnd.first; clusterIter2 != clusterEnd.second; clusterIter2++)
0408     {
0409       if (clusterIter == clusterIter2)
0410       {
0411         continue;
0412       }
0413       RawCluster* recoCluster2 = clusterIter2->second;
0414
0415       CLHEP::Hep3Vector E_vec_cluster2 = RawClusterUtility::GetECoreVec(*recoCluster2, vertex);
0416
0417       float clus2E = E_vec_cluster2.mag();
0418       float clus2_eta = E_vec_cluster2.pseudoRapidity();
0419       float clus2_phi = E_vec_cluster2.phi();
0420       float clus2_pt = E_vec_cluster2.perp();
0421       //clus2_pt *= rnd->Gaus(1,smear);
0422       float clus2_chisq = recoCluster2->get_chi2();
0423
0424       if (clus2_pt < pt2ClusCut || clus2_pt > ptMaxCut) continue;
0425       if (clus2_chisq > clus_chisq_cut) continue;
0426
0427
0428       TLorentzVector photon2;
0429       photon2.SetPtEtaPhiE(clus2_pt, clus2_eta, clus2_phi, clus2E);
0430
0431       if (fabs(clusE - clus2E) / (clusE + clus2E) > maxAlpha) continue;
0432
0433       if (photon1.DeltaR(photon2) > maxDr) continue;
0434
0435       TLorentzVector pi0 = photon1 + photon2;
0436       if (pi0.Pt() < pi0ptcut) continue;
0437
0438
0439       h_pt1->Fill(photon1.Pt());
0440       h_pt2->Fill(photon2.Pt());
0441       h_InvMass->Fill(pi0.M());
0442       h_InvMass_w->Fill(pi0.M(),weight);
0443       h_mass_eta_lt[lt_eta]->Fill(pi0.M());
0444       h_mass_eta_lt_rw[lt_eta]->Fill(pi0.M());
0445     }
0446   }  // clus1 loop
0447
0448
0449   ht_phi.clear();
0450   ht_eta.clear();
0451
0452   return Fun4AllReturnCodes::EVENT_OK;
0453 }
0454
0455
0456 int CaloAna::End(PHCompositeNode* /*topNode*/)
0457 {
0458   outfile->cd();
0459
0460   outfile->Write();
0461   outfile->Close();
0462   delete outfile;
0463   hm->dumpHistos(outfilename, "UPDATE");
0464   return 0;
0465 }
0466
0467 float CaloAna::getWeight(int ieta, float pt){
0468   if (ieta < 0 || ieta > 95) return 0;
0469   float val = h_pt_rw[ieta]->GetBinContent(h_pt_rw[ieta]->FindBin(pt));
0470   if (val==0) return 0;
0471   return 1/val;
0472 }
0473
0474
0475 TF1* CaloAna::fitHistogram(TH1* h)
0476 {
0477   TF1* fitFunc = new TF1("fitFunc", "[0]/[2]/2.5*exp(-0.5*((x-[1])/[2])^2) + [3] + [4]*x + [5]*x^2 + [6]*x^3", h->GetXaxis()->GetXmin(), h->GetXaxis()->GetXmax());
0478
0479   fitFunc->SetParameter(0, 5);
0480   fitFunc->SetParameter(1, target_pi0_mass);
0481   fitFunc->SetParameter(2, 0.01);
0482   fitFunc->SetParameter(3, 0.0);
0483   fitFunc->SetParameter(4, 0.0);
0484   fitFunc->SetParameter(5, 0.0);
0485   fitFunc->SetParameter(6, 100);
0486
0487   fitFunc->SetParLimits(0, 0,10);
0488   fitFunc->SetParLimits(1, 0.113, 0.25);
0489   fitFunc->SetParLimits(2, 0.01, 0.04);
0490   fitFunc->SetParLimits(3,-2 ,1 );
0491   fitFunc->SetParLimits(4,0 ,40 );
0492   fitFunc->SetParLimits(5, -150,50 );
0493   fitFunc->SetParLimits(6, 0,200 );
0494
0495   fitFunc->SetRange(0.05, 0.7);
0496
0497   // Perform the fit
0498   h->Fit("fitFunc", "QN");
0499
0500   return fitFunc;
0501 }
0502
0503
0504 void CaloAna::fitEtaSlices(const std::string& infile, const std::string& fitOutFile, const std::string& cdbFile)
0505 {
0506   TFile* fin = new TFile(infile.c_str());
0507   std::cout << "getting hists" << std::endl;
0508   TH1F* h_peak_eta = new TH1F("h_peak_eta", "", 96, 0, 96);
0509   TH1F* h_sigma_eta = new TH1F("h_sigma_eta", "", 96, 0, 96);
0510   TH1F* h_p3_eta = new TH1F("h_p3_eta", "", 96, 0, 96);
0511   TH1F* h_p4_eta = new TH1F("h_p4_eta", "", 96, 0, 96);
0512   TH1F* h_p5_eta = new TH1F("h_p5_eta", "", 96, 0, 96);
0513   TH1F* h_p6_eta = new TH1F("h_p6_eta", "", 96, 0, 96);
0514   TH1F* h_p0_eta = new TH1F("h_p0_eta", "", 96, 0, 96);
0515   if (!fin)
0516   {
0517     std::cout << "pi0EtaByEta::fitEtaSlices null fin" << std::endl;
0518     exit(1);
0519   }
0520   TH1F* h_M_eta[96];
0521   for (int i = 0; i < 96; i++)
0522   {
0523     h_M_eta[i] = (TH1F*) fin->Get(Form("h_mass_eta_lt_rw%d", i));
0524     h_M_eta[i]->Scale(1./h_M_eta[i]->Integral(),"width");
0525   }
0526
0527   TF1* fitFunOut[96];
0528   for (int i = 0; i < 96; i++)
0529   {
0530     if (!h_M_eta[i])
0531     {
0532       std::cout << "pi0EtaByEta::fitEtaSlices null hist" << std::endl;
0533     }
0534
0535     fitFunOut[i] = fitHistogram(h_M_eta[i]);
0536     fitFunOut[i]->SetName(Form("f_pi0_eta%d",i));
0537     float mass_val_out = fitFunOut[i]->GetParameter(1);
0538     float mass_err_out = fitFunOut[i]->GetParError(1);
0539     h_peak_eta->SetBinContent(i + 1, mass_val_out);
0540     if (isnan(h_M_eta[i]->GetEntries())){
0541        h_peak_eta->SetBinError(i + 1, 0);
0542        continue;
0543     }
0544     h_peak_eta->SetBinError(i + 1, mass_err_out);
0545     h_sigma_eta->SetBinContent(i+1,fitFunOut[i]->GetParameter(2));
0546     h_sigma_eta->SetBinError(i+1,fitFunOut[i]->GetParError(2));
0547     h_p3_eta->SetBinContent(i+1,fitFunOut[i]->GetParameter(3));
0548     h_p3_eta->SetBinError(i+1,fitFunOut[i]->GetParError(3));
0549     h_p4_eta->SetBinContent(i+1,fitFunOut[i]->GetParameter(4));
0550     h_p4_eta->SetBinError(i+1,fitFunOut[i]->GetParError(4));
0551     h_p5_eta->SetBinContent(i+1,fitFunOut[i]->GetParameter(5));
0552     h_p5_eta->SetBinError(i+1,fitFunOut[i]->GetParError(5));
0553     h_p6_eta->SetBinContent(i+1,fitFunOut[i]->GetParameter(6));
0554     h_p6_eta->SetBinError(i+1,fitFunOut[i]->GetParError(6));
0555     h_p0_eta->SetBinContent(i+1,fitFunOut[i]->GetParameter(0));
0556     h_p0_eta->SetBinError(i+1,fitFunOut[i]->GetParError(0));
0557   }
0558
0559   CDBTTree* cdbttree1 = new CDBTTree(cdbFile.c_str());
0560   CDBTTree* cdbttree2 = new CDBTTree(cdbFile.c_str());
0561
0562   std::string m_fieldname = "Femc_datadriven_qm1_correction";
0563
0564   for (int i = 0; i < 96; i++)
0565   {
0566     for (int j = 0; j < 256; j++)
0567     {
0568       float correction = target_pi0_mass / h_peak_eta->GetBinContent(i + 1);
0569       unsigned int key = TowerInfoDefs::encode_emcal(i, j);
0570       float val1 = cdbttree1->GetFloatValue(key, m_fieldname);
0571       cdbttree2->SetFloatValue(key, m_fieldname, val1 * correction);
0572     }
0573   }
0574
0575   cdbttree2->Commit();
0576   cdbttree2->WriteCDBTTree();
0577   delete cdbttree2;
0578   delete cdbttree1;
0579
0580   TFile* fit_out = new TFile(fitOutFile.c_str(), "recreate");
0581   fit_out->cd();
0582   for (auto& i : h_M_eta)
0583   {
0584     i->Write();
0585     delete i;
0586   }
0587   for (auto& i : fitFunOut)
0588   {
0589     i->Write();
0590     delete i;
0591   }
0592
0593   h_p3_eta->Write();
0594   h_p4_eta->Write();
0595   h_p5_eta->Write();
0596   h_p6_eta->Write();
0597   h_p0_eta->Write();
0598   h_sigma_eta->Write();
0599   h_peak_eta->Write();
0600   fin->Close();
0601
0602   std::cout << "finish fitting suc" << std::endl;
0603
0604   return;
0605 }