9de1b27471bbce5fbaf35d4925c6be8b45b7a73c
[u/mrichter/AliRoot.git] / 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  Double_t pi=acos(-1.);
165
166  IceEvent* evt=new IceEvent();
167
168  evt->SetTrackCopy(1);
169  evt->SetDevCopy(1);
170
171  // Branch in the tree for the event structure
172  otree->Branch("IceEvent","IceEvent",&evt,fBsize,fSplit); 
173
174  // Create the particle database and extend it with some F2000 specific definitions
175  if (!fPdg) fPdg=new TDatabasePDG();
176  Double_t me=fPdg->GetParticle(11)->Mass();
177  fPdg->AddParticle("brems"   ,"brems"   ,0,1,0,0,"none",10001001,0,0);
178  fPdg->AddParticle("deltae"  ,"deltae"  ,me,1,0,-3,"Lepton",10001002,0,0);
179  fPdg->AddParticle("pairprod","pairprod",0,1,0,0,"none",10001003,0,0);
180  fPdg->AddParticle("nucl_int","nucl_Int",0,1,0,0,"none",10001004,0,0);
181  fPdg->AddParticle("mu_pair" ,"mu_pair" ,0,1,0,0,"none",10001005,0,0);
182  fPdg->AddParticle("hadrons" ,"hadrons" ,0,1,0,0,"none",10001006,0,0);
183  fPdg->AddParticle("fiberlaser","fiberlaser",0,1,0,0,"none",10002100,0,0);
184  fPdg->AddParticle("n2laser"   ,"n2laser"   ,0,1,0,0,"none",10002101,0,0);
185  fPdg->AddParticle("yaglaser"  ,"yaglaser"  ,0,1,0,0,"none",10002201,0,0);
186  fPdg->AddParticle("z_primary","z_primary",0,1,0,0,"none",10003000,0,0);
187  fPdg->AddParticle("a_primary","a_primary",0,1,0,0,"none",10003500,0,0);
188
189  // Fill the database with geometry, calib. etc... parameters
190  // for all the devices
191  FillOMdbase();
192
193  // Set the fit definitions according to the F2000 header info
194  SetFitdefs();
195
196  for (Int_t jentry=0; jentry<nentries; jentry++)
197  {
198   if (rdmc_revt(fInput,&fHeader,&fEvent) != 0) break;
199
200   // Reset the complete Event structure
201   evt->Reset();
202
203   evt->SetRunNumber(fEvent.nrun);
204   evt->SetEventNumber(fEvent.enr);
205   evt->SetMJD(fEvent.mjd,fEvent.secs,fEvent.nsecs);
206
207   PutMcTracks(evt);
208
209   PutRecoTracks(evt);
210
211   PutHits(evt);
212
213   if (!(jentry%printfreq))
214   {
215    evt->HeaderData();
216   }
217
218   // Write the complete structure to the output Tree
219   otree->Fill();
220  }
221
222  if (evt) delete evt;
223 }
224 ///////////////////////////////////////////////////////////////////////////
225 void IceF2k::FillOMdbase()
226 {
227 // Fill the database with geometry, calib. etc... parameters 
228 // for all the devices.
229
230  if (fHeader.nch<=0) return;
231
232  if (fOmdb)
233  {
234   fOmdb->Reset();
235  }
236  else
237  {
238   fOmdb=new AliObjMatrix();
239   fOmdb->SetNameTitle("OMDBASE","The OM geometry, calib. etc... database");
240   fOmdb->SetOwner();
241  }
242
243  IceAOM* dev=0;
244  Double_t pos[3]={0,0,0};
245  for (Int_t i=0; i<fHeader.nch; i++)
246  {
247   dev=new IceAOM();
248   dev->SetUniqueID(i+1);
249   dev->SetSlotName("TYPE",1);
250   dev->SetSlotName("ORIENT",2);
251   dev->SetSlotName("T0",3);
252   dev->SetSlotName("ALPHA",4);
253   dev->SetSlotName("KADC",5);
254   dev->SetSlotName("KTOT",6);
255   dev->SetSlotName("KTDC",7);
256
257   pos[0]=fHeader.x[i];
258   pos[1]=fHeader.y[i];
259   pos[2]=fHeader.z[i];
260   dev->SetPosition(pos,"car");
261   dev->SetSignal(fHeader.type[i],1);
262   dev->SetSignal((Float_t)fHeader.costh[i],2);
263   dev->SetSignal(fHeader.cal[i].t_0,3);
264   dev->SetSignal(fHeader.cal[i].alpha_t,4);
265   dev->SetSignal(fHeader.cal[i].beta_a,5);
266   dev->SetSignal(fHeader.cal[i].beta_tot,6);
267   dev->SetSignal(fHeader.cal[i].beta_t,7);
268   fOmdb->EnterObject(i+1,1,dev);
269  }
270 }
271 ///////////////////////////////////////////////////////////////////////////
272 void IceF2k::SetFitdefs()
273 {
274 // Obtain the names of the variables for each fit procedure from the
275 // f2000 header. Each different fit procedure is then stored as a separate
276 // hit of an AliDevice object and the various fit variables are stored
277 // as separate signal slots of the corresponding hit.
278 // Via the GetFitdefs() memberfunction this AliDevice object can be
279 // retrieved and stored in the ROOT output file if wanted.
280 // The name of the object is FitDefinitions and the stored data can be
281 // inspected via the AliDevice::Data() memberfunction and looks as follows :
282 //
283 //  *AliDevice::Data* Id :0 Name : FitDefinitions
284 //    Position Vector in car coordinates : 0 0 0
285 //    Err. in car coordinates : 0 0 0
286 //  The following 8 hits are registered :
287 //  *AliSignal::Data* Id :0
288 //    Position Vector in car coordinates : 0 0 0
289 //    Err. in car coordinates : 0 0 0
290 //    Owned by device : AliDevice Name : FitDefinitions
291 //    Slot : 1 Signal value : 1 name : id
292 //    Slot : 2 Signal value : 2 name : rchi2
293 //    Slot : 3 Signal value : 3 name : prob
294 //    Slot : 4 Signal value : 4 name : sigth
295 //    Slot : 5 Signal value : 5 name : covmin
296 //    Slot : 6 Signal value : 6 name : covmax
297 //    Slot : 7 Signal value : 7 name : cutflag
298 //    Slot : 8 Signal value : 8 name : chi2
299 //  *AliSignal::Data* Id :1
300 //    Position Vector in car coordinates : 0 0 0
301 //    Err. in car coordinates : 0 0 0
302 //    Owned by device : AliDevice Name : FitDefinitions
303 //    Slot : 1 Signal value : 1 name : id
304 //    Slot : 2 Signal value : 2 name : rchi2
305 //    Slot : 3 Signal value : 3 name : prob
306 // etc....  
307 //
308 // This memberfunction is based on the original idea/code by Adam Bouchta.
309
310  if (fHeader.n_fit<=0) return;
311
312  if (fFitdefs)
313  {
314   fFitdefs->Reset(1);
315  }
316  else
317  {
318   fFitdefs=new AliDevice();
319  }
320
321  fFitdefs->SetName("FitDefinitions");
322  fFitdefs->SetHitCopy (1);
323
324  AliSignal s;
325  s.Reset();
326
327  for (Int_t i=0; i<fHeader.n_fit; i++)
328  {
329   s.SetUniqueID(fHeader.def_fit[i].id);
330
331   for (Int_t j=0; j<fHeader.def_fit[i].nwords; j++)
332   {
333    s.SetSlotName(TString(fHeader.def_fit[i].words[j]),j+1);
334    s.SetSignal(j+1,j+1);
335   }
336
337   fFitdefs->AddHit(s);
338   s.Reset(1);
339  }
340 }
341 ///////////////////////////////////////////////////////////////////////////
342 void IceF2k::PutMcTracks(IceEvent* evt)
343 {
344 // Get the MC tracks from the F2000 file into the IcePack structure.
345 // Note : MC tracks are given negative track id's in the event structure.
346 // This memberfunction is based on the original code by Adam Bouchta.
347
348  if (!evt || fEvent.ntrack<=0) return;
349
350  // Loop over all the tracks and add them to the current event
351  AliTrack t;
352  Double_t vec[3];
353  AliPosition r;
354  Ali3Vector p;
355  Int_t tid=0;
356  Int_t idpdg=0;
357  Int_t idf2k=0;
358  for (Int_t i=0; i<fEvent.ntrack; i++)
359  {
360   t.Reset ();
361
362   // Beginpoint of the track
363   vec[0]=fEvent.gen[i].x;
364   vec[1]=fEvent.gen[i].y;
365   vec[2]=fEvent.gen[i].z;
366   r.SetPosition(vec,"car");
367   t.SetBeginPoint(r);
368
369   // Endpoint of the track
370   vec[0]+=fEvent.gen[i].length*fEvent.gen[i].px;
371   vec[1]+=fEvent.gen[i].length*fEvent.gen[i].py;
372   vec[2]+=fEvent.gen[i].length*fEvent.gen[i].pz;
373   r.SetPosition(vec,"car");
374   t.SetEndPoint(r);
375
376   // Momentum in GeV/c
377   vec[0]=fEvent.gen[i].e*fEvent.gen[i].px*1e-3;
378   vec[1]=fEvent.gen[i].e*fEvent.gen[i].py*1e-3;
379   vec[2]=fEvent.gen[i].e*fEvent.gen[i].pz*1e-3;
380   p.SetVector (vec,"car");
381   t.Set3Momentum(p);
382
383   // MC tracks are indicated by negative track id's
384   tid=fEvent.gen[i].tag;
385   t.SetId(-abs(tid));
386
387   idf2k=fEvent.gen[i].id;
388   idpdg=0;
389   if (idf2k>1000)
390   {
391    idpdg=idf2k+10000000;
392   }
393   else if (idf2k <= 48)
394   {
395    idpdg=fPdg->ConvertGeant3ToPdg(idf2k);
396   }
397   else
398   {
399    if (idf2k==201) idpdg=12;
400    if (idf2k==202) idpdg=14;
401    if (idf2k==203) idpdg=16;
402    if (idf2k==204) idpdg=-12;
403    if (idf2k==205) idpdg=-14;
404    if (idf2k==206) idpdg=-16;
405   }
406
407   t.SetParticleCode(idpdg);
408   t.SetName(fPdg->GetParticle(idpdg)->GetName());
409   t.SetTitle("MC track");
410   t.SetMass(fPdg->GetParticle(idpdg)->Mass());
411   t.SetCharge(fPdg->GetParticle(idpdg)->Charge()/3.);
412
413   evt->AddTrack(t);
414  }
415
416  // Create the pointers to each particle's parent particle.
417  Int_t txid=0;
418  Int_t parid=0;
419  for (Int_t itk=1; itk<=evt->GetNtracks (); itk++)
420  {
421   AliTrack* tx=evt->GetTrack(itk);
422
423   if (!tx) continue;
424
425   txid=tx->GetId();
426
427   parid=-1;
428   for (Int_t j=0; j<fEvent.ntrack; j++)
429   {
430    tid=fEvent.gen[j].tag;
431    if (-abs(tid) == txid) parid=fEvent.gen[j].parent;
432   }
433
434   if (parid<0) continue;
435
436   AliTrack* tpar=evt->GetIdTrack(-abs(parid));
437
438   if (!tpar) continue;
439
440   tx->SetParentTrack(tpar);
441  }
442 }
443 ///////////////////////////////////////////////////////////////////////////
444 void IceF2k::PutRecoTracks(IceEvent* evt)
445 {
446 // Get the reconstructed tracks from the F2000 file into the IcePack structure.
447 // Note : Reco tracks are given positive track id's in the event structure.
448 // This memberfunction is based on the original code by Adam Bouchta.
449
450  if (!evt || fEvent.nfit<=0) return;
451
452  // Loop over all the tracks and add them to the current event
453  AliTrack t;
454  Double_t vec[3];
455  AliPosition r;
456  Ali3Vector p;
457  Int_t tid=0;
458  Int_t idpdg=0;
459  Int_t idf2k=0;
460  for (Int_t i=0; i<fEvent.nfit; i++)
461  {
462   t.Reset ();
463
464   // Beginpoint of the track
465   vec[0]=fEvent.rec[i].x;
466   vec[1]=fEvent.rec[i].y;
467   vec[2]=fEvent.rec[i].z;
468   r.SetPosition(vec,"car");
469   t.SetBeginPoint(r);
470
471   // Endpoint of the track
472   vec[0]+=fEvent.rec[i].length*fEvent.rec[i].px;
473   vec[1]+=fEvent.rec[i].length*fEvent.rec[i].py;
474   vec[2]+=fEvent.rec[i].length*fEvent.rec[i].pz;
475   r.SetPosition(vec,"car");
476   t.SetEndPoint(r);
477
478   // Momentum in GeV/c
479   if (fEvent.rec[i].e > 0)
480   {
481    vec[0]=fEvent.rec[i].e*fEvent.rec[i].px*1e-3;
482    vec[1]=fEvent.rec[i].e*fEvent.rec[i].py*1e-3;
483    vec[2]=fEvent.rec[i].e*fEvent.rec[i].pz*1e-3;
484   }
485   else // Give the track a nominal momentum of 1 GeV/c
486   {
487    vec[0]=fEvent.rec[i].px;
488    vec[1]=fEvent.rec[i].py;
489    vec[2]=fEvent.rec[i].pz;
490   }
491   p.SetVector (vec,"car");
492   t.Set3Momentum(p);
493
494   // Use the fit number as track id
495   tid=fEvent.rec[i].tag;
496   t.SetId(abs(tid));
497
498   idf2k=fEvent.rec[i].id;
499   idpdg=0;
500   if (idf2k>1000)
501   {
502    idpdg=idf2k+10000000;
503   }
504   else if (idf2k <= 48)
505   {
506    idpdg=fPdg->ConvertGeant3ToPdg(idf2k);
507   }
508   else
509   {
510    if (idf2k==201) idpdg=12;
511    if (idf2k==202) idpdg=14;
512    if (idf2k==203) idpdg=16;
513    if (idf2k==204) idpdg=-12;
514    if (idf2k==205) idpdg=-14;
515    if (idf2k==206) idpdg=-16;
516   }
517
518   t.SetParticleCode(idpdg);
519   t.SetName(fPdg->GetParticle(idpdg)->GetName());
520   t.SetTitle("RECO track");
521   t.SetMass(fPdg->GetParticle(idpdg)->Mass());
522   t.SetCharge(fPdg->GetParticle(idpdg)->Charge()/3.);
523
524   // Retrieve the various fit parameters for this track
525   AliSignal* fitdata=fFitdefs->GetIdHit(i);
526   for (Int_t jval=0; jval<fEvent.fresult[i].nval; jval++)
527   {
528    fitdata->SetSignal(fEvent.fresult[i].val[jval],jval+1);
529   }
530
531   // Store the various fit parameters for this track
532   t.SetFitDetails(fitdata);
533
534   // Store the various reco tracks as track hypotheses.
535   // A copy of the first reco track is entered as a new track instance
536   // into the event and all reco tracks (incl. the first one) are
537   // stored as hypotheses linked to this new reco track.
538   if (i==0)
539   {
540    evt->AddTrack(t);
541    AliTrack* tx=evt->GetTrack(evt->GetNtracks());
542    Int_t nrec=evt->GetNtracks(1);
543    tx->SetId(nrec+1);
544   }
545   AliTrack* tx=evt->GetTrack(evt->GetNtracks());
546   if (tx) tx->AddTrackHypothesis(t);
547  }
548 }
549 ///////////////////////////////////////////////////////////////////////////
550 void IceF2k::PutHits(IceEvent* evt)
551 {
552 // Get the hit and waveform info from the F2000 file into the IcePack structure.
553 // This memberfunction is based on the original code by Adam Bouchta.
554
555  if (!evt) return;
556
557  // Loop over all the hits and add them to the current event
558  IceAOM om;
559  AliSignal s;
560  s.SetSlotName("ADC",1);
561  s.SetSlotName("LE",2);
562  s.SetSlotName("TOT",3);
563  Int_t chan=0;
564  Int_t maxchan=800;
565  if (fOmdb) maxchan=fHeader.nch;
566  IceAOM* omx=0;
567  AliSignal* sx=0;
568  Int_t tid=0;
569  AliTrack* tx=0;
570  for (Int_t i=0; i<fEvent.nhits; i++)
571  {
572   chan=fEvent.h[i].ch+1;
573   if (chan>maxchan) continue; // Channels 9001, 9002 etc are trigger channels
574
575   // Get corresponding device from the current event structure  
576   omx=(IceAOM*)evt->GetIdDevice(chan);
577   if (!omx)
578   {
579    if (fOmdb)
580    {
581     omx=(IceAOM*)fOmdb->GetObject(chan,1);
582     evt->AddDevice(omx);
583    }
584    else
585    {
586     om.Reset(1);
587     om.SetUniqueID(chan);
588     evt->AddDevice(om);
589    }
590    omx=(IceAOM*)evt->GetIdDevice(chan);
591   }
592
593   if (!omx) continue;
594
595   s.Reset();
596   s.SetUniqueID(fEvent.h[i].id);
597   s.SetSignal(fEvent.h[i].amp,1);
598   s.SetSignal(fEvent.h[i].t,2);
599   s.SetSignal(fEvent.h[i].tot,3);
600
601   omx->AddHit(s);
602
603   sx=omx->GetHit(omx->GetNhits());
604   if (!sx) continue;
605
606   // Bi-directional link between this hit and the track that caused the ADC value.
607   // This F2K info is probably only present for MC tracks.
608   tid=fEvent.h[i].ma;
609   if (tid > 0)
610   {
611    tx=evt->GetIdTrack(tid); // Reco tracks
612    if (!tx) tx=evt->GetIdTrack(-tid); // MC tracks
613    if (tx) sx->AddLink(tx);
614   }
615   else
616   {
617    if (tid == -2) s.SetNameTitle("N","Noise");
618    if (tid == -3) s.SetNameTitle("A","Afterpulse");
619   }
620  }
621
622  // Loop over all the waveforms and add the histo(s) to the corresponding OM's
623  TH1F histo;
624  Int_t nbins=0;
625  Float_t xlow=0;
626  Float_t xup=0;
627  TString hname;
628  for (Int_t iwf=0; iwf<fEvent.nwf; iwf++)
629  {
630   chan=fEvent.wf[iwf].om;
631   if (chan<=0 || chan>maxchan) continue; // Skip trigger channels
632
633   // Get corresponding device from the current event structure  
634   omx=(IceAOM*)evt->GetIdDevice(chan);
635   if (!omx)
636   {
637    if (fOmdb)
638    {
639     omx=(IceAOM*)fOmdb->GetObject(chan,1);
640     evt->AddDevice(omx);
641    }
642    else
643    {
644     om.Reset(1);
645     om.SetUniqueID(chan);
646     evt->AddDevice(om);
647    }
648    omx=(IceAOM*)evt->GetIdDevice(chan);
649   }
650
651   if (!omx) continue;
652
653   omx->SetSlotName("BASELINE",omx->GetNnames()+1);
654   omx->SetSignal(fEvent.wf[iwf].baseline,"BASELINE");
655
656   // Fill the waveform histogram
657   hname="OM";
658   hname+=chan;
659   hname+="-WF";
660   hname+=omx->GetNwaveforms()+1;
661
662   histo.Reset();
663   histo.SetName(hname.Data());
664   nbins=fEvent.wf[iwf].ndigi;
665   xlow=fEvent.wf[iwf].t_start;
666   xup=xlow+float(nbins)*fEvent.wf[iwf].t_bin;
667   histo.SetBins(nbins,xlow,xup);
668
669   for (Int_t jbin=1; jbin<=fEvent.wf[iwf].ndigi; jbin++)
670   {
671    histo.SetBinContent(jbin,fEvent.wf[iwf].digi[jbin]);
672   }
673
674   omx->SetWaveform(&histo,omx->GetNwaveforms()+1);
675  }
676
677  // Set bi-directional links between hits and reco track hypotheses.
678  // Note : Reco tracks are recognised by their positive id.
679  Int_t hid=0;
680  TObjArray* rectracks=evt->GetTracks(1);
681  for (Int_t jtk=0; jtk<rectracks->GetEntries(); jtk++)
682  {
683   tx=(AliTrack*)rectracks->At(jtk);
684   if (!tx) continue;
685   
686   for (Int_t jhyp=1; jhyp<=tx->GetNhypotheses(); jhyp++)
687   {
688    AliTrack* hypx=tx->GetTrackHypothesis(jhyp);
689    if (!hypx) continue;
690
691    // Loop over all combinations of F2K fits and used OM hits
692    for (Int_t k=0; k<fEvent.nfit_uses; k++)
693    {
694     if (fEvent.fit_uses[k].useid != hypx->GetId()) continue;
695     hid=fEvent.fit_uses[k].hitid;
696     sx=evt->GetIdHit(hid,"IceAOM");
697     if (sx) sx->AddLink(hypx);
698    }
699   }
700  }
701 }
702 ///////////////////////////////////////////////////////////////////////////