Commit afedd59b authored by dauncey's avatar dauncey

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

parents 4a65be59 6a80a1fa
......@@ -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";
}
}
}
......
......@@ -51,7 +51,7 @@ string hName(string sName, string endcap, string specificName){
class AnalysisTrgC2d : public AnalysisBase {
public:
AnalysisTrgC2d(const std::string &sRoot="", const std::string &sOut="")
AnalysisTrgC2d(const std::string &sRoot="", const std::string &sOut="" )
: AnalysisBase("AnalysisTrgC2d",sRoot,sOut) {
gEnlapsedTime = new TGraph();
......@@ -60,11 +60,16 @@ public:
hC3d[e] = new TH2F(hName(fName, sEndcapLabel[e], "C3d").c_str(),
hName(fName, sEndcapLabel[e], "C3d").c_str(),
100, 0, 1,
100, 0, 1);
200, -.5, .5,
200, -.5, .5);
for(int l=0; l<Geometry::kNumberOfLayers; l++){
gC3d[e][l] = new TGraph();
gC3d[e][l]->SetName( hName(fName, sEndcapLabel[e], sLayerLabel[l], "graphC3d").c_str() );
#ifdef JOHAN
hChargeVsRadius[e][l] = new TH2F( hName(fName, sEndcapLabel[e], sLayerLabel[l], "QvsR").c_str(),
hName(fName, sEndcapLabel[e], sLayerLabel[l], "QvsR").c_str(),
......@@ -237,7 +242,7 @@ public:
}
}
// cluShapes[e][l] = new TList;
// projShapes[e][l] = new TList;
projShapes[e][l] = new TList;
// cluShapesAll[e][l] = new TList;
}
}
......@@ -251,9 +256,15 @@ public:
// // clearing the lists
for(unsigned e(0); e<Geometry::kNumberOfEndcaps; e++) {
hC3d[e]->Reset("ICE");
for(unsigned l(0); l<Geometry::kNumberOfLayers; l++) {
gC3d[e][l]->Set(0);
// cluShapes[e][l]->Clear();
// projShapes[e][l]->Clear();
projShapes[e][l]->Clear();
// cluShapesAll[e][l]->Clear();
// hOccupancy_simHit[e][l]->Reset("ICE");
#ifdef C2D_OCCUPANCY
......@@ -262,6 +273,9 @@ public:
hOccupancyC2dEne[e][l]->Reset("ICE");
hOccupancyTrgHitEne[e][l]->Reset("ICE");
#endif
// hOccupancy_cilCoo[e][l]->Reset("ICE");
// vProjEneMips[e][l].clear();
}
......@@ -376,13 +390,13 @@ public:
}
}
// projections
//for(unsigned int i_par=0; i_par<vPhoton.size(); i_par++){
// Particle proj = vPhoton.at(i_par).project( Geometry::layerZ(e,l) );
// TEllipse* projShape = new TEllipse( proj.position().x(), proj.position().y(), cluRadius);
//projShapes[e][l]->Add( projShape );
/* store signal particle projections */
for(unsigned int i_par=0; i_par<vSignal.size(); i_par++){
Particle proj = vSignal.at(i_par).project( Geometry::layerZ(e,l) );
TEllipse* projShape = new TEllipse( proj.position().x(), proj.position().y(), cluRadius);
projShapes[e][l]->Add( projShape );
//vProjEneMips[e][l].push_back(0.);
//}
}
/**** Getting Clusters ****/
std::vector<TrgC2d> &vTrgC2d( event.trgC2ds(e,l) );
......@@ -432,8 +446,13 @@ public:
i_clu++;
/* for 3d clustering analysis */
hC3d[e]->Fill( clu->center().first/Geometry::layerZ(e, l), clu->center().second/Geometry::layerZ(e, l) );
hC3d[e]->Fill( clu->center().first/Geometry::layerZ(e, l),
clu->center().second/Geometry::layerZ(e, l)
);
gC3d[e][l]->SetPoint( gC3d[e][l]->GetN(),
clu->center().first/Geometry::layerZ(e, l),
clu->center().second/Geometry::layerZ(e, l)
);
/* only clusters with sim sig ene */
if( clu->containsSigEne() ){
nCluWithSimEne++;
......@@ -493,7 +512,7 @@ public:
else{
bool found=false;
for( unsigned i=0; i<waferIds.size(); i++ ){
if(clu->trgHits().at(i_tc)->wafer() == waferIds.at(i) )
if(clu->trgHits().at(i_tc)->wafer() == (int)waferIds.at(i) )
found=true;
}
if( !found )
......@@ -573,7 +592,7 @@ public:
//
/* write single event file */
writeEventOnDisk=true;
writeEventOnDisk=false;
if(writeEventOnDisk)
this->write( fEventNumber );
......@@ -593,13 +612,20 @@ public:
std::string fileName("AnalysisTrgC2dSpecific");
TFile fs( (fileName+"_"+sEvt.c_str()+".root").c_str(), "RECREATE" );
gEnlapsedTime->Write("TimeEnlapsed");
hC3d[0]->Write("endcap0_C3d");
hC3d[1]->Write("endcap1_C3d");
//gC3d[0]->Write("endcap0_C3dgraph");
//gC3d[1]->Write("endcap1_C3dgraph");
for(unsigned e(0); e<Geometry::kNumberOfEndcaps; e++) {
for(unsigned l(0); l<Geometry::kNumberOfLayers; l++) {
std::string label("");
// TFile f( (label+"Evt"+sEvt+sEndcapLabel[e]+sLayerLabel[l]+".root").c_str(), "RECREATE" );
TDirectory* dir = fs.mkdir( (label+"Evt"+sEvt+sEndcapLabel[e]+sLayerLabel[l]).c_str() );
dir->cd();
gC3d[e][l]->Write();
#ifdef C2D_OCCUPANCY
hOccupancyEmpty[e][l]->Write();
hOccupancyC2d[e][l]->Write();
......@@ -611,8 +637,7 @@ public:
// hOccupancy_simHit[e][l]->Write();
// cluShapesAll[e][l]->Write("cluShapesAll", TObject::kSingleKey);
// cluShapes[e][l]->Write("cluShapes", TObject::kSingleKey);
// projShapes[e][l]->Write("projShapes", TObject::kSingleKey);
// f.Close();
projShapes[e][l]->Write("projShapes", TObject::kSingleKey);
}
}
......@@ -631,7 +656,6 @@ private:
TH2Poly* hOccupancy_cilCoo[Geometry::kNumberOfEndcaps][Geometry::kNumberOfLayers];
TH2Poly* hOccupancy_simHit[Geometry::kNumberOfEndcaps][Geometry::kNumberOfLayers];
TH1F *hNcluPerLayer[Geometry::kNumberOfEndcaps][Geometry::kNumberOfLayers];
/* energy relaterd histos */
......@@ -659,8 +683,8 @@ private:
TH2F *hCluRealCluEneRatioVsDistance[Geometry::kNumberOfEndcaps][Geometry::kNumberOfLayers];
TH2F *hC3d[Geometry::kNumberOfEndcaps];
TGraph *gC3d[Geometry::kNumberOfEndcaps][Geometry::kNumberOfLayers];
TTree pippa;
#ifdef JOHAN
TH2F *hChargeVsRadius[Geometry::kNumberOfEndcaps][Geometry::kNumberOfLayers];
......
This diff is collapsed.
......@@ -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) {
......
......@@ -1236,6 +1236,15 @@ static double fTrgHitSelectorMipsCut;
static unsigned fTrgC2dCreatorMethod;
static unsigned fTrgC3dCreatorMethod;
// 2dC
static double fENERGY_CUT_Et;
static double fSEED_ENERGY_CUT_Et;
static double fENERGY_CUT;
static double fSEED_ENERGY_CUT;
// 3dC
static double fNormCooCut;
protected:
//const Geometry &fGeometry;
......@@ -1286,8 +1295,19 @@ double Event::fTrgHitSelectorMipsCut=5.0;
unsigned Event::fTrgC2dCreatorMethod=4; /*near neighbours (version 3)*/
unsigned Event::fTrgC3dCreatorMethod=0; /* eta phi clustering */
unsigned Event::fTrgC3dCreatorMethod=2; /* eta phi clustering */
/* energy thr for 2d clustering (in mips) */
double Event::fENERGY_CUT_Et = 15.;
double Event::fSEED_ENERGY_CUT_Et = 5.;
double Event::fENERGY_CUT = 25.;
double Event::fSEED_ENERGY_CUT = 10.;
/* cut in normalized coordinate for 3d clustering */
double Event::fNormCooCut = .01;
//double Event::fNormCooCut=0.2; /* norm clustering cu in x/z y/z plane */
#endif
......
......@@ -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;
......
......@@ -287,6 +287,21 @@ public:
}
}
pair<float,float> centerNorm(){
pair<float,float> p(0, 0);
float layerZ = Geometry::layerZ( this->endcap(), this->layer() );
if( vTrgHit.size() > 0 ){
return pair<float,float>( fCenter.first/layerZ, fCenter.second/layerZ );
}
else{
return p;
}
}
pair<float,float> centerRhoPhi(){
pair<float,float> rphi(0, 0);
......@@ -352,6 +367,10 @@ public:
}
/* DISTANCES */
double distXY(TrgC2d& trgObj){
......
......@@ -12,6 +12,8 @@
class TrgC3d {
public:
TrgC3d() : fParticle(0) {
fCenterNorm.first = 0;
fCenterNorm.second = 0;
}
~TrgC3d() {
......@@ -24,6 +26,14 @@ public:
bool addTrgC2d(TrgC2d *s) {
assert(s->endcap()==fEndcap);
vTrgC2d.push_back(s);
if( vTrgC2d.size() == 1 ){
fCenterNorm.first = ( s->centerNorm().first ) ;
fCenterNorm.second = ( s->centerNorm().second ) ;
}
else{
fCenterNorm.first = ( fCenterNorm.first + s->centerNorm().first ) / 2;
fCenterNorm.second = ( fCenterNorm.second + s->centerNorm().second ) / 2;
}
return true;
}
......@@ -169,6 +179,15 @@ public:
}
*/
pair<float,float> centerNorm(){
return fCenterNorm;
}
float C2dDistanceNorm( pair<float,float> C2dCenter ){
double dist = sqrt( (fCenterNorm.first-C2dCenter.first)*(fCenterNorm.first-C2dCenter.first) + (fCenterNorm.second-C2dCenter.second)*(fCenterNorm.second-C2dCenter.second) ) ;
return dist;
}
Point position() const {
return fPosition;
}
......@@ -189,6 +208,14 @@ public:
return e;
}
bool containsSigEne(){
for(unsigned i(0);i<vTrgC2d.size();i++)
if( vTrgC2d[i]->containsSigEne() )
return true;
return false;
}
double mips() const {
double e(0.0);
for(unsigned i(0);i<vTrgC2d.size();i++) {
......@@ -197,6 +224,14 @@ public:
return e;
}
double transverseMips() const {
double e(0.0);
for(unsigned i(0);i<vTrgC2d.size();i++) {
e+=vTrgC2d[i]->transverseMips();
}
return e;
}
unsigned numberOfTrgC2ds() const {
return vTrgC2d.size();
}
......@@ -239,6 +274,24 @@ public:
return e;
}
double simMips(bool signal=false) const {
double e(0.0);
for(unsigned i(0);i<vTrgC2d.size();i++) {
e+=vTrgC2d[i]->simMips(signal);
}
return e;
}
double simTransverseMips(bool signal=false) const {
double e(0.0);
for(unsigned i(0);i<vTrgC2d.size();i++) {
e+=vTrgC2d[i]->simTransverseMips(signal);
}
return e;
}
/*
Point2D scale2D() const {
assert(vTrgC2d.size()>0);
......@@ -316,6 +369,8 @@ private:
double fChiSquared;
unsigned fNumberOfDof;
Particle *fParticle;
pair<float,float> fCenterNorm;
};
#endif
void trgC2dCreator() {
void trgC2dCreator(
)