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