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