Commit cd0771c4 authored by Luca Mastrolorenzo's avatar Luca Mastrolorenzo

problem in phi-response solved in the CoefficientDetClu.hh

parent afa80453
......@@ -22,10 +22,10 @@ public:
fEventSelect=0;
fPrintLevel=0;
for(unsigned e(0);e<Geometry::kNumberOfEndcaps;e++){
for(unsigned e(0);e<=Geometry::kNumberOfEndcaps;e++){
if(e<Geometry::kNumberOfEndcaps) {
std::ostringstream sl,st;
sl << "Endcap" << std::setw(2) << std::setfill('0') << e;
sl << "Endcap" << std::setw(1) << std::setfill('0') << e;
st << "Endcap " << e;
sEndcapLabel[e]=sl.str();
sEndcapTitle[e]=st.str();
......@@ -43,7 +43,7 @@ public:
sLayerTitle[l]=st.str();
} else {
sLayerLabel[l]="LayerAll";
sLayerTitle[l]="All layers";
sLayerTitle[l]="AllLayers";
}
}
}
......
......@@ -61,10 +61,14 @@ public:
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++){
h_CluEtaResponse[e] = new TH1F(("h_CluEtaResponse_"+sEndcapLabel[e]).c_str(), ("h_CluEtaResponse_"+sEndcapTitle[e]).c_str(), 200, -0.5, 0.5);
h_CluPhiResponse[e] = new TH1F(("h_CluPhiResponse_"+sEndcapLabel[e]).c_str(), ("h_CluPhiResponse_"+sEndcapTitle[e]).c_str(), 200, -0.5, 0.5);
h_CluEtaResponseVsEtaTrue[e] = new TH2F(("h_CluEtaResponseVsEtaTrue_"+sEndcapLabel[e]).c_str(),("h_CluEtaResponseVsEtaTrue_"+sEndcapTitle[e]).c_str(),200,1.4,3.2, 200, -0.5, 0.5);
h_CluPhiResponseVsPhiTrue[e] = new TH2F(("h_CluPhiResponseVsPhiTrue_"+sEndcapLabel[e]).c_str(),("h_CluPhiResponseVsPhiTrue_"+sEndcapTitle[e]).c_str(),200,-3.5,3.5, 200, -0.5, 0.5);
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);
......@@ -80,7 +84,8 @@ public:
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(),
......@@ -121,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;
......@@ -145,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 ;
}
}
......@@ -163,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
......@@ -185,7 +193,7 @@ 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 ****/
......@@ -234,32 +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;
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].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_CluEtaResponse[e]->Fill( deltaEta( Eta_C2d_Max, vSigProj[i].FourP().Eta() ) );
h_CluPhiResponse[e]->Fill( deltaPhi( Phi_C2d_Max, vSigProj[i].FourP().Phi() ) );
h_CluEtaResponseVsEtaTrue[e]->Fill( fabs(vSigProj[i].FourP().Eta()), deltaPhi(Eta_C2d_Max, vSigProj[i].FourP().Eta() ) / vSigProj[i].FourP().Eta());
h_CluPhiResponseVsPhiTrue[e]->Fill( vSigProj[i].FourP().Phi(), deltaPhi(Phi_C2d_Max, vSigProj[i].FourP().Phi() ) / vSigProj[i].FourP().Phi() );
h_CluEtaResponse[Geometry::kNumberOfEndcaps]->Fill( deltaEta( Eta_C2d_Max, vSigProj[i].FourP().Eta() ) );
h_CluPhiResponse[Geometry::kNumberOfEndcaps]->Fill( deltaPhi( Phi_C2d_Max, vSigProj[i].FourP().Phi() ) );
h_CluEtaResponseVsEtaTrue[Geometry::kNumberOfEndcaps]->Fill( fabs(vSigProj[i].FourP().Eta()), deltaPhi(Eta_C2d_Max, vSigProj[i].FourP().Eta() ) / vSigProj[i].FourP().Eta());
h_CluPhiResponseVsPhiTrue[Geometry::kNumberOfEndcaps]->Fill( vSigProj[i].FourP().Phi(), deltaPhi(Phi_C2d_Max, vSigProj[i].FourP().Phi() ) / vSigProj[i].FourP().Phi() );
//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].FourP().Eta() ) / vSigProj[i].FourP().Eta() );
h_CluPhiResponseVsPhiC2d ->Fill( Phi_C2d_Max, 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].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){
......@@ -363,10 +386,10 @@ private:
// TH1F *h_CluPhiResponse;
// TH2F *h_CluEtaResponseVsEtaTrue;
// TH2F *h_CluPhiResponseVsPhiTrue;
TH1F *h_CluEtaResponse[Geometry::kNumberOfEndcaps+1];
TH1F *h_CluPhiResponse[Geometry::kNumberOfEndcaps+1];
TH2F *h_CluEtaResponseVsEtaTrue[Geometry::kNumberOfEndcaps+1];
TH2F *h_CluPhiResponseVsPhiTrue[Geometry::kNumberOfEndcaps+1];
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;
......@@ -387,6 +410,8 @@ private:
TH1F *h_NTrgHitPerC2d;
TH1F *h_NTrgHitPerLayer;
TH1F *h_NC2dPerLayer;
TH1F *h_PhiTrueFromXY;
TH1F *h_PhiTrueFrom4P;
double deltaPhi( double phi1, double phi2) {
......
......@@ -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;
......
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