CaloDataAnaRun24pp/macro/makePhotonPairs.C

0001 #include <sPhenixStyle.C>
0002
0003 void makePhotonPairs(const std::string &fin = "commissioning.root")
0004 {
0005   SetsPhenixStyle();
0006   TFile *f = new TFile(fin.c_str());
0007
0008   TTree *T = (TTree *) f->Get("T");
0009
0010   std::vector<float> *clusterE{0};
0011   std::vector<float> *clusterPhi{0};
0012   std::vector<float> *clusterEta{0};
0013   std::vector<float> *clusterPt{0};
0014   std::vector<float> *clusterChi2{0};
0015   std::vector<float> *clusterNtow{0};
0016   std::vector<float> *clusterTowMaxE{0};
0017   std::vector<float> *clusterECore{0};
0018
0019   T->SetBranchAddress("clusterE", &clusterE);
0020   T->SetBranchAddress("clusterPhi", &clusterPhi);
0021   T->SetBranchAddress("clusterEta", &clusterEta);
0022   T->SetBranchAddress("clusterPt", &clusterPt);
0023   T->SetBranchAddress("clusterChi2", &clusterChi2);
0024   T->SetBranchAddress("clusterNtow", &clusterNtow);
0025   T->SetBranchAddress("clusterTowMaxE", &clusterTowMaxE);
0026   T->SetBranchAddress("clusterECore", &clusterECore);
0027
0028   TH1 *mass = new TH1F("invMass", "invMass", 100, 0, 1);
0029   mass->SetTitle(";M_{#gamma#gamma};");
0030
0031   float maxChi2 = 4;
0032   float minClusE = 1;
0033   float maxAsym = 0.7;
0034
0035   for (int i = 0; i < T->GetEntries(); i++)
0036   {
0037     T->GetEntry(i);
0038
0039     for (int clus1 = 0; clus1 < clusterE->size(); clus1++)
0040     {
0041       if (clusterChi2->at(clus1) > maxChi2) continue;
0042
0043       if (clusterE->at(clus1) < minClusE) continue;
0044
0045       TLorentzVector photon1;
0046       photon1.SetPtEtaPhiE(clusterPt->at(clus1), clusterEta->at(clus1), clusterPhi->at(clus1), clusterE->at(clus1));
0047
0048       for (int clus2 = 0; clus2 < clusterE->size(); clus2++)
0049       {
0050         if (clus2 <= clus1) continue;
0051
0052         if (clusterE->at(clus2) < minClusE) continue;
0053
0054         if (clusterChi2->at(clus2) > maxChi2) continue;
0055
0056         float asym = sqrt(pow(clusterE->at(clus1) - clusterE->at(clus2), 2)) / (clusterE->at(clus1) + clusterE->at(clus2));
0057         if (asym > maxAsym) continue;
0058
0059         TLorentzVector photon2;
0060         photon2.SetPtEtaPhiE(clusterPt->at(clus2), clusterEta->at(clus2), clusterPhi->at(clus2), clusterE->at(clus2));
0061
0062         TLorentzVector meson = photon1 + photon2;
0063
0064         mass->Fill(meson.M());
0065       }
0066     }
0067   }
0068   TCanvas *cMass = new TCanvas("cMass", "cMass");
0069   mass->Draw("ep");
0070 }