analysis_INTT/src/Hit.h

0001 #ifndef HIT_H
0002 #define HIT_H
0003
0004 #include <algorithm>
0005 #include <cmath>
0006 #include <iostream>
0007 #include <stdio.h>
0008
0009 #include <TF1.h>
0010 #include <TH1F.h>
0011 #include <TObject.h>
0012 #include <TRandom3.h>
0013 #include <TVector3.h>
0014
0015 using namespace std;
0016
0017 class Hit : public TObject
0018 {
0019   public:
0020     Hit();
0021     Hit(float, float, float, float, float, float, int);
0022     Hit(float, float, float, float, float, float, int, float);
0023     Hit(float, float, float, float, float, float, int, float, unsigned int);
0024     Hit(float, float); // simple construct
0025     ~Hit();
0026
0027     float posX();
0028     float posY();
0029     float posZ();
0030     float rho();
0031     float vtxX();
0032     float vtxY();
0033     float vtxZ();
0034     float Eta();
0035     float Phi();
0036     float R();
0037     int Layer();
0038     float PhiSize() { return _phisize; }
0039     unsigned int ClusADC() { return _clusadc; }
0040     pair<float, float> Edge();
0041     int MatchedG4P_trackID() { return matchedG4P_trackID; };
0042     int MatchedG4P_ancestor_trackID() { return matchedG4P_ancestor_trackID; };
0043     float MatchedG4P_Pt() { return matchedG4P_Pt; };
0044     float MatchedG4P_Eta() { return matchedG4P_Eta; };
0045     float MatchedG4P_Phi() { return matchedG4P_Phi; };
0046
0047     void Update();
0048     void SetPos(float, float, float);
0049     void SetVtx(float, float, float);
0050     void SetEdge(float, float);
0051     void SetMatchedG4P(int, int, int, int, int);
0052     void SetMatchedTkl();
0053     bool IsMatchedTkl();
0054     void Print();
0055
0056     TVector3 VecPos();
0057     TVector3 VecVtx();
0058     TVector3 VecRel();
0059
0060   private:
0061     float _x;
0062     float _y;
0063     float _z;
0064     float _vtxX;
0065     float _vtxY;
0066     float _vtxZ;
0067     float _eta;
0068     float _phi;
0069     float _R;
0070     int _layer;
0071     float _phisize;
0072     unsigned int _clusadc;
0073     pair<float, float> _edge;
0074     bool _matched_tkl;
0075     TVector3 vechit;
0076     TVector3 vecvtx;
0077     TVector3 vecrel;
0078     int matchedG4P_trackID; // only for simulation (matching)
0079     int matchedG4P_ancestor_trackID;
0080     int matchedG4P_Pt;
0081     int matchedG4P_Eta;
0082     int matchedG4P_Phi;
0083 };
0084
0085 Hit::Hit()
0086 {
0087     vechit.SetXYZ(0, 0, 0);
0088     vecvtx.SetXYZ(0, 0, 0);
0089     _layer = -999;
0090     vecrel = vechit - vecvtx;
0091     _eta = vecrel.Eta();
0092     _phi = vecrel.Phi();
0093     _R = vecrel.Mag();
0094     _matched_tkl = false;
0095 }
0096
0097 Hit::Hit(float x, float y, float z, float vtxX, float vtxY, float vtxZ, int layer)
0098 {
0099     vechit.SetXYZ(x, y, z);
0100     vecvtx.SetXYZ(vtxX, vtxY, vtxZ);
0101     _layer = layer;
0102     vecrel = vechit - vecvtx;
0103     _eta = vecrel.Eta();
0104     _phi = vecrel.Phi();
0105     _R = vecrel.Mag();
0106     _matched_tkl = false;
0107 }
0108
0109 Hit::Hit(float x, float y, float z, float vtxX, float vtxY, float vtxZ, int layer, float phisize)
0110 {
0111     vechit.SetXYZ(x, y, z);
0112     vecvtx.SetXYZ(vtxX, vtxY, vtxZ);
0113     _layer = layer;
0114     _phisize = phisize;
0115     vecrel = vechit - vecvtx;
0116     _eta = vecrel.Eta();
0117     _phi = vecrel.Phi();
0118     _R = vecrel.Mag();
0119     _matched_tkl = false;
0120 }
0121
0122 Hit::Hit(float x, float y, float z, float vtxX, float vtxY, float vtxZ, int layer, float phisize, unsigned int clusadc)
0123 {
0124     vechit.SetXYZ(x, y, z);
0125     vecvtx.SetXYZ(vtxX, vtxY, vtxZ);
0126     _layer = layer;
0127     _phisize = phisize;
0128     _clusadc = clusadc;
0129     vecrel = vechit - vecvtx;
0130     _eta = vecrel.Eta();
0131     _phi = vecrel.Phi();
0132     _R = vecrel.Mag();
0133     _matched_tkl = false;
0134 }
0135
0136 Hit::Hit(float eta, float phi)
0137 {
0138     _eta = eta;
0139     _phi = phi;
0140     _matched_tkl = false;
0141 }
0142
0143 Hit::~Hit() {}
0144
0145 float Hit::posX() { return (vechit.X()); }
0146
0147 float Hit::posY() { return (vechit.Y()); }
0148
0149 float Hit::posZ() { return (vechit.Z()); }
0150
0151 float Hit::rho() { return (sqrt(vechit.X() * vechit.X() + vechit.Y() * vechit.Y())); }
0152
0153 float Hit::vtxX() { return (vecvtx.X()); }
0154
0155 float Hit::vtxY() { return (vecvtx.Y()); }
0156
0157 float Hit::vtxZ() { return (vecvtx.Z()); }
0158
0159 float Hit::Eta() { return (_eta); } // with respect to the vertex that it associates to
0160
0161 float Hit::Phi() { return (_phi); } // with respect to the vertex that it associates to
0162
0163 float Hit::R() { return (_R); } // with respect to the vertex that it associates to
0164
0165 int Hit::Layer() { return (_layer); }
0166
0167 pair<float, float> Hit::Edge() { return (_edge); }
0168
0169 void Hit::SetPos(float x, float y, float z) { vechit.SetXYZ(x, y, z); }
0170
0171 void Hit::SetVtx(float vtxX, float vtxY, float vtxZ) { vecvtx.SetXYZ(vtxX, vtxY, vtxZ); }
0172
0173 void Hit::SetEdge(float edge1, float edge2) { _edge = make_pair(edge1, edge2); }
0174
0175 void Hit::Update()
0176 {
0177     vecrel = vechit - vecvtx;
0178     _eta = vecrel.Eta();
0179     _phi = vecrel.Phi();
0180     _R = vecrel.Mag();
0181 }
0182
0183 void Hit::SetMatchedTkl() { _matched_tkl = true; }
0184
0185 bool Hit::IsMatchedTkl() { return _matched_tkl; }
0186
0187 TVector3 Hit::VecPos() { return (vechit); }
0188
0189 TVector3 Hit::VecVtx() { return (vecvtx); }
0190
0191 TVector3 Hit::VecRel() { return (vecrel); }
0192
0193 void Hit::SetMatchedG4P(int trackID, int ancestor_trackID, int Pt, int Eta, int Phi)
0194 {
0195     matchedG4P_trackID = trackID;
0196     matchedG4P_ancestor_trackID = ancestor_trackID;
0197     matchedG4P_Pt = Pt;
0198     matchedG4P_Eta = Eta;
0199     matchedG4P_Phi = Phi;
0200 }
0201
0202 void Hit::Print()
0203 {
0204     printf("[Hit::Print()] (posX, posY, posZ) = (%f, %f, %f), (vtxX, vtxY, vtxZ) = (%f, %f, %f), (eta, phi) = (%f, %f) \n", vechit.X(), vechit.Y(), vechit.Z(), vecvtx.X(), vecvtx.Y(), vecvtx.Z(),
0205            vecrel.Eta(), vecrel.Phi());
0206 }
0207
0208 void UpdateHits(vector<Hit *> &Hits, vector<float> PV)
0209 {
0210     for (auto &hit : Hits)
0211     {
0212         hit->SetVtx(PV[0], PV[1], PV[2]);
0213         hit->Update();
0214     }
0215 }
0216
0217 int RandomHit(int set)
0218 {
0219     return static_cast<int>(35 * set);
0220 }
0221
0222 Hit *RandomHit(float vx, float vy, float vz, int layer)
0223 {
0224     gRandom->SetSeed(0);
0225     // The 26 unique z positions for INTT
0226     vector<float> zpos = {-22.57245, -20.57245, -18.57245, -16.57245, -14.57245, -12.57245, -10.97245, -9.372450, -7.772450, -6.172450, -4.572450, -2.972450, -1.372450,
0227                           0.4275496, 2.0275495, 3.6275494, 5.2275495, 6.8275494, 8.4275493, 10.027549, 11.627549, 13.627549, 15.627549, 17.627550, 19.627550, 21.627550};
0228     int zpos_idx = gRandom->Integer(26);
0229
0230     float layer_radius[4] = {7.453, 8.046, 9.934, 10.569};
0231     float phiMin = -TMath::Pi();
0232     float phiMax = TMath::Pi();
0233     float phi = phiMin + (phiMax - phiMin) * gRandom->Rndm();
0234     float x = layer_radius[layer] * cos(phi);
0235     float y = layer_radius[layer] * sin(phi);
0236     float z = zpos[zpos_idx];                                 // layer_radius[layer] / tan(2 * atan(exp(-eta)));
0237     Hit *randhit = new Hit(x, y, z, vx, vy, vz, layer, 1, 1); // assign the phisize and clusadc to 1 for random clusters, which should be ok
0238
0239     return randhit;
0240 }
0241
0242 float theta2pseudorapidity(float theta) { return -1. * TMath::Log(TMath::Tan(theta / 2)); }
0243
0244 TF1 *ClusADCCut(int constscale, float etascale)
0245 {
0246     TF1 *f = new TF1("f", Form("%d*TMath::CosH(%f*x)", constscale, etascale), -10, 10);
0247     return f;
0248 }
0249
0250 TH1F *ClusADCCut_StepFunc_INTTPrivate()
0251 {
0252     std::vector<float> thetastep = {0.001, 15, 20, 25, 30, 35, 45, 55, 125, 135, 145, 150, 155, 160, 165, 179.999};
0253     std::reverse(thetastep.begin(), thetastep.end());
0254     std::vector<float> adccut_theta = {225, 165, 135, 120, 105, 90, 75, 60, 75, 90, 105, 120, 135, 165, 225};
0255     float etastep_array[thetastep.size()];
0256     for (int i = 0; i < thetastep.size(); i++)
0257     {
0258         etastep_array[i] = theta2pseudorapidity(thetastep[i] * TMath::Pi() / 180);
0259     }
0260
0261     TH1F *hm_cut_inttprivate = new TH1F("hm_cut_inttprivate", "hm_cut_inttprivate", thetastep.size() - 1, etastep_array);
0262     for (int j = 0; j < hm_cut_inttprivate->GetNbinsX(); j++)
0263     {
0264         hm_cut_inttprivate->SetBinContent(j + 1, adccut_theta[j]);
0265     }
0266
0267     return hm_cut_inttprivate;
0268 }
0269
0270 TH1F *ClusADCCut_StepFunc(int constscale, float etascale)
0271 {
0272     TF1 *f_cut = ClusADCCut(constscale, etascale);
0273
0274     std::vector<float> adcstep;
0275     for (int i = 0; i < 20; i++)
0276     {
0277         adcstep.push_back(20 + i * 30);
0278     }
0279     std::vector<float> etastep;
0280     for (int i = 0; i < 20; i++)
0281     {
0282         etastep.insert(etastep.begin(), f_cut->GetX(adcstep[i], -10, 0));
0283         etastep.push_back(f_cut->GetX(adcstep[i], 0, 10));
0284     }
0285
0286     etastep.erase(std::remove_if(etastep.begin(), etastep.end(), [](float x) { return std::isnan(x); }), etastep.end());
0287
0288     float etastep_array[etastep.size()];
0289     for (int i = 0; i < etastep.size(); i++)
0290     {
0291         etastep_array[i] = etastep[i];
0292     }
0293     TH1F *hm_cut = new TH1F("hm_cut", "hm_cut", etastep.size() - 1, etastep_array);
0294     for (int j = 0; j < hm_cut->GetNbinsX(); j++)
0295     {
0296         if (hm_cut->GetBinLowEdge(j + 1) < 0)
0297         {
0298             hm_cut->SetBinContent(j + 1, f_cut->Eval(hm_cut->GetBinLowEdge(j + 1)));
0299         }
0300         else
0301         {
0302             hm_cut->SetBinContent(j + 1, f_cut->Eval(hm_cut->GetBinCenter(j + 1) + hm_cut->GetBinWidth(j + 1) / 2));
0303         }
0304     }
0305
0306     return hm_cut;
0307 }
0308
0309 int ConstADCCut(int set)
0310 {
0311     float cut = 0;
0312     switch (set)
0313     {
0314     case 0:
0315         cut = 35;
0316         break;
0317     case 1:
0318         cut = 0;
0319         break;
0320     case 2:
0321         cut = 50;
0322         break;
0323     }
0324     return cut;
0325 }
0326
0327 int ConstClusPhiCut(int set)
0328 {
0329     float cut = 0;
0330     switch (set)
0331     {
0332     case 0:
0333         cut = 40;
0334         break;
0335     case 1:
0336         cut = 1E6;
0337         break;
0338     }
0339     return cut;
0340 }
0341
0342 #endif