RALICE/icepack/iceconvert/IceF2k.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 IceF2k
  20 // Conversion of Amanda F2K data into IceEvent physics event structures.
  21 //
  22 // Usage example :
  23 // ---------------
  24 //
  25 // gSystem->Load("ralice");
  26 // gSystem->Load("icepack");
  27 // gSystem->Load("iceconvert");
  28 //
  29 // // Output file for the event structures
  30 // TFile* ofile=new TFile("events.root","RECREATE","F2K data in IceEvent structure");
  31 // TTree* otree=new TTree("T","Data of an Amanda run");
  32 //
  33 // // Limit the number of entries for testing
  34 // Int_t nentries=300;
  35 //
  36 // // Print frequency to produce a short summary print every printfreq events
  37 // Int_t printfreq=10;
  38 //
  39 // // Split level for the output structures
  40 // Int_t split=2;
  41 //
  42 // // Buffer size for the output structures
  43 // Int_t bsize=32000;
  44 //
  45 // IceF2k q("run8000.f2k",split,bsize);
  46 // q.Loop(otree,nentries,printfreq);
  47 //
  48 // // Select various objects to be added to the output file
  49 //
  50 // AliObjMatrix* omdb=q.GetOMdbase();
  51 // if (omdb) omdb->Write();
  52 //
  53 // AliDevice* fitdefs=q.GetFitdefs();
  54 // if (fitdefs) fitdefs->Write();
  55 //
  56 // TDatabasePDG* pdg=q.GetPDG();
  57 // if (pdg) pdg->Write();
  58 //
  59 // // Close output file
  60 // ofile->Write();
  61 // ofile->Close();
  62 //
  63 //--- Author: Nick van Eijndhoven 11-mar-2005 Utrecht University
  64 //- Modified: NvE $Date$ Utrecht University
  65 ///////////////////////////////////////////////////////////////////////////
  66
  67 #include "IceF2k.h"
  68 #include "Riostream.h"
  69
  70 ClassImp(IceF2k) // Class implementation to enable ROOT I/O
  71
  72 IceF2k::IceF2k(char* fname,Int_t split,Int_t bsize)
  73 {
  74 // Default constructor.
  75 // Initialise the input file and data structres to be converted.
  76 // Also the required split level and buffer size of the output tree
  77 // can be specified in this constructor.
  78 // By default tree=0, split=0 and bsize=32000.
  79
  80  fSplit=split;
  81  fBsize=bsize;
  82
  83  fPdg=0;
  84  fOmdb=0;
  85  fFitdefs=0;
  86
  87  if (!fname)
  88  {
  89   cout << " *IceF2k ctor* No data input file specified." << endl;
  90   return;
  91  }
  92
  93  // Open the input file in the default ascii format (autodetection) for reading
  94  fInput=rdmc_mcopen(fname,"r",RDMC_DEFAULT_ASCII_F);
  95
  96  if (!fInput)
  97  {
  98   cout << " *IceF2k ctor* No input file found with name : " << fname << endl;
  99   return;
 100  }
 101
 102  // Initialise the event structure
 103  rdmc_init_mevt(&fEvent);
 104
 105  // Read the file header information
 106  rdmc_rarr(fInput,&fHeader);
 107 }
 108 ///////////////////////////////////////////////////////////////////////////
 109 IceF2k::~IceF2k()
 110 {
 111 // Default destructor.
 112  if (fPdg)
 113  {
 114   delete fPdg;
 115   fPdg=0;
 116  }
 117
 118  if (fOmdb)
 119  {
 120   delete fOmdb;
 121   fOmdb=0;
 122  }
 123
 124  if (fFitdefs)
 125  {
 126   delete fFitdefs;
 127   fFitdefs=0;
 128  }
 129 }
 130 ///////////////////////////////////////////////////////////////////////////
 131 TDatabasePDG* IceF2k::GetPDG()
 132 {
 133 // Provide pointer to the PDG database
 134  return fPdg;
 135 }
 136 ///////////////////////////////////////////////////////////////////////////
 137 AliObjMatrix* IceF2k::GetOMdbase()
 138 {
 139 // Provide pointer to the OM geometry, calib. etc... database
 140  return fOmdb;
 141 }
 142 ///////////////////////////////////////////////////////////////////////////
 143 AliDevice* IceF2k::GetFitdefs()
 144 {
 145 // Provide pointer to the fit definitions
 146  return fFitdefs;
 147 }
 148 ///////////////////////////////////////////////////////////////////////////
 149 void IceF2k::Loop(TTree* otree,Int_t nentries,Int_t printfreq)
 150 {
 151 // Loop over the specified number of entries and convert the
 152 // F2K data into the IceEvent structure.
 153 // The output will be written on the output tree specified as "otree".
 154 // If otree=0, a default standard output tree will be created.
 155 // If nentries<0 (default) all the entries of the input file
 156 // will be processed.
 157 // Every "printfreq" events a short event summary will be printed.
 158 // The default value is printfreq=1.
 159
 160  if (!fInput || fSplit<0) return;
 161
 162  if (!otree) otree=new TTree("T","F2K Data");
 163
 164  IceEvent* evt=new IceEvent();
 165
 166  evt->SetTrackCopy(1);
 167  evt->SetDevCopy(1);
 168
 169  // Branch in the tree for the event structure
 170  otree->Branch("IceEvent","IceEvent",&evt,fBsize,fSplit);
 171
 172  // Create the particle database and extend it with some F2000 specific definitions
 173  if (!fPdg) fPdg=new TDatabasePDG();
 174  Double_t me=fPdg->GetParticle(11)->Mass();
 175  fPdg->AddParticle("brems"   ,"brems"   ,0,1,0,0,"none",10001001,0,0);
 176  fPdg->AddParticle("deltae"  ,"deltae"  ,me,1,0,-3,"Lepton",10001002,0,0);
 177  fPdg->AddParticle("pairprod","pairprod",0,1,0,0,"none",10001003,0,0);
 178  fPdg->AddParticle("nucl_int","nucl_Int",0,1,0,0,"none",10001004,0,0);
 179  fPdg->AddParticle("mu_pair" ,"mu_pair" ,0,1,0,0,"none",10001005,0,0);
 180  fPdg->AddParticle("hadrons" ,"hadrons" ,0,1,0,0,"none",10001006,0,0);
 181  fPdg->AddParticle("fiberlaser","fiberlaser",0,1,0,0,"none",10002100,0,0);
 182  fPdg->AddParticle("n2laser"   ,"n2laser"   ,0,1,0,0,"none",10002101,0,0);
 183  fPdg->AddParticle("yaglaser"  ,"yaglaser"  ,0,1,0,0,"none",10002201,0,0);
 184  fPdg->AddParticle("z_primary","z_primary",0,1,0,0,"none",10003000,0,0);
 185  fPdg->AddParticle("a_primary","a_primary",0,1,0,0,"none",10003500,0,0);
 186
 187  // Fill the database with geometry, calib. etc... parameters
 188  // for all the devices
 189  FillOMdbase();
 190
 191  // Set the fit definitions according to the F2000 header info
 192  SetFitdefs();
 193
 194  for (Int_t jentry=0; jentry<nentries; jentry++)
 195  {
 196   if (rdmc_revt(fInput,&fHeader,&fEvent) != 0) break;
 197
 198   // Reset the complete Event structure
 199   evt->Reset();
 200
 201   evt->SetRunNumber(fEvent.nrun);
 202   evt->SetEventNumber(fEvent.enr);
 203   evt->SetMJD(fEvent.mjd,fEvent.secs,fEvent.nsecs);
 204
 205   PutMcTracks(evt);
 206
 207   PutRecoTracks(evt);
 208
 209   PutHits(evt);
 210
 211   if (!(jentry%printfreq))
 212   {
 213    evt->HeaderData();
 214   }
 215
 216   // Write the complete structure to the output Tree
 217   otree->Fill();
 218  }
 219
 220  if (evt) delete evt;
 221 }
 222 ///////////////////////////////////////////////////////////////////////////
 223 void IceF2k::FillOMdbase()
 224 {
 225 // Fill the database with geometry, calib. etc... parameters
 226 // for all the devices.
 227
 228  if (fHeader.nch<=0) return;
 229
 230  if (fOmdb)
 231  {
 232   fOmdb->Reset();
 233  }
 234  else
 235  {
 236   fOmdb=new AliObjMatrix();
 237   fOmdb->SetNameTitle("OMDBASE","The OM geometry, calib. etc... database");
 238   fOmdb->SetOwner();
 239  }
 240
 241  IceAOM* dev=0;
 242  Double_t pos[3]={0,0,0};
 243  for (Int_t i=0; i<fHeader.nch; i++)
 244  {
 245   dev=new IceAOM();
 246   dev->SetUniqueID(i+1);
 247   dev->SetSlotName("TYPE",1);
 248   dev->SetSlotName("ORIENT",2);
 249   dev->SetSlotName("T0",3);
 250   dev->SetSlotName("ALPHA",4);
 251   dev->SetSlotName("KADC",5);
 252   dev->SetSlotName("KTOT",6);
 253   dev->SetSlotName("KTDC",7);
 254
 255   pos[0]=fHeader.x[i];
 256   pos[1]=fHeader.y[i];
 257   pos[2]=fHeader.z[i];
 258   dev->SetPosition(pos,"car");
 259   dev->SetSignal(fHeader.type[i],1);
 260   dev->SetSignal((Float_t)fHeader.costh[i],2);
 261   dev->SetSignal(fHeader.cal[i].t_0,3);
 262   dev->SetSignal(fHeader.cal[i].alpha_t,4);
 263   dev->SetSignal(fHeader.cal[i].beta_a,5);
 264   dev->SetSignal(fHeader.cal[i].beta_tot,6);
 265   dev->SetSignal(fHeader.cal[i].beta_t,7);
 266   fOmdb->EnterObject(i+1,1,dev);
 267  }
 268 }
 269 ///////////////////////////////////////////////////////////////////////////
 270 void IceF2k::SetFitdefs()
 271 {
 272 // Obtain the names of the variables for each fit procedure from the
 273 // f2000 header. Each different fit procedure is then stored as a separate
 274 // hit of an AliDevice object and the various fit variables are stored
 275 // as separate signal slots of the corresponding hit.
 276 // Via the GetFitdefs() memberfunction this AliDevice object can be
 277 // retrieved and stored in the ROOT output file if wanted.
 278 // The name of the object is FitDefinitions and the stored data can be
 279 // inspected via the AliDevice::Data() memberfunction and looks as follows :
 280 //
 281 //  *AliDevice::Data* Id :0 Name : FitDefinitions
 282 //    Position Vector in car coordinates : 0 0 0
 283 //    Err. in car coordinates : 0 0 0
 284 //  The following 8 hits are registered :
 285 //  *AliSignal::Data* Id :0
 286 //    Position Vector in car coordinates : 0 0 0
 287 //    Err. in car coordinates : 0 0 0
 288 //    Owned by device : AliDevice Name : FitDefinitions
 289 //    Slot : 1 Signal value : 1 name : id
 290 //    Slot : 2 Signal value : 2 name : rchi2
 291 //    Slot : 3 Signal value : 3 name : prob
 292 //    Slot : 4 Signal value : 4 name : sigth
 293 //    Slot : 5 Signal value : 5 name : covmin
 294 //    Slot : 6 Signal value : 6 name : covmax
 295 //    Slot : 7 Signal value : 7 name : cutflag
 296 //    Slot : 8 Signal value : 8 name : chi2
 297 //  *AliSignal::Data* Id :1
 298 //    Position Vector in car coordinates : 0 0 0
 299 //    Err. in car coordinates : 0 0 0
 300 //    Owned by device : AliDevice Name : FitDefinitions
 301 //    Slot : 1 Signal value : 1 name : id
 302 //    Slot : 2 Signal value : 2 name : rchi2
 303 //    Slot : 3 Signal value : 3 name : prob
 304 // etc....
 305 //
 306 // This memberfunction is based on the original idea/code by Adam Bouchta.
 307
 308  if (fHeader.n_fit<=0) return;
 309
 310  if (fFitdefs)
 311  {
 312   fFitdefs->Reset(1);
 313  }
 314  else
 315  {
 316   fFitdefs=new AliDevice();
 317  }
 318
 319  fFitdefs->SetName("FitDefinitions");
 320  fFitdefs->SetHitCopy (1);
 321
 322  AliSignal s;
 323  s.Reset();
 324
 325  for (Int_t i=0; i<fHeader.n_fit; i++)
 326  {
 327   s.SetUniqueID(fHeader.def_fit[i].id);
 328
 329   for (Int_t j=0; j<fHeader.def_fit[i].nwords; j++)
 330   {
 331    s.SetSlotName(TString(fHeader.def_fit[i].words[j]),j+1);
 332    s.SetSignal(j+1,j+1);
 333   }
 334
 335   fFitdefs->AddHit(s);
 336   s.Reset(1);
 337  }
 338 }
 339 ///////////////////////////////////////////////////////////////////////////
 340 void IceF2k::PutMcTracks(IceEvent* evt)
 341 {
 342 // Get the MC tracks from the F2000 file into the IcePack structure.
 343 // Note : MC tracks are given negative track id's in the event structure.
 344 // This memberfunction is based on the original code by Adam Bouchta.
 345
 346  if (!evt || fEvent.ntrack<=0) return;
 347
 348  // Loop over all the tracks and add them to the current event
 349  AliTrack t;
 350  Double_t vec[3];
 351  AliPosition r;
 352  Ali3Vector p;
 353  Int_t tid=0;
 354  Int_t idpdg=0;
 355  Int_t idf2k=0;
 356  for (Int_t i=0; i<fEvent.ntrack; i++)
 357  {
 358   t.Reset ();
 359
 360   // Beginpoint of the track
 361   vec[0]=fEvent.gen[i].x;
 362   vec[1]=fEvent.gen[i].y;
 363   vec[2]=fEvent.gen[i].z;
 364   r.SetPosition(vec,"car");
 365   t.SetBeginPoint(r);
 366
 367   // Endpoint of the track
 368   vec[0]+=fEvent.gen[i].length*fEvent.gen[i].px;
 369   vec[1]+=fEvent.gen[i].length*fEvent.gen[i].py;
 370   vec[2]+=fEvent.gen[i].length*fEvent.gen[i].pz;
 371   r.SetPosition(vec,"car");
 372   t.SetEndPoint(r);
 373
 374   // Momentum in GeV/c
 375   vec[0]=fEvent.gen[i].e*fEvent.gen[i].px*1e-3;
 376   vec[1]=fEvent.gen[i].e*fEvent.gen[i].py*1e-3;
 377   vec[2]=fEvent.gen[i].e*fEvent.gen[i].pz*1e-3;
 378   p.SetVector (vec,"car");
 379   t.Set3Momentum(p);
 380
 381   // MC tracks are indicated by negative track id's
 382   tid=fEvent.gen[i].tag;
 383   t.SetId(-abs(tid));
 384
 385   idf2k=fEvent.gen[i].id;
 386   idpdg=0;
 387   if (idf2k>1000)
 388   {
 389    idpdg=idf2k+10000000;
 390   }
 391   else if (idf2k <= 48)
 392   {
 393    idpdg=fPdg->ConvertGeant3ToPdg(idf2k);
 394   }
 395   else
 396   {
 397    if (idf2k==201) idpdg=12;
 398    if (idf2k==202) idpdg=14;
 399    if (idf2k==203) idpdg=16;
 400    if (idf2k==204) idpdg=-12;
 401    if (idf2k==205) idpdg=-14;
 402    if (idf2k==206) idpdg=-16;
 403   }
 404
 405   t.SetParticleCode(idpdg);
 406   t.SetName(fPdg->GetParticle(idpdg)->GetName());
 407   t.SetTitle("MC track");
 408   t.SetMass(fPdg->GetParticle(idpdg)->Mass());
 409   t.SetCharge(fPdg->GetParticle(idpdg)->Charge()/3.);
 410
 411   evt->AddTrack(t);
 412  }
 413
 414  // Create the pointers to each particle's parent particle.
 415  Int_t txid=0;
 416  Int_t parid=0;
 417  for (Int_t itk=1; itk<=evt->GetNtracks (); itk++)
 418  {
 419   AliTrack* tx=evt->GetTrack(itk);
 420
 421   if (!tx) continue;
 422
 423   txid=tx->GetId();
 424
 425   parid=-1;
 426   for (Int_t j=0; j<fEvent.ntrack; j++)
 427   {
 428    tid=fEvent.gen[j].tag;
 429    if (-abs(tid) == txid) parid=fEvent.gen[j].parent;
 430   }
 431
 432   if (parid<0) continue;
 433
 434   AliTrack* tpar=evt->GetIdTrack(-abs(parid));
 435
 436   if (!tpar) continue;
 437
 438   tx->SetParentTrack(tpar);
 439  }
 440 }
 441 ///////////////////////////////////////////////////////////////////////////
 442 void IceF2k::PutRecoTracks(IceEvent* evt)
 443 {
 444 // Get the reconstructed tracks from the F2000 file into the IcePack structure.
 445 // Note : Reco tracks are given positive track id's in the event structure.
 446 // This memberfunction is based on the original code by Adam Bouchta.
 447
 448  if (!evt || fEvent.nfit<=0) return;
 449
 450  // Loop over all the tracks and add them to the current event
 451  AliTrack t;
 452  Double_t vec[3];
 453  AliPosition r;
 454  Ali3Vector p;
 455  Int_t tid=0;
 456  Int_t idpdg=0;
 457  Int_t idf2k=0;
 458  for (Int_t i=0; i<fEvent.nfit; i++)
 459  {
 460   t.Reset ();
 461
 462   // Beginpoint of the track
 463   vec[0]=fEvent.rec[i].x;
 464   vec[1]=fEvent.rec[i].y;
 465   vec[2]=fEvent.rec[i].z;
 466   r.SetPosition(vec,"car");
 467   t.SetBeginPoint(r);
 468
 469   // Endpoint of the track
 470   vec[0]+=fEvent.rec[i].length*fEvent.rec[i].px;
 471   vec[1]+=fEvent.rec[i].length*fEvent.rec[i].py;
 472   vec[2]+=fEvent.rec[i].length*fEvent.rec[i].pz;
 473   r.SetPosition(vec,"car");
 474   t.SetEndPoint(r);
 475
 476   // Momentum in GeV/c
 477   if (fEvent.rec[i].e > 0)
 478   {
 479    vec[0]=fEvent.rec[i].e*fEvent.rec[i].px*1e-3;
 480    vec[1]=fEvent.rec[i].e*fEvent.rec[i].py*1e-3;
 481    vec[2]=fEvent.rec[i].e*fEvent.rec[i].pz*1e-3;
 482   }
 483   else // Give the track a nominal momentum of 1 GeV/c
 484   {
 485    vec[0]=fEvent.rec[i].px;
 486    vec[1]=fEvent.rec[i].py;
 487    vec[2]=fEvent.rec[i].pz;
 488   }
 489   p.SetVector (vec,"car");
 490   t.Set3Momentum(p);
 491
 492   // Use the fit number as track id
 493   tid=fEvent.rec[i].tag;
 494   t.SetId(abs(tid));
 495
 496   idf2k=fEvent.rec[i].id;
 497   idpdg=0;
 498   if (idf2k>1000)
 499   {
 500    idpdg=idf2k+10000000;
 501   }
 502   else if (idf2k <= 48)
 503   {
 504    idpdg=fPdg->ConvertGeant3ToPdg(idf2k);
 505   }
 506   else
 507   {
 508    if (idf2k==201) idpdg=12;
 509    if (idf2k==202) idpdg=14;
 510    if (idf2k==203) idpdg=16;
 511    if (idf2k==204) idpdg=-12;
 512    if (idf2k==205) idpdg=-14;
 513    if (idf2k==206) idpdg=-16;
 514   }
 515
 516   t.SetParticleCode(idpdg);
 517   t.SetName(fPdg->GetParticle(idpdg)->GetName());
 518   t.SetTitle("RECO track");
 519   t.SetMass(fPdg->GetParticle(idpdg)->Mass());
 520   t.SetCharge(fPdg->GetParticle(idpdg)->Charge()/3.);
 521
 522   // Retrieve the various fit parameters for this track
 523   AliSignal* fitdata=fFitdefs->GetIdHit(i);
 524   for (Int_t jval=0; jval<fEvent.fresult[i].nval; jval++)
 525   {
 526    fitdata->SetSignal(fEvent.fresult[i].val[jval],jval+1);
 527   }
 528
 529   // Store the various fit parameters for this track
 530   t.SetFitDetails(fitdata);
 531
 532   // Store the various reco tracks as track hypotheses.
 533   // A copy of the first reco track is entered as a new track instance
 534   // into the event and all reco tracks (incl. the first one) are
 535   // stored as hypotheses linked to this new reco track.
 536   if (i==0)
 537   {
 538    evt->AddTrack(t);
 539    AliTrack* tx=evt->GetTrack(evt->GetNtracks());
 540    Int_t nrec=evt->GetNtracks(1);
 541    tx->SetId(nrec+1);
 542   }
 543   AliTrack* tx=evt->GetTrack(evt->GetNtracks());
 544   if (tx) tx->AddTrackHypothesis(t);
 545  }
 546 }
 547 ///////////////////////////////////////////////////////////////////////////
 548 void IceF2k::PutHits(IceEvent* evt)
 549 {
 550 // Get the hit and waveform info from the F2000 file into the IcePack structure.
 551 // This memberfunction is based on the original code by Adam Bouchta.
 552
 553  if (!evt) return;
 554
 555  // Loop over all the hits and add them to the current event
 556  IceAOM om;
 557  AliSignal s;
 558  s.SetSlotName("ADC",1);
 559  s.SetSlotName("LE",2);
 560  s.SetSlotName("TOT",3);
 561  Int_t chan=0;
 562  Int_t maxchan=800;
 563  if (fOmdb) maxchan=fHeader.nch;
 564  IceAOM* omx=0;
 565  AliSignal* sx=0;
 566  Int_t tid=0;
 567  AliTrack* tx=0;
 568  for (Int_t i=0; i<fEvent.nhits; i++)
 569  {
 570   chan=fEvent.h[i].ch+1;
 571   if (chan>maxchan) continue; // Channels 9001, 9002 etc are trigger channels
 572
 573   // Get corresponding device from the current event structure
 574   omx=(IceAOM*)evt->GetIdDevice(chan);
 575   if (!omx)
 576   {
 577    if (fOmdb)
 578    {
 579     omx=(IceAOM*)fOmdb->GetObject(chan,1);
 580     evt->AddDevice(omx);
 581    }
 582    else
 583    {
 584     om.Reset(1);
 585     om.SetUniqueID(chan);
 586     evt->AddDevice(om);
 587    }
 588    omx=(IceAOM*)evt->GetIdDevice(chan);
 589   }
 590
 591   if (!omx) continue;
 592
 593   s.Reset();
 594   s.SetUniqueID(fEvent.h[i].id);
 595   s.SetSignal(fEvent.h[i].amp,1);
 596   s.SetSignal(fEvent.h[i].t,2);
 597   s.SetSignal(fEvent.h[i].tot,3);
 598
 599   omx->AddHit(s);
 600
 601   sx=omx->GetHit(omx->GetNhits());
 602   if (!sx) continue;
 603
 604   // Bi-directional link between this hit and the track that caused the ADC value.
 605   // This F2K info is probably only present for MC tracks.
 606   tid=fEvent.h[i].ma;
 607   if (tid > 0)
 608   {
 609    tx=evt->GetIdTrack(tid); // Reco tracks
 610    if (!tx) tx=evt->GetIdTrack(-tid); // MC tracks
 611    if (tx) sx->AddLink(tx);
 612   }
 613   else
 614   {
 615    if (tid == -2) s.SetNameTitle("N","Noise");
 616    if (tid == -3) s.SetNameTitle("A","Afterpulse");
 617   }
 618  }
 619
 620  // Loop over all the waveforms and add the histo(s) to the corresponding OM's
 621  TH1F histo;
 622  Int_t nbins=0;
 623  Float_t xlow=0;
 624  Float_t xup=0;
 625  TString hname;
 626  for (Int_t iwf=0; iwf<fEvent.nwf; iwf++)
 627  {
 628   chan=fEvent.wf[iwf].om;
 629   if (chan<=0 || chan>maxchan) continue; // Skip trigger channels
 630
 631   // Get corresponding device from the current event structure
 632   omx=(IceAOM*)evt->GetIdDevice(chan);
 633   if (!omx)
 634   {
 635    if (fOmdb)
 636    {
 637     omx=(IceAOM*)fOmdb->GetObject(chan,1);
 638     evt->AddDevice(omx);
 639    }
 640    else
 641    {
 642     om.Reset(1);
 643     om.SetUniqueID(chan);
 644     evt->AddDevice(om);
 645    }
 646    omx=(IceAOM*)evt->GetIdDevice(chan);
 647   }
 648
 649   if (!omx) continue;
 650
 651   omx->SetSlotName("BASELINE",omx->GetNnames()+1);
 652   omx->SetSignal(fEvent.wf[iwf].baseline,"BASELINE");
 653
 654   // Fill the waveform histogram
 655   hname="OM";
 656   hname+=chan;
 657   hname+="-WF";
 658   hname+=omx->GetNwaveforms()+1;
 659
 660   histo.Reset();
 661   histo.SetName(hname.Data());
 662   nbins=fEvent.wf[iwf].ndigi;
 663   xlow=fEvent.wf[iwf].t_start;
 664   xup=xlow+float(nbins)*fEvent.wf[iwf].t_bin;
 665   histo.SetBins(nbins,xlow,xup);
 666
 667   for (Int_t jbin=1; jbin<=fEvent.wf[iwf].ndigi; jbin++)
 668   {
 669    histo.SetBinContent(jbin,fEvent.wf[iwf].digi[jbin]);
 670   }
 671
 672   omx->SetWaveform(&histo,omx->GetNwaveforms()+1);
 673  }
 674
 675  // Set bi-directional links between hits and reco track hypotheses.
 676  // Note : Reco tracks are recognised by their positive id.
 677  Int_t hid=0;
 678  TObjArray* rectracks=evt->GetTracks(1);
 679  for (Int_t jtk=0; jtk<rectracks->GetEntries(); jtk++)
 680  {
 681   tx=(AliTrack*)rectracks->At(jtk);
 682   if (!tx) continue;
 683
 684   for (Int_t jhyp=1; jhyp<=tx->GetNhypotheses(); jhyp++)
 685   {
 686    AliTrack* hypx=tx->GetTrackHypothesis(jhyp);
 687    if (!hypx) continue;
 688
 689    // Loop over all combinations of F2K fits and used OM hits
 690    for (Int_t k=0; k<fEvent.nfit_uses; k++)
 691    {
 692     if (fEvent.fit_uses[k].useid != hypx->GetId()) continue;
 693     hid=fEvent.fit_uses[k].hitid;
 694     sx=evt->GetIdHit(hid,"IceAOM");
 695     if (sx) sx->AddLink(hypx);
 696    }
 697   }
 698  }
 699 }
 700 ///////////////////////////////////////////////////////////////////////////