c2a86b6ebc86e8015366bd443f98a897c186e95f
[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 the output file.
93 // // The output file is not explicitly closed here
94 // // to allow ineractive investigation of the data tree
95 // // when this macro is run in an interactive ROOT/CINT session.
96 // ofile->Write();
97 //
98 //--- Author: Nick van Eijndhoven 11-mar-2005 Utrecht University
99 //- Modified: NvE $Date$ Utrecht University
100 ///////////////////////////////////////////////////////////////////////////
101  
102 #include "IceF2k.h"
103 #include "Riostream.h"
104
105 ClassImp(IceF2k) // Class implementation to enable ROOT I/O
106
107 IceF2k::IceF2k(const char* name,const char* title) : AliJob(name,title)
108 {
109 // Default constructor.
110 // By default maxevent=-1, split=99, bsize=32000, printfreq=1.
111
112  fSplit=99;
113  fBsize=32000;
114  fMaxevt=-1;
115  fPrintfreq=1;
116  fInfile="";
117  fOutfile=0;
118
119  fPdg=0;
120  fOmdb=0;
121  fFitdefs=0;
122 }
123 ///////////////////////////////////////////////////////////////////////////
124 IceF2k::~IceF2k()
125 {
126 // Default destructor.
127
128  if (fPdg)
129  {
130   delete fPdg;
131   fPdg=0;
132  }
133
134  if (fOmdb)
135  {
136   delete fOmdb;
137   fOmdb=0;
138  }
139
140  if (fFitdefs)
141  {
142   delete fFitdefs;
143   fFitdefs=0;
144  }
145 }
146 ///////////////////////////////////////////////////////////////////////////
147 void IceF2k::SetMaxEvents(Int_t n)
148 {
149 // Set the maximum number of events to be processed.
150 // n=-1 implies processing of the complete input file, which is the default
151 // initialisation in the constructor.
152  fMaxevt=n;
153 }
154 ///////////////////////////////////////////////////////////////////////////
155 void IceF2k::SetPrintFreq(Int_t f)
156 {
157 // Set the printfrequency to produce info every f events.
158 // f=1 is the default initialisation in the constructor.
159  if (f>=0) fPrintfreq=f;
160 }
161 ///////////////////////////////////////////////////////////////////////////
162 void IceF2k::SetSplitLevel(Int_t split)
163 {
164 // Set the split level for the ROOT data file.
165 // split=99 is the default initialisation in the constructor.
166  if (split>=0) fSplit=split;
167 }
168 ///////////////////////////////////////////////////////////////////////////
169 void IceF2k::SetBufferSize(Int_t bsize)
170 {
171 // Set the buffer size for the ROOT data file.
172 // bsize=32000 is the default initialisation in the constructor.
173  if (bsize>=0) fBsize=bsize;
174 }
175 ///////////////////////////////////////////////////////////////////////////
176 void IceF2k::SetInputFile(TString name)
177 {
178 // Set the name of the F2K input file.
179  fInfile=name;
180 }
181 ///////////////////////////////////////////////////////////////////////////
182 void IceF2k::SetOutputFile(TFile* ofile)
183 {
184 // Set the output file for the ROOT data.
185  fOutfile=ofile;
186 }
187 ///////////////////////////////////////////////////////////////////////////
188 TDatabasePDG* IceF2k::GetPDG()
189 {
190 // Provide pointer to the PDG database
191  return fPdg;
192 }
193 ///////////////////////////////////////////////////////////////////////////
194 AliObjMatrix* IceF2k::GetOMdbase()
195 {
196 // Provide pointer to the OM geometry, calib. etc... database
197  return fOmdb;
198 }
199 ///////////////////////////////////////////////////////////////////////////
200 AliDevice* IceF2k::GetFitdefs()
201 {
202 // Provide pointer to the fit definitions
203  return fFitdefs;
204 }
205 ///////////////////////////////////////////////////////////////////////////
206 void IceF2k::Exec(Option_t* opt)
207 {
208 // Job to loop over the specified number of events and convert the 
209 // F2K data into the IceEvent structure.
210 // If maxevents<0 (default) all the entries of the input file
211 // will be processed.
212 // Every "printfreq" events a short event summary will be printed.
213 // The default value is printfreq=1.
214 // The output will be written on a standard output tree named "T".
215 //
216 // Notes :
217 // -------
218 // 1) This class is derived from AliJob, allowing a task based processing.
219 //    After the conversion of an F2K event into an IceEvent structure,
220 //    the processing of all available sub-tasks (if any) is invoked.
221 //    This provides an event-by-event (sub)task processing before the
222 //    final data structures are written out.
223 // 2) The main object in this job environment is an IceEvent* pointer.
224
225  if (fInfile=="")
226  {
227   cout << " *IceF2k Exec* No data input file specified." << endl;
228   return;
229  }
230
231  // Open the input file in the default ascii format (autodetection) for reading 
232  fInput=rdmc_mcopen(fInfile.Data(),"r",RDMC_DEFAULT_ASCII_F);
233
234  if (!fInput)
235  {
236   cout << " *IceF2k Exec* No input file found with name : " << fInfile.Data() << endl;
237   return;
238  }
239
240  // Initialise the event structure 
241  rdmc_init_mevt(&fEvent);
242
243  // Read the file header information
244  rdmc_rarr(fInput,&fHeader);
245
246  TTree* otree=0;
247  if (fOutfile)
248  {
249   otree=new TTree("T","F2K Data converted to IceEvent structures");
250   otree->SetDirectory(fOutfile);
251  }
252
253  IceEvent* evt=new IceEvent();
254  evt->SetTrackCopy(1);
255  evt->SetDevCopy(1);
256
257  // Branch in the tree for the event structure
258  if (otree) otree->Branch("IceEvent","IceEvent",&evt,fBsize,fSplit); 
259
260  // Create the particle database and extend it with some F2000 specific definitions
261  if (!fPdg) fPdg=new TDatabasePDG();
262  Double_t me=fPdg->GetParticle(11)->Mass();
263  fPdg->AddParticle("brems"   ,"brems"   ,0,1,0,0,"none",10001001,0,0);
264  fPdg->AddParticle("deltae"  ,"deltae"  ,me,1,0,-3,"Lepton",10001002,0,0);
265  fPdg->AddParticle("pairprod","pairprod",0,1,0,0,"none",10001003,0,0);
266  fPdg->AddParticle("nucl_int","nucl_Int",0,1,0,0,"none",10001004,0,0);
267  fPdg->AddParticle("mu_pair" ,"mu_pair" ,0,1,0,0,"none",10001005,0,0);
268  fPdg->AddParticle("hadrons" ,"hadrons" ,0,1,0,0,"none",10001006,0,0);
269  fPdg->AddParticle("fiberlaser","fiberlaser",0,1,0,0,"none",10002100,0,0);
270  fPdg->AddParticle("n2laser"   ,"n2laser"   ,0,1,0,0,"none",10002101,0,0);
271  fPdg->AddParticle("yaglaser"  ,"yaglaser"  ,0,1,0,0,"none",10002201,0,0);
272  fPdg->AddParticle("z_primary","z_primary",0,1,0,0,"none",10003000,0,0);
273  fPdg->AddParticle("a_primary","a_primary",0,1,0,0,"none",10003500,0,0);
274
275  // Fill the database with geometry, calib. etc... parameters
276  // for all the devices
277  FillOMdbase();
278
279  // Set the fit definitions according to the F2000 header info
280  SetFitdefs();
281
282  // Initialise the job working environment
283  SetMainObject(evt);
284  if (fOutfile)
285  {
286   AddObject(fOutfile);
287   AddObject(otree);
288  }
289
290  cout << " ***" << endl;
291  cout << " *** Start processing of job " << GetName() << " ***" << endl;
292  cout << " ***" << endl;
293  cout << " F2K input file : " << fInfile.Data() << endl;
294  cout << " Maximum number of events to be processed : " << fMaxevt << endl;
295  cout << " Print frequency : " << fPrintfreq << endl;
296  if (fOutfile)
297  {
298   cout << " ROOT output file : " << fOutfile->GetName() << endl;
299   cout << " Output characteristics : splitlevel = " << fSplit << " buffersize = " << fBsize << endl;
300  }
301
302  ListEnvironment();
303  
304  Int_t nevt=0;
305  while (!rdmc_revt(fInput,&fHeader,&fEvent))
306  {
307   if (fMaxevt>-1 && nevt>=fMaxevt) break;
308
309   // Reset the complete Event structure
310   evt->Reset();
311
312   evt->SetRunNumber(fEvent.nrun);
313   evt->SetEventNumber(fEvent.enr);
314   evt->SetMJD(fEvent.mjd,fEvent.secs,fEvent.nsecs);
315
316   PutMcTracks();
317
318   PutRecoTracks();
319
320   PutHits();
321
322   // Invoke all available sub-tasks (if any)
323   CleanTasks();
324   ExecuteTasks(opt);
325
326   if (fPrintfreq)
327   {
328    if (!(nevt%fPrintfreq)) evt->HeaderData();
329   }
330
331   // Write the complete structure to the output Tree
332   if (otree) otree->Fill();
333
334   // Update event counter
335   nevt++;
336  }
337
338  // Remove the IceEvent object from the environment
339  // and delete it as well
340  if (evt)
341  {
342   RemoveObject(evt);
343   delete evt;
344  }
345 }
346 ///////////////////////////////////////////////////////////////////////////
347 void IceF2k::FillOMdbase()
348 {
349 // Fill the database with geometry, calib. etc... parameters 
350 // for all the devices.
351
352  if (fHeader.nch<=0) return;
353
354  Int_t geocal=fHeader.is_calib.geo;
355  Int_t adccal=fHeader.is_calib.adc;
356  Int_t tdccal=fHeader.is_calib.tdc;
357  Int_t totcal=fHeader.is_calib.tot;
358  Int_t utccal=fHeader.is_calib.utc;
359
360  TF1 fadccal("fadccal","(x-[1])*[0]");
361  TF1 fadcdecal("fadcdecal","(x/[0])+[1]");
362  fadccal.SetParName(0,"BETA-ADC");
363  fadccal.SetParName(1,"PED-ADC");
364  fadcdecal.SetParName(0,"BETA-ADC");
365  fadcdecal.SetParName(1,"PED-ADC");
366
367  TF1 ftdccal("ftdccal","(x*[0])-[1]-([0]-1.)*32767.-[2]/sqrt([3])");
368  TF1 ftdcdecal("ftdcdecal","(x+([0]-1.)*32767.+[1]+[2]/sqrt([3]))/[0]");
369  ftdccal.SetParName(0,"BETA-TDC");
370  ftdccal.SetParName(1,"T0");
371  ftdccal.SetParName(2,"ALPHA-TDC");
372  ftdccal.SetParName(3,"ADC-SLEW");
373  ftdcdecal.SetParName(0,"BETA-TDC");
374  ftdcdecal.SetParName(1,"T0");
375  ftdcdecal.SetParName(2,"ALPHA-TDC");
376  ftdcdecal.SetParName(3,"ADC-SLEW");
377
378  TF1 ftotcal("ftotcal","x*[0]");
379  TF1 ftotdecal("ftotdecal","x/[0]");
380  ftotcal.SetParName(0,"BETA-TOT");
381  ftotdecal.SetParName(0,"BETA-TOT");
382
383  if (fOmdb)
384  {
385   fOmdb->Reset();
386  }
387  else
388  {
389   fOmdb=new AliObjMatrix();
390   fOmdb->SetNameTitle("OMDBASE","The OM geometry, calib. etc... database");
391   fOmdb->SetOwner();
392  }
393
394  IceAOM* dev=0;
395  Double_t pos[3]={0,0,0};
396  for (Int_t i=0; i<fHeader.nch; i++)
397  {
398   dev=new IceAOM();
399   dev->SetUniqueID(i+1);
400
401   dev->SetSlotName("ADC",1);
402   dev->SetSlotName("LE",2);
403   dev->SetSlotName("TOT",3);
404
405   dev->SetSlotName("TYPE",4);
406   dev->SetSlotName("ORIENT",5);
407   dev->SetSlotName("THRESH",6);
408   dev->SetSlotName("SENSIT",7);
409   dev->SetSlotName("BETA-TDC",8);
410   dev->SetSlotName("T0",9);
411   dev->SetSlotName("ALPHA-TDC",10);
412   dev->SetSlotName("PED-ADC",11);
413   dev->SetSlotName("BETA-ADC",12);
414   dev->SetSlotName("KAPPA-ADC",13);
415   dev->SetSlotName("PED-TOT",14);
416   dev->SetSlotName("BETA-TOT",15);
417   dev->SetSlotName("KAPPA-TOT",16);
418
419   pos[0]=fHeader.x[i];
420   pos[1]=fHeader.y[i];
421   pos[2]=fHeader.z[i];
422   dev->SetPosition(pos,"car");
423
424   fadccal.SetParameter(0,fHeader.cal[i].beta_a);
425   fadccal.SetParameter(1,fHeader.cal[i].ped);
426   fadcdecal.SetParameter(0,fHeader.cal[i].beta_a);
427   if (!fHeader.cal[i].beta_a) fadcdecal.SetParameter(0,1);
428   fadcdecal.SetParameter(1,fHeader.cal[i].ped);
429
430   ftdccal.SetParameter(0,fHeader.cal[i].beta_t);
431   ftdccal.SetParameter(1,fHeader.cal[i].t_0);
432   ftdccal.SetParameter(2,fHeader.cal[i].alpha_t);
433   ftdccal.SetParameter(3,1.e20);
434   ftdcdecal.SetParameter(0,fHeader.cal[i].beta_t);
435   if (!fHeader.cal[i].beta_t) ftdcdecal.SetParameter(0,1);
436   ftdcdecal.SetParameter(1,fHeader.cal[i].t_0);
437   ftdcdecal.SetParameter(2,fHeader.cal[i].alpha_t);
438   ftdcdecal.SetParameter(3,1.e20);
439
440   ftotcal.SetParameter(0,fHeader.cal[i].beta_tot);
441   ftotdecal.SetParameter(0,fHeader.cal[i].beta_tot);
442   if (!fHeader.cal[i].beta_tot) ftotdecal.SetParameter(0,1);
443
444   if (adccal)
445   {
446    dev->SetDecalFunction(&fadcdecal,1);
447   }
448   else
449   {
450    dev->SetCalFunction(&fadccal,1);
451   }
452
453   if (tdccal)
454   {
455    dev->SetDecalFunction(&ftdcdecal,2);
456   }
457   else
458   {
459    dev->SetCalFunction(&ftdccal,2);
460   }
461
462   if (totcal)
463   {
464    dev->SetDecalFunction(&ftotdecal,3);
465   }
466   else
467   {
468    dev->SetCalFunction(&ftotcal,3);
469   }
470
471   dev->SetSignal(fHeader.type[i],4);
472   dev->SetSignal((Float_t)fHeader.costh[i],5);
473   dev->SetSignal(fHeader.thresh[i],6);
474   dev->SetSignal(fHeader.sensit[i],7);
475   dev->SetSignal(fHeader.cal[i].beta_t,8);
476   dev->SetSignal(fHeader.cal[i].t_0,9);
477   dev->SetSignal(fHeader.cal[i].alpha_t,10);
478   dev->SetSignal(fHeader.cal[i].ped,11);
479   dev->SetSignal(fHeader.cal[i].beta_a,12);
480   dev->SetSignal(fHeader.cal[i].kappa,13);
481   dev->SetSignal(fHeader.cal[i].ped_tot,14);
482   dev->SetSignal(fHeader.cal[i].beta_tot,15);
483   dev->SetSignal(fHeader.cal[i].kappa_tot,16);
484
485   fOmdb->EnterObject(i+1,1,dev);
486  }
487 }
488 ///////////////////////////////////////////////////////////////////////////
489 void IceF2k::SetFitdefs()
490 {
491 // Obtain the names of the variables for each fit procedure from the
492 // f2000 header. Each different fit procedure is then stored as a separate
493 // hit of an AliDevice object and the various fit variables are stored
494 // as separate signal slots of the corresponding hit.
495 // Via the GetFitdefs() memberfunction this AliDevice object can be
496 // retrieved and stored in the ROOT output file if wanted.
497 // The name of the object is FitDefinitions and the stored data can be
498 // inspected via the AliDevice::Data() memberfunction and looks as follows :
499 //
500 //  *AliDevice::Data* Id :0 Name : FitDefinitions
501 //    Position Vector in car coordinates : 0 0 0
502 //    Err. in car coordinates : 0 0 0
503 //  The following 8 hits are registered :
504 //  *AliSignal::Data* Id :0
505 //    Position Vector in car coordinates : 0 0 0
506 //    Err. in car coordinates : 0 0 0
507 //    Owned by device : AliDevice Name : FitDefinitions
508 //    Slot : 1 Signal value : 1 name : id
509 //    Slot : 2 Signal value : 2 name : rchi2
510 //    Slot : 3 Signal value : 3 name : prob
511 //    Slot : 4 Signal value : 4 name : sigth
512 //    Slot : 5 Signal value : 5 name : covmin
513 //    Slot : 6 Signal value : 6 name : covmax
514 //    Slot : 7 Signal value : 7 name : cutflag
515 //    Slot : 8 Signal value : 8 name : chi2
516 //  *AliSignal::Data* Id :1
517 //    Position Vector in car coordinates : 0 0 0
518 //    Err. in car coordinates : 0 0 0
519 //    Owned by device : AliDevice Name : FitDefinitions
520 //    Slot : 1 Signal value : 1 name : id
521 //    Slot : 2 Signal value : 2 name : rchi2
522 //    Slot : 3 Signal value : 3 name : prob
523 // etc....  
524 //
525 // This memberfunction is based on the original idea/code by Adam Bouchta.
526
527  if (fHeader.n_fit<=0) return;
528
529  if (fFitdefs)
530  {
531   fFitdefs->Reset(1);
532  }
533  else
534  {
535   fFitdefs=new AliDevice();
536  }
537
538  fFitdefs->SetName("FitDefinitions");
539  fFitdefs->SetHitCopy (1);
540
541  AliSignal s;
542  s.Reset();
543
544  for (Int_t i=0; i<fHeader.n_fit; i++)
545  {
546   s.SetUniqueID(fHeader.def_fit[i].id);
547
548   for (Int_t j=0; j<fHeader.def_fit[i].nwords; j++)
549   {
550    s.SetSlotName(TString(fHeader.def_fit[i].words[j]),j+1);
551    s.SetSignal(j+1,j+1);
552   }
553
554   fFitdefs->AddHit(s);
555   s.Reset(1);
556  }
557 }
558 ///////////////////////////////////////////////////////////////////////////
559 void IceF2k::PutMcTracks()
560 {
561 // Get the MC tracks from the F2000 file into the IcePack structure.
562 // Note : MC tracks are given negative track id's in the event structure.
563 // This memberfunction is based on the original code by Adam Bouchta.
564
565  IceEvent* evt=(IceEvent*)GetMainObject();
566  if (!evt || fEvent.ntrack<=0) return;
567
568  // Loop over all the tracks and add them to the current event
569  AliTrack t;
570  Double_t vec[3];
571  AliPosition r;
572  Ali3Vector p;
573  Int_t tid=0;
574  Int_t idpdg=0;
575  Int_t idf2k=0;
576  for (Int_t i=0; i<fEvent.ntrack; i++)
577  {
578   t.Reset ();
579
580   // Beginpoint of the track
581   vec[0]=fEvent.gen[i].x;
582   vec[1]=fEvent.gen[i].y;
583   vec[2]=fEvent.gen[i].z;
584   r.SetPosition(vec,"car");
585   t.SetBeginPoint(r);
586
587   // Endpoint of the track
588   vec[0]+=fEvent.gen[i].length*fEvent.gen[i].px;
589   vec[1]+=fEvent.gen[i].length*fEvent.gen[i].py;
590   vec[2]+=fEvent.gen[i].length*fEvent.gen[i].pz;
591   r.SetPosition(vec,"car");
592   t.SetEndPoint(r);
593
594   // Momentum in GeV/c
595   vec[0]=fEvent.gen[i].e*fEvent.gen[i].px*1e-3;
596   vec[1]=fEvent.gen[i].e*fEvent.gen[i].py*1e-3;
597   vec[2]=fEvent.gen[i].e*fEvent.gen[i].pz*1e-3;
598   p.SetVector (vec,"car");
599   t.Set3Momentum(p);
600
601   // MC tracks are indicated by negative track id's
602   tid=fEvent.gen[i].tag;
603   t.SetId(-abs(tid));
604
605   idf2k=fEvent.gen[i].id;
606   idpdg=0;
607   if (idf2k>1000)
608   {
609    idpdg=idf2k+10000000;
610   }
611   else if (idf2k <= 48)
612   {
613    idpdg=fPdg->ConvertGeant3ToPdg(idf2k);
614   }
615   else
616   {
617    if (idf2k==201) idpdg=12;
618    if (idf2k==202) idpdg=14;
619    if (idf2k==203) idpdg=16;
620    if (idf2k==204) idpdg=-12;
621    if (idf2k==205) idpdg=-14;
622    if (idf2k==206) idpdg=-16;
623   }
624
625   t.SetParticleCode(idpdg);
626   t.SetName(fPdg->GetParticle(idpdg)->GetName());
627   t.SetTitle("MC track");
628   t.SetMass(fPdg->GetParticle(idpdg)->Mass());
629   t.SetCharge(fPdg->GetParticle(idpdg)->Charge()/3.);
630
631   evt->AddTrack(t);
632  }
633
634  // Create the pointers to each particle's parent particle.
635  Int_t txid=0;
636  Int_t parid=0;
637  for (Int_t itk=1; itk<=evt->GetNtracks (); itk++)
638  {
639   AliTrack* tx=evt->GetTrack(itk);
640
641   if (!tx) continue;
642
643   txid=tx->GetId();
644
645   parid=-1;
646   for (Int_t j=0; j<fEvent.ntrack; j++)
647   {
648    tid=fEvent.gen[j].tag;
649    if (-abs(tid) == txid) parid=fEvent.gen[j].parent;
650   }
651
652   if (parid<0) continue;
653
654   AliTrack* tpar=evt->GetIdTrack(-abs(parid));
655
656   if (!tpar) continue;
657
658   tx->SetParentTrack(tpar);
659  }
660 }
661 ///////////////////////////////////////////////////////////////////////////
662 void IceF2k::PutRecoTracks()
663 {
664 // Get the reconstructed tracks from the F2000 file into the IcePack structure.
665 // Note : Reco tracks are given positive track id's in the event structure.
666 // This memberfunction is based on the original code by Adam Bouchta.
667
668  IceEvent* evt=(IceEvent*)GetMainObject();
669  if (!evt || fEvent.nfit<=0) return;
670
671  // Loop over all the tracks and add them to the current event
672  AliTrack t;
673  Double_t vec[3];
674  AliPosition r;
675  Ali3Vector p;
676  Int_t tid=0;
677  Int_t idpdg=0;
678  Int_t idf2k=0;
679  for (Int_t i=0; i<fEvent.nfit; i++)
680  {
681   t.Reset ();
682
683   // Beginpoint of the track
684   vec[0]=fEvent.rec[i].x;
685   vec[1]=fEvent.rec[i].y;
686   vec[2]=fEvent.rec[i].z;
687   r.SetPosition(vec,"car");
688   t.SetBeginPoint(r);
689
690   // Endpoint of the track
691   vec[0]+=fEvent.rec[i].length*fEvent.rec[i].px;
692   vec[1]+=fEvent.rec[i].length*fEvent.rec[i].py;
693   vec[2]+=fEvent.rec[i].length*fEvent.rec[i].pz;
694   r.SetPosition(vec,"car");
695   t.SetEndPoint(r);
696
697   // Momentum in GeV/c
698   if (fEvent.rec[i].e > 0)
699   {
700    vec[0]=fEvent.rec[i].e*fEvent.rec[i].px*1e-3;
701    vec[1]=fEvent.rec[i].e*fEvent.rec[i].py*1e-3;
702    vec[2]=fEvent.rec[i].e*fEvent.rec[i].pz*1e-3;
703   }
704   else // Give the track a nominal momentum of 1 GeV/c
705   {
706    vec[0]=fEvent.rec[i].px;
707    vec[1]=fEvent.rec[i].py;
708    vec[2]=fEvent.rec[i].pz;
709   }
710   p.SetVector (vec,"car");
711   t.Set3Momentum(p);
712
713   // Use the fit number as track id
714   tid=fEvent.rec[i].tag;
715   t.SetId(abs(tid));
716
717   idf2k=fEvent.rec[i].id;
718   idpdg=0;
719   if (idf2k>1000)
720   {
721    idpdg=idf2k+10000000;
722   }
723   else if (idf2k <= 48)
724   {
725    idpdg=fPdg->ConvertGeant3ToPdg(idf2k);
726   }
727   else
728   {
729    if (idf2k==201) idpdg=12;
730    if (idf2k==202) idpdg=14;
731    if (idf2k==203) idpdg=16;
732    if (idf2k==204) idpdg=-12;
733    if (idf2k==205) idpdg=-14;
734    if (idf2k==206) idpdg=-16;
735   }
736
737   t.SetParticleCode(idpdg);
738   t.SetName(fPdg->GetParticle(idpdg)->GetName());
739   t.SetTitle("RECO track");
740   t.SetMass(fPdg->GetParticle(idpdg)->Mass());
741   t.SetCharge(fPdg->GetParticle(idpdg)->Charge()/3.);
742
743   // Retrieve the various fit parameters for this track
744   AliSignal* fitdata=fFitdefs->GetIdHit(i);
745   for (Int_t jval=0; jval<fEvent.fresult[i].nval; jval++)
746   {
747    fitdata->SetSignal(fEvent.fresult[i].val[jval],jval+1);
748   }
749
750   // Store the various fit parameters for this track
751   t.SetFitDetails(fitdata);
752
753   // Store the various reco tracks as track hypotheses.
754   // A copy of the first reco track is entered as a new track instance
755   // into the event and all reco tracks (incl. the first one) are
756   // stored as hypotheses linked to this new reco track.
757   if (i==0)
758   {
759    evt->AddTrack(t);
760    AliTrack* tx=evt->GetTrack(evt->GetNtracks());
761    Int_t nrec=evt->GetNtracks(1);
762    tx->SetId(nrec+1);
763   }
764   AliTrack* tx=evt->GetTrack(evt->GetNtracks());
765   if (tx) tx->AddTrackHypothesis(t);
766  }
767 }
768 ///////////////////////////////////////////////////////////////////////////
769 void IceF2k::PutHits()
770 {
771 // Get the hit and waveform info from the F2000 file into the IcePack structure.
772 // This memberfunction is based on the original code by Adam Bouchta.
773
774  IceEvent* evt=(IceEvent*)GetMainObject();
775  if (!evt) return;
776
777  // Loop over all the hits and add them to the current event
778  IceAOM om;
779  AliSignal s;
780  s.SetSlotName("ADC",1);
781  s.SetSlotName("LE",2);
782  s.SetSlotName("TOT",3);
783  Int_t chan=0;
784  Int_t maxchan=800;
785  if (fOmdb) maxchan=fHeader.nch;
786  IceAOM* omx=0;
787  AliSignal* sx=0;
788  Int_t tid=0;
789  AliTrack* tx=0;
790  Float_t adc=0;
791  for (Int_t i=0; i<fEvent.nhits; i++)
792  {
793   chan=fEvent.h[i].ch+1;
794   if (chan>maxchan) continue; // Channels 9001, 9002 etc are trigger channels
795
796   // Get corresponding device from the current event structure  
797   omx=(IceAOM*)evt->GetIdDevice(chan);
798   if (!omx)
799   {
800    if (fOmdb)
801    {
802     omx=(IceAOM*)fOmdb->GetObject(chan,1);
803     evt->AddDevice(omx);
804    }
805    else
806    {
807     om.Reset(1);
808     om.SetUniqueID(chan);
809     evt->AddDevice(om);
810    }
811    omx=(IceAOM*)evt->GetIdDevice(chan);
812   }
813
814   if (!omx) continue;
815
816   s.Reset();
817   s.SetUniqueID(fEvent.h[i].id);
818   s.SetSignal(fEvent.h[i].amp,1);
819   s.SetSignal(fEvent.h[i].t,2);
820   s.SetSignal(fEvent.h[i].tot,3);
821
822   omx->AddHit(s);
823
824   sx=omx->GetHit(omx->GetNhits());
825   if (!sx) continue;
826
827   // ADC dependent TDC (de)calibration function for this hit
828   TF1* fcal=omx->GetCalFunction("LE");
829   TF1* fdecal=omx->GetDecalFunction("LE");
830   if (fcal) sx->SetCalFunction(fcal,2);
831   if (fdecal) sx->SetDecalFunction(fdecal,2);
832   fcal=sx->GetCalFunction(2);
833   fdecal=sx->GetDecalFunction(2);
834   adc=sx->GetSignal(1,-4);
835   if (adc>0)
836   {
837    if (fcal) fcal->SetParameter(3,adc);
838    if (fdecal) fdecal->SetParameter(3,adc);
839   }
840   else
841   {
842    if (fcal) fcal->SetParameter(3,1.e20);
843    if (fdecal) fdecal->SetParameter(3,1.e20);
844   }
845
846   // Bi-directional link between this hit and the track that caused the ADC value.
847   // This F2K info is probably only present for MC tracks.
848   tid=fEvent.h[i].ma;
849   if (tid > 0)
850   {
851    tx=evt->GetIdTrack(tid); // Reco tracks
852    if (!tx) tx=evt->GetIdTrack(-tid); // MC tracks
853    if (tx) sx->AddLink(tx);
854   }
855   else
856   {
857    if (tid == -2) s.SetNameTitle("N","Noise");
858    if (tid == -3) s.SetNameTitle("A","Afterpulse");
859   }
860  }
861
862  // Loop over all the waveforms and add the histo(s) to the corresponding OM's
863  TH1F histo;
864  Int_t nbins=0;
865  Float_t xlow=0;
866  Float_t xup=0;
867  TString hname;
868  for (Int_t iwf=0; iwf<fEvent.nwf; iwf++)
869  {
870   chan=fEvent.wf[iwf].om;
871   if (chan<=0 || chan>maxchan) continue; // Skip trigger channels
872
873   // Get corresponding device from the current event structure  
874   omx=(IceAOM*)evt->GetIdDevice(chan);
875   if (!omx)
876   {
877    if (fOmdb)
878    {
879     omx=(IceAOM*)fOmdb->GetObject(chan,1);
880     evt->AddDevice(omx);
881    }
882    else
883    {
884     om.Reset(1);
885     om.SetUniqueID(chan);
886     evt->AddDevice(om);
887    }
888    omx=(IceAOM*)evt->GetIdDevice(chan);
889   }
890
891   if (!omx) continue;
892
893   omx->SetSlotName("BASELINE",omx->GetNnames()+1);
894   omx->SetSignal(fEvent.wf[iwf].baseline,"BASELINE");
895
896   // Fill the waveform histogram
897   hname="OM";
898   hname+=chan;
899   hname+="-WF";
900   hname+=omx->GetNwaveforms()+1;
901
902   histo.Reset();
903   histo.SetName(hname.Data());
904   nbins=fEvent.wf[iwf].ndigi;
905   xlow=fEvent.wf[iwf].t_start;
906   xup=xlow+float(nbins)*fEvent.wf[iwf].t_bin;
907   histo.SetBins(nbins,xlow,xup);
908
909   for (Int_t jbin=1; jbin<=fEvent.wf[iwf].ndigi; jbin++)
910   {
911    histo.SetBinContent(jbin,fEvent.wf[iwf].digi[jbin-1]);
912   }
913
914   omx->SetWaveform(&histo,omx->GetNwaveforms()+1);
915  }
916
917  // Set bi-directional links between hits and reco track hypotheses.
918  // Note : Reco tracks are recognised by their positive id.
919  Int_t hid=0;
920  TObjArray* rectracks=evt->GetTracks(1);
921  for (Int_t jtk=0; jtk<rectracks->GetEntries(); jtk++)
922  {
923   tx=(AliTrack*)rectracks->At(jtk);
924   if (!tx) continue;
925   
926   for (Int_t jhyp=1; jhyp<=tx->GetNhypotheses(); jhyp++)
927   {
928    AliTrack* hypx=tx->GetTrackHypothesis(jhyp);
929    if (!hypx) continue;
930
931    // Loop over all combinations of F2K fits and used OM hits
932    for (Int_t k=0; k<fEvent.nfit_uses; k++)
933    {
934     if (fEvent.fit_uses[k].useid != hypx->GetId()) continue;
935     hid=fEvent.fit_uses[k].hitid;
936     sx=evt->GetIdHit(hid,"IceAOM");
937     if (sx) sx->AddLink(hypx);
938    }
939   }
940  }
941 }
942 ///////////////////////////////////////////////////////////////////////////