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 // The maximum number of good hits per Amanda OM to be used for the reconstruction
  31 // can be specified via the memberfunction SetMaxHitsA().
  32 // By default all good hits of each Amanda OM are used but the user may want
  33 // to restrict this number to the first n hits of each Amanda OM to account
  34 // for possible noise and/or afterpulse signals that are not recognised by the
  35 // hit cleaning procedure.
  36 //
  37 // Information about the actual parameter settings can be found in the event
  38 // structure itself via the device named "IceLinefit".
  39 //
  40 // The reconstructed track is stored in the IceEvent structure with as
  41 // default "IceLinefit" as the name of the track.
  42 // This track name identifier can be modified by the user via the
  43 // SetTrackName() memberfunction. This will allow unique identification
  44 // of tracks which are produced when re-processing existing data with
  45 // different criteria.
  46 // The track 3-momentum is set to the reconstructed velocity, normalised
  47 // to 1 GeV. The mass and charge of the track are left 0.
  48 // The r0 and t0 can be obtained from the reference point of the track,
  49 // whereas the t0 ia also available from the track timestamp .
  50 // By default the charge of the produced tracks is set to 0, since
  51 // no distinction can be made between positive or negative tracks.
  52 // However, the user can define the track charge by invokation
  53 // of the memberfunction SetCharge().
  54 // This facility may be used to distinguish tracks produced by the
  55 // various reconstruction algorithms in a (3D) colour display
  56 // (see the class AliHelix for further details).
  57 // The value of beta=v/c for the reconstructed velocity is available
  58 // from the fitdetails as stored for the reconstructed track.
  59 //
  60 // For further details the user is referred to NIM A524 (2004) 169.
  61 //
  62 // Note : This algorithm works best on data which has been calibrated
  63 //        (IceCalibrate), cross talk corrected (IceXtalk) and cleaned
  64 //        from noise hits etc. (IceCleanHits).
  65 //
  66 //--- Author: Nick van Eijndhoven 10-mar-2006 Utrecht University
  67 //- Modified: NvE $Date$ Utrecht University
  68 ///////////////////////////////////////////////////////////////////////////
  69
  70 #include "IceLinefit.h"
  71 #include "Riostream.h"
  72
  73 ClassImp(IceLinefit) // Class implementation to enable ROOT I/O
  74
  75 IceLinefit::IceLinefit(const char* name,const char* title) : TTask(name,title)
  76 {
  77 // Default constructor.
  78  fMaxmodA=999;
  79  fMinmodA=0;
  80  fMaxhitsA=0;
  81  fTrackname="IceLinefit";
  82  fCharge=0;
  83 }
  84 ///////////////////////////////////////////////////////////////////////////
  85 IceLinefit::~IceLinefit()
  86 {
  87 // Default destructor.
  88 }
  89 ///////////////////////////////////////////////////////////////////////////
  90 void IceLinefit::SetMaxModA(Int_t nmax)
  91 {
  92 // Set the maximum number of good Amanda modules that may have fired
  93 // in order to process this event.
  94 // This allows suppression of processing (high-energy) cascade events
  95 // with this linefit tracking to prevent waisting cpu time for cases
  96 // in which tracking doesn't make sense anyhow.
  97 // Furthermore it allows selection of low multiplicity events for processing.
  98 // By default the maximum number of Amanda modules is set to 999 in the ctor,
  99 // which implies no selection on maximum module multiplicity.
 100 // See also the memberfunction SetMinModA().
 101  fMaxmodA=nmax;
 102 }
 103 ///////////////////////////////////////////////////////////////////////////
 104 void IceLinefit::SetMinModA(Int_t nmin)
 105 {
 106 // Set the minimum number of good Amanda modules that must have fired
 107 // in order to process this event.
 108 // This allows selection of a minimal multiplicity for events to be processed.
 109 // By default the minimum number of Amanda modules is set to 0 in the ctor,
 110 // which implies no selection on minimum module multiplicity.
 111 // See also the memberfunction SetMaxModA().
 112  fMinmodA=nmin;
 113 }
 114 ///////////////////////////////////////////////////////////////////////////
 115 void IceLinefit::SetMaxHitsA(Int_t nmax)
 116 {
 117 // Set the maximum number of good hits per Amanda module to be processed.
 118 //
 119 // Special values :
 120 // nmax = 0 : No maximum limit set; all good hits will be processed
 121 //      < 0 : No hits will be processed
 122 //
 123 // In case the user selects a maximum number of good hits per module, all the
 124 // hits of each module will be ordered w.r.t. increasing hit time (LE).
 125 // This allows selection of processing e.g. only the first good hits etc...
 126 // By default the maximum number of hits per Amanda modules is set to 0 in the ctor,
 127 // which implies just processing all good hits without any maximum limit.
 128  fMaxhitsA=nmax;
 129 }
 130 ///////////////////////////////////////////////////////////////////////////
 131 void IceLinefit::SetTrackName(TString s)
 132 {
 133 // Set (alternative) name identifier for the produced first guess tracks.
 134 // This allows unique identification of (newly) produced linefit tracks
 135 // in case of re-processing of existing data with different criteria.
 136 // By default the produced first guess tracks have the name "IceLinefit"
 137 // which is set in the constructor of this class.
 138  fTrackname=s;
 139 }
 140 ///////////////////////////////////////////////////////////////////////////
 141 void IceLinefit::SetCharge(Float_t charge)
 142 {
 143 // Set user defined charge for the produced first guess tracks.
 144 // This allows identification of these tracks on color displays.
 145 // By default the produced first guess tracks have charge=0
 146 // which is set in the constructor of this class.
 147  fCharge=charge;
 148 }
 149 ///////////////////////////////////////////////////////////////////////////
 150 void IceLinefit::Exec(Option_t* opt)
 151 {
 152 // Implementation of the linefit reconstruction.
 153
 154  TString name=opt;
 155  AliJob* parent=(AliJob*)(gROOT->GetListOfTasks()->FindObject(name.Data()));
 156
 157  if (!parent) return;
 158
 159  IceEvent* evt=(IceEvent*)parent->GetObject("IceEvent");
 160  if (!evt) return;
 161
 162  // Enter the reco parameters as a device in the event
 163  AliSignal params;
 164  params.SetNameTitle("IceLinefit","IceLinefit reco parameters");
 165  params.SetSlotName("MaxmodA",1);
 166  params.SetSlotName("MinmodA",2);
 167  params.SetSlotName("MaxhitsA",3);
 168
 169  params.SetSignal(fMaxmodA,1);
 170  params.SetSignal(fMinmodA,2);
 171  params.SetSignal(fMaxhitsA,3);
 172
 173  evt->AddDevice(params);
 174
 175  if (fMaxhitsA<0) return;
 176
 177  // Fetch all fired Amanda OMs for this event
 178  TObjArray* aoms=evt->GetDevices("IceAOM");
 179  if (!aoms) return;
 180  Int_t naoms=aoms->GetEntries();
 181  if (!naoms) return;
 182
 183  // Check for the minimum and/or maximum number of good fired Amanda OMs
 184  Int_t ngood=0;
 185  for (Int_t iom=0; iom<naoms; iom++)
 186  {
 187   IceGOM* omx=(IceGOM*)aoms->At(iom);
 188   if (!omx) continue;
 189   if (omx->GetDeadValue("ADC") || omx->GetDeadValue("LE") || omx->GetDeadValue("TOT")) continue;
 190   ngood++;
 191  }
 192  if (ngood<fMinmodA || ngood>fMaxmodA) return;
 193
 194  const Float_t c=0.299792; // Light speed in vacuum in meters per ns
 195
 196  AliSignal* sx=0;
 197  Ali3Vector rom,sumr;
 198  Ali3Vector rt,sumrt;
 199  Float_t thit;
 200  Float_t sumt=0,sumt2=0;
 201  TObjArray hits;
 202  TObjArray* ordered;
 203  Int_t nh;
 204
 205  // Loop over all OMs and hits to determine the linefit parameters.
 206  // Also all the used hits are recorded for association with the track.
 207  for (Int_t iom=0; iom<naoms; iom++)
 208  {
 209   IceGOM* omx=(IceGOM*)aoms->At(iom);
 210   if (!omx) continue;
 211   if (omx->GetDeadValue("LE")) continue;
 212   rom=(Ali3Vector)omx->GetPosition();
 213   // Use the specified good hits of this OM
 214   ordered=0;
 215   if (fMaxhitsA>0 && omx->GetNhits()>fMaxhitsA) ordered=omx->SortHits("LE",1,0,7);
 216   nh=0;
 217   for (Int_t ih=1; ih<=omx->GetNhits(); ih++)
 218   {
 219    if (ordered)
 220    {
 221     if (nh>=fMaxhitsA) break;
 222     sx=(AliSignal*)ordered->At(ih-1);
 223    }
 224    else
 225    {
 226     sx=omx->GetHit(ih);
 227    }
 228    if (!sx) continue;
 229    if (sx->GetDeadValue("ADC") || sx->GetDeadValue("LE") || sx->GetDeadValue("TOT")) continue;
 230
 231    thit=sx->GetSignal("LE",7);
 232    rt=rom*thit;
 233    sumr+=rom;
 234    sumrt+=rt;
 235    sumt+=thit;
 236    sumt2+=thit*thit;
 237
 238    // Record this hit for association with the track
 239    hits.Add(sx);
 240    nh++;
 241   }
 242  }
 243
 244  Int_t nused=hits.GetEntries();
 245  if (!nused) return;
 246
 247  sumr/=float(nused);
 248  sumrt/=float(nused);
 249  sumt/=float(nused);
 250  sumt2/=float(nused);
 251
 252  Ali3Vector v;
 253  Ali3Vector temp;
 254  temp=sumr*sumt;
 255  v=sumrt-temp;
 256  Float_t dum=sumt2-(sumt*sumt);
 257  if (dum) v/=dum;
 258
 259  Float_t beta=v.GetNorm()/c;
 260  AliSignal fitstats;
 261  fitstats.SetNameTitle("Fitstats","Fit stats for IceLinefit");
 262  fitstats.SetSlotName("Beta",1);
 263  fitstats.SetSignal(beta,1);
 264
 265  Ali3Vector r;
 266  temp=v*sumt;
 267  r=sumr-temp;
 268
 269  AliTrack t;
 270  t.SetNameTitle(fTrackname.Data(),"IceLinefit linefit track");
 271  t.SetCharge(fCharge);
 272  evt->AddTrack(t);
 273  AliTrack* trk=evt->GetTrack(evt->GetNtracks());
 274  if (!trk) return;
 275
 276  trk->SetId(evt->GetNtracks(1)+1);
 277
 278  Ali3Vector p;
 279  Float_t vec[3];
 280  v.GetVector(vec,"sph");
 281  vec[0]=1;
 282  p.SetVector(vec,"sph");
 283
 284  AliPosition r0;
 285  r0.SetPosition(r);
 286  r0.SetTimestamp((AliTimestamp&)*evt);
 287  AliTimestamp* t0=r0.GetTimestamp();
 288  t0->Add(0,0,(int)sumt);
 289
 290  trk->Set3Momentum(p);
 291  trk->SetReferencePoint(r0);
 292  trk->SetTimestamp(*t0);
 293  trk->SetFitDetails(fitstats);
 294
 295  // Link the used hits to the track (and vice versa)
 296  for (Int_t i=0; i<nused; i++)
 297  {
 298   sx=(AliSignal*)hits.At(i);
 299   if (sx) sx->AddTrack(*trk);
 300  }
 301 }
 302 ///////////////////////////////////////////////////////////////////////////