12-sep-2006 NvE Memberfunctions GetNslots and AddNamedSlot introduced and various
[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 // Note : Sometimes the filtering/reco process which produced the F2K file
33 //        may have introduced a shift (i.e. offset) in the hit times w.r.t.
34 //        the actual trigger time. The aim of this is to obtain the hit times
35 //        centered more or less around zero.
36 //        In case of real data, this is recorded in the F2K data itself and
37 //        as such will be taken automatically into account by this IceF2k
38 //        processor such that all times will be provided again unshifted.
39 //        In other words, all times will be w.r.t. the actual trigger time
40 //        as recorded in the trigger data device named "Trigger" in the IceEvent
41 //        structure.
42 //        In case of simulated data this shift is not available in the F2K data.
43 //        The offset denoted in the F2K record is related to the time of the
44 //        primary interaction to put it well ahead of the detector trigger.
45 //        This primary interaction time, however, is irrelevant for the
46 //        reconstruction of the recorded hit patterns. 
47 //        If a user had introduced a shift in producing the MC data,
48 //        very frequently (but not always) a value of -19000 is used.
49 //        For the IceF2k processing, the user can manually introduce a
50 //        time offset in case of MC data via the memberfunction SetMcToffset().
51 //        This user defined offset value will then be used to correct all
52 //        the hit times such that they will be provided again unshifted
53 //        w.r.t. the actual trigger time as recorded in the device named 
54 //        "Trigger" in the IceEvent structure.
55 //        By default the MC time offset is set to 0 in the constructor
56 //        of this class.
57 //
58 // Usage example :
59 // ---------------
60 //
61 // Note : This example creates automatically the ROOT output file, which
62 //        is the most user friendly way of running the conversion job. 
63 //        In the subdirectory /macros the example macro icef2k.cc provides
64 //        an example of how to create a ROOT output file yourself and passing
65 //        this file via a pointer to IceF2k. 
66 //
67 // gSystem->Load("ralice");
68 // gSystem->Load("icepack");
69 // gSystem->Load("iceconvert");
70 //
71 // IceF2k q("IceF2k","F2K to IcePack data structure conversion");
72 //
73 // // Limit the number of entries for testing
74 // q.SetMaxEvents(10);
75 //
76 // // Print frequency to produce a short summary print every printfreq events
77 // q.SetPrintFreq(1);
78 //
79 // // Split level for the output structures
80 // q.SetSplitLevel(0);
81 //
82 // // Buffer size for the output structures
83 // q.SetBufferSize(32000);
84 //
85 // // The F2K input filename(s)
86 // q.AddInputFile("run7825.f2k");
87 //
88 // // Output file for the event structures
89 // q.SetOutputFile("events.root");
90 //
91 // ///////////////////////////////////////////////////////////////////
92 // // Here the user can specify his/her sub-tasks to be executed
93 // // on an event-by-event basis after the IceEvent structure
94 // // has been filled and before the data is written out.
95 // // Sub-tasks (i.e. a user classes derived from TTask) are entered
96 // // as follows :
97 // //
98 // //    MyXtalk task1("task1","Cross talk correction");
99 // //    MyClean task2("task2","Hit cleaning");
100 // //    q.Add(&task1);
101 // //    q.Add(&task2);
102 // //
103 // // The sub-tasks will be executed in the order as they are entered.
104 // ///////////////////////////////////////////////////////////////////
105 //
106 // // Perform the conversion and execute subtasks (if any)
107 // // on an event-by-event basis
108 // q.ExecuteJob();
109 //
110 // // Select various objects to be added to the output file
111 //
112 // TFile* ofile=q.GetOutputFile();
113 // 
114 // if (ofile)
115 // {
116 //  ofile->cd(); // Switch to the output file directory
117 //
118 //  AliObjMatrix* omdb=q.GetOMdbase();
119 //  if (omdb) omdb->Write();
120 //
121 //  AliDevice* fitdefs=q.GetFitdefs();
122 //  if (fitdefs) fitdefs->Write();
123 //
124 //  TDatabasePDG* pdg=q.GetPDG();
125 //  if (pdg) pdg->Write();
126 //
127 //  // Flush additional objects to the output file.
128 //  // The output file is not explicitly closed here
129 //  // to allow interactive investigation of the data tree
130 //  // when this macro is run in an interactive ROOT/CINT session.
131 //  ofile->Write();
132 // }
133 //
134 //--- Author: Nick van Eijndhoven 11-mar-2005 Utrecht University
135 //- Modified: NvE $Date$ Utrecht University
136 ///////////////////////////////////////////////////////////////////////////
137  
138 #include "IceF2k.h"
139 #include "Riostream.h"
140
141 ClassImp(IceF2k) // Class implementation to enable ROOT I/O
142
143 IceF2k::IceF2k(const char* name,const char* title) : AliJob(name,title)
144 {
145 // Default constructor.
146 // By default maxevent=-1, split=0, bsize=32000, printfreq=1.
147
148  fSplit=0;
149  fBsize=32000;
150  fMaxevt=-1;
151  fPrintfreq=1;
152  fInfiles=0;
153  fOutfile=0;
154
155  fPdg=0;
156  fOmdb=0;
157  fFitdefs=0;
158  fTrigdefs=0;
159  fToffset=0;
160  fMctoffset=0;
161  fMctracks=3;
162 }
163 ///////////////////////////////////////////////////////////////////////////
164 IceF2k::~IceF2k()
165 {
166 // Default destructor.
167
168  if (fInfiles)
169  {
170   delete fInfiles;
171   fInfiles=0;
172  }
173
174  if (fPdg)
175  {
176   delete fPdg;
177   fPdg=0;
178  }
179
180  if (fOmdb)
181  {
182   delete fOmdb;
183   fOmdb=0;
184  }
185
186  if (fFitdefs)
187  {
188   delete fFitdefs;
189   fFitdefs=0;
190  }
191
192  if (fTrigdefs)
193  {
194   delete fTrigdefs;
195   fTrigdefs=0;
196  }
197 }
198 ///////////////////////////////////////////////////////////////////////////
199 void IceF2k::SetMaxEvents(Int_t n)
200 {
201 // Set the maximum number of events to be processed.
202 // n=-1 implies processing of the complete input file, which is the default
203 // initialisation in the constructor.
204  fMaxevt=n;
205 }
206 ///////////////////////////////////////////////////////////////////////////
207 void IceF2k::SetPrintFreq(Int_t f)
208 {
209 // Set the printfrequency to produce info every f events.
210 // f=1 is the default initialisation in the constructor.
211  if (f>=0) fPrintfreq=f;
212 }
213 ///////////////////////////////////////////////////////////////////////////
214 void IceF2k::SetSplitLevel(Int_t split)
215 {
216 // Set the split level for the ROOT data file.
217 // split=0 is the default initialisation in the constructor.
218  if (split>=0) fSplit=split;
219 }
220 ///////////////////////////////////////////////////////////////////////////
221 void IceF2k::SetBufferSize(Int_t bsize)
222 {
223 // Set the buffer size for the ROOT data file.
224 // bsize=32000 is the default initialisation in the constructor.
225  if (bsize>=0) fBsize=bsize;
226 }
227 ///////////////////////////////////////////////////////////////////////////
228 void IceF2k::SetMcToffset(Float_t toffset)
229 {
230 // Set a user defined time offset for Monte Carlo data.
231 // A very frequently (but not always) used value is -19000.
232 // See the introductory docs of this class for further details.
233  fMctoffset=toffset;
234 }
235 ///////////////////////////////////////////////////////////////////////////
236 void IceF2k::SelectMcTracks(Int_t mode)
237 {
238 // User selection of MC tracks to be stored in the event structure.
239 //
240 // mode = 0 : No MC tracks are stored
241 //        1 : Only muon and muon-neutrino MC tracks are stored
242 //        2 : All lepton MC tracks are stored
243 //        3 : All MC tracks (incl. brems, pairprod etc...) are stored
244 //
245 // By default mode=3 is set in the constructor of this class.
246
247  if (mode<0 || mode >3) return;
248  fMctracks=mode;
249 }
250 ///////////////////////////////////////////////////////////////////////////
251 void IceF2k::SetInputFile(TString name)
252 {
253 // Set the name of the F2K input file.
254 // This function has become obsolete but is kept for backward compatibility.
255 // The user is advised to use AddInputFile() instead, which allows processing
256 // of multiple F2K input files.
257 // This function will reset the list of all F2K input files and put the specified
258 // filename at the first position.
259 // Additional F2K input files can be specified via AddInputFile().
260
261  if (fInfiles) delete fInfiles;
262
263  fInfiles=new TObjArray();
264  fInfiles->SetOwner();
265
266  TObjString* s=new TObjString();
267  s->SetString(name);
268  fInfiles->Add(s);
269 }
270 ///////////////////////////////////////////////////////////////////////////
271 void IceF2k::AddInputFile(TString name)
272 {
273 // Add the name of this F2K input file to the list to be processed.
274
275  if (!fInfiles)
276  {
277   fInfiles=new TObjArray();
278   fInfiles->SetOwner();
279  }
280
281  TObjString* s=new TObjString();
282  s->SetString(name);
283  fInfiles->Add(s);
284 }
285 ///////////////////////////////////////////////////////////////////////////
286 void IceF2k::SetOutputFile(TFile* ofile)
287 {
288 // Set the output file for the ROOT data.
289  if (fOutfile) delete fOutfile;
290  fOutfile=ofile;
291 }
292 ///////////////////////////////////////////////////////////////////////////
293 void IceF2k::SetOutputFile(TString name)
294 {
295 // Create the output file for the ROOT data.
296  if (fOutfile) delete fOutfile;
297  fOutfile=new TFile(name.Data(),"RECREATE","F2K data in IceEvent structure");
298 }
299 ///////////////////////////////////////////////////////////////////////////
300 TFile* IceF2k::GetOutputFile()
301 {
302 // Provide pointer to the ROOT output file.
303  return fOutfile;
304 }
305 ///////////////////////////////////////////////////////////////////////////
306 TDatabasePDG* IceF2k::GetPDG()
307 {
308 // Provide pointer to the PDG database
309  return fPdg;
310 }
311 ///////////////////////////////////////////////////////////////////////////
312 AliObjMatrix* IceF2k::GetOMdbase()
313 {
314 // Provide pointer to the OM geometry, calib. etc... database
315  return fOmdb;
316 }
317 ///////////////////////////////////////////////////////////////////////////
318 AliDevice* IceF2k::GetFitdefs()
319 {
320 // Provide pointer to the fit definitions
321  return fFitdefs;
322 }
323 ///////////////////////////////////////////////////////////////////////////
324 AliDevice* IceF2k::GetTrigdefs()
325 {
326 // Provide pointer to the trigger definitions
327  return fTrigdefs;
328 }
329 ///////////////////////////////////////////////////////////////////////////
330 void IceF2k::Exec(Option_t* opt)
331 {
332 // Job to loop over the specified number of events and convert the 
333 // F2K data into the IceEvent structure.
334 // If maxevents<0 (default) all the entries of the input file
335 // will be processed.
336 // Every "printfreq" events a short event summary will be printed.
337 // The default value is printfreq=1.
338 // The output will be written on a standard output tree named "T".
339 //
340 // Notes :
341 // -------
342 // 1) This class is derived from AliJob, allowing a task based processing.
343 //    After the conversion of an F2K event into an IceEvent structure,
344 //    the processing of all available sub-tasks (if any) is invoked.
345 //    This provides an event-by-event (sub)task processing before the
346 //    final data structures are written out.
347 // 2) The main object in this job environment is an IceEvent* pointer.
348
349  if (!fInfiles)
350  {
351   cout << " *IceF2k Exec* No data input file(s) specified." << endl;
352   return;
353  }
354
355  Int_t ninfiles=fInfiles->GetEntries();
356  if (!ninfiles)
357  {
358   cout << " *IceF2k Exec* No data input file(s) specified." << endl;
359   return;
360  }
361
362  TTree* otree=0;
363  if (fOutfile)
364  {
365   otree=new TTree("T","F2K Data converted to IceEvent structures");
366   otree->SetDirectory(fOutfile);
367  }
368
369  IceEvent* evt=new IceEvent();
370  evt->SetTrackCopy(1);
371  evt->SetDevCopy(1);
372
373  // Branch in the tree for the event structure
374  if (otree) otree->Branch("IceEvent","IceEvent",&evt,fBsize,fSplit); 
375
376  // Create the particle database and extend it with some F2000 specific definitions
377  if (!fPdg) fPdg=new TDatabasePDG();
378  Double_t me=fPdg->GetParticle(11)->Mass();
379  fPdg->AddParticle("brems"   ,"brems"   ,0,1,0,0,"none",10001001,0,0);
380  fPdg->AddParticle("deltae"  ,"deltae"  ,me,1,0,-3,"Lepton",10001002,0,0);
381  fPdg->AddParticle("pairprod","pairprod",0,1,0,0,"none",10001003,0,0);
382  fPdg->AddParticle("nucl_int","nucl_Int",0,1,0,0,"none",10001004,0,0);
383  fPdg->AddParticle("mu_pair" ,"mu_pair" ,0,1,0,0,"none",10001005,0,0);
384  fPdg->AddParticle("hadrons" ,"hadrons" ,0,1,0,0,"none",10001006,0,0);
385  fPdg->AddParticle("fiberlaser","fiberlaser",0,1,0,0,"none",10002100,0,0);
386  fPdg->AddParticle("n2laser"   ,"n2laser"   ,0,1,0,0,"none",10002101,0,0);
387  fPdg->AddParticle("yaglaser"  ,"yaglaser"  ,0,1,0,0,"none",10002201,0,0);
388  fPdg->AddParticle("z_primary","z_primary",0,1,0,0,"none",10003000,0,0);
389  fPdg->AddParticle("a_primary","a_primary",0,1,0,0,"none",10003500,0,0);
390
391  // Initialise the job working environment
392  SetMainObject(evt);
393  if (fOutfile)
394  {
395   AddObject(fOutfile);
396   AddObject(otree);
397  }
398
399  TString inputfile;
400
401  cout << " ***" << endl;
402  cout << " *** Start processing of job " << GetName() << " ***" << endl;
403  cout << " ***" << endl;
404  for (Int_t i=0; i<ninfiles; i++)
405  {
406   TObjString* sx=(TObjString*)fInfiles->At(i);
407   if (!sx) continue;
408   inputfile=sx->GetString(); 
409   cout << " F2K input file : " << inputfile.Data() << endl;
410  }
411  cout << " Maximum number of events to be processed : " << fMaxevt << endl;
412  cout << " Print frequency : " << fPrintfreq << endl;
413  if (fOutfile)
414  {
415   cout << " ROOT output file : " << fOutfile->GetName() << endl;
416   cout << " Output characteristics : splitlevel = " << fSplit << " buffersize = " << fBsize << endl;
417  }
418
419  ListEnvironment();
420
421  Int_t nevt=0;
422  for (Int_t ifile=0; ifile<ninfiles; ifile++)
423  {
424   TObjString* sx=(TObjString*)fInfiles->At(ifile);
425   if (!sx) continue;
426
427   inputfile=sx->GetString(); 
428   if (inputfile=="") continue;
429
430   // Open the input file in the default ascii format (autodetection) for reading 
431   fInput=rdmc_mcopen(inputfile.Data(),"r",RDMC_DEFAULT_ASCII_F);
432
433   if (!fInput)
434   {
435    cout << " *IceF2k Exec* No input file found with name : " << inputfile.Data() << endl;
436    continue;
437   }
438
439   // Initialise the event structure 
440   rdmc_init_mevt(&fEvent);
441
442   // Read the file header information
443   rdmc_rarr(fInput,&fHeader);
444
445   // Fill the database with geometry, calib. etc... parameters
446   // for all the devices
447   FillOMdbase();
448
449   // Set the fit definitions according to the F2000 header info
450   SetFitdefs();
451
452   // Set the trigger definitions according to the F2000 header info
453   SetTrigdefs();
454  
455   while (!rdmc_revt(fInput,&fHeader,&fEvent))
456   {
457    if (fMaxevt>-1 && nevt>=fMaxevt) break;
458
459    // Reset the complete Event structure
460    evt->Reset();
461
462    evt->SetRunNumber(fEvent.nrun);
463    evt->SetEventNumber(fEvent.enr);
464    evt->SetMJD(fEvent.mjd,fEvent.secs,fEvent.nsecs);
465
466    // Take trigger offset into account which might have been
467    // introduced during the filtering process.
468    // For simulated data this will be treated separately in PutMcTracks().
469    fToffset=fEvent.t_offset;
470
471    PutTrigger();
472
473    PutMcTracks();
474
475    PutRecoTracks();
476
477    PutHits();
478
479    // Invoke all available sub-tasks (if any)
480    CleanTasks();
481    ExecuteTasks(opt);
482
483    if (fPrintfreq)
484    {
485     if (!(nevt%fPrintfreq)) evt->HeaderData();
486    }
487
488    // Write the complete structure to the output Tree
489    if (otree) otree->Fill();
490
491    // Update event counter
492    nevt++;
493   }
494   if (fMaxevt>-1 && nevt>=fMaxevt) break;
495  }
496
497  // Flush possible memory resident data to the output file
498  if (fOutfile) fOutfile->Write();
499
500  // Remove the IceEvent object from the environment
501  // and delete it as well
502  if (evt)
503  {
504   RemoveObject(evt);
505   delete evt;
506  }
507 }
508 ///////////////////////////////////////////////////////////////////////////
509 void IceF2k::FillOMdbase()
510 {
511 // Fill the database with geometry, calib. etc... parameters 
512 // for all the devices.
513
514  if (fHeader.nch<=0)
515  {
516   if (fOmdb)
517   {
518    delete fOmdb;
519    fOmdb=0;
520   }
521   return;
522  }
523
524  Int_t adccal=fHeader.is_calib.adc;
525  Int_t tdccal=fHeader.is_calib.tdc;
526  Int_t totcal=fHeader.is_calib.tot;
527
528  TF1 fadccal("fadccal","(x-[1])*[0]");
529  TF1 fadcdecal("fadcdecal","(x/[0])+[1]");
530  fadccal.SetParName(0,"BETA-ADC");
531  fadccal.SetParName(1,"PED-ADC");
532  fadcdecal.SetParName(0,"BETA-ADC");
533  fadcdecal.SetParName(1,"PED-ADC");
534
535  TF1 ftdccal("ftdccal","(x*[0])-[1]-([0]-1.)*32767.-[2]/sqrt([3])");
536  TF1 ftdcdecal("ftdcdecal","(x+([0]-1.)*32767.+[1]+[2]/sqrt([3]))/[0]");
537  ftdccal.SetParName(0,"BETA-TDC");
538  ftdccal.SetParName(1,"T0");
539  ftdccal.SetParName(2,"ALPHA-TDC");
540  ftdccal.SetParName(3,"ADC-SLEW");
541  ftdcdecal.SetParName(0,"BETA-TDC");
542  ftdcdecal.SetParName(1,"T0");
543  ftdcdecal.SetParName(2,"ALPHA-TDC");
544  ftdcdecal.SetParName(3,"ADC-SLEW");
545
546  TF1 ftotcal("ftotcal","x*[0]");
547  TF1 ftotdecal("ftotdecal","x/[0]");
548  ftotcal.SetParName(0,"BETA-TOT");
549  ftotdecal.SetParName(0,"BETA-TOT");
550
551  if (fOmdb)
552  {
553   fOmdb->Reset();
554  }
555  else
556  {
557   fOmdb=new AliObjMatrix();
558   fOmdb->SetNameTitle("OMDBASE","The OM geometry, calib. etc... database");
559   fOmdb->SetOwner();
560  }
561
562  IceAOM* dev=0;
563  Double_t pos[3]={0,0,0};
564  for (Int_t i=0; i<fHeader.nch; i++)
565  {
566   dev=new IceAOM();
567   dev->SetUniqueID(i+1);
568   // Slots to hold the various (de)calibration functions  
569   dev->SetSlotName("ADC",1);
570   dev->SetSlotName("LE",2);
571   dev->SetSlotName("TOT",3);
572   // Slots to hold hardware parameters
573   dev->SetSlotName("TYPE",4);
574   dev->SetSlotName("ORIENT",5);
575   dev->SetSlotName("THRESH",6);
576   dev->SetSlotName("SENSIT",7);
577
578   pos[0]=fHeader.x[i];
579   pos[1]=fHeader.y[i];
580   pos[2]=fHeader.z[i];
581   dev->SetPosition(pos,"car");
582
583   fadccal.SetParameter(0,fHeader.cal[i].beta_a);
584   fadccal.SetParameter(1,fHeader.cal[i].ped);
585   fadcdecal.SetParameter(0,fHeader.cal[i].beta_a);
586   if (!fHeader.cal[i].beta_a) fadcdecal.SetParameter(0,1);
587   fadcdecal.SetParameter(1,fHeader.cal[i].ped);
588
589   ftdccal.SetParameter(0,fHeader.cal[i].beta_t);
590   ftdccal.SetParameter(1,fHeader.cal[i].t_0);
591   ftdccal.SetParameter(2,fHeader.cal[i].alpha_t);
592   ftdccal.SetParameter(3,1.e20);
593   ftdcdecal.SetParameter(0,fHeader.cal[i].beta_t);
594   if (!fHeader.cal[i].beta_t) ftdcdecal.SetParameter(0,1);
595   ftdcdecal.SetParameter(1,fHeader.cal[i].t_0);
596   ftdcdecal.SetParameter(2,fHeader.cal[i].alpha_t);
597   ftdcdecal.SetParameter(3,1.e20);
598
599   ftotcal.SetParameter(0,fHeader.cal[i].beta_tot);
600   ftotdecal.SetParameter(0,fHeader.cal[i].beta_tot);
601   if (!fHeader.cal[i].beta_tot) ftotdecal.SetParameter(0,1);
602
603   if (adccal)
604   {
605    dev->SetDecalFunction(&fadcdecal,1);
606   }
607   else
608   {
609    dev->SetCalFunction(&fadccal,1);
610   }
611
612   if (tdccal)
613   {
614    dev->SetDecalFunction(&ftdcdecal,2);
615   }
616   else
617   {
618    dev->SetCalFunction(&ftdccal,2);
619   }
620
621   if (totcal)
622   {
623    dev->SetDecalFunction(&ftotdecal,3);
624   }
625   else
626   {
627    dev->SetCalFunction(&ftotcal,3);
628   }
629
630   dev->SetSignal(fHeader.type[i],4);
631   dev->SetSignal((Float_t)fHeader.costh[i],5);
632   dev->SetSignal(fHeader.thresh[i],6);
633   dev->SetSignal(fHeader.sensit[i],7);
634
635   fOmdb->EnterObject(i+1,1,dev);
636  }
637 }
638 ///////////////////////////////////////////////////////////////////////////
639 void IceF2k::SetFitdefs()
640 {
641 // Obtain the names of the variables for each fit procedure from the
642 // f2000 header. Each different fit procedure is then stored as a separate
643 // "hit" of an AliDevice object and the various fit variables are stored
644 // as separate signal slots of the corresponding "hit".
645 // Via the GetFitdefs() memberfunction this AliDevice object can be
646 // retrieved and stored in the ROOT output file if wanted.
647 // The name of the object is FitDefinitions and the stored data can be
648 // inspected via the AliDevice::Data() memberfunction and looks as follows :
649 //
650 //  *AliDevice::Data* Id :0 Name : FitDefinitions
651 //    Position Vector in car coordinates : 0 0 0
652 //    Err. in car coordinates : 0 0 0
653 //  The following 8 hits are registered :
654 //  *AliSignal::Data* Id :0
655 //    Position Vector in car coordinates : 0 0 0
656 //    Err. in car coordinates : 0 0 0
657 //    Owned by device : AliDevice Name : FitDefinitions
658 //    Slot : 1 Signal value : 0 name : id
659 //    Slot : 2 Signal value : 0 name : rchi2
660 //    Slot : 3 Signal value : 0 name : prob
661 //    Slot : 4 Signal value : 0 name : sigth
662 //    Slot : 5 Signal value : 0 name : covmin
663 //    Slot : 6 Signal value : 0 name : covmax
664 //    Slot : 7 Signal value : 0 name : cutflag
665 //    Slot : 8 Signal value : 0 name : chi2
666 //  *AliSignal::Data* Id :1
667 //    Position Vector in car coordinates : 0 0 0
668 //    Err. in car coordinates : 0 0 0
669 //    Owned by device : AliDevice Name : FitDefinitions
670 //    Slot : 1 Signal value : 0 name : id
671 //    Slot : 2 Signal value : 0 name : rchi2
672 //    Slot : 3 Signal value : 0 name : prob
673 // etc....  
674 //
675 // This memberfunction is based on the original idea/code by Adam Bouchta.
676
677  if (fHeader.n_fit<=0)
678  {
679   if (fFitdefs)
680   {
681    delete fFitdefs;
682    fFitdefs=0;
683   }
684   return;
685  }
686
687  if (fFitdefs)
688  {
689   fFitdefs->Reset(1);
690  }
691  else
692  {
693   fFitdefs=new AliDevice();
694  }
695
696  fFitdefs->SetName("FitDefinitions");
697  fFitdefs->SetHitCopy (1);
698
699  AliSignal s;
700  s.Reset();
701
702  for (Int_t i=0; i<fHeader.n_fit; i++)
703  {
704   s.SetUniqueID(fHeader.def_fit[i].id);
705   s.SetName(TString(fHeader.def_fit[i].tag));
706
707   for (Int_t j=0; j<fHeader.def_fit[i].nwords; j++)
708   {
709    s.SetSlotName(TString(fHeader.def_fit[i].words[j]),j+1);
710    s.SetSignal(0,j+1);
711   }
712
713   fFitdefs->AddHit(s);
714   s.Reset(1);
715  }
716 }
717 ///////////////////////////////////////////////////////////////////////////
718 void IceF2k::SetTrigdefs()
719 {
720 // Obtain the names of the variables for each trigger procedure from the
721 // f2000 header. Each different trigger procedure is then stored as a separate
722 // "hit" of an AliDevice object and the various trigger variables are stored
723 // as separate signal slots of the corresponding "hit".
724 // Via the GetFitdefs() memberfunction this AliDevice object can be
725 // retrieved and stored in the ROOT output file if wanted.
726 // The name of the object is TrigDefinitions and the stored data can be
727 // inspected via the AliDevice::Data() memberfunction and looks as follows :
728 //
729 //  *AliDevice::Data* Id : 0 Name : TrigDefinitions
730 //    Position Vector in car (rad) coordinates : 0 0 0
731 //    Err. in car (rad) coordinates : 0 0 0
732 //  The following 9 hits are registered : 
733 //  *AliSignal::Data* Id : 1 Name : main
734 //    Position Vector in car (rad) coordinates : 0 0 0
735 //    Err. in car (rad) coordinates : 0 0 0
736 //    Owned by device : AliDevice Id : 0 Name : TrigDefinitions
737 //     Slot : 1 Signal value : 0 name : trig_pulse_le
738 //     Slot : 2 Signal value : 0 name : trig_pulse_tot
739 //     Slot : 3 Signal value : 0 name : regi_flag
740 //   *AliSignal::Data* Id : 2 Name : amaa
741 //     Position Vector in car (rad) coordinates : 0 0 0
742 //     Err. in car (rad) coordinates : 0 0 0
743 //     Owned by device : AliDevice Id : 0 Name : TrigDefinitions
744 //     Slot : 1 Signal value : 0 name : trig_pulse_le
745 //     Slot : 2 Signal value : 0 name : trig_pulse_tot
746 //     Slot : 3 Signal value : 0 name : regi_flag
747 //   *AliSignal::Data* Id : 3 Name : amab10
748 //     Position Vector in car (rad) coordinates : 0 0 0
749 //     Err. in car (rad) coordinates : 0 0 0
750 //     Owned by device : AliDevice Id : 0 Name : TrigDefinitions
751 //     Slot : 1 Signal value : 0 name : trig_pulse_le
752 //     Slot : 2 Signal value : 0 name : trig_pulse_tot
753 //     Slot : 3 Signal value : 0 name : regi_flag
754 // etc....  
755
756  if (fHeader.n_trigger<=0)
757  {
758   if (fTrigdefs)
759   {
760    delete fTrigdefs;
761    fTrigdefs=0;
762   }
763   return;
764  }
765
766  if (fTrigdefs)
767  {
768   fTrigdefs->Reset(1);
769  }
770  else
771  {
772   fTrigdefs=new AliDevice();
773  }
774
775  fTrigdefs->SetName("TrigDefinitions");
776  fTrigdefs->SetHitCopy (1);
777
778  AliSignal s;
779  s.Reset();
780
781  for (Int_t i=0; i<fHeader.n_trigger; i++)
782  {
783   s.SetUniqueID(fHeader.def_trig[i].id);
784   s.SetName(TString(fHeader.def_trig[i].tag));
785
786   for (Int_t j=0; j<fHeader.def_trig[i].nwords; j++)
787   {
788    s.SetSlotName(TString(fHeader.def_trig[i].words[j]),j+1);
789    s.SetSignal(0,j+1);
790   }
791
792   fTrigdefs->AddHit(s);
793   s.Reset(1);
794  }
795 }
796 ///////////////////////////////////////////////////////////////////////////
797 void IceF2k::PutMcTracks()
798 {
799 // Get the MC tracks from the F2000 file into the IcePack structure.
800 // Note : MC tracks are given negative track id's in the event structure.
801 // This memberfunction is based on the original code by Adam Bouchta.
802
803  IceEvent* evt=(IceEvent*)GetMainObject();
804  if (!evt || fEvent.ntrack<=0) return;
805
806  // User defined trigger offset in case of simulated data.
807  // The offset in the F2K file is meant to put the primary interaction
808  // well ahead of the detector trigger.
809  // See the introductory docs of this IceF2k class for further details.
810  fToffset=fMctoffset;
811
812  if (!fMctracks) return;
813
814  // Loop over all the tracks and add them to the current event
815  AliTrack t;
816  Double_t vec[3];
817  AliPosition r;
818  Ali3Vector p;
819  Int_t tid=0;
820  Int_t idpdg=0;
821  Int_t idf2k=0;
822  for (Int_t i=0; i<fEvent.ntrack; i++)
823  {
824   t.Reset ();
825
826   // Beginpoint of the track
827   vec[0]=fEvent.gen[i].x;
828   vec[1]=fEvent.gen[i].y;
829   vec[2]=fEvent.gen[i].z;
830   r.SetPosition(vec,"car");
831   t.SetBeginPoint(r);
832
833   // Endpoint of the track
834   vec[0]+=fEvent.gen[i].length*fEvent.gen[i].px;
835   vec[1]+=fEvent.gen[i].length*fEvent.gen[i].py;
836   vec[2]+=fEvent.gen[i].length*fEvent.gen[i].pz;
837   r.SetPosition(vec,"car");
838   t.SetEndPoint(r);
839
840   // Momentum in GeV/c
841   vec[0]=fEvent.gen[i].e*fEvent.gen[i].px*1e-3;
842   vec[1]=fEvent.gen[i].e*fEvent.gen[i].py*1e-3;
843   vec[2]=fEvent.gen[i].e*fEvent.gen[i].pz*1e-3;
844   p.SetVector (vec,"car");
845   t.Set3Momentum(p);
846
847   // MC tracks are indicated by negative track id's
848   tid=fEvent.gen[i].tag;
849   t.SetId(-abs(tid));
850
851   idf2k=fEvent.gen[i].id;
852   idpdg=0;
853   if (idf2k>1000)
854   {
855    idpdg=idf2k+10000000;
856   }
857   else if (idf2k <= 48)
858   {
859    idpdg=fPdg->ConvertGeant3ToPdg(idf2k);
860   }
861   else
862   {
863    if (idf2k==201) idpdg=12;
864    if (idf2k==202) idpdg=14;
865    if (idf2k==203) idpdg=16;
866    if (idf2k==204) idpdg=-12;
867    if (idf2k==205) idpdg=-14;
868    if (idf2k==206) idpdg=-16;
869   }
870
871   // Check for the user selected MC track storage
872   if (fMctracks==1) // Store only muon and muon-neutrino tracks
873   {
874    if (abs(idpdg)!=13 && abs(idpdg)!=14) continue;
875   }
876   else if (fMctracks==2) // Store all lepton tracks
877   {
878    if (abs(idpdg)<11 || abs(idpdg)>16) continue;
879   }
880
881   t.SetParticleCode(idpdg);
882   t.SetName(fPdg->GetParticle(idpdg)->GetName());
883   t.SetTitle("MC track");
884   t.SetMass(fPdg->GetParticle(idpdg)->Mass());
885   t.SetCharge(fPdg->GetParticle(idpdg)->Charge()/3.);
886
887   evt->AddTrack(t);
888  }
889
890  // Create the pointers to each particle's parent particle.
891  Int_t txid=0;
892  Int_t parid=0;
893  for (Int_t itk=1; itk<=evt->GetNtracks (); itk++)
894  {
895   AliTrack* tx=evt->GetTrack(itk);
896
897   if (!tx) continue;
898
899   txid=tx->GetId();
900
901   parid=-1;
902   for (Int_t j=0; j<fEvent.ntrack; j++)
903   {
904    tid=fEvent.gen[j].tag;
905    if (-abs(tid) == txid) parid=fEvent.gen[j].parent;
906   }
907
908   if (parid<0) continue;
909
910   AliTrack* tpar=evt->GetIdTrack(-abs(parid));
911
912   if (!tpar) continue;
913
914   tx->SetParentTrack(tpar);
915  }
916 }
917 ///////////////////////////////////////////////////////////////////////////
918 void IceF2k::PutRecoTracks()
919 {
920 // Get the reconstructed tracks from the F2000 file into the IcePack structure.
921 // Note : Reco tracks are given positive track id's in the event structure.
922 // This memberfunction is based on the original code by Adam Bouchta.
923
924  IceEvent* evt=(IceEvent*)GetMainObject();
925  if (!evt || fEvent.nfit<=0) return;
926
927  // Loop over all the tracks and add them to the current event
928  AliTrack t;
929  Double_t vec[3];
930  AliPosition r;
931  Ali3Vector p;
932  Int_t tid=0;
933  Int_t idpdg=0;
934  Int_t idf2k=0;
935  for (Int_t i=0; i<fEvent.nfit; i++)
936  {
937   t.Reset ();
938
939   // Beginpoint of the track
940   vec[0]=fEvent.rec[i].x;
941   vec[1]=fEvent.rec[i].y;
942   vec[2]=fEvent.rec[i].z;
943   r.SetPosition(vec,"car");
944   t.SetBeginPoint(r);
945
946   // Endpoint of the track
947   vec[0]+=fEvent.rec[i].length*fEvent.rec[i].px;
948   vec[1]+=fEvent.rec[i].length*fEvent.rec[i].py;
949   vec[2]+=fEvent.rec[i].length*fEvent.rec[i].pz;
950   r.SetPosition(vec,"car");
951   t.SetEndPoint(r);
952
953   // Momentum in GeV/c
954   if (fEvent.rec[i].e > 0)
955   {
956    vec[0]=fEvent.rec[i].e*fEvent.rec[i].px*1e-3;
957    vec[1]=fEvent.rec[i].e*fEvent.rec[i].py*1e-3;
958    vec[2]=fEvent.rec[i].e*fEvent.rec[i].pz*1e-3;
959   }
960   else // Give the track a nominal momentum of 1 GeV/c
961   {
962    vec[0]=fEvent.rec[i].px;
963    vec[1]=fEvent.rec[i].py;
964    vec[2]=fEvent.rec[i].pz;
965   }
966   p.SetVector (vec,"car");
967   t.Set3Momentum(p);
968
969   // Use the fit number as track id
970   tid=fEvent.rec[i].tag;
971   t.SetId(abs(tid));
972
973   idf2k=fEvent.rec[i].id;
974   idpdg=0;
975   if (idf2k>1000)
976   {
977    idpdg=idf2k+10000000;
978   }
979   else if (idf2k <= 48)
980   {
981    idpdg=fPdg->ConvertGeant3ToPdg(idf2k);
982   }
983   else
984   {
985    if (idf2k==201) idpdg=12;
986    if (idf2k==202) idpdg=14;
987    if (idf2k==203) idpdg=16;
988    if (idf2k==204) idpdg=-12;
989    if (idf2k==205) idpdg=-14;
990    if (idf2k==206) idpdg=-16;
991   }
992
993   t.SetParticleCode(idpdg);
994   t.SetNameTitle("Sieglinde","RECO track");
995   t.SetMass(fPdg->GetParticle(idpdg)->Mass());
996   t.SetCharge(fPdg->GetParticle(idpdg)->Charge()/3.);
997
998   // Retrieve the various fit parameters for this track
999   AliSignal* fitdata=fFitdefs->GetIdHit(i);
1000   for (Int_t jval=0; jval<fEvent.fresult[i].nval; jval++)
1001   {
1002    fitdata->SetSignal(fEvent.fresult[i].val[jval],jval+1);
1003   }
1004
1005   // Store the various fit parameters for this track
1006   t.SetFitDetails(fitdata);
1007
1008   // Store the various reco tracks as track hypotheses.
1009   // A copy of the first reco track is entered as a new track instance
1010   // into the event and all reco tracks (incl. the first one) are
1011   // stored as hypotheses linked to this new reco track.
1012   if (i==0)
1013   {
1014    evt->AddTrack(t);
1015    AliTrack* tx=evt->GetTrack(evt->GetNtracks());
1016    Int_t nrec=evt->GetNtracks(1);
1017    tx->SetId(nrec+1);
1018   }
1019   AliTrack* tx=evt->GetTrack(evt->GetNtracks());
1020   if (tx) tx->AddTrackHypothesis(t);
1021  }
1022 }
1023 ///////////////////////////////////////////////////////////////////////////
1024 void IceF2k::PutHits()
1025 {
1026 // Get the hit and waveform info from the F2000 file into the IcePack structure.
1027 // This memberfunction is based on the original code by Adam Bouchta.
1028
1029  IceEvent* evt=(IceEvent*)GetMainObject();
1030  if (!evt) return;
1031
1032  // Loop over all the hits and add them to the current event
1033  IceAOM om;
1034  AliSignal s;
1035  s.SetSlotName("ADC",1);
1036  s.SetSlotName("LE",2);
1037  s.SetSlotName("TOT",3);
1038  Int_t chan=0;
1039  Int_t maxchan=800;
1040  if (fOmdb) maxchan=fHeader.nch;
1041  IceAOM* omx=0;
1042  AliSignal* sx=0;
1043  Int_t tid=0;
1044  AliTrack* tx=0;
1045  Float_t adc=0;
1046  Float_t adcfirst=0; // Adc value of the first hit of an OM
1047  for (Int_t i=0; i<fEvent.nhits; i++)
1048  {
1049   chan=fEvent.h[i].ch+1;
1050   if (chan>maxchan) continue; // Channels 9001, 9002 etc are trigger channels
1051
1052   // Get corresponding device from the current event structure  
1053   omx=(IceAOM*)evt->GetIdDevice(chan);
1054   if (!omx)
1055   {
1056    if (fOmdb)
1057    {
1058     omx=(IceAOM*)fOmdb->GetObject(chan,1);
1059     evt->AddDevice(omx);
1060    }
1061    else
1062    {
1063     om.Reset(1);
1064     om.SetUniqueID(chan);
1065     evt->AddDevice(om);
1066    }
1067    omx=(IceAOM*)evt->GetIdDevice(chan);
1068   }
1069
1070   if (!omx) continue;
1071
1072   adc=fEvent.h[i].amp;
1073
1074   // Multiple hits in the same OM with the same ADC value
1075   // are indicated by "*" in the F2K file.
1076   // This corresponds to a value of -2 in the data structure.
1077   if (int(adc) == -2)
1078   {
1079    adc=adcfirst;
1080   }
1081   else
1082   {
1083    adcfirst=adc;
1084   }
1085   s.Reset();
1086   s.SetUniqueID(fEvent.h[i].id);
1087   s.SetSignal(adc,1);
1088   s.SetSignal((fEvent.h[i].t-fToffset),2);
1089   s.SetSignal(fEvent.h[i].tot,3);
1090
1091   omx->AddHit(s);
1092
1093   sx=omx->GetHit(omx->GetNhits());
1094   if (!sx) continue;
1095
1096   // ADC dependent TDC (de)calibration function for this hit
1097   TF1* fcal=omx->GetCalFunction("LE");
1098   TF1* fdecal=omx->GetDecalFunction("LE");
1099   if (fcal) sx->SetCalFunction(fcal,2);
1100   if (fdecal) sx->SetDecalFunction(fdecal,2);
1101   fcal=sx->GetCalFunction(2);
1102   fdecal=sx->GetDecalFunction(2);
1103   adc=sx->GetSignal(1,-4);
1104   if (adc>0)
1105   {
1106    if (fcal) fcal->SetParameter(3,adc);
1107    if (fdecal) fdecal->SetParameter(3,adc);
1108   }
1109   else
1110   {
1111    if (fcal) fcal->SetParameter(3,1.e20);
1112    if (fdecal) fdecal->SetParameter(3,1.e20);
1113   }
1114
1115   // Bi-directional link between this hit and the track that caused the ADC value.
1116   // This F2K info is probably only present for MC tracks.
1117   tid=fEvent.h[i].ma;
1118   if (tid > 0)
1119   {
1120    tx=evt->GetIdTrack(tid); // Reco tracks
1121    if (!tx) tx=evt->GetIdTrack(-tid); // MC tracks
1122    if (tx) sx->AddTrack(*tx);
1123   }
1124   else
1125   {
1126    if (tid == -2) sx->SetNameTitle("N","Noise");
1127    if (tid == -3) sx->SetNameTitle("A","Afterpulse");
1128   }
1129  }
1130
1131  // Loop over all the waveforms and add the histo(s) to the corresponding OM's
1132  TH1F histo;
1133  Int_t nbins=0;
1134  Float_t xlow=0;
1135  Float_t xup=0;
1136  TString hname;
1137  for (Int_t iwf=0; iwf<fEvent.nwf; iwf++)
1138  {
1139   chan=fEvent.wf[iwf].om;
1140   if (chan<=0 || chan>maxchan) continue; // Skip trigger channels
1141
1142   // Get corresponding device from the current event structure  
1143   omx=(IceAOM*)evt->GetIdDevice(chan);
1144   if (!omx)
1145   {
1146    if (fOmdb)
1147    {
1148     omx=(IceAOM*)fOmdb->GetObject(chan,1);
1149     evt->AddDevice(omx);
1150    }
1151    else
1152    {
1153     om.Reset(1);
1154     om.SetUniqueID(chan);
1155     evt->AddDevice(om);
1156    }
1157    omx=(IceAOM*)evt->GetIdDevice(chan);
1158   }
1159
1160   if (!omx) continue;
1161
1162   hname="BASELINE-WF";
1163   hname+=omx->GetNwaveforms()+1;
1164   omx->AddNamedSlot(hname);
1165   omx->SetSignal(fEvent.wf[iwf].baseline,hname);
1166
1167   // Fill the waveform histogram
1168   hname="OM";
1169   hname+=chan;
1170   hname+="-WF";
1171   hname+=omx->GetNwaveforms()+1;
1172
1173   histo.Reset();
1174   histo.SetName(hname.Data());
1175   nbins=fEvent.wf[iwf].ndigi;
1176   xlow=fEvent.wf[iwf].t_start;
1177   xup=xlow+float(nbins)*fEvent.wf[iwf].t_bin;
1178   histo.SetBins(nbins,xlow,xup);
1179
1180   for (Int_t jbin=1; jbin<=fEvent.wf[iwf].ndigi; jbin++)
1181   {
1182    histo.SetBinContent(jbin,fEvent.wf[iwf].baseline-fEvent.wf[iwf].digi[jbin-1]);
1183   }
1184
1185   omx->SetWaveform(&histo,omx->GetNwaveforms()+1);
1186  }
1187
1188  // Set bi-directional links between hits and reco track hypotheses.
1189  // Note : Reco tracks are recognised by their positive id.
1190  Int_t hid=0;
1191  TObjArray* rectracks=evt->GetTracks(1);
1192  for (Int_t jtk=0; jtk<rectracks->GetEntries(); jtk++)
1193  {
1194   tx=(AliTrack*)rectracks->At(jtk);
1195   if (!tx) continue;
1196   
1197   for (Int_t jhyp=1; jhyp<=tx->GetNhypotheses(); jhyp++)
1198   {
1199    AliTrack* hypx=tx->GetTrackHypothesis(jhyp);
1200    if (!hypx) continue;
1201
1202    // Loop over all combinations of F2K fits and used OM hits
1203    for (Int_t k=0; k<fEvent.nfit_uses; k++)
1204    {
1205     if (fEvent.fit_uses[k].useid != hypx->GetId()) continue;
1206     hid=fEvent.fit_uses[k].hitid;
1207     sx=evt->GetIdHit(hid,"IceAOM");
1208     if (sx) sx->AddTrack(*hypx);
1209    }
1210   }
1211  }
1212 }
1213 ///////////////////////////////////////////////////////////////////////////
1214 void IceF2k::PutTrigger()
1215 {
1216 // Get the trigger info from the F2000 file into the IcePack structure.
1217
1218  if (!fTrigdefs) return;
1219
1220  IceEvent* evt=(IceEvent*)GetMainObject();
1221  if (!evt || fEvent.ntrig<=0) return;
1222
1223  AliDevice trig;
1224  trig.SetNameTitle("Trigger","Amanda/IceCube event triggers");
1225  AliSignal s;
1226  TString trigname;
1227  TString slotname;
1228  Int_t id=0;
1229  Int_t nval=0;
1230  for (Int_t i=0; i<fEvent.ntrig; i++)
1231  {
1232   id=fEvent.ptrig[i].id;
1233   nval=fEvent.ptrig[i].nval;
1234   if (!nval) continue;
1235   AliSignal* tdef=fTrigdefs->GetIdHit(id+1);
1236   if (!tdef) continue;
1237   trigname=tdef->GetName();
1238   s.Reset(1);
1239   s.SetName(trigname);
1240   s.SetUniqueID(id+1);
1241   for (Int_t jval=0; jval<fEvent.ptrig[i].nval; jval++)
1242   {
1243    slotname=tdef->GetSlotName(jval+1);
1244    s.SetSlotName(slotname,jval+1);
1245    s.SetSignal(fEvent.ptrig[i].val[jval],jval+1);
1246   }
1247   trig.AddHit(s);
1248  }
1249
1250  // Store the trigger data into the IceEvent structure
1251  evt->AddDevice(trig);
1252 }
1253 ///////////////////////////////////////////////////////////////////////////