RALICE/icepack/IcePandel.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 IcePandel
  20 // TTask derived class to perform Pandel fitting.
  21 //
  22 // The code in this processor is based on the algorithms as developed
  23 // by Oladipo Fadiran and George Japaridze (Clark Atlanta University, USA).
  24 //
  25 // Use the UseTracks memberfunction to specify the first guess tracks
  26 // to be processed by the minimiser.
  27 // By default only the first encountered IceDwalk track will be processed.
  28 //
  29 // Use the SelectHits memberfunction to specify the hits to be used.
  30 // By default only the hits associated to the first guess track are used.
  31 //
  32 // The fit processor (TFitter which is basically Minuit) printlevel
  33 // can be selected via the memberfunction SetPrintLevel.
  34 // By default all printout is suppressed (i.e. level=-2).
  35 //
  36 // An example of how to invoke this processor after Xtalk, hit cleaning
  37 // and a direct walk first guess estimate can be found in the ROOT example
  38 // macro icepandel.cc which resides in the /macros subdirectory.
  39 //
  40 // The minimisation results are stored in the IceEvent structure as
  41 // tracks with as default the name "IcePandel" (just like the first guess
  42 // results of e.g. IceDwalk).
  43 // This track name identifier can be modified by the user via the
  44 // SetTrackName() memberfunction. This will allow unique identification
  45 // of tracks which are produced when re-processing existing data with
  46 // different criteria.
  47 // By default the charge of the produced tracks is set to 0, since
  48 // no distinction can be made between positive or negative tracks.
  49 // However, the user can define the track charge by invokation
  50 // of the memberfunction SetCharge().
  51 // This facility may be used to distinguish tracks produced by the
  52 // various reconstruction algorithms in a (3D) colour display
  53 // (see the class AliHelix for further details).
  54 // A pointer to the first guess track which was used as input is available
  55 // via the GetParentTrack facility of these "IcePandel" tracks.
  56 // Furthermore, all the hits that were used in the minisation are available
  57 // via the GetSignal facility of a certain track.
  58 // The statistics of the TFitter result are stored as an AliSignal object
  59 // in the track, which can be obtained via the GetFitDetails memberfunction.
  60 // Currently no overall probability has yet been defined to indicate the
  61 // quality of a certain fit result. So, for the time being the user has
  62 // to judge the fit quality his/herself by means of the various TFitter stats
  63 // as accessible via the GetFitDetails facility.
  64 //
  65 // An example of how the various data can be accessed is given below,
  66 // where "evt" indicates the pointer to the IceEvent structure.
  67 //
  68 // Example for accessing data :
  69 // ----------------------------
  70 // TObjArray* tracks=evt->GetTracks("IcePandel");
  71 // if (!tracks) return;
  72 // AliPosition* r0=0;
  73 // Float_t fcn=0;
  74 // for (Int_t jtk=0; jtk<tracks->GetEntries(); jtk++)
  75 // {
  76 //  AliTrack* tx=(AliTrack*)tracks->At(jtk);
  77 //  if (!tx) continue;
  78 //  tx->Data("sph");
  79 //  r0=tx->GetReferencePoint();
  80 //  if (r0) r0->Data();
  81 //  sx=(AliSignal*)tx->GetFitDetails();
  82 //  if (sx) fcn=sx->GetSignal("FCN");
  83 //  AliTrack* tx2=tx->GetParentTrack();
  84 //  if (!tx2) continue;
  85 //  tx2->Data("sph");
  86 //  r0=tx2->GetReferencePoint();
  87 //  if (r0) r0->Data();
  88 // }
  89 //
  90 // Notes :
  91 // -------
  92 // 1) This processor only works properly on data which are Time and ADC
  93 //    calibrated and contain tracks from first guess algorithms like
  94 //    e.g. IceDwalk.
  95 // 2) In view of the usage of TFitter/Minuit minimisation, a global pointer
  96 //    to the instance of this class (gIcePandel) and a global static
  97 //    wrapper function (IcePandelFCN) have been introduced, to allow the
  98 //    actual minimisation to be performed via the memberfunction FitFCN.
  99 //    This implies that in a certain processing job only 1 instance of
 100 //    this IcePandel class may occur.
 101 //
 102 //--- Author: Nick van Eijndhoven 09-feb-2006 Utrecht University
 103 //- Modified: NvE $Date$ Utrecht University
 104 ///////////////////////////////////////////////////////////////////////////
 105
 106 #include "IcePandel.h"
 107 #include "Riostream.h"
 108
 109 // Global pointer to the instance of this object
 110  IcePandel* gIcePandel=0;
 111
 112 // TFitter/Minuit interface to IcePandel::FitFCN
 113  void IcePandelFCN(Int_t& npar,Double_t* gin,Double_t& f,Double_t* u,Int_t flag)
 114  {
 115   if (gIcePandel) gIcePandel->FitFCN(npar,gin,f,u,flag);
 116  }
 117
 118 ClassImp(IcePandel) // Class implementation to enable ROOT I/O
 119
 120 IcePandel::IcePandel(const char* name,const char* title) : TTask(name,title)
 121 {
 122 // Default constructor.
 123  fFirst=1;
 124  fPrint=-2;
 125  fSelhits=1;
 126  fEvt=0;
 127  fUseNames=0;
 128  fUseNtk=0;
 129  fTrack=0;
 130  fHits=0;
 131  fFitter=0;
 132  fTrackname="IcePandel";
 133  fCharge=0;
 134
 135  // Set the global pointer to this instance
 136  gIcePandel=this;
 137 }
 138 ///////////////////////////////////////////////////////////////////////////
 139 IcePandel::~IcePandel()
 140 {
 141 // Default destructor.
 142  if (fUseNames)
 143  {
 144   delete fUseNames;
 145   fUseNames=0;
 146  }
 147  if (fUseNtk)
 148  {
 149   delete fUseNtk;
 150   fUseNtk=0;
 151  }
 152  if (fHits)
 153  {
 154   delete fHits;
 155   fHits=0;
 156  }
 157  if (fFitter)
 158  {
 159   delete fFitter;
 160   fFitter=0;
 161  }
 162 }
 163 ///////////////////////////////////////////////////////////////////////////
 164 void IcePandel::Exec(Option_t* opt)
 165 {
 166 // Implementation of the hit fitting procedure.
 167
 168  TString name=opt;
 169  AliJob* parent=(AliJob*)(gROOT->GetListOfTasks()->FindObject(name.Data()));
 170
 171  if (!parent) return;
 172
 173  fEvt=(IceEvent*)parent->GetObject("IceEvent");
 174  if (!fEvt) return;
 175
 176  if (!fUseNames) UseTracks("IceDwalk",1);
 177
 178  Int_t nclasses=fUseNames->GetEntries(); // Number of first guess classes to be processed
 179  Int_t ntkmax=0; // Max. number of tracks for a certain class
 180  TObjString* strx=0;
 181  TString str;
 182
 183  if (fFirst)
 184  {
 185   cout << " *IcePandel* First guess selections to be processed (-1=all)." << endl;
 186   for (Int_t i=0; i<nclasses; i++)
 187   {
 188    strx=(TObjString*)fUseNames->At(i);
 189    if (!strx) continue;
 190    str=strx->GetString();
 191    ntkmax=fUseNtk->At(i);
 192    cout << " Maximally " << ntkmax << " track(s) per event for procedure : " << str.Data() << endl;
 193   }
 194   cout << " *IcePandel* Hit selection mode : " << fSelhits << endl;
 195   cout << endl;
 196   fFirst=0;
 197  }
 198
 199  Double_t pi=acos(-1.);
 200
 201  // Initialisation of the minimisation processor
 202  Double_t arglist[100];
 203  if (!fFitter) fFitter=new TFitter();
 204
 205  // The number of reconstructed tracks already present in the event
 206  Int_t ntkreco=fEvt->GetNtracks(1);
 207
 208  if (!fHits)
 209  {
 210   fHits=new TObjArray();
 211  }
 212  else
 213  {
 214   fHits->Clear();
 215  }
 216
 217  // If selected, use all the good quality hits of the complete event
 218  if (fSelhits==0)
 219  {
 220   TObjArray* hits=fEvt->GetHits("IceGOM");
 221   for (Int_t ih=0; ih<hits->GetEntries(); ih++)
 222   {
 223    AliSignal* sx=(AliSignal*)hits->At(ih);
 224    if (!sx) continue;
 225    if (sx->GetDeadValue("ADC") || sx->GetDeadValue("LE") || sx->GetDeadValue("TOT")) continue;
 226    fHits->Add(sx);
 227   }
 228  }
 229
 230  // Track by track processing of the selected first guess classes
 231  Float_t vec[3];
 232  Float_t err[3];
 233  Float_t margin,minmarg=25; // (minimal) margin for the fitter r0 search area
 234  Ali3Vector p;
 235  AliPosition pos;
 236  AliTrack tkfit;
 237  tkfit.SetNameTitle(fTrackname.Data(),"IcePandel fit result");
 238  AliSignal fitstats;
 239  fitstats.SetNameTitle("Fitstats","TFitter stats for Pandel fit");
 240  fitstats.SetSlotName("IER",1);
 241  fitstats.SetSlotName("FCN",2);
 242  fitstats.SetSlotName("EDM",3);
 243  fitstats.SetSlotName("NVARS",4);
 244  Float_t x,y,z,theta,phi;
 245  Float_t xmin,xmax,ymin,ymax,zmin,zmax,thetamin,thetamax,phimin,phimax;
 246  Double_t amin,edm,errdef; // Minimisation stats
 247  Int_t ier,nvpar,nparx;    // Minimisation stats
 248  Int_t ntk=0;
 249  Int_t nsig=0;
 250  for (Int_t iclass=0; iclass<nclasses; iclass++) // Loop over first guess classes
 251  {
 252   strx=(TObjString*)fUseNames->At(iclass);
 253   if (!strx) continue;
 254   str=strx->GetString();
 255   ntkmax=fUseNtk->At(iclass);
 256   TObjArray* tracks=fEvt->GetTracks(str);
 257   ntk=tracks->GetEntries();
 258   if (ntkmax>0 && ntk>ntkmax) ntk=ntkmax;
 259
 260   for (Int_t jtk=0; jtk<ntk; jtk++) // Loop over tracks of a certain class
 261   {
 262    fTrack=(AliTrack*)tracks->At(jtk);
 263    if (!fTrack) continue;
 264
 265    AliPosition* r0=fTrack->GetReferencePoint();
 266    if (!r0) continue;
 267
 268    // If selected, use only the first guess track associated hits
 269    if (fSelhits==1)
 270    {
 271     fHits->Clear();
 272     nsig=fTrack->GetNsignals();
 273     for (Int_t is=1; is<=nsig; is++)
 274     {
 275      AliSignal* sx=fTrack->GetSignal(is);
 276      if (!sx) continue;
 277      if (!sx->GetDevice()->InheritsFrom("IceGOM")) continue;
 278      if (sx->GetDeadValue("ADC") || sx->GetDeadValue("LE") || sx->GetDeadValue("TOT")) continue;
 279      fHits->Add(sx);
 280     }
 281    }
 282
 283    if (!fHits->GetEntries()) continue;
 284
 285    r0->GetVector(vec,"car");
 286    r0->GetErrors(err,"car");
 287
 288    x=vec[0];
 289    y=vec[1];
 290    z=vec[2];
 291    margin=5.*err[0];
 292    if (margin<minmarg) margin=minmarg;
 293    xmin=x-margin;
 294    xmax=x+margin;
 295    margin=5.*err[1];
 296    if (margin<minmarg) margin=minmarg;
 297    ymin=y-margin;
 298    ymax=y+margin;
 299    margin=5.*err[2];
 300    if (margin<minmarg) margin=minmarg;
 301    zmin=z-margin;
 302    zmax=z+margin;
 303
 304    p=fTrack->Get3Momentum();
 305    p.GetVector(vec,"sph");
 306
 307    theta=vec[1];
 308    phi=vec[2];
 309    thetamin=theta-(pi/3.);
 310    if (thetamin<0) thetamin=0;
 311    thetamax=theta+(pi/3.);
 312    if (thetamax>pi) thetamax=pi;
 313    phimin=phi-(pi/2.);
 314    phimax=phi+(pi/2.);
 315
 316    // Process this first guess track with its associated hits
 317    fFitter->Clear();
 318
 319    arglist[0]=fPrint;
 320    if (fPrint==-2) arglist[0]=-1;
 321    fFitter->ExecuteCommand("SET PRINT",arglist,1);
 322    if (fPrint==-2) fFitter->ExecuteCommand("SET NOWARNINGS",arglist,0);
 323
 324    fFitter->SetFitMethod("loglikelihood");
 325
 326    fFitter->SetParameter(0,"r0x",x,0.001,xmin,xmax);
 327    fFitter->SetParameter(1,"r0y",y,0.001,ymin,ymax);
 328    fFitter->SetParameter(2,"r0z",z,0.001,zmin,zmax);
 329    fFitter->SetParameter(3,"theta",theta,0.001,thetamin,thetamax);
 330    fFitter->SetParameter(4,"phi",phi,0.001,phimin,phimax);
 331
 332    fFitter->SetFCN(IcePandelFCN);
 333
 334    arglist[0]=0;
 335    ier=fFitter->ExecuteCommand("MINIMIZE",arglist,0);
 336
 337    fFitter->GetStats(amin,edm,errdef,nvpar,nparx);
 338    fitstats.Reset();
 339    fitstats.SetSignal(ier,1);
 340    fitstats.SetSignal(amin,2);
 341    fitstats.SetSignal(edm,3);
 342    fitstats.SetSignal(nvpar,4);
 343
 344    // Resulting parameters after minimisation
 345    vec[0]=fFitter->GetParameter(0);
 346    vec[1]=fFitter->GetParameter(1);
 347    vec[2]=fFitter->GetParameter(2);
 348    err[0]=fFitter->GetParError(0);
 349    err[1]=fFitter->GetParError(1);
 350    err[2]=fFitter->GetParError(2);
 351    pos.SetPosition(vec,"car");
 352    pos.SetPositionErrors(err,"car");
 353
 354    vec[0]=1.;
 355    vec[1]=fFitter->GetParameter(3);
 356    vec[2]=fFitter->GetParameter(4);
 357    err[0]=0.;
 358    err[1]=fFitter->GetParError(3);
 359    err[2]=fFitter->GetParError(4);
 360    p.SetVector(vec,"sph");
 361    p.SetErrors(err,"sph");
 362
 363    // Enter the fit result as a track in the event structure
 364    tkfit.Reset();
 365    ntkreco++;
 366    tkfit.SetId(ntkreco);
 367    tkfit.SetCharge(fCharge);
 368    tkfit.SetParentTrack(fTrack);
 369    AliTimestamp* tt0=r0->GetTimestamp();
 370    pos.SetTimestamp(*tt0);
 371    tkfit.SetTimestamp(*tt0);
 372    tkfit.SetReferencePoint(pos);
 373    tkfit.Set3Momentum(p);
 374    tkfit.SetFitDetails(fitstats);
 375    for (Int_t ihit=0; ihit<fHits->GetEntries(); ihit++)
 376    {
 377     AliSignal* sx=(AliSignal*)fHits->At(ihit);
 378     if (sx) tkfit.AddSignal(*sx);
 379    }
 380    fEvt->AddTrack(tkfit);
 381   } // End loop over tracks
 382  } // End loop over first guess classes
 383
 384 }
 385 ///////////////////////////////////////////////////////////////////////////
 386 void IcePandel::SetPrintLevel(Int_t level)
 387 {
 388 // Set the fitter (Minuit) print level.
 389 // See the TFitter and TMinuit docs for details.
 390 //
 391 // Note : level=-2 suppresses also all fit processor warnings.
 392 //
 393 // The default in the constructor is level=-2.
 394
 395  fPrint=level;
 396 }
 397 ///////////////////////////////////////////////////////////////////////////
 398 void IcePandel::UseTracks(TString classname,Int_t n)
 399 {
 400 // Specification of the first guess tracks to be used.
 401 //
 402 // classname : Specifies the first guess algorithm (e.g. "IceDwalk");
 403 // n : Specifies the max. number of these tracks to be used
 404 //
 405 // Note : n<0 will use all the existing tracks of the specified classname
 406 //
 407 // The default is n=-1.
 408 //
 409 // Consecutive invokations of this memberfunction with different classnames
 410 // will result in an incremental effect.
 411 //
 412 // Example :
 413 // ---------
 414 // UseTracks("IceDwalk",5);
 415 // UseTracks("IceLfit",2);
 416 // UseTracks("IceJams");
 417 //
 418 // This will use the first 5 IceDwalk, the first 2 IceLfit and all the
 419 // IceJams tracks which are encountered in the event structure.
 420
 421  if (!fUseNames)
 422  {
 423   fUseNames=new TObjArray();
 424   fUseNames->SetOwner();
 425  }
 426
 427  if (!fUseNtk) fUseNtk=new TArrayI();
 428
 429  // Check if this classname has already been specified before
 430  TString s;
 431  Int_t nen=fUseNames->GetEntries();
 432  for (Int_t i=0; i<nen; i++)
 433  {
 434   TObjString* sx=(TObjString*)fUseNames->At(i);
 435   if (!sx) continue;
 436   s=sx->GetString();
 437   if (s==classname) return;
 438  }
 439
 440  // New classname to be added into the storage
 441  if (nen >= fUseNames->GetSize()) fUseNames->Expand(nen+1);
 442  if (nen >= fUseNtk->GetSize()) fUseNtk->Set(nen+1);
 443
 444  TObjString* name=new TObjString();
 445  name->SetString(classname);
 446  fUseNames->Add(name);
 447  fUseNtk->AddAt(n,nen);
 448 }
 449 ///////////////////////////////////////////////////////////////////////////
 450 void IcePandel::SelectHits(Int_t mode)
 451 {
 452 // Specification of the hits to be used in the minimisation.
 453 //
 454 // mode = 0 : All hit cleaning survived hits of the complete event are used
 455 //        1 : Only the associated hits are used for each first guess track
 456 //
 457 // The default is mode=1.
 458
 459  if (mode==0 || mode==1) fSelhits=mode;
 460 }
 461 ///////////////////////////////////////////////////////////////////////////
 462 void IcePandel::SetTrackName(TString s)
 463 {
 464 // Set (alternative) name identifier for the produced tracks.
 465 // This allows unique identification of (newly) produced pandel tracks
 466 // in case of re-processing of existing data with different criteria.
 467 // By default the produced tracks have the name "IcePandel" which is
 468 // set in the constructor of this class.
 469  fTrackname=s;
 470 }
 471 ///////////////////////////////////////////////////////////////////////////
 472 void IcePandel::SetCharge(Float_t charge)
 473 {
 474 // Set user defined charge for the produced tracks.
 475 // This allows identification of these tracks on color displays.
 476 // By default the produced tracks have charge=0 which is set in the
 477 // constructor of this class.
 478  fCharge=charge;
 479 }
 480 ///////////////////////////////////////////////////////////////////////////
 481 void IcePandel::FitFCN(Int_t&,Double_t*,Double_t& f,Double_t* x,Int_t)
 482 {
 483 // The Pandel function used for the minimisation process.
 484
 485  const Float_t c=0.299792;           // Light speed in vacuum in meters per ns
 486  const Float_t nice=1.35634;         // Refractive index of ice
 487  const Float_t thetac=acos(1./nice); // Cherenkov angle (in radians)
 488  const Float_t lambda=33.3;          // Light scattering length in ice
 489  const Float_t labs=98;              // Light absorbtion length in ice
 490  const Float_t cice=c/nice;          // Light speed in ice in meters per ns
 491  const Float_t tau=557;
 492  const Double_t rho=((1./tau)+(cice/labs));
 493
 494  f=0;
 495
 496  // The first original guess track parameters
 497  AliPosition* tr0=fTrack->GetReferencePoint();
 498  AliTimestamp* tt0=tr0->GetTimestamp();
 499  Float_t t0=fEvt->GetDifference(tt0,"ns");
 500  Ali3Vector tp=fTrack->Get3Momentum();
 501
 502  // The new r0 and p vectors from the minimisation
 503  Float_t vec[3];
 504
 505  AliPosition r0;
 506  vec[0]=x[0];
 507  vec[1]=x[1];
 508  vec[2]=x[2];
 509  r0.SetPosition(vec,"car");
 510
 511  Ali3Vector p;
 512  vec[0]=1;
 513  vec[1]=x[3];
 514  vec[2]=x[4];
 515  p.SetVector(vec,"sph");
 516
 517  Int_t nhits=fHits->GetEntries();
 518  AliPosition rhit;
 519  Ali3Vector r12;
 520  Float_t d,dist,thit,tgeo;
 521  Double_t tres,zeta,pandel;
 522  for (Int_t i=0; i<nhits; i++)
 523  {
 524   AliSignal* sx=(AliSignal*)fHits->At(i);
 525   if (!sx) continue;
 526   IceGOM* omx=(IceGOM*)sx->GetDevice();
 527   if (!omx) continue;
 528   rhit=omx->GetPosition();
 529   d=fTrack->GetDistance(rhit);
 530   zeta=d/lambda;
 531   d*=sin(thetac);
 532   r12=rhit-r0;
 533   dist=p.Dot(r12)+d*tan(thetac);
 534   tgeo=t0+dist/c;
 535   thit=sx->GetSignal("LE",7);
 536   tres=thit-tgeo;
 537
 538   // The Pandel function evaluation
 539   // Avoid minimiser problem for infinite derivative on Pandel surface
 540   // This problem will be absent when using a smooth convoluted Pandel
 541   if (tres>0 && zeta>=1)
 542   {
 543    pandel=pow(rho,zeta)*pow(tres,(zeta-1.))*exp(-rho*tres)/TMath::Gamma(zeta);
 544   }
 545   else
 546   {
 547    pandel=1.e-6;
 548   }
 549
 550   // Use 10*log10 expression to obtain intuitive decibel scale
 551   f-=10.*log10(pandel);
 552  }
 553 }
 554 ///////////////////////////////////////////////////////////////////////////