Commit b5cf91c7 authored by vpalladi's avatar vpalladi

Merge branch 'master' of gitlab.doc.ic.ac.uk:vpalladi/HgcTpgSim

parents ea5b6cab 88df6798
......@@ -22,17 +22,28 @@ public:
fEventSelect=0;
fPrintLevel=0;
for(unsigned e(0);e<=Geometry::kNumberOfEndcaps;e++){
if(e<Geometry::kNumberOfEndcaps) {
std::ostringstream sl,st;
sl << "Endcap" << std::setw(1) << std::setfill('0') << e;
st << "Endcap " << e;
sEndcapLabel[e]=sl.str();
sEndcapTitle[e]=st.str();
}else{
sEndcapLabel[e]="AllEndcaps";
sEndcapTitle[e]="AllEndcaps";
}
}
for(unsigned l(0);l<=Geometry::kNumberOfLayers;l++) {
if(l<Geometry::kNumberOfLayers) {
std::ostringstream sl,st;
sl << "Layer" << std::setw(2) << std::setfill('0') << l;
st << "Layer " << l;
sLayerLabel[l]=sl.str();
sLayerTitle[l]=st.str();
std::ostringstream sl,st;
sl << "Layer" << std::setw(2) << std::setfill('0') << l;
st << "Layer " << l;
sLayerLabel[l]=sl.str();
sLayerTitle[l]=st.str();
} else {
sLayerLabel[l]="LayerAll";
sLayerTitle[l]="All layers";
sLayerLabel[l]="LayerAll";
sLayerTitle[l]="AllLayers";
}
}
}
......@@ -105,6 +116,8 @@ protected:
std::string sLayerLabel[Geometry::kNumberOfLayers+1];
std::string sLayerTitle[Geometry::kNumberOfLayers+1];
std::string sEndcapLabel[Geometry::kNumberOfEndcaps+1];
std::string sEndcapTitle[Geometry::kNumberOfEndcaps+1];
unsigned fEventNumber;
unsigned fEventSelect;
......
......@@ -50,17 +50,26 @@ public:
h_N = new TH1F("h_N","N", 53, 0, 53);
h_CluResponse = new TH1F("h_CluResponse","h_CluResponse",100,0,2);
h_CluResponse2D = new TH2F("h_CluResponse2D","h_CluResponse2D",600,0,300,600,0,300);
h_CluEtaResponse = new TH1F("h_CluEtaResponse", "h_CluEtaResponse", 200, -0.5, 0.5);
h_CluPhiResponse = new TH1F("h_CluPhiResponse", "h_CluPhiResponse", 200, -0.5, 0.5);
h_CluEtaResponseVsEtaTrue = new TH2F("h_CluEtaResponseVsEtaTrue","h_CluEtaResponseVsEtaTrue",200,1.4,3.2, 200, -0.5, 0.5);
h_CluPhiResponseVsPhiTrue = new TH2F("h_CluPhiResponseVsPhiTrue","h_CluPhiResponseVsPhiTrue",200,-3.5,3.5, 200, -0.5, 0.5);
// h_CluEtaResponse = new TH1F("h_CluEtaResponse", "h_CluEtaResponse", 200, -0.5, 0.5);
// h_CluPhiResponse = new TH1F("h_CluPhiResponse", "h_CluPhiResponse", 200, -0.5, 0.5);
// h_CluEtaResponseVsEtaTrue = new TH2F("h_CluEtaResponseVsEtaTrue","h_CluEtaResponseVsEtaTrue",200,1.4,3.2, 200, -0.5, 0.5);
// h_CluPhiResponseVsPhiTrue = new TH2F("h_CluPhiResponseVsPhiTrue","h_CluPhiResponseVsPhiTrue",200,-3.5,3.5, 200, -0.5, 0.5);
h_CluEtaResponseVsEtaC2d = new TH2F("h_CluEtaResponseVsEtaC2d","h_CluEtaResponseVsEtaC2d",200,1.4,3.2, 200, -0.5, 0.5);
h_CluPhiResponseVsPhiC2d = new TH2F("h_CluPhiResponseVsPhiC2d","h_CluPhiResponseVsPhiC2d",200,-3.5,3.5, 200, -0.5, 0.5);
h_CluEtaResponseVsEtTrue = new TH2F("h_CluEtaResponseVsEtTrue","h_CluEtaResponseVsEtTrue",150,0,150, 200, -0.5, 0.5);
h_CluPhiResponseVsEtTrue = new TH2F("h_CluPhiResponseVsEtTrue","h_CluPhiResponseVsEtTrue",150,0,150, 200, -0.5, 0.5);
h_CluEtaC2dVsEtaTrue = new TH2F("h_CluEtaC2dVsEtaTrue","h_CluEtaC2dVsEtaTrue", 200, 1.4, 3.2, 200, 1.4, 3.2);
h_CluPhiC2dVsPhiTrue = new TH2F("h_CluPhiC2dVsPhiTrue","h_CluPhiC2dVsPhiTrue", 200, -3.5, 3.5, 200, -3.5, 3.5);
for(unsigned e(0); e<=Geometry::kNumberOfEndcaps; e++){
std::cout << ("endcap name "+sEndcapLabel[e]).c_str() << std::endl;
for(unsigned l(0); l<=Geometry::kNumberOfLayers; l++){
h_CluEtaResponse[e][l] = new TH1F(("h_CluEtaResponse_"+sEndcapLabel[e]+"_"+sLayerLabel[l]).c_str(), ("h_CluEtaResponse_"+sEndcapTitle[e]+"_"+sLayerTitle[l]).c_str(), 200, -0.5, 0.5);
h_CluPhiResponse[e][l] = new TH1F(("h_CluPhiResponse_"+sEndcapLabel[e]+"_"+sLayerLabel[l]).c_str(), ("h_CluPhiResponse_"+sEndcapTitle[e]+"_"+sLayerTitle[l]).c_str(), 200, -0.5, 0.5);
h_CluEtaResponseVsEtaTrue[e][l] = new TH2F(("h_CluEtaResponseVsEtaTrue_"+sEndcapLabel[e]+"_"+sLayerLabel[l]).c_str(),("h_CluEtaResponseVsEtaTrue_"+sEndcapTitle[e]+"_"+sLayerTitle[l]).c_str(),200,1.4,3.2, 200, -0.5, 0.5);
h_CluPhiResponseVsPhiTrue[e][l] = new TH2F(("h_CluPhiResponseVsPhiTrue_"+sEndcapLabel[e]+"_"+sLayerLabel[l]).c_str(),("h_CluPhiResponseVsPhiTrue_"+sEndcapTitle[e]+"_"+sLayerTitle[l]).c_str(),200,-3.5,3.5, 200, -0.5, 0.5);
}
}
h_CluXResponse = new TH1F("h_CluXResponse", "h_CluXResponse", 200, -10, 10);
h_CluYResponse = new TH1F("h_CluYResponse", "h_CluYResponse", 200, -10, 10);
......@@ -69,11 +78,14 @@ public:
h_MipTmaxC2d_SimMipTmaxC2d = new TH1F("h_MipTmaxC2d_SimMipTmaxC2d","h_MipTmaxC2d_SimMipTmaxC2d",100,0,2);
h_dAB = new TH2F("h_dAB","h_dAB",600, 300,360,100,0,25);
h_dAB_vs_SigEt = new TH2F("h_dAB_vs_SigEt","h_dAB_vs_SigEt",500,0,500,100,0,25);
h_RatoMipTmaxC2d_MipTsecmaxC2d = new TH1F("h_RatoMipTmaxC2d_MipTsecmaxC2d","h_RatoMipTmaxC2d_MipTsecmaxC2d",20,0,10);
h_dmaxC2dSig_vs_dsecmaxC2dSig = new TH2F("h_dmaxC2dSig_vs_dsecmaxC2dSig","h_dmaxC2dSig_vs_dsecmaxC2dSig",100,0,50,100,0,50);
h_C2dMipTdensity_vs_C2dMipT = new TH2F("h_C2dMipTdensity_vs_C2dMipT","h_C2dMipTdensity_vs_C2dMipT",100,0,200,100,0,200);
h_NTrgHitPerC2d = new TH1F("h_NTrgHitPerC2d", "h_NTrgHitPerC2d", 50, 0, 50);
h_NTrgHitPerLayer = new TH1F("h_NTrgHitPerLayer", "h_NTrgHitPerLayer", 75, 0, 75);
h_NC2dPerLayer = new TH1F("h_NC2dPerLayer", "h_NC2dPerLayer", 15, 0, 15);
h_PhiTrueFromXY = new TH1F("h_PhiTrueFromXY","h_PhiTrueFromXY",100,-3.5,7);
h_PhiTrueFrom4P = new TH1F("h_PhiTrueFrom4P","h_PhiTrueFrom4P",100,-3.5,7);
//for(unsigned l(0);l<=Geometry::kNumberOfLayers;l++) {
// h_ELayer_vs_R[l]=new TH2F((fName+sLayerLabel[l]+"h_ELayer_vs_R").c_str(),
// (sLayerTitle[l]+";MIP vs Resolution").c_str(),
......@@ -81,6 +93,7 @@ public:
//}
double deltaEta(double eta1, double eta2);
double deltaPhi(double phi1, double phi2);
double dist2D(double x1, double y1, double x2, double y2);
}
......@@ -113,6 +126,9 @@ public:
}else if(particleType=="photon"){
std::cout << "Using Photons! \n";
vSig=event.photons(e);
}else if(particleType=="electron"){
std::cout << "Using Electrons! \n";
vSig=event.electrons(e);
}
if(debug)std::cout << "e = " << e << " signal size " << vSig.size() <<std::endl;
......@@ -137,7 +153,7 @@ public:
vSigProjTest.push_back(vSig[i].project(Geometry::layerZ(e,l)));
}
for(unsigned i(0);i<vSigProjTest.size();i++){
if( fabs( vSigProjTest[i].FourP().Eta() ) > eta_min && fabs( vSigProjTest[i].FourP().Eta() ) < eta_max ) goodLayer = goodLayer + 1 ;
if( fabs( vSigProjTest[i].position().eta() ) > eta_min && fabs( vSigProjTest[i].position().eta() ) < eta_max ) goodLayer = goodLayer + 1 ;
}
}
......@@ -155,13 +171,13 @@ public:
double SimMIPt_C2d_All = 0.;
double SimMIPt_TOT_layer = 0.;
pair<float,float> CenterEtaPhi(0., 0.);
float Eta_C2d_Max = 0.;
float Phi_C2d_Max = 0.;
double Eta_C2d_Max = 0.;
double Phi_C2d_Max = 0.;
pair<float,float> Center(0.,0.);
float X_C2d_Max = 0.;
float Y_C2d_Max = 0.;
float X_C2d_SecMax = 0.;
float Y_C2d_SecMax = 0.;
double X_C2d_Max = 0.;
double Y_C2d_Max = 0.;
double X_C2d_SecMax = 0.;
double Y_C2d_SecMax = 0.;
int NtrgHit_layer=0;
// Find projection of pion position into layer, and fill a vector at each layer
......@@ -177,12 +193,13 @@ public:
for(unsigned i(0);i<vSigProj.size();i++){
if(debug)std::cout << " ----PiProj Pt "<< vSigProj[i].FourP().Eta() << " rho = " << vSigProj[i].position().rho() << std::endl;
if( !(vSigProj.size()>0 && vSigProj[i].FourP().Pt()>0 && fabs( vSigProj[i].FourP().Eta() ) > 1.6 && fabs( vSigProj[i].FourP().Eta() ) < 2.8 ) ) continue;
if( !(vSigProj.size()>0 && vSigProj[i].FourP().Pt()>0 && fabs( vSigProj[i].position().eta() ) > 1.6 && fabs( vSigProj[i].position().eta() ) < 2.8 ) ) continue;
useful=true;
/**** Getting Clusters ****/
std::vector<TrgC2d> &vTrgC2d( event.trgC2ds(e,l) );
std::sort( vTrgC2d.begin(), vTrgC2d.end(), eneOrderC2d );
//std::sort( vTrgC2d.begin(), vTrgC2d.end(), eneOrderC2d );
std::sort( vTrgC2d.begin(), vTrgC2d.end(), eneTranOrderC2d );
std::cout << "Number of reconstructed clusters = "<< vTrgC2d.size() << std::endl;
h_NC2dPerLayer->Fill(vTrgC2d.size());
......@@ -200,9 +217,12 @@ public:
Center = vTrgC2d[0].center();
X_C2d_Max = Center.first;
Y_C2d_Max = Center.second;
h_C2dMipTdensity_vs_C2dMipT->Fill(vTrgC2d[0].transverseMips(), vTrgC2d[0].transverseMips()/vTrgC2d[0].numberOfTrgHits());
if(vTrgC2d.size()>1){
h_RatoMipTmaxC2d_MipTsecmaxC2d->Fill(vTrgC2d[0].transverseMips()/vTrgC2d[1].transverseMips());
X_C2d_SecMax = (vTrgC2d[1].center()).first;
Y_C2d_SecMax = (vTrgC2d[1].center()).second;
h_dmaxC2dSig_vs_dsecmaxC2dSig->Fill( dist2D(X_C2d_Max, Y_C2d_Max, vSigProj[i].position().x(), vSigProj[i].position().y()), dist2D(X_C2d_SecMax, Y_C2d_SecMax, vSigProj[i].position().x(), vSigProj[i].position().y()) );
}
h_NTrgHitPerC2d->Fill(vTrgC2d[c].numberOfTrgHits());
NtrgHit_layer += vTrgC2d[c].numberOfTrgHits();
......@@ -222,23 +242,47 @@ public:
//it works because there is one pion per endcap!!
t_weightDep_T[l] = MIPt_C2d_Max * vSigProj[i].FourP().Et();
h_MipTmaxC2d_SimMipTmaxC2d->Fill( MIPt_C2d_Max / SimMIPt_C2d_Max );
std::cout << " Pion position --- eta,phi =" << vSigProj[i].FourP().Eta() << ", "<< vSigProj[i].FourP().Phi() << " (x,y) = "<< vSigProj[i].position().x() << ", " << vSigProj[i].position().y() << std::endl;
std::cout << " position eta,phi = " << Eta_C2d_Max << ", "<< Phi_C2d_Max << std::endl;
h_CluEtaResponse->Fill( deltaEta( Eta_C2d_Max, vSigProj[i].FourP().Eta() ) );
h_CluPhiResponse->Fill( deltaPhi( Phi_C2d_Max, vSigProj[i].FourP().Phi() ) );
h_CluEtaResponseVsEtaTrue ->Fill( fabs(vSigProj[i].FourP().Eta()), deltaPhi(Eta_C2d_Max, vSigProj[i].FourP().Eta() ) / vSigProj[i].FourP().Eta());
h_CluPhiResponseVsPhiTrue ->Fill( vSigProj[i].FourP().Phi(), deltaPhi(Phi_C2d_Max, vSigProj[i].FourP().Phi() ) / vSigProj[i].FourP().Phi() );
h_CluEtaResponseVsEtaC2d ->Fill( fabs(Eta_C2d_Max), deltaPhi(Eta_C2d_Max, vSigProj[i].FourP().Eta() ) / vSigProj[i].FourP().Eta() );
h_CluPhiResponseVsPhiC2d ->Fill( Phi_C2d_Max, deltaPhi(Phi_C2d_Max, vSigProj[i].FourP().Phi() ) / vSigProj[i].FourP().Phi() );
if(debug) std::cout << " Pion position --- eta,phi =" << vSigProj[i].FourP().Eta() << ", "<< vSigProj[i].position().phi() << " atan2(y,x) = " << TMath::ATan2(vSigProj[i].position().y(), vSigProj[i].position().x() )<<" (x,y) = "<< vSigProj[i].position().x() << ", " << vSigProj[i].position().y() << std::endl;
if(debug) std::cout << " position eta,phi = " << Eta_C2d_Max << ", "<< Phi_C2d_Max << std::endl;
double PtTrue = TMath::Sqrt(vSigProj[i].momentum().x()*vSigProj[i].momentum().x() + vSigProj[i].momentum().y()*vSigProj[i].momentum().y());
double PhiTrue= vSigProj[i].position().phi() ;
double EtaTrue= vSigProj[i].position().eta() ;
h_PhiTrueFromXY->Fill(PhiTrue);
h_PhiTrueFrom4P->Fill(vSigProj[i].position().phi());
//h_CluEtaResponse->Fill( deltaEta( Eta_C2d_Max, vSigProj[i].FourP().Eta() ) );
//h_CluPhiResponse->Fill( deltaPhi( Phi_C2d_Max, vSigProj[i].position().phi() ) );
//h_CluEtaResponseVsEtaTrue ->Fill( fabs(vSigProj[i].position().eta()), deltaPhi(Eta_C2d_Max, vSigProj[i].position().eta() ) / vSigProj[i].position().eta());
//h_CluPhiResponseVsPhiTrue ->Fill( vSigProj[i].position().phi(), deltaPhi(Phi_C2d_Max, vSigProj[i].position().phi() ) / vSigProj[i].position().phi() );
h_CluEtaResponse[e][l]->Fill( deltaEta( Eta_C2d_Max, vSigProj[i].position().eta() ) );
h_CluPhiResponse[e][l]->Fill( deltaPhi( Phi_C2d_Max, vSigProj[i].position().phi() ) );
h_CluEtaResponseVsEtaTrue[e][l]->Fill( fabs(vSigProj[i].position().eta()), deltaEta(Eta_C2d_Max, vSigProj[i].position().eta() ) );
h_CluPhiResponseVsPhiTrue[e][l]->Fill( vSigProj[i].position().phi(), deltaPhi(Phi_C2d_Max, vSigProj[i].position().phi() ) );
h_CluEtaResponse[Geometry::kNumberOfEndcaps][l]->Fill( deltaEta( Eta_C2d_Max, vSigProj[i].position().eta() ) );
h_CluPhiResponse[Geometry::kNumberOfEndcaps][l]->Fill( deltaPhi( Phi_C2d_Max, vSigProj[i].position().phi() ) );
h_CluEtaResponseVsEtaTrue[Geometry::kNumberOfEndcaps][l]->Fill( fabs(vSigProj[i].position().eta()), deltaEta(Eta_C2d_Max, vSigProj[i].position().eta() ) );
h_CluPhiResponseVsPhiTrue[Geometry::kNumberOfEndcaps][l]->Fill( vSigProj[i].position().phi(), deltaPhi(Phi_C2d_Max, vSigProj[i].position().phi() ) );
h_CluEtaResponse[e][Geometry::kNumberOfLayers]->Fill( deltaEta( Eta_C2d_Max, vSigProj[i].position().eta() ) );
h_CluPhiResponse[e][Geometry::kNumberOfLayers]->Fill( deltaPhi( Phi_C2d_Max, vSigProj[i].position().phi() ) );
h_CluEtaResponseVsEtaTrue[e][Geometry::kNumberOfLayers]->Fill( fabs(vSigProj[i].position().eta()), deltaEta(Eta_C2d_Max, vSigProj[i].position().eta() ) );
h_CluPhiResponseVsPhiTrue[e][Geometry::kNumberOfLayers]->Fill( vSigProj[i].position().phi(), deltaPhi(Phi_C2d_Max, vSigProj[i].position().phi() ) );
h_CluEtaResponse[Geometry::kNumberOfEndcaps][Geometry::kNumberOfLayers]->Fill( deltaEta( Eta_C2d_Max, vSigProj[i].position().eta() ) );
h_CluPhiResponse[Geometry::kNumberOfEndcaps][Geometry::kNumberOfLayers]->Fill( deltaPhi( Phi_C2d_Max, PhiTrue ) );
h_CluEtaResponseVsEtaTrue[Geometry::kNumberOfEndcaps][Geometry::kNumberOfLayers]->Fill( fabs(vSigProj[i].position().eta()), deltaEta(Eta_C2d_Max, vSigProj[i].position().eta() ));
h_CluPhiResponseVsPhiTrue[Geometry::kNumberOfEndcaps][Geometry::kNumberOfLayers]->Fill( vSigProj[i].position().phi(), deltaPhi(Phi_C2d_Max, vSigProj[i].position().phi() ) );
h_CluEtaResponseVsEtaC2d ->Fill( fabs(Eta_C2d_Max), deltaPhi(Eta_C2d_Max, vSigProj[i].position().eta() ) / vSigProj[i].position().eta() );
h_CluPhiResponseVsPhiC2d ->Fill( Phi_C2d_Max, deltaPhi(Phi_C2d_Max, vSigProj[i].position().phi() ) / vSigProj[i].position().phi() );
if(particleType=="pion"){
h_CluEtaResponseVsEtTrue ->Fill( vSigProj[i].FourP().Et(), deltaPhi(Eta_C2d_Max, vSigProj[i].FourP().Eta() ) / vSigProj[i].FourP().Eta(), vSig[0].FourP().Theta() );
h_CluPhiResponseVsEtTrue ->Fill( vSigProj[i].FourP().Et(), deltaPhi(Phi_C2d_Max, vSigProj[i].FourP().Phi() ) / vSigProj[i].FourP().Phi(), vSig[0].FourP().Theta() );
}else if(particleType=="photon"){
h_CluEtaResponseVsEtTrue ->Fill( vSigProj[i].FourP().Et(), deltaPhi(Eta_C2d_Max, vSigProj[i].FourP().Eta() ) / vSigProj[i].FourP().Eta() );
h_CluPhiResponseVsEtTrue ->Fill( vSigProj[i].FourP().Et(), deltaPhi(Phi_C2d_Max, vSigProj[i].FourP().Phi() ) / vSigProj[i].FourP().Phi() );
h_CluEtaResponseVsEtTrue ->Fill( vSigProj[i].FourP().Et(), deltaPhi(Eta_C2d_Max, vSigProj[i].position().eta() ) / vSigProj[i].position().eta(), vSig[0].FourP().Theta() );
h_CluPhiResponseVsEtTrue ->Fill( vSigProj[i].FourP().Et(), deltaPhi(Phi_C2d_Max, vSigProj[i].position().phi() ) / vSigProj[i].position().phi(), vSig[0].FourP().Theta() );
}else if(particleType=="photon" || particleType=="electron" ){
h_CluEtaResponseVsEtTrue ->Fill( vSigProj[i].FourP().Et(), deltaPhi(Eta_C2d_Max, vSigProj[i].position().eta() ) / vSigProj[i].position().eta() );
h_CluPhiResponseVsEtTrue ->Fill( vSigProj[i].FourP().Et(), deltaPhi(Phi_C2d_Max, vSigProj[i].position().phi() ) / vSigProj[i].position().phi() );
}
h_CluEtaC2dVsEtaTrue ->Fill( fabs(vSigProj[i].FourP().Eta()), Eta_C2d_Max );
h_CluPhiC2dVsPhiTrue ->Fill( vSigProj[i].FourP().Phi(), Phi_C2d_Max );
h_CluEtaC2dVsEtaTrue ->Fill( fabs(vSigProj[i].position().eta()), Eta_C2d_Max );
h_CluPhiC2dVsPhiTrue ->Fill( vSigProj[i].position().phi(), Phi_C2d_Max );
h_CluXResponse->Fill( X_C2d_Max - vSigProj[i].position().x() ) ;
h_CluYResponse->Fill( Y_C2d_Max - vSigProj[i].position().y() ) ;
if(vTrgC2d.size()>1){
......@@ -338,10 +382,14 @@ private:
TH1F *h_CluResponse;
TH2F *h_CluResponse2D;
TH1F *h_CluEtaResponse;
TH1F *h_CluPhiResponse;
TH2F *h_CluEtaResponseVsEtaTrue;
TH2F *h_CluPhiResponseVsPhiTrue;
// TH1F *h_CluEtaResponse;
// TH1F *h_CluPhiResponse;
// TH2F *h_CluEtaResponseVsEtaTrue;
// TH2F *h_CluPhiResponseVsPhiTrue;
TH1F *h_CluEtaResponse[Geometry::kNumberOfEndcaps+1][Geometry::kNumberOfLayers+1];
TH1F *h_CluPhiResponse[Geometry::kNumberOfEndcaps+1][Geometry::kNumberOfLayers+1];
TH2F *h_CluEtaResponseVsEtaTrue[Geometry::kNumberOfEndcaps+1][Geometry::kNumberOfLayers+1];
TH2F *h_CluPhiResponseVsPhiTrue[Geometry::kNumberOfEndcaps+1][Geometry::kNumberOfLayers+1];
TH2F *h_CluEtaResponseVsEtaC2d;
TH2F *h_CluPhiResponseVsPhiC2d;
TH2F *h_CluEtaResponseVsEtTrue;
......@@ -356,9 +404,14 @@ private:
TH1F *h_MipTmaxC2d_SimMipTmaxC2d;
TH2F *h_dAB;
TH2F *h_dAB_vs_SigEt;
TH1F *h_RatoMipTmaxC2d_MipTsecmaxC2d;
TH2F *h_dmaxC2dSig_vs_dsecmaxC2dSig;
TH2F *h_C2dMipTdensity_vs_C2dMipT;
TH1F *h_NTrgHitPerC2d;
TH1F *h_NTrgHitPerLayer;
TH1F *h_NC2dPerLayer;
TH1F *h_PhiTrueFromXY;
TH1F *h_PhiTrueFrom4P;
double deltaPhi( double phi1, double phi2) {
......@@ -381,6 +434,11 @@ private:
double dPhi = deltaPhi(phi1, phi2);
return sqrt(dEta*dEta+dPhi*dPhi);
}
double dist2D(double x1, double y1, double x2, double y2){
double dx = ( x1 - x2 );
double dy = ( y1 - y2 );
return sqrt( dx*dx + dy*dy );
}
};
......
......@@ -48,6 +48,11 @@ bool eneOrderC2d( TrgC2d A, TrgC2d B)
return A.mips() > B.mips();
}
bool eneTranOrderC2d( TrgC2d A, TrgC2d B)
{
return A.transverseMips() > B.transverseMips();
}
class Event {
public:
......
......@@ -88,7 +88,7 @@ public:
HepMC::FourVector p4(p.momentum());
fMomentum=Point(p4.x(),p4.y(),p4.z());
fFourP=TLorentzVector( p4.px(), p4.py(), p4.pz(), p4.e() );
fFourP=TLorentzVector( p4.x(), p4.y(), p4.z(), p4.e() );
}
......@@ -155,18 +155,22 @@ public:
p.fFourP=fFourP;
} else {
double dPhi(dz*q*fBField/fMomentum.z());
x+=(fMomentum.x()*sin(dPhi)+fMomentum.y()*(1.0-cos(dPhi)))/(q*fBField);
y+=(fMomentum.y()*sin(dPhi)-fMomentum.x()*(1.0-cos(dPhi)))/(q*fBField);
p.fPosition=Point(x,y,z);
double px(fMomentum.x()*cos(dPhi)+fMomentum.y()*sin(dPhi));
double py(fMomentum.y()*cos(dPhi)-fMomentum.x()*sin(dPhi));
p.fMomentum=Point(px,py,fMomentum.z());
p.fFourP=fFourP;
p.fFourP.SetZ(fMomentum.z());
double dPhi(dz*q*fBField/fMomentum.z());
x+=(fMomentum.x()*sin(dPhi)+fMomentum.y()*(1.0-cos(dPhi)))/(q*fBField);
y+=(fMomentum.y()*sin(dPhi)-fMomentum.x()*(1.0-cos(dPhi)))/(q*fBField);
p.fPosition=Point(x,y,z);
double px(fMomentum.x()*cos(dPhi)+fMomentum.y()*sin(dPhi));
double py(fMomentum.y()*cos(dPhi)-fMomentum.x()*sin(dPhi));
p.fMomentum=Point(px,py,fMomentum.z());
//p.fFourP=fFourP;
p.fFourP.SetX(p.fMomentum.x());
p.fFourP.SetY(p.fMomentum.y());
p.fFourP.SetZ(p.fMomentum.z());
}
return p;
......
......@@ -19,8 +19,8 @@ class SubAnalysisPhotonTrgHit {
public:
enum {
nMCut=10,
nRCut=10
nMCut=9,
nRCut=9
};
SubAnalysisPhotonTrgHit(const std::string &sRoot="",
......@@ -140,9 +140,10 @@ public:
assert(vPhoton.size()==1);
for(unsigned j(0);j<nMCut;j++) {
double mCut(5.0+5.0*j);
double mCut(8.0+8.0*j);
for(unsigned k(0);k<nRCut;k++) {
double rCut(2.0+2.0*k);
double rCutSq(rCut*rCut);
TVectorD vd[2]={TVectorD(Coefficients::kNumberOfCoefficients),
TVectorD(Coefficients::kNumberOfCoefficients)};
......@@ -159,11 +160,18 @@ public:
for(unsigned c(0);c<Geometry::numberOfTriggerCells(l,w);c++) {
Point pCell(Geometry::point(vTrgHit[c]));
Point dp(pCell-proj.position());
if(dp.rho()<rCut) {
if(vTrgHit[c].mips()>=mCut) ml+=vTrgHit[c].mips();
if(vTrgHit[c].transverseMips()>=0.25*mCut) mlt+=vTrgHit[c].transverseMips();
}
//Point dp(pCell-proj.position());
//if(dp.rho()<rCut)
double dRho(pCell.deltaRho(proj.position()));
double dRhoPhi(pCell.rhoDeltaPhi(proj.position()));
//std::cout << "dp.rho() = " << dp.rho() << " =? " << sqrt(dRho*dRho+dRhoPhi*dRhoPhi) << std::endl;
if(3.0*3.0*dRho*dRho+dRhoPhi*dRhoPhi<rCutSq)
{
if(vTrgHit[c].mips()>=mCut) ml+=vTrgHit[c].mips();
if(vTrgHit[c].transverseMips()>=0.25*mCut) mlt+=vTrgHit[c].transverseMips();
}
}
}
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment