Long-lived Neutralinos
In some Gauge-Mediated Supersymmetry Breaking (GMSB) scenarios, the ILC could produce
e+ e- -> neutralino_1 neutralino_1 -> photon Gravitino photon Gravitino
In this case the neutralino_1 mass and lifetime can only be reconstructed from the observed photon shower.
Creating Long-lived neutral particles
How to create the SUSY spectrum with SPheno, simulate the interactions with WHiZard and track them through the detector with Geant4/Mokka, is described on the respective Wiki pages.
I have created samples with a neutralino mass of m_x = 151 GeV and a neutralino lifetime of t_x = 0.2ns; 2.0 ns and 11.0 ns. They are located in the DESY dCache
/pnfs/desy.de/flc/users/nwattime/neutralinoAnalysis/simulatedFiles/ILD00/neutralinos/
Event Reconstruction
I have used the RootTreeWriter for the event reconstruction. My ClusterWriteEngine(ClusterWriteEngine.cc,ClusterWriteEngine.hh) fills a root tree with
int iEvt; // event number int nCluster; // # cluster per event int nGamma; // # photons per event float etMiss; // missing transverse energy [GeV] float esumCalos; // energy sum in calorimeters [GeV] int nHitsCalos; // # hits in calorimeters float esumECAL; // energy sum in ECAL [GeV] int nHitsECAL; // # hits in ECAL float esumHCAL; // energy sum in HCAL [GeV] int nHitsHCAL; // # hits in HCAL float energy[MAXCLUSTER]; // energy per cluster [GeV] float clusterPos[MAXCLUSTER][3]; // cluster center of gravity float thetaDir[MAXCLUSTER]; // cluster angle theta [degree] float phiDir[MAXCLUSTER]; // cluster angle phi [degree] float cosThetaPos[MAXCLUSTER]; // cosine of cluster position theta float phiPos[MAXCLUSTER]; // cluster position phi [degree] float dca[MAXCLUSTER]; // distance of closest approach to IP [mm] float clusterAxis[MAXCLUSTER][3]; // cluster axis of inertia int nHitsPerCluster[MAXCLUSTER]; // #hits per cluster float trueDist[MAXCLUSTER]; // true distance of MCparticle to IP [mm] float trueGammaEnergy[MAXCLUSTER]; // true energy of MCparticle gamma [GeV] float trueGammaDirX[MAXCLUSTER]; // true x direction of MCparticle gamma [GeV] float trueGammaDirY[MAXCLUSTER]; // true y direction of MCparticle gamma [GeV] float trueGammaDirZ[MAXCLUSTER]; // true z direction of MCparticle gamma [GeV] float trueGammaPosX[MAXCLUSTER]; // true x position of MCparticle gamma [GeV] float trueGammaPosY[MAXCLUSTER]; // true y position of MCparticle gamma [GeV] float trueGammaPosZ[MAXCLUSTER]; // true z position of MCparticle gamma [GeV] float trueNeutralinoEnergy[MAXCLUSTER]; // true energy of MCparticle neutralino [GeV] float trueNeutralinoDirX[MAXCLUSTER]; // true x direction of MCparticle neutralino [GeV] float trueNeutralinoDirY[MAXCLUSTER]; // true y direction of MCparticle neutralino [GeV] float trueNeutralinoDirZ[MAXCLUSTER]; // true z direction of MCparticle neutralino [GeV] float trueThetaDir[MAXCLUSTER]; // true gamma angle theta [degree] float truePhiDir[MAXCLUSTER]; // true gamma angle phi [degree] int trueNcluster; // true # cluster / event
Background
I have only considered background that yields exactly the same detector signature as my signal, namely two highly energetic photons, missing energy and nothing else. This turned out to be e+ e- -> neutralino neutralino gamma gamma events. I have taken them from the official LoI production. The production is called Slac_SM_LDCPrime_02Sc_ppr002, and the chosen files are named
- M06-06-p03_ppr002_n1n1aa_w11751_500_LDCPrime_02Sc_LCP_ep+1.0_em-1.0_Slac_SM
- M06-06-p03_ppr002_n1n1aa_w11752_500_LDCPrime_02Sc_LCP_ep-1.0_em+1.0_Slac_SM
- M06-06-p03_ppr002_n2n2aa_w11755_500_LDCPrime_02Sc_LCP_ep+1.0_em-1.0_Slac_SM
- M06-06-p03_ppr002_n2n2aa_w11756_500_LDCPrime_02Sc_LCP_ep-1.0_em+1.0_Slac_SM
- M06-06-p03_ppr002_n3n3aa_w11759_500_LDCPrime_02Sc_LCP_ep+1.0_em-1.0_Slac_SM
- M06-06-p03_ppr002_n3n3aa_w11760_500_LDCPrime_02Sc_LCP_ep-1.0_em+1.0_Slac_SM
To yield the proper beam polarization of P(e+,e-) = (+80%,-60%), the samples were mixed as 2% of the ep+1.0 em-1.0 with 72% of the ep-1.0 em+1.0 sample. For Root trees this can be done with
float statWeight[nFiles]; //= nExpected*polWeight/nTrue outTree->SetWeight(statWeight);
Event Selection
Signal and background can be seperated by some simple cuts
//selection criteria const int _minNhitsEcal = 1500; //minimum # hits in ECAL: 1500 const int _maxNhitsEcal = 6000; //maximum # hits in ECAL: 6000 const float _EminEcal = 80; //minimum energy in ECAL: 80 GeV const float _EmaxEcal = 450; //maximum energy in ECAL: 450 GeV const float _minEt_miss = 30; //minimum transverse energy: 30 GeV const float _minNgamma = 2; //minimum # of cluster: 2 const float _minEcluster = 20; //minimum energy of cluster: 20 GeV const float _maxCosTheta = 0.75; //maximum cluster theta position: 0.75
the attached event selection root script shows how I did the event selection.
Mass Reconstruction
The neutralino mass (m_x) can be reconstructed from the endpoint of the photon energy spectrum, according to
E_max = 1/4 [ sqrt(s) + sqrt(s - 4*m_x^2) ] -> m_x = 1/2 * sqrt( s + [4*E_max - sqrt(s)]^2 )
In the mass reconstruction root script I have fitted the upper edge with
E(A;m_x^2;S;B) = A/( 1 + exp[ { E + 1/4( sqrt(s) + sqrt(s-4*m_x^2 ) } / S ] ) + B
where "A" is the amplitude, "S" represents the slope, "B" corresponds to underlying background and sqrt(s) is the center of mass energy.
The fit is relatively unstable, thus binning and starting conditions should be chosen carefully.
Lifetime Reconstruction
The neutralino decay length (l_x) can be reconstructed from the vector "<x>" between the interaction point and the photon shower center of gravity, the angle "Psi" between this vector and the photon shower main principal axis, and the decay angle "Delta" between the neutralino and the photon
cos(Delta) = E_x / p_x - m_x^2/(2*p_x*E_gamma) -> l_x = <x> * sin(Psi)/sin(Delta)
where the neutralino energy E_x=sqrt(s)/2 corresponds to half of the center of mass energy and the neutralino momentum is given as p_x = sqrt( E_x - m_x).
The neutralino decay length than translates into the neutralino lifetime (t_x) with
t_x = c / ( l_x *sqrt[ E_x^2/ m_x^2 - 1 ] )
As shown in the lifetime reconstruction root script, I have fitted the resulting spectrum with the rooFit package. I have chosen the convolution of a Gaussian, to account for remaining background, and an exponential function:
f(t) = f_sig ( exp[ -t/t_x] otimes 1/[sigma*2*pi]*exp[-t'^2/2] ) otimes f_bkg (1/[sigma*2*pi]*exp[-t'^2/2]