RALICE/icepack/IceLinefit.cxx

   1 /*******************************************************************************
   2  * Copyright(c) 2003, IceCube Experiment at the South Pole. All rights reserved.
   3  *
   4  * Author: The IceCube RALICE-based Offline Project.
   5  * Contributors are mentioned in the code where appropriate.
   6  *
   7  * Permission to use, copy, modify and distribute this software and its
   8  * documentation strictly for non-commercial purposes is hereby granted
   9  * without fee, provided that the above copyright notice appears in all
  10  * copies and that both the copyright notice and this permission notice
  11  * appear in the supporting documentation.
  12  * The authors make no claims about the suitability of this software for
  13  * any purpose. It is provided "as is" without express or implied warranty.
  14  *******************************************************************************/
  15
  16 // $Id$
  17
  18 ///////////////////////////////////////////////////////////////////////////
  19 // Class IceLinefit
  20 // TTask derived class to perform a linefit track reconstruction.
  21 // The procedure is based on the method described in the Amanda publication
  22 // in Nuclear Instruments and Methods A524 (2004) 179-180.
  23 // To prevent waisting CPU time in trying to reconstruct (high-energy) cascade
  24 // events, or to select specifically reconstruction of low multiplicity events,
  25 // the user may invoke the memberfunctions SetMaxModA() and SetMinModA.
  26 // This allows selection of events for processing with a certain maximum
  27 // and/or minimum number of good Amanda OMs firing.
  28 // By default the minimum and maximum are set to 0 and 999, respectively,
  29 // in the constructor, which implies no multiplicity selection.
  30 //
  31 // The reconstructed track is stored in the IceEvent structure with as
  32 // default "IceLinefit" as the name of the track.
  33 // This track name identifier can be modified by the user via the
  34 // SetTrackName() memberfunction. This will allow unique identification
  35 // of tracks which are produced when re-processing existing data with
  36 // different criteria.
  37 // The track 3-momentum is set to the reconstructed velocity, normalised
  38 // to 1 GeV. The mass and charge of the track are left 0.
  39 // The r0 and t0 can be obtained from the reference point of the track,
  40 // whereas the t0 ia also available from the track timestamp .
  41 // By default the charge of the produced tracks is set to 0, since
  42 // no distinction can be made between positive or negative tracks.
  43 // However, the user can define the track charge by invokation
  44 // of the memberfunction SetCharge().
  45 // This facility may be used to distinguish tracks produced by the
  46 // various reconstruction algorithms in a (3D) colour display
  47 // (see the class AliHelix for further details).
  48 //
  49 // For further details the user is referred to NIM A524 (2004) 169.
  50 //
  51 // Note : This algorithm works best on data which has been calibrated
  52 //        (IceCalibrate), cross talk corrected (IceXtalk) and cleaned
  53 //        from noise hits etc. (IceCleanHits).
  54 //
  55 //--- Author: Nick van Eijndhoven 10-mar-2006 Utrecht University
  56 //- Modified: NvE $Date$ Utrecht University
  57 ///////////////////////////////////////////////////////////////////////////
  58
  59 #include "IceLinefit.h"
  60 #include "Riostream.h"
  61
  62 ClassImp(IceLinefit) // Class implementation to enable ROOT I/O
  63
  64 IceLinefit::IceLinefit(const char* name,const char* title) : TTask(name,title)
  65 {
  66 // Default constructor.
  67  fMaxmodA=999;
  68  fMinmodA=0;
  69  fTrackname="IceLinefit";
  70  fCharge=0;
  71 }
  72 ///////////////////////////////////////////////////////////////////////////
  73 IceLinefit::~IceLinefit()
  74 {
  75 // Default destructor.
  76 }
  77 ///////////////////////////////////////////////////////////////////////////
  78 void IceLinefit::SetMaxModA(Int_t nmax)
  79 {
  80 // Set the maximum number of good Amanda modules that may have fired
  81 // in order to process this event.
  82 // This allows suppression of processing (high-energy) cascade events
  83 // with this linefit tracking to prevent waisting cpu time for cases
  84 // in which tracking doesn't make sense anyhow.
  85 // Furthermore it allows selection of low multiplicity events for processing.
  86 // By default the maximum number of Amanda modules is set to 999 in the ctor,
  87 // which implies no selection on maximum module multiplicity.
  88 // See also the memberfunction SetMinModA().
  89  fMaxmodA=nmax;
  90 }
  91 ///////////////////////////////////////////////////////////////////////////
  92 void IceLinefit::SetMinModA(Int_t nmin)
  93 {
  94 // Set the minimum number of good Amanda modules that must have fired
  95 // in order to process this event.
  96 // This allows selection of a minimal multiplicity for events to be processed.
  97 // By default the minimum number of Amanda modules is set to 0 in the ctor,
  98 // which implies no selection on minimum module multiplicity.
  99 // See also the memberfunction SetMaxModA().
 100  fMinmodA=nmin;
 101 }
 102 ///////////////////////////////////////////////////////////////////////////
 103 void IceLinefit::SetTrackName(TString s)
 104 {
 105 // Set (alternative) name identifier for the produced first guess tracks.
 106 // This allows unique identification of (newly) produced linefit tracks
 107 // in case of re-processing of existing data with different criteria.
 108 // By default the produced first guess tracks have the name "IceLinefit"
 109 // which is set in the constructor of this class.
 110  fTrackname=s;
 111 }
 112 ///////////////////////////////////////////////////////////////////////////
 113 void IceLinefit::SetCharge(Float_t charge)
 114 {
 115 // Set user defined charge for the produced first guess tracks.
 116 // This allows identification of these tracks on color displays.
 117 // By default the produced first guess tracks have charge=0
 118 // which is set in the constructor of this class.
 119  fCharge=charge;
 120 }
 121 ///////////////////////////////////////////////////////////////////////////
 122 void IceLinefit::Exec(Option_t* opt)
 123 {
 124 // Implementation of the linefit reconstruction.
 125
 126  TString name=opt;
 127  AliJob* parent=(AliJob*)(gROOT->GetListOfTasks()->FindObject(name.Data()));
 128
 129  if (!parent) return;
 130
 131  IceEvent* evt=(IceEvent*)parent->GetObject("IceEvent");
 132  if (!evt) return;
 133
 134  // Fetch all fired Amanda OMs for this event
 135  TObjArray* aoms=evt->GetDevices("IceAOM");
 136  Int_t naoms=aoms->GetEntries();
 137  if (!naoms) return;
 138
 139  // Check for the minimum and/or maximum number of good fired Amanda OMs
 140  Int_t ngood=0;
 141  for (Int_t iom=0; iom<naoms; iom++)
 142  {
 143   IceGOM* omx=(IceGOM*)aoms->At(iom);
 144   if (!omx) continue;
 145   if (omx->GetDeadValue("ADC") || omx->GetDeadValue("LE") || omx->GetDeadValue("TOT")) continue;
 146   ngood++;
 147  }
 148  if (ngood<fMinmodA || ngood>fMaxmodA) return;
 149
 150  AliSignal* sx=0;
 151  Ali3Vector rom,sumr;
 152  Ali3Vector rt,sumrt;
 153  Float_t thit;
 154  Float_t sumt=0,sumt2=0;
 155  TObjArray hits;
 156
 157  // Loop over all OMs and hits to determine the linefit parameters.
 158  // Also all the used hits are recorded for association with the track.
 159  for (Int_t iom=0; iom<naoms; iom++)
 160  {
 161   IceGOM* omx=(IceGOM*)aoms->At(iom);
 162   if (!omx) continue;
 163   if (omx->GetDeadValue("LE")) continue;
 164   rom=(Ali3Vector)omx->GetPosition();
 165   // Use all the good hits of this OM
 166   for (Int_t ih=1; ih<=omx->GetNhits(); ih++)
 167   {
 168    sx=omx->GetHit(ih);
 169    if (!sx) continue;
 170    if (sx->GetDeadValue("ADC") || sx->GetDeadValue("LE") || sx->GetDeadValue("TOT")) continue;
 171
 172    thit=sx->GetSignal("LE",7);
 173    rt=rom*thit;
 174    sumr+=rom;
 175    sumrt+=rt;
 176    sumt+=thit;
 177    sumt2+=thit*thit;
 178
 179    // Record this hit for association with the track
 180    hits.Add(sx);
 181   }
 182  }
 183
 184  Int_t nused=hits.GetEntries();
 185  if (!nused) return;
 186
 187  sumr/=float(nused);
 188  sumrt/=float(nused);
 189  sumt/=float(nused);
 190  sumt2/=float(nused);
 191
 192  Ali3Vector v;
 193  Ali3Vector temp;
 194  temp=sumr*sumt;
 195  v=sumrt-temp;
 196  Float_t dum=sumt2-(sumt*sumt);
 197  if (dum) v/=dum;
 198
 199  Ali3Vector r;
 200  temp=v*sumt;
 201  r=sumr-temp;
 202
 203  AliTrack t;
 204  t.SetNameTitle(fTrackname.Data(),"IceLinefit linefit track");
 205  t.SetCharge(fCharge);
 206  evt->AddTrack(t);
 207  AliTrack* trk=evt->GetTrack(evt->GetNtracks());
 208  if (!trk) return;
 209
 210  Ali3Vector p;
 211  Float_t vec[3];
 212  v.GetVector(vec,"sph");
 213  vec[0]=1;
 214  p.SetVector(vec,"sph");
 215
 216  AliPosition r0;
 217  r0.SetPosition(r);
 218  r0.SetTimestamp((AliTimestamp&)*evt);
 219  AliTimestamp* t0=r0.GetTimestamp();
 220  t0->Add(0,0,(int)sumt);
 221
 222  trk->Set3Momentum(p);
 223  trk->SetReferencePoint(r0);
 224  trk->SetTimestamp(*t0);
 225
 226  // Link the used hits to the track (and vice versa)
 227  for (Int_t i=0; i<nused; i++)
 228  {
 229   sx=(AliSignal*)hits.At(i);
 230   if (sx) sx->AddLink(trk);
 231  }
 232 }
 233 ///////////////////////////////////////////////////////////////////////////