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  Int_t naoms=aoms->GetEntries();
 180  if (!naoms) return;
 181
 182  // Check for the minimum and/or maximum number of good fired Amanda OMs
 183  Int_t ngood=0;
 184  for (Int_t iom=0; iom<naoms; iom++)
 185  {
 186   IceGOM* omx=(IceGOM*)aoms->At(iom);
 187   if (!omx) continue;
 188   if (omx->GetDeadValue("ADC") || omx->GetDeadValue("LE") || omx->GetDeadValue("TOT")) continue;
 189   ngood++;
 190  }
 191  if (ngood<fMinmodA || ngood>fMaxmodA) return;
 192
 193  const Float_t c=0.299792; // Light speed in vacuum in meters per ns
 194
 195  AliSignal* sx=0;
 196  Ali3Vector rom,sumr;
 197  Ali3Vector rt,sumrt;
 198  Float_t thit;
 199  Float_t sumt=0,sumt2=0;
 200  TObjArray hits;
 201  TObjArray* ordered;
 202  Int_t nh;
 203
 204  // Loop over all OMs and hits to determine the linefit parameters.
 205  // Also all the used hits are recorded for association with the track.
 206  for (Int_t iom=0; iom<naoms; iom++)
 207  {
 208   IceGOM* omx=(IceGOM*)aoms->At(iom);
 209   if (!omx) continue;
 210   if (omx->GetDeadValue("LE")) continue;
 211   rom=(Ali3Vector)omx->GetPosition();
 212   // Use the specified good hits of this OM
 213   ordered=0;
 214   if (fMaxhitsA>0 && omx->GetNhits()>fMaxhitsA) ordered=omx->SortHits("LE",1,0,7);
 215   nh=0;
 216   for (Int_t ih=1; ih<=omx->GetNhits(); ih++)
 217   {
 218    if (ordered)
 219    {
 220     if (nh>=fMaxhitsA) break;
 221     sx=(AliSignal*)ordered->At(ih-1);
 222    }
 223    else
 224    {
 225     sx=omx->GetHit(ih);
 226    }
 227    if (!sx) continue;
 228    if (sx->GetDeadValue("ADC") || sx->GetDeadValue("LE") || sx->GetDeadValue("TOT")) continue;
 229
 230    thit=sx->GetSignal("LE",7);
 231    rt=rom*thit;
 232    sumr+=rom;
 233    sumrt+=rt;
 234    sumt+=thit;
 235    sumt2+=thit*thit;
 236
 237    // Record this hit for association with the track
 238    hits.Add(sx);
 239    nh++;
 240   }
 241  }
 242
 243  Int_t nused=hits.GetEntries();
 244  if (!nused) return;
 245
 246  sumr/=float(nused);
 247  sumrt/=float(nused);
 248  sumt/=float(nused);
 249  sumt2/=float(nused);
 250
 251  Ali3Vector v;
 252  Ali3Vector temp;
 253  temp=sumr*sumt;
 254  v=sumrt-temp;
 255  Float_t dum=sumt2-(sumt*sumt);
 256  if (dum) v/=dum;
 257
 258  Float_t beta=v.GetNorm()/c;
 259  AliSignal fitstats;
 260  fitstats.SetNameTitle("Fitstats","Fit stats for IceLinefit");
 261  fitstats.SetSlotName("Beta",1);
 262  fitstats.SetSignal(beta,1);
 263
 264  Ali3Vector r;
 265  temp=v*sumt;
 266  r=sumr-temp;
 267
 268  AliTrack t;
 269  t.SetNameTitle(fTrackname.Data(),"IceLinefit linefit track");
 270  t.SetCharge(fCharge);
 271  evt->AddTrack(t);
 272  AliTrack* trk=evt->GetTrack(evt->GetNtracks());
 273  if (!trk) return;
 274
 275  trk->SetId(evt->GetNtracks(1)+1);
 276
 277  Ali3Vector p;
 278  Float_t vec[3];
 279  v.GetVector(vec,"sph");
 280  vec[0]=1;
 281  p.SetVector(vec,"sph");
 282
 283  AliPosition r0;
 284  r0.SetPosition(r);
 285  r0.SetTimestamp((AliTimestamp&)*evt);
 286  AliTimestamp* t0=r0.GetTimestamp();
 287  t0->Add(0,0,(int)sumt);
 288
 289  trk->Set3Momentum(p);
 290  trk->SetReferencePoint(r0);
 291  trk->SetTimestamp(*t0);
 292  trk->SetFitDetails(fitstats);
 293
 294  // Link the used hits to the track (and vice versa)
 295  for (Int_t i=0; i<nused; i++)
 296  {
 297   sx=(AliSignal*)hits.At(i);
 298   if (sx) sx->AddTrack(*trk);
 299  }
 300 }
 301 ///////////////////////////////////////////////////////////////////////////