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