MVTXTracklet/src/Tracklet.h

0001 #ifndef TRACKLET_H
0002 #define TRACKLET_H
0003
0004 #include <TTree.h>
0005
0006 #include "GenHadron.h"
0007 #include "Hit.h"
0008 #include "Utilities.h"
0009
0010 using namespace std;
0011
0012 class Tracklet : public TObject
0013 {
0014   public:
0015     Tracklet();
0016     Tracklet(Hit *hit1, Hit *hit2);
0017     ~Tracklet();
0018
0019     void SetHits(Hit *hit1, Hit *hit2);
0020     Hit *Hit1();
0021     Hit *Hit2();
0022     float dEta();
0023     float dPhi();
0024     float dR();
0025     float Eta();
0026     float Phi();
0027     float tklVtxZ();
0028     int LayerComb(); // Layer combination: can only be 12, 23, 13
0029     void SetMatchedGenHardon();
0030     bool IsMatchedGenHadron();
0031     void SetGenHadron(GenHadron *genhad);
0032     GenHadron *MatchedGenHadron();
0033
0034   private:
0035     Hit *_hit1;
0036     Hit *_hit2;
0037     float _deta;
0038     float _dphi;
0039     bool _matched_genhadron;
0040     GenHadron *_genhad;
0041 };
0042
0043 Tracklet::Tracklet()
0044 {
0045     _hit1 = nullptr;
0046     _hit2 = nullptr;
0047     _matched_genhadron = false;
0048     _genhad = nullptr;
0049 }
0050
0051 Tracklet::Tracklet(Hit *hit1, Hit *hit2)
0052 {
0053     _hit1 = hit1;
0054     _hit2 = hit2;
0055     _deta = _hit2->Eta() - _hit1->Eta();
0056     _dphi = deltaPhi(_hit2->Phi(), _hit1->Phi());
0057     _matched_genhadron = false;
0058     _genhad = nullptr;
0059 }
0060
0061 Tracklet::~Tracklet() {}
0062
0063 void Tracklet::SetHits(Hit *hit1, Hit *hit2)
0064 {
0065     _hit1 = hit1;
0066     _hit2 = hit2;
0067     _deta = _hit2->Eta() - _hit1->Eta();
0068     _dphi = deltaPhi(_hit2->Phi(), _hit1->Phi());
0069     _matched_genhadron = false;
0070     _genhad = nullptr;
0071 }
0072
0073 Hit *Tracklet::Hit1() { return (_hit1); }
0074
0075 Hit *Tracklet::Hit2() { return (_hit2); }
0076
0077 float Tracklet::dEta() { return (_deta); }
0078
0079 float Tracklet::dPhi() { return (_dphi); }
0080
0081 float Tracklet::dR() { return (sqrt(_deta * _deta + _dphi * _dphi)); }
0082
0083 float Tracklet::Eta() { return (0.5 * (_hit1->Eta() + _hit2->Eta())); }
0084
0085 float Tracklet::Phi() { return (0.5 * (_hit1->Phi() + _hit2->Phi())); }
0086
0087 float Tracklet::tklVtxZ() { return (_hit1->posZ() - _hit1->rho() * (_hit2->posZ() - _hit1->posZ()) / (_hit2->rho() - _hit1->rho())); }
0088
0089 int Tracklet::LayerComb() { return (_hit1->Layer()) * 10 + (_hit2->Layer()); }
0090
0091 void Tracklet::SetMatchedGenHardon() { _matched_genhadron = true; }
0092
0093 bool Tracklet::IsMatchedGenHadron() { return _matched_genhadron; }
0094
0095 void Tracklet::SetGenHadron(GenHadron *genhad)
0096 {
0097     _matched_genhadron = true;
0098     _genhad = genhad;
0099 }
0100
0101 GenHadron *Tracklet::MatchedGenHadron() { return _genhad; }
0102
0103 /*------------------------------------------------------------------------------------------*/
0104
0105 class TrackletData
0106 {
0107   public:
0108     vector<vector<Hit *>> MVTXlayers = {{}, {}, {}};
0109     vector<Tracklet *> ProtoTkls, RecoTkls, RecoTkls_GenMatched;
0110     vector<GenHadron *> GenHadrons;
0111     int MVTXl1, MVTXl2;
0112
0113     int event, NhitsLayer1, NPrototkl, NRecotkl_Raw, NRecotkl_GenMatched, NGenHadron, CentralityBin;
0114     float PV_z, TruthPV_trig_z, TruthPV_mostNpart_z;
0115     // bool pu0_sel, vz_sel, o_sel, proc_sel, evt_sel, gen_sel;
0116     bool pu0_sel, trig;
0117     int process; // single diffractive process
0118     vector<float> recoclus_eta_l1, recoclus_eta_l2, recoclus_eta_l3, recoclus_phi_l1, recoclus_phi_l2, recoclus_phi_l3;
0119     vector<float> prototkl_eta, prototkl_phi, prototkl_clus2eta, prototkl_clus2phi, prototkl_deta, prototkl_dphi, prototkl_dR;
0120     vector<float> recotklraw_eta, recotklraw_phi, recotklraw_clus2eta, recotklraw_clus2phi, recotklraw_deta, recotklraw_dphi, recotklraw_dR;
0121     vector<float> recotklgm_eta, recotklgm_phi, recotklgm_clus2eta, recotklgm_clus2phi, recotklgm_deta, recotklgm_dphi, recotklgm_dR;
0122     vector<float> GenHadron_Pt, GenHadron_eta, GenHadron_phi, GenHadron_E, GenHadron_matched_Pt, GenHadron_matched_eta, GenHadron_matched_phi, GenHadron_matched_E;
0123 };
0124
0125 void RecoclusToMinitree(TrackletData &tkldata)
0126 {
0127     for (size_t i = 0; i < tkldata.MVTXlayers[0].size(); i++)
0128     {
0129         tkldata.recoclus_eta_l1.push_back(tkldata.MVTXlayers[0][i]->Eta());
0130         tkldata.recoclus_phi_l1.push_back(tkldata.MVTXlayers[0][i]->Phi());
0131     }
0132
0133     for (size_t i = 0; i < tkldata.MVTXlayers[1].size(); i++)
0134     {
0135         tkldata.recoclus_eta_l2.push_back(tkldata.MVTXlayers[1][i]->Eta());
0136         tkldata.recoclus_phi_l2.push_back(tkldata.MVTXlayers[1][i]->Phi());
0137     }
0138
0139     for (size_t i = 0; i < tkldata.MVTXlayers[2].size(); i++)
0140     {
0141         tkldata.recoclus_eta_l3.push_back(tkldata.MVTXlayers[2][i]->Eta());
0142         tkldata.recoclus_phi_l3.push_back(tkldata.MVTXlayers[2][i]->Phi());
0143     }
0144 }
0145
0146 void SetMinitree(TTree *outTree, TrackletData &tkldata)
0147 {
0148     outTree->Branch("event", &tkldata.event);
0149     outTree->Branch("NhitsLayer1", &tkldata.NhitsLayer1);
0150     outTree->Branch("NPrototkl", &tkldata.NPrototkl);
0151     outTree->Branch("NRecotkl_Raw", &tkldata.NRecotkl_Raw);
0152     outTree->Branch("NRecotkl_GenMatched", &tkldata.NRecotkl_GenMatched);
0153     outTree->Branch("NGenHadron", &tkldata.NGenHadron);
0154     outTree->Branch("CentralityBin", &tkldata.CentralityBin);
0155     outTree->Branch("PV_z", &tkldata.PV_z);
0156     outTree->Branch("TruthPV_trig_z", &tkldata.TruthPV_trig_z);
0157     outTree->Branch("TruthPV_mostNpart_z", &tkldata.TruthPV_mostNpart_z);
0158     outTree->Branch("pu0_sel", &tkldata.pu0_sel);
0159     outTree->Branch("trig", &tkldata.trig);
0160     outTree->Branch("process", &tkldata.process);
0161     // outTree->Branch("vz_sel", &tkldata.vz_sel);
0162     // outTree->Branch("o_sel", &tkldata.o_sel);
0163     // outTree->Branch("proc_sel", &tkldata.proc_sel);
0164     // outTree->Branch("evt_sel", &tkldata.evt_sel);
0165     // outTree->Branch("gen_sel", &tkldata.gen_sel);
0166     outTree->Branch("recoclus_eta_l1", &tkldata.recoclus_eta_l1);
0167     outTree->Branch("recoclus_phi_l1", &tkldata.recoclus_phi_l1);
0168     outTree->Branch("recoclus_eta_l2", &tkldata.recoclus_eta_l2);
0169     outTree->Branch("recoclus_phi_l2", &tkldata.recoclus_phi_l2);
0170     outTree->Branch("recoclus_eta_l3", &tkldata.recoclus_eta_l3);
0171     outTree->Branch("recoclus_phi_l3", &tkldata.recoclus_phi_l3);
0172     outTree->Branch("prototkl_eta", &tkldata.prototkl_eta);
0173     outTree->Branch("prototkl_phi", &tkldata.prototkl_phi);
0174     outTree->Branch("prototkl_clus2eta", &tkldata.prototkl_clus2eta);
0175     outTree->Branch("prototkl_clus2phi", &tkldata.prototkl_clus2phi);
0176     outTree->Branch("prototkl_deta", &tkldata.prototkl_deta);
0177     outTree->Branch("prototkl_dphi", &tkldata.prototkl_dphi);
0178     outTree->Branch("prototkl_dR", &tkldata.prototkl_dR);
0179     outTree->Branch("recotklraw_eta", &tkldata.recotklraw_eta);
0180     outTree->Branch("recotklraw_phi", &tkldata.recotklraw_phi);
0181     outTree->Branch("recotklraw_clus2eta", &tkldata.recotklraw_clus2eta);
0182     outTree->Branch("recotklraw_clus2phi", &tkldata.recotklraw_clus2phi);
0183     outTree->Branch("recotklraw_deta", &tkldata.recotklraw_deta);
0184     outTree->Branch("recotklraw_dphi", &tkldata.recotklraw_dphi);
0185     outTree->Branch("recotklraw_dR", &tkldata.recotklraw_dR);
0186     outTree->Branch("recotklgm_eta", &tkldata.recotklgm_eta);
0187     outTree->Branch("recotklgm_phi", &tkldata.recotklgm_phi);
0188     outTree->Branch("recotklgm_clus2eta", &tkldata.recotklgm_clus2eta);
0189     outTree->Branch("recotklgm_clus2phi", &tkldata.recotklgm_clus2phi);
0190     outTree->Branch("recotklgm_deta", &tkldata.recotklgm_deta);
0191     outTree->Branch("recotklgm_dphi", &tkldata.recotklgm_dphi);
0192     outTree->Branch("recotklgm_dR", &tkldata.recotklgm_dR);
0193     outTree->Branch("GenHadron_Pt", &tkldata.GenHadron_Pt);
0194     outTree->Branch("GenHadron_eta", &tkldata.GenHadron_eta);
0195     outTree->Branch("GenHadron_phi", &tkldata.GenHadron_phi);
0196     outTree->Branch("GenHadron_E", &tkldata.GenHadron_E);
0197     outTree->Branch("GenHadron_matched_Pt", &tkldata.GenHadron_matched_Pt);
0198     outTree->Branch("GenHadron_matched_eta", &tkldata.GenHadron_matched_eta);
0199     outTree->Branch("GenHadron_matched_phi", &tkldata.GenHadron_matched_phi);
0200     outTree->Branch("GenHadron_matched_E", &tkldata.GenHadron_matched_E);
0201 }
0202
0203 void ResetVec(TrackletData &tkldata)
0204 {
0205     for (size_t i = 0; i < tkldata.MVTXlayers.size(); i++)
0206     {
0207         CleanVec(tkldata.MVTXlayers[i]);
0208     }
0209     CleanVec(tkldata.recoclus_eta_l1);
0210     CleanVec(tkldata.recoclus_phi_l1);
0211     CleanVec(tkldata.recoclus_eta_l2);
0212     CleanVec(tkldata.recoclus_phi_l2);
0213     CleanVec(tkldata.recoclus_eta_l3);
0214     CleanVec(tkldata.recoclus_phi_l3);
0215     CleanVec(tkldata.ProtoTkls);
0216     CleanVec(tkldata.RecoTkls);
0217     CleanVec(tkldata.RecoTkls_GenMatched);
0218     CleanVec(tkldata.GenHadrons);
0219     CleanVec(tkldata.prototkl_eta);
0220     CleanVec(tkldata.prototkl_phi);
0221     CleanVec(tkldata.prototkl_clus2eta);
0222     CleanVec(tkldata.prototkl_clus2phi);
0223     CleanVec(tkldata.prototkl_deta);
0224     CleanVec(tkldata.prototkl_dphi);
0225     CleanVec(tkldata.prototkl_dR);
0226     CleanVec(tkldata.recotklraw_eta);
0227     CleanVec(tkldata.recotklraw_phi);
0228     CleanVec(tkldata.recotklraw_clus2eta);
0229     CleanVec(tkldata.recotklraw_clus2phi);
0230     CleanVec(tkldata.recotklraw_deta);
0231     CleanVec(tkldata.recotklraw_dphi);
0232     CleanVec(tkldata.recotklraw_dR);
0233     CleanVec(tkldata.recotklgm_eta);
0234     CleanVec(tkldata.recotklgm_phi);
0235     CleanVec(tkldata.recotklgm_clus2eta);
0236     CleanVec(tkldata.recotklgm_clus2phi);
0237     CleanVec(tkldata.recotklgm_deta);
0238     CleanVec(tkldata.recotklgm_dphi);
0239     CleanVec(tkldata.recotklgm_dR);
0240     CleanVec(tkldata.GenHadron_Pt);
0241     CleanVec(tkldata.GenHadron_eta);
0242     CleanVec(tkldata.GenHadron_phi);
0243     CleanVec(tkldata.GenHadron_E);
0244     CleanVec(tkldata.GenHadron_matched_Pt);
0245     CleanVec(tkldata.GenHadron_matched_eta);
0246     CleanVec(tkldata.GenHadron_matched_phi);
0247     CleanVec(tkldata.GenHadron_matched_E);
0248 }
0249
0250 bool compare_dR(Tracklet *a, Tracklet *b) { return a->dR() < b->dR(); }
0251
0252 void ProtoTracklets(TrackletData &tkldata, float cutdr)
0253 {
0254     float Cut_dR = cutdr; // Nominal: 0.5
0255
0256     int iComb = 0;
0257     // Get all combinations
0258     for (auto ihitl1 : tkldata.MVTXlayers[tkldata.MVTXl1])
0259     {
0260         for (auto ihitl2 : tkldata.MVTXlayers[tkldata.MVTXl2])
0261         {
0262             // iComb++;
0263             // if (iComb % 500000 == 0)
0264             //     fprintf(stderr, "processing %d of %d combinations (%.3f %%)\n", iComb, nhits1 * nhits2, (float)iComb / (nhits1 * nhits2) * 100.);
0265
0266             float dEta = ihitl2->Eta() - ihitl1->Eta();
0267             float dPhi = deltaPhi(ihitl2->Phi(), ihitl1->Phi());
0268             float dR = sqrt(dEta * dEta + dPhi * dPhi);
0269
0270             if (dR < Cut_dR)
0271             {
0272                 Tracklet *tmptkl = new Tracklet(ihitl1, ihitl2);
0273                 tkldata.ProtoTkls.push_back(tmptkl);
0274
0275                 tkldata.prototkl_eta.push_back(tmptkl->Hit1()->Eta());
0276                 tkldata.prototkl_phi.push_back(tmptkl->Hit1()->Phi());
0277                 tkldata.prototkl_clus2eta.push_back(tmptkl->Hit2()->Eta());
0278                 tkldata.prototkl_clus2phi.push_back(tmptkl->Hit2()->Phi());
0279                 tkldata.prototkl_deta.push_back(tmptkl->dEta());
0280                 tkldata.prototkl_dphi.push_back(tmptkl->dPhi());
0281                 tkldata.prototkl_dR.push_back(tmptkl->dR());
0282             }
0283             else
0284                 continue;
0285         }
0286     }
0287
0288     tkldata.NPrototkl = tkldata.ProtoTkls.size();
0289 }
0290
0291 void RecoTracklets(TrackletData &tkldata)
0292 {
0293     sort(tkldata.ProtoTkls.begin(), tkldata.ProtoTkls.end(), compare_dR);
0294
0295     int itkl = 0;
0296     for (auto &tkl : tkldata.ProtoTkls)
0297     {
0298         if (tkl->Hit1()->IsMatchedTkl() || tkl->Hit2()->IsMatchedTkl())
0299         {
0300             continue;
0301         }
0302         else
0303         {
0304             tkldata.RecoTkls.push_back(tkl);
0305             tkl->Hit1()->SetMatchedTkl();
0306             tkl->Hit2()->SetMatchedTkl();
0307
0308             tkldata.recotklraw_eta.push_back(tkl->Hit1()->Eta());
0309             tkldata.recotklraw_phi.push_back(tkl->Hit1()->Phi());
0310             tkldata.recotklraw_clus2eta.push_back(tkl->Hit2()->Eta());
0311             tkldata.recotklraw_clus2phi.push_back(tkl->Hit2()->Phi());
0312             tkldata.recotklraw_deta.push_back(tkl->dEta());
0313             tkldata.recotklraw_dphi.push_back(tkl->dPhi());
0314             tkldata.recotklraw_dR.push_back(tkl->dR());
0315         }
0316     }
0317
0318     tkldata.NRecotkl_Raw = tkldata.RecoTkls.size();
0319 }
0320
0321 void GenMatch_Recotkl(TrackletData &tkldata)
0322 {
0323     for (auto &tkl : tkldata.RecoTkls)
0324     {
0325         if (tkl->IsMatchedGenHadron())
0326             continue;
0327         for (auto &ghadron : tkldata.GenHadrons)
0328         {
0329             if (ghadron->IsMatchedToRecotkl() || tkl->IsMatchedGenHadron())
0330                 continue;
0331             // Matching criteria
0332             if (deltaR(tkl->Eta(), tkl->Phi(), ghadron->Eta(), ghadron->Phi()) > 0.05)
0333                 continue;
0334             else
0335             {
0336                 tkl->SetMatchedGenHardon();
0337                 ghadron->SetMatchedToRecotkl();
0338                 tkl->SetGenHadron(ghadron);
0339                 tkldata.RecoTkls_GenMatched.push_back(tkl);
0340                 tkldata.recotklgm_eta.push_back(tkl->Hit1()->Eta());
0341                 tkldata.recotklgm_phi.push_back(tkl->Hit1()->Phi());
0342                 tkldata.recotklgm_clus2eta.push_back(tkl->Hit2()->Eta());
0343                 tkldata.recotklgm_clus2phi.push_back(tkl->Hit2()->Phi());
0344                 tkldata.recotklgm_deta.push_back(tkl->dEta());
0345                 tkldata.recotklgm_dphi.push_back(tkl->dPhi());
0346                 tkldata.recotklgm_dR.push_back(tkl->dR());
0347                 tkldata.GenHadron_matched_Pt.push_back(ghadron->Pt());
0348                 tkldata.GenHadron_matched_eta.push_back(ghadron->Eta());
0349                 tkldata.GenHadron_matched_phi.push_back(ghadron->Phi());
0350                 tkldata.GenHadron_matched_E.push_back(ghadron->E());
0351             }
0352         }
0353     }
0354
0355     tkldata.NRecotkl_GenMatched = tkldata.RecoTkls_GenMatched.size();
0356 }
0357
0358 #endif