///////////////////////////////////////////////////////////////////////
//                                                                   //
//  T49XiFinder applies the xifinder algorithm to mini-DST data.     //
//                                                                   //
///////////////////////////////////////////////////////////////////////

#include <stdio.h>

#include "TObjArray.h"

#include "T49ParticleRoot.h"
#include "T49VertexRoot.h"

#include "T49XiFinder.h"
#include "T49Trkstep.h"

ClassImp(T49XiFinder)

//______________________________________________________________________________
 T49XiFinder::T49XiFinder()
{
  //
  // T49XiFinder default constructor
  //

  fNeutralInputArray    = 0;
  fChargedInputArray    = 0;

  fXiArray              = 0;
  fNeutralDaughterArray = 0;
  fChargedDaughterArray = 0;

  fTrkstepCharged       = 0;
  fTrkstepXi            = 0;
  
  fZtarget              = 0;
  fZmin                 = 0;
  fZstepsize            = 0;
  fDcaXcut              = 0;
  fDcaYcut              = 0;
  fDcaRcut              = 0;
  fSideCut              = 0;
  fPtarmCut             = 0;
  fPtMaxCut             = 0;
  fPtotMaxCut           = 0;
  fXtargCut             = 0;
  fYtargCut             = 0;
  fV0MassMin            = 0;
  fV0MassMax            = 0;
  fV0Zmin               = 0;
  fV0ArmPtMax           = 0;
  fChargedXtargMin      = 0;
  fChargedYtargMin      = 0;
  fChargedNPointMin     = 0;

}

//______________________________________________________________________________
 T49XiFinder::~T49XiFinder()
{
  //
  // T49XiFinder destructor
  //

  if (fNeutralInputArray) {
    fNeutralInputArray->Delete();
    delete fNeutralInputArray;
  }
  if (fChargedInputArray) {
    fChargedInputArray->Delete();
    delete fChargedInputArray;
  }
  if (fXiArray) {
    fXiArray->Delete();
    delete fXiArray;
  }
  if (fNeutralDaughterArray) {
    fNeutralDaughterArray->Delete();
    delete fNeutralDaughterArray;
  }
  if (fChargedDaughterArray) {
    fChargedDaughterArray->Delete();
    delete fChargedDaughterArray;
  }

}

//______________________________________________________________________________
 void T49XiFinder::InitGlobalParameter()
{
  //
  // Initializes the T49XiFinder with the standard parameter for the global 
  // reconstruction and 256tb data.
  //

  fZtarget          = -581.1;
  fZmin             = -560.0; 
  fZstepsize        =    2.0;
  fDcaXcut          =    1.0;  
  fDcaYcut          =    1.0;  
  fDcaRcut          =  999.0;
  fSideCut          =    0.0;
  fPtarmCut         =    0.30;
  fPtMaxCut         =   50.0;
  fPtotMaxCut       =  500.0;
  fXtargCut         =    3.0;
  fYtargCut         =    3.0;
  fV0MassMin        =    1.101;
  fV0MassMax        =    1.131;
  fV0Zmin           = -510.0;
  fV0ArmPtMax       =    0.25;
  fChargedXtargMin  =    1.0;
  fChargedYtargMin  =    0.0;
  fChargedNPointMin =   10; 

  if (fTrkstepCharged) delete fTrkstepCharged;
  fTrkstepCharged = new T49Trkstep();
  fTrkstepCharged->SetKeys("3905","RU00","V00A","COMB/STD+","FULL","R95A");
  fTrkstepCharged->InitDS("trkstep_charged");  

  if (fTrkstepXi)      delete fTrkstepXi;
  fTrkstepXi      = new T49Trkstep();
  fTrkstepXi->SetKeys("3905","RU00","V00A","COMB/STD+","FULL","R95A");
  fTrkstepXi->InitDS("trkstep_xi");  

}

//______________________________________________________________________________
 Bool_t T49XiFinder::FindXis()
{
  //
  // Looks for Xi/Omega candidates employing the same procedure than the
  // xifinder of the reconstruction chain.
  //

  Int_t trkstepFlag = 0;

  if ((!fTrkstepXi) || (!fTrkstepCharged)) {
    printf("<T49XiFinder::FindXis> Define the T49Trkstep objects first.n");
    return kFALSE;
  }    

  // Create output arrays if neccessary
  if (!fXiArray) {
    fXiArray = new TObjArray();
  }
  else {
    fXiArray->Delete();
  }
  if (!fNeutralDaughterArray) {
    fNeutralDaughterArray = new TObjArray();
  }
  else {
    fNeutralDaughterArray->Delete();
  }
  if (!fChargedDaughterArray) {
    fChargedDaughterArray = new TObjArray();
  }
  else {
    fChargedDaughterArray->Delete();
  }

  // Pair combinations possible?
  if ((fChargedInputArray->GetEntries() > 0) &&
      (fNeutralInputArray->GetEntries() > 0)) {

    // Determine the largest vertez z of the lambda
    Float_t zLamMax = -9999.0;
    Float_t zLabMax = -9999.0;
    for (Int_t iNeutral = 0;
	 iNeutral < fNeutralInputArray->GetEntries();
         iNeutral++) {
      T49ParticleRoot *neutral = (T49ParticleRoot *) 
                                 fNeutralInputArray->UncheckedAt(iNeutral);
      if ((neutral->GetIflag() & 1)) {
        if (neutral->GetZFirst() > zLamMax) {
          zLamMax = neutral->GetZFirst();
        }
      }
      if ((neutral->GetIflag() & 2)) {
        if (neutral->GetZFirst() > zLabMax) {
          zLabMax = neutral->GetZFirst();
        }
      }
    }

    for (Int_t iCharged = 0; 
         iCharged < fChargedInputArray->GetEntries(); 
         iCharged++) {

      T49ParticleRoot *charged = (T49ParticleRoot *) 
                                 fChargedInputArray->UncheckedAt(iCharged);

      // Antiparticle ?
      Bool_t anti = kFALSE;
      if (charged->GetCharge() > 0) {
        anti = kTRUE;
      }

      // Reset the found flag and dcas
      for (Int_t iNeutral = 0; 
           iNeutral < fNeutralInputArray->GetEntries(); 
           iNeutral++) {
        T49ParticleRoot *neutral = (T49ParticleRoot *) 
                                   fNeutralInputArray->UncheckedAt(iNeutral);
        neutral->SetIdDet(0);
        neutral->SetTofX(999.0);
        neutral->SetTofY(999.0);
      }

      // Initialize trkstep. Tracking has to start from track position at target 
      // or the start vertex position (stored in the last point coordinates).
      fTrkstepCharged->SetCharge(charged->GetCharge());
      fTrkstepCharged->SetP(charged->GetPx()
                           ,charged->GetPy()
                           ,charged->GetPz());
      fTrkstepCharged->SetPos(charged->GetXLast()
                             ,charged->GetYLast()
                             ,charged->GetZLast());
      trkstepFlag = fTrkstepCharged->TrackTo(charged->GetZFirst()); 
      if (trkstepFlag) {
        printf("<T49XiFinder::FindXis> Trkstep error charged 0 (%d)n",trkstepFlag);
        break;
      }

      // Determine starting position 
      Float_t zstart = charged->GetZFirst();
      if (anti) {
        if (zstart > zLabMax) {
          zstart = zLabMax;
	}
      }
      else {
        if (zstart > zLamMax) {
          zstart = zLamMax;
        }
      }

      // Step backwards
      for (Float_t zstep = zstart; zstep > fZmin; zstep -= fZstepsize) {

        trkstepFlag = fTrkstepCharged->TrackTo(zstep); 
        if (trkstepFlag) {
          printf("<T49XiFinder::FindXis> Trkstep error charged 1 (%d)n",trkstepFlag);
          break;
	}

        Float_t x_c = fTrkstepCharged->GetX();
        Float_t y_c = fTrkstepCharged->GetY();

        for (Int_t iNeutral = 0; 
             iNeutral < fNeutralInputArray->GetEntries(); 
             iNeutral++) {

          T49ParticleRoot *neutral = (T49ParticleRoot *) 
                                     fNeutralInputArray->UncheckedAt(iNeutral);

          Float_t x_n     = neutral->GetXFirst();
          Float_t y_n     = neutral->GetYFirst();
          Float_t z_n     = neutral->GetZFirst();
          Float_t px_n    = neutral->GetPx();
          Float_t py_n    = neutral->GetPy();
          Float_t pz_n    = neutral->GetPz();
          Int_t   iflag_n = neutral->GetIflag();
          Int_t   found_n = neutral->GetIdDet();

          // Take only combinations where the (anti)lambda is in the right mass window
          if (( anti &&  (iflag_n & 2) ) ||
              (!anti &&  (iflag_n & 1))) {

            // (Anti)lambda vertex must be downstream
            if (zstep < (z_n + fZstepsize)) {

              if (fVerbose > 2) {
                printf("<T49XiFinder::FindXis>");
                printf(" combination charged %d + neutral %d, step = %fn"
                      ,iCharged,iNeutral,zstep); 
	      }

	      // Check whether combination was already found
              if (found_n == 0) {

		// Calculate dcas
                Float_t dcaX = TMath::Abs((x_n + (px_n/pz_n * (zstep - z_n))) -  x_c);
                Float_t dcaY = TMath::Abs((y_n + (py_n/pz_n * (zstep - z_n))) -  y_c);
                Float_t dcaR = TMath::Sqrt(dcaX*dcaX + dcaY*dcaY);

                if (fVerbose > 2) {
                  printf("<T49XiFinder::FindXis>");
                  printf(" dcaX = %f, dcaY = %f, dcaR = %fn",dcaX,dcaY,dcaR);
		}

		// Minimum in dca?
                Float_t vdca1 = neutral->GetTofX();
                Float_t vdca2 = neutral->GetTofY();
                if ((dcaX  <= fDcaXcut) && 
                    (dcaY  <= fDcaYcut) &&
                    (vdca1 <      dcaR) &&
                    (vdca1 <     vdca2) &&
                    (vdca2 <     999.0)) {

                  if (fVerbose > 1) {
                    printf("<T49XiFinder::FindXis>");
                    printf(" minimum in DCA foundn");
		  }

		  // Mark this combination as already found
                  neutral->SetIdDet(1);

		  // Find vertex z
                  Float_t a  = zstep;
                  Float_t b  = zstep +       fZstepsize;
                  Float_t c  = zstep + 2.0 * fZstepsize;
                  Float_t fa = dcaR;
                  Float_t fb = vdca1;
                  Float_t fc = vdca2;
                  Float_t dn = (((b-a)*(fb-fc)) - ((b-c)*(fb-fa))) * 2.0;
                  if (dn == 0) {
                    printf("<T49XiFinder::FindXis> Interpolation errorn");
                    goto failed_r_cut;
		  }
                  Float_t decayz = b - ((((b-a)*(b-a)*(fb-fc)) 
                                       - ((b-c)*(b-c)*(fb-fa))) / dn); 

		  // Track to the found vertex position
                  trkstepFlag = fTrkstepCharged->TrackTo(decayz);
                  if (trkstepFlag) {
                    printf("<T49XiFinder::FindXis> Trkstep error charged 2 (%d)n",trkstepFlag);
                    goto failed_r_cut;
	          }

                  x_c = fTrkstepCharged->GetX();
                  y_c = fTrkstepCharged->GetY();

		  // Calculate dcas again
                  dcaX = TMath::Abs((x_n + (px_n/pz_n * (zstep - z_n))) -  x_c);
                  dcaY = TMath::Abs((y_n + (py_n/pz_n * (zstep - z_n))) -  y_c);
                  dcaR = TMath::Sqrt(dcaX*dcaX + dcaY*dcaY);

                  // Momentum of charged particle at the Xi vertex
                  Float_t px_c = fTrkstepCharged->GetPx();
                  Float_t py_c = fTrkstepCharged->GetPy();
                  Float_t pz_c = fTrkstepCharged->GetPz();

		  // Good Xi candidate
                  if (dcaR > fDcaRcut) goto failed_r_cut;
                  if (fVerbose > 1) {
                    printf("<T49XiFinder::FindXis>");
                    printf(" candidate passed DCA cutn");
		  }

                  // Calculate the Xi decay vertex position
                  Float_t x_xi = x_n + (px_n/pz_n * (decayz - z_n));
                  Float_t y_xi = y_n + (py_n/pz_n * (decayz - z_n));
                  Float_t z_xi = decayz;

		  // The Xi invariant mass
                  Float_t e_pion_c = charged->GetE(kPionMass);
                  Float_t e_kaon_c = charged->GetE(kKaonMass);
                  Float_t e_n      = neutral->GetE(kLambdaMass);
                  Float_t p_n_c    = (px_c*px_n) + (py_c*py_n) + (pz_c*pz_n); 
                  Float_t diff_xi  = e_pion_c*e_n - p_n_c; 
                  Float_t diff_om  = e_kaon_c*e_n - p_n_c;
                  Float_t minv_xi  = kPionMass*kPionMass 
                                   + kLambdaMass*kLambdaMass 
                                   + 2.0 * diff_xi;
                  Float_t minv_om  = kKaonMass*kKaonMass 
                                   + kLambdaMass*kLambdaMass 
                                   + 2.0 * diff_om;
                  if (minv_xi > 0.0) {
                    minv_xi = TMath::Sqrt(minv_xi);
		  }
                  else {
                    minv_xi = -999.0;
		  }
                  if (minv_om > 0.0) {
                    minv_om = TMath::Sqrt(minv_om);
		  }
                  else {
                    minv_om = -999.0;
		  }

                  if ((minv_xi > 0.0) && (minv_om > 0.0)) {

  		    // The Xi momentum components at the decay vertex
                    Float_t px_xi   = px_c + px_n;
                    Float_t py_xi   = py_c + py_n;
                    Float_t pz_xi   = pz_c + pz_n;
                    Float_t ptot_xi = TMath::Sqrt(px_xi*px_xi 
                                                + py_xi*py_xi 
                                                + pz_xi*pz_xi);

		    // The armenteros variables
                    Float_t ptotsq_n = px_n*px_n + py_n*py_n + pz_n*pz_n;
                    Float_t ppar_n   = (px_xi*px_n + py_xi*py_n + pz_xi*pz_n)
                                     / ptot_xi;
                    Float_t ptarm    = ptotsq_n - (ppar_n*ppar_n);
                    if (ptarm >= 0.0) {
                      ptarm = TMath::Sqrt(ptarm);
		    }
                    else {
                      ptarm = 9999.0;
		    }

		    // The Xi parameters in the target plane
                    Float_t xtarg_xi = 999.0;
                    Float_t ytarg_xi = 999.0;
                    Float_t ztarg_xi = 999.0;
                    fTrkstepXi->SetCharge(charged->GetCharge());
                    fTrkstepXi->SetPos(x_xi,y_xi,z_xi);
                    fTrkstepXi->SetP(px_xi,py_xi,pz_xi);
                    trkstepFlag = fTrkstepXi->TrackTo(fZtarget);
                    if (trkstepFlag) {
                      printf("<T49XiFinder::FindXis> Trkstep error xi (%d)n",trkstepFlag);
                      break;
	            }
                    else {
                      xtarg_xi = fTrkstepXi->GetX();
                      ytarg_xi = fTrkstepXi->GetY();
                      ztarg_xi = fTrkstepXi->GetZ();
		    }
                    Float_t pxtarg_xi  = fTrkstepXi->GetPx();
                    Float_t pytarg_xi  = fTrkstepXi->GetPy();
                    Float_t pztarg_xi  = fTrkstepXi->GetPz();
                    Float_t pttarg_xi  = TMath::Sqrt(pxtarg_xi*pxtarg_xi 
                                                   + pytarg_xi*pytarg_xi);
                    Float_t x1minx2_xi = x_n + (px_n/pz_n * (fZtarget - z_n)) 
                                       - charged->GetTofX();

		    if (fVerbose) {
                      printf("<T49XiFinder::FindXis>");
                      printf(" found a candidate before cutsn");
		    }

		    // Apply the xifinder pair cuts
                    if (ptarm                 > fPtarmCut  ) goto found_one;
		    if (pttarg_xi             > fPtMaxCut  ) goto found_one;
                    if (ptot_xi               > fPtotMaxCut) goto found_one;
                    if (TMath::Abs(xtarg_xi)  > fXtargCut  ) goto found_one;
                    if (TMath::Abs(ytarg_xi)  > fYtargCut  ) goto found_one;
                    if (fSideCut * x_n * x_xi < 0.0        ) goto found_one;

		    if (fVerbose) {
                      printf("<T49XiFinder::FindXis>");
                      printf(" found a candidate after cutsn");
		    }

		    // Create a new particle structure for the neutral
		    // daughter and add it the to the output list
                    T49ParticleRoot *n_dght = new T49ParticleRoot(neutral);
                    fNeutralDaughterArray->Add(n_dght);
                       
		    // Create a new particle structure for the charged
		    // daughter and add it the to the output list
                    T49ParticleRoot *c_dght = new T49ParticleRoot(charged);
                    c_dght->SetPx(px_c);
                    c_dght->SetPy(py_c);
                    c_dght->SetPz(pz_c);
                    fChargedDaughterArray->Add(c_dght);

                    // Create a new vertex structure for this Xi/Omega 
		    // candidate and add it to the output list
                    T49VertexRoot *xi = new T49VertexRoot();
                    xi->SetX(x_xi);
                    xi->SetY(y_xi);
                    xi->SetZ(z_xi);
                    xi->SetXTarg(xtarg_xi);
                    xi->SetYTarg(ytarg_xi);
                    xi->SetPxXi(pxtarg_xi);
                    xi->SetPyXi(pytarg_xi);
                    xi->SetPzXi(pztarg_xi);
                    xi->SetX1minX2(x1minx2_xi);
                    xi->SetDaughters(c_dght,n_dght);
		    // For testing purposes
                    xi->SetPchi2(minv_xi);
                    fXiArray->Add(xi);

		  } // if: xi and omega mass ok

		} // if: minimum in dca

		// Store last two dcas
	        failed_r_cut:
                neutral->SetTofX(dcaR);
                neutral->SetTofY(vdca1);

	      } // if: lambda already used

	    } // if: lambda downstream

	  } // if: combination (anti)particle <-> (anti)particle

	  found_one:
	  continue;

	} // loop: neutral

      } //  loop: trkstep backwards for charged

    } // loop: charged 

  } // if: pair combinations

  return kTRUE;

}

//______________________________________________________________________________
 void T49XiFinder::CreateChargedInput(TObjArray *inputArray, Bool_t primary)
{
  //
  // Defines the charge input array by applying the xifinder cuts to the 
  // pion/kaon candidate.
  //

  if (!fTrkstepCharged) {
    printf("<T49XiFinder::SetChargedInput> Define the T49Trkstep objects first.n");
  }
  else {

    if (fChargedInputArray) {
      fChargedInputArray->Delete();
    }
    else {
      fChargedInputArray = new TObjArray();
    }
  
    for (Int_t idx = 0; idx < inputArray->GetEntries(); idx++) {

      T49ParticleRoot *charged = (T49ParticleRoot *) inputArray->UncheckedAt(idx);

      // Check whether track is trackable
      fTrkstepCharged->SetCharge(charged->GetCharge());
      if (primary) {
        charged->SetXLast(charged->GetBx());
        charged->SetYLast(charged->GetBy());
        charged->SetZLast(fZtarget);
      } 
      fTrkstepCharged->SetPos(charged->GetXLast()
                             ,charged->GetYLast()
                             ,charged->GetZLast());
      fTrkstepCharged->SetP(charged->GetPx()
                           ,charged->GetPy()
                           ,charged->GetPz());

      if (!(fTrkstepCharged->TrackTo(charged->GetZFirst()))) {

        if (!(fTrkstepCharged->TrackTo(fZtarget))) {

          Float_t xtarg_charged = fTrkstepCharged->GetX();
          Float_t ytarg_charged = fTrkstepCharged->GetY();

          // Check impact parameter and number of points only for primary tracks
          Bool_t accept = kTRUE;
          if (primary) {
            if ((TMath::Abs(ytarg_charged) < fChargedYtargMin ) &&
                (TMath::Abs(xtarg_charged) < fChargedXtargMin ) && 
                (charged->GetNPoint(3)     < fChargedNPointMin)) {
              accept = kFALSE;
	    }
	  }

          if (accept) {

   	    // Add to the charge daughter list
            T49ParticleRoot *chargedPart = new T49ParticleRoot();
            chargedPart->CopyParticle(charged);
            chargedPart->SetTofX(xtarg_charged);
            chargedPart->SetTofY(ytarg_charged);
            fChargedInputArray->Add(chargedPart);

	  }

        }

      }

    }

    if (fVerbose) {
      printf("<T49XiFinder::CreateChargedInput>");
      printf(" accepted %d charged input particlesn"
            ,fChargedInputArray->GetEntries());
    }

  }

}

//______________________________________________________________________________
 void T49XiFinder::CreateNeutralInput(TObjArray *inputArray, Bool_t primary)
{
  //
  // Defines the neutral input array by applying the xifinder cuts to the 
  // V0 candidate.
  //

  if (fNeutralInputArray) {
    fNeutralInputArray->Delete();
  }
  else {
    fNeutralInputArray = new TObjArray();
  }
  
  for (Int_t idx = 0; idx < inputArray->GetEntries(); idx++) {

    T49VertexRoot *neutral = (T49VertexRoot *) inputArray->UncheckedAt(idx);

    Int_t flag = 0;

    // Apply V0 cuts only for primary tracks
    Bool_t accept = kTRUE;
    if (primary) {
      if ((neutral->GetZ()            < fV0Zmin    ) &&
          (neutral->GetArmenterosPt() > fV0ArmPtMax)) {
        accept = kFALSE;
      }
    }

    if (accept) {

      // Check for a lambda candidate
      Float_t minvLambda     = neutral->GetInvariantMassLambda();
      if ((minvLambda     >= fV0MassMin) &&
          (minvLambda     <= fV0MassMax)) {
        flag += 1;
      }

      // Check for a antilambda candidate
      Float_t minvAntiLambda = neutral->GetInvariantMassAntiLambda();
      if ((minvAntiLambda >= fV0MassMin) &&
          (minvAntiLambda <= fV0MassMax)) {
        flag += 2;
      }

      // Add to the daughter list if neutral or antineutral
      if (flag) {
        T49ParticleRoot *neutralPart = new T49ParticleRoot();
        neutralPart->SetPx(neutral->GetPx());
        neutralPart->SetPy(neutral->GetPy());
        neutralPart->SetPz(neutral->GetPz());
        neutralPart->SetXFirst(neutral->GetX());
        neutralPart->SetYFirst(neutral->GetY());
        neutralPart->SetZFirst(neutral->GetZ());
        neutralPart->SetIflag(flag);
        neutralPart->SetIdDet(0);
        neutralPart->SetBx(neutral->GetXTarg());
        neutralPart->SetBy(neutral->GetYTarg());
        neutralPart->SetLabel(neutral->GetNTrkFit());
        fNeutralInputArray->Add(neutralPart);
      }

    }

  }

  if (fVerbose) {
    printf("<T49XiFinder::CreateNeutralInput>");
    printf(" accepted %d neutral input particlesn"
          ,fNeutralInputArray->GetEntries());
  }

}

//______________________________________________________________________________
 void T49XiFinder::SetTrkstepCharged(T49Trkstep *trkstep)
{
  //
  // Define T49Trkstep for the charged tracks.
  //

  if (fTrkstepCharged) delete fTrkstepCharged;
  fTrkstepCharged = trkstep;

}

//______________________________________________________________________________
 void T49XiFinder::SetTrkstepXi(T49Trkstep *trkstep)
{
  //
  // Define T49Trkstep for the Xi/Omega.
  //

  if (fTrkstepXi) delete fTrkstepXi;
  fTrkstepXi = trkstep;

}

//______________________________________________________________________________
 void T49XiFinder::DeleteChargedInput()
{
  //
  // Deletes the charged input list
  //

  if (fChargedInputArray) {
    fChargedInputArray->Delete();
    delete fChargedInputArray;
    fChargedInputArray = 0;
  }

}

//______________________________________________________________________________
 void T49XiFinder::DeleteNeutralInput()
{
  //
  // Deletes the neutral input list
  //

  if (fNeutralInputArray) {
    fNeutralInputArray->Delete();
    delete fNeutralInputArray;
    fNeutralInputArray = 0;
  }

}








ROOT page - Class index - Class Hierarchy - Top of the page

This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.