3 #include <calobase/RawCluster.h>
4 #include <calobase/RawClusterContainer.h>
5 #include <calobase/RawClusterUtility.h>
6 #include <calobase/RawTower.h>
7 #include <calobase/RawTowerContainer.h>
8 #include <calobase/RawTowerGeom.h>
9 #include <calobase/RawTowerGeomContainer.h>
29 #include <TLorentzVector.h>
65 for (
int nj=0; nj < 96; nj++)
68 for (
int nk=0; nk < 256; nk++)
80 cout <<
"LiteCaloEval::Init(PHCompositeNode *topNode) Initializing" << endl;
88 pairInvMassTotal =
new TH1F(
"pairInvMassTotal",
"Pair Mass Histo", 70, 0.0, 0.7);
89 mass_eta =
new TH2F(
"mass_eta",
"2d Pair Mass Histo", 70, 0.0, 0.7, 400,-1.5,1.5);
90 mass_eta_phi =
new TH3F(
"mass_eta_phi",
"3d Pair Mass Histo", 70, 0.0, 0.7, 150,-1.5,1.5, 256, -3.142, 3.142);
95 for (
int i = 0; i<96; i++)
97 for (
int j = 0; j<258; j++)
103 TString hist_name =
"emc_ieta" + i1 +
"_phi"+ j1;
105 cemc_hist_eta_phi[i][j] =
new TH1F(hist_name.Data(),
"Hist_ieta_phi_",70,0.0,0.7);
111 for (
int i = 0; i < 96; i++)
113 gStyle->SetOptFit(1);
116 TString
b =
"eta_" +
a;
118 eta_hist[i] =
new TH1F (b.Data(),
"eta and all phi's",70,0.0,0.7);
125 _eventTree =
new TTree(
"_eventTree",
"An event level info Tree");
147 cout <<
"Beginning of the event " <<
_ievent << endl;
148 cout <<
"====================================" << endl;
154 RawClusterContainer *recal_clusters = findNode::getClass<RawClusterContainer>(topNode,
"CLUSTER_POS_COR_CEMC");
157 string towernode =
"TOWER_CALIB_" +
_caloname;
158 RawTowerContainer* _towers = findNode::getClass<RawTowerContainer>(topNode, towernode.c_str());
161 cerr <<
PHWHERE <<
" ERROR: Can't find " << towernode << endl;
165 string towergeomnode =
"TOWERGEOM_" +
_caloname;
166 RawTowerGeomContainer* towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, towergeomnode.c_str());
169 cout <<
PHWHERE <<
": Could not find node " << towergeomnode << endl;
179 GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode,
"GlobalVertexMap");
182 if (!vertexmap->
empty())
208 for (t_rclusiter = t_rbegin_end.first; t_rclusiter != t_rbegin_end.second; ++t_rclusiter)
210 RawCluster *t_recalcluster = t_rclusiter->second;
213 if (cluse > 0.1) inCs++;
215 if (cluse > 1.0) savCs[iCs++] = t_recalcluster;
226 for(
int jCs=0; jCs<iCs; jCs++)
238 float tt_clus_energy = E_vec_cluster.
mag();
239 if (tt_clus_energy > 1.0 )
242 float tt_clus_pt = E_vec_cluster.
perp();
243 float tt_clus_phi = E_vec_cluster.
getPhi();
252 for(
int kCs=0; kCs<iCs; kCs++)
254 if(jCs==kCs)
continue;
258 float tt2_clus_energy = E_vec_cluster2.
mag();
259 if (tt2_clus_energy > 1.0)
262 alphaCut =
abs(tt_clus_energy - tt2_clus_energy)/(tt_clus_energy + tt_clus_energy);
266 float tt2_clus_pt = E_vec_cluster2.
perp();
267 float tt2_clus_phi = E_vec_cluster2.
getPhi();
269 TLorentzVector pho1, pho2, pi0lv;
271 pho1.SetPtEtaPhiE(tt_clus_pt, tt_clus_eta,tt_clus_phi,tt_clus_energy);
272 pho2.SetPtEtaPhiE(tt2_clus_pt, tt2_clus_eta,tt2_clus_phi,tt2_clus_energy);
274 if (pho1.DeltaR(pho2) > 0.8)
continue;
277 float pairInvMass=pi0lv.M();
283 std::vector<int> toweretas;
284 std::vector<int> towerphis;
285 std::vector<float> towerenergies;
295 for (toweriter = towers.first; toweriter != towers.second; ++toweriter)
304 toweretas.push_back(towereta);
305 towerphis.push_back(towerphi);
306 towerenergies.push_back(towerenergy);
315 int maxTowerIndex = max_element(towerenergies.begin(), towerenergies.end()) - towerenergies.begin();
324 if (tt_clus_energy > 2.5 && tt2_clus_energy > 1.5)
327 mass_eta->Fill(pairInvMass,tt_clus_eta);
mass_eta_phi->Fill(pairInvMass,tt_clus_eta, tt_clus_phi);
372 TFile *
f =
new TFile(_filename);
374 TTree *
t1 = (TTree*)f->Get(
"_eventTree");
388 t1->SetBranchAddress(
"_eventNumber", &_eventNumber);
389 t1->SetBranchAddress(
"_nClusters", &_nClusters);
390 t1->SetBranchAddress(
"_clusterIDs", _clusterIDs);
391 t1->SetBranchAddress(
"_clusterEnergies", _clusterEnergies);
392 t1->SetBranchAddress(
"_clusterPts", _clusterPts);
393 t1->SetBranchAddress(
"_clusterEtas", _clusterEtas);
394 t1->SetBranchAddress(
"_clusterPhis", _clusterPhis);
395 t1->SetBranchAddress(
"_maxTowerEtas", _maxTowerEtas);
396 t1->SetBranchAddress(
"_maxTowerPhis", _maxTowerPhis);
403 TLorentzVector *savClusLV[10000];
405 int nEntries = (
int)t1->GetEntriesFast();
408 if (nevts < 0 || nEntries < nevts)
411 for (
int i=0; i < nevts2; i++)
419 for (
int j=0; j < nClusters; j++)
424 eta = _clusterEtas[j];
425 phi = _clusterPhis[j];
426 E = _clusterEnergies[j];
430 savClusLV[j] =
new TLorentzVector();
431 savClusLV[j]->SetPtEtaPhiE(pt, eta, phi, E);
438 TLorentzVector *pho1, *pho2;
440 for(
int jCs=0; jCs<iCs; jCs++)
442 pho1 = savClusLV[jCs];
445 for(
int kCs=0; kCs<iCs; kCs++)
447 if(jCs==kCs)
continue;
449 pho2 = savClusLV[kCs];
451 TLorentzVector pi0lv;
453 if (pho1->DeltaR(*pho2) > 0.8)
continue;
456 float pairInvMass=pi0lv.M();
463 eta_hist[_maxTowerEtas[jCs]]->Fill(pairInvMass);
477 TFile *parFile =
new TFile(
"parFile",
"RECREATE");
479 TTree *parTree =
new TTree(
"parTree",
"Tree with fit mean saved");
482 float corrval[100000];
484 parTree->Branch(
"ieta", &ieta,
"ieta/I");
485 parTree->Branch(
"iphi", &iphi,
"iphi/I");
486 parTree->Branch(
"corrval", corrval,
"corrval[ieta*1000+iphi]");
490 for (
int i=0; i<96; i++)
492 for (
int j=0; j<=258; j++)
500 corrval[i*1000+j] = fit_result->GetParameter(1);
504 eta_hist[i]->Fit(
"gaus",
"",
"", 0.095, 0.175);
505 fit_result =
eta_hist[i]->GetFunction(
"gaus");
508 corrval[i*1000+j] = fit_result->GetParameter(1);