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


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


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

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

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]


NandaWattimena/Neutralinos (last edited 2010-03-26 09:31:02 by NandaWattimena)