Forgot to add the new IceRoot to CVS registry.
[u/mrichter/AliRoot.git] / RALICE / icepack / iceconvert / IceRoot.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 IceRoot
20 // Conversion of simple Root data into IceEvent data structures.
21 // This class reads data from the simple Root files as output by Martijn
22 // Duvoort's Walnut analyser.
23 // An indication of the active DAQ system is available in the IceEvent structure
24 // via a device named "Daq". Here the various daq systems (TWR, Muon, ...)
25 // from which the actual hits (ADC, LE, TOT) eventually will be composed
26 // are indicated as "signals" of the device itself. 
27 // This class is derived from AliJob providing a task-based processing
28 // structure on an event-by-event basis.
29 // The main object in the job environment is an IceEvent* pointer.
30 // In case the user has provided sub-tasks, these will be executed
31 // on an event-by-event basis after the IceEvent structure has been filled
32 // with the simple Root data and before the final structures are written out.
33 // Note that the data structures are only written out if an outputfile has
34 // been specified via the SetOutputFile memberfunction.
35 // In case no outputfile has been specified, this class provides a facility
36 // to investigate/analyse simple Root data using the Ralice/IcePack analysis tools.
37 //
38 // Usage example :
39 // ---------------
40 //
41 // gSystem->Load("ralice");
42 // gSystem->Load("icepack");
43 // gSystem->Load("iceconvert");
44 //
45 // IceRoot q("IceRoot","Simple Root data to IcePack data structure conversion");
46 //
47 // // Limit the number of entries for testing
48 // q.SetMaxEvents(10);
49 //
50 // // Print frequency to produce a short summary print every printfreq events
51 // q.SetPrintFreq(1);
52 //
53 // // The simple Root data input filename(s)
54 // q.AddInputFile("test-i3.root");
55 //
56 // // Output file for the event structures
57 // q.SetOutputFile("events.root");
58 //
59 // ///////////////////////////////////////////////////////////////////
60 // // Here the user can specify his/her sub-tasks to be executed
61 // // on an event-by-event basis after the IceEvent structure
62 // // has been filled and before the data is written out.
63 // // Sub-tasks (i.e. a user classes derived from TTask) are entered
64 // // as follows :
65 // //
66 // //    MyXtalk task1("task1","Cross talk correction");
67 // //    MyClean task2("task2","Hit cleaning");
68 // //    q.Add(&task1);
69 // //    q.Add(&task2);
70 // //
71 // // The sub-tasks will be executed in the order as they are entered.
72 // ///////////////////////////////////////////////////////////////////
73 //
74 // // Perform the conversion and execute subtasks (if any)
75 // // on an event-by-event basis
76 // q.ExecuteJob();
77 //
78 //--- Author: Garmt de Vries-Uiterweerd 13-Mar-2007 Utrecht University
79 //- Modified: GdV $Date$ Utrecht University
80 ///////////////////////////////////////////////////////////////////////////
81  
82 #include "IceRoot.h"
83 #include "Riostream.h"
84
85 ClassImp(IceRoot) // Class implementation to enable ROOT I/O
86
87 IceRoot::IceRoot(const char* name,const char* title) : AliJob(name,title)
88 {
89 // Default constructor.
90 // By default maxevent=-1, split=0, bsize=32000, printfreq=1.
91
92  fSplit=0;
93  fBsize=32000;
94  fMaxevt=-1;
95  fPrintfreq=1;
96  fInfiles=0;
97  fOutfile=0;
98 }
99 ///////////////////////////////////////////////////////////////////////////
100 IceRoot::~IceRoot()
101 {
102 // Default destructor.
103
104  if (fInfiles)
105  {
106   delete fInfiles;
107   fInfiles=0;
108  }
109 }
110 ///////////////////////////////////////////////////////////////////////////
111 void IceRoot::SetMaxEvents(Int_t n)
112 {
113 // Set the maximum number of events to be processed.
114 // n=-1 implies processing of the complete input file, which is the default
115 // initialisation in the constructor.
116  fMaxevt=n;
117 }
118 ///////////////////////////////////////////////////////////////////////////
119 void IceRoot::SetPrintFreq(Int_t f)
120 {
121 // Set the printfrequency to produce info every f events.
122 // f=1 is the default initialisation in the constructor.
123  if (f>=0) fPrintfreq=f;
124 }
125 ///////////////////////////////////////////////////////////////////////////
126 void IceRoot::SetSplitLevel(Int_t split)
127 {
128 // Set the split level for the ROOT data file.
129 // split=0 is the default initialisation in the constructor.
130  if (split>=0) fSplit=split;
131 }
132 ///////////////////////////////////////////////////////////////////////////
133 void IceRoot::SetBufferSize(Int_t bsize)
134 {
135 // Set the buffer size for the ROOT data file.
136 // bsize=32000 is the default initialisation in the constructor.
137  if (bsize>=0) fBsize=bsize;
138 }
139 ///////////////////////////////////////////////////////////////////////////
140 void IceRoot::AddInputFile(TString name)
141 {
142 // Add the name of this simple Root data input file to the list to be processed.
143
144  if (!fInfiles)
145  {
146   fInfiles=new TObjArray();
147   fInfiles->SetOwner();
148  }
149
150  TObjString* s=new TObjString();
151  s->SetString(name);
152  fInfiles->Add(s);
153 }
154 ///////////////////////////////////////////////////////////////////////////
155 void IceRoot::SetOutputFile(TFile* ofile)
156 {
157 // Set the output file for the ROOT data.
158  if (fOutfile) delete fOutfile;
159  fOutfile=ofile;
160 }
161 ///////////////////////////////////////////////////////////////////////////
162 void IceRoot::SetOutputFile(TString name)
163 {
164 // Create the output file for the ROOT data.
165  if (fOutfile) delete fOutfile;
166  fOutfile=new TFile(name.Data(),"RECREATE","Simple Root data in IceEvent structure");
167 }
168 ///////////////////////////////////////////////////////////////////////////
169 TFile* IceRoot::GetOutputFile()
170 {
171 // Provide pointer to the ROOT output file.
172  return fOutfile;
173 }
174 ///////////////////////////////////////////////////////////////////////////
175 void IceRoot::Exec(Option_t* opt)
176 {
177 // Job to loop over the specified number of events and convert the 
178 // simple Root data into the IceEvent structure.
179 // If maxevents<0 (default) all the entries of the input file
180 // will be processed.
181 // Every "printfreq" events a short event summary will be printed.
182 // The default value is printfreq=1.
183 // The output will be written on a standard output tree named "T".
184 //
185 // Notes :
186 // -------
187 // 1) This class is derived from AliJob, allowing a task based processing.
188 //    After the conversion of a simple Root data event into an IceEvent structure,
189 //    the processing of all available sub-tasks (if any) is invoked.
190 //    This provides an event-by-event (sub)task processing before the
191 //    final data structures are written out.
192 // 2) The main object in this job environment is an IceEvent* pointer.
193
194  if (!fInfiles)
195  {
196   cout << " *IceRoot Exec* No data input file(s) specified." << endl;
197   return;
198  }
199
200  Int_t ninfiles=fInfiles->GetEntries();
201  if (!ninfiles)
202  {
203   cout << " *IceRoot Exec* No data input file(s) specified." << endl;
204   return;
205  }
206
207  // Create output tree if necessary
208  TTree* otree=0;
209  if (fOutfile)
210  {
211   otree=new TTree("T","Simple Root data converted to IceEvent structures");
212   otree->SetDirectory(fOutfile);
213  }
214
215  // Create IceEvent structure
216  IceEvent* evt=new IceEvent();
217  evt->SetTrackCopy(1);
218  evt->SetDevCopy(1);
219
220  // Branch in the tree for the event structure
221  if (otree) otree->Branch("IceEvent","IceEvent",&evt,fBsize,fSplit); 
222
223  // Initialise the job working environment
224  SetMainObject(evt);
225  if (fOutfile)
226  {
227   AddObject(fOutfile);
228   AddObject(otree);
229  }
230
231  TString inputfile;
232
233  cout << " ***" << endl;
234  cout << " *** Start processing of job " << GetName() << " ***" << endl;
235  cout << " ***" << endl;
236  for (Int_t i=0; i<ninfiles; i++)
237  {
238   TObjString* sx=(TObjString*)fInfiles->At(i);
239   if (!sx) continue;
240   inputfile=sx->GetString(); 
241   cout << " Simple Root data input file : " << inputfile.Data() << endl;
242  }
243  cout << " Maximum number of events to be processed : " << fMaxevt << endl;
244  cout << " Print frequency : " << fPrintfreq << endl;
245  if (fOutfile)
246  {
247   cout << " ROOT output file : " << fOutfile->GetName() << endl;
248   cout << " Output characteristics : splitlevel = " << fSplit << " buffersize = " << fBsize << endl;
249  }
250
251  ListEnvironment();
252
253  // Set DAQ device info
254  AliDevice daq;
255  daq.SetName("Daq");
256  daq.SetSlotName("TWR",1);
257  daq.SetSignal(1,1);
258
259  // Some variables
260  Double_t triggertime=0;
261  Int_t    eventtimemjd=0;
262  Int_t    eventtimemjds=0;
263  Int_t    eventtimemjdns=0;
264  Int_t    eventid=0;
265  Int_t    runid=0;
266  Int_t    String=0;
267  Int_t    OM=0;
268  Float_t  baseline=0;
269  Short_t  type=0;
270  Int_t    numberbins=0;
271  Float_t  starttime=0;
272  Int_t    ntracks=0;
273
274  Int_t nevt=0;
275  Int_t lastevent=0;
276  Int_t omid=0;
277
278  TString hname;
279  TH1F histo;
280  IceAOM om;
281  IceAOM* omx=0;
282  AliTrack t;
283  Double_t vec[3];
284  AliPosition r;
285  Ali3Vector p;
286
287  Float_t pi=acos(-1.);
288
289  Float_t binwidth=10;
290  Int_t firstomonstring[20];
291  firstomonstring[0]=0;
292  firstomonstring[1]=1;
293  firstomonstring[2]=21;
294  firstomonstring[3]=41;
295  firstomonstring[4]=61;
296  firstomonstring[5]=87;
297  firstomonstring[6]=123;
298  firstomonstring[7]=159;
299  firstomonstring[8]=195;
300  firstomonstring[9]=231;
301  firstomonstring[10]=267;
302  firstomonstring[11]=303;
303  firstomonstring[12]=345;
304  firstomonstring[13]=387;
305  firstomonstring[14]=429;
306  firstomonstring[15]=471;
307  firstomonstring[16]=513;
308  firstomonstring[17]=555;
309  firstomonstring[18]=597;
310  firstomonstring[19]=639;
311
312  for (Int_t ifile=0; ifile<ninfiles; ifile++)
313  {
314   TObjString* sx=(TObjString*)fInfiles->At(ifile);
315   if (!sx) continue;
316
317   inputfile=sx->GetString(); 
318   if (inputfile=="") continue;
319
320   // Open the simple Root data input file
321   fInput=new TFile(inputfile);
322
323   if (!fInput)
324   {
325    cout << " *IceRoot Exec* No input file found with name : " << inputfile.Data() << endl;
326    continue;
327   }
328
329   // Get simple Root tree
330   fTree=(TTree*)fInput->Get("T");
331   fTree->SetBranchAddress("triggertime",&triggertime);
332   fTree->SetBranchAddress("eventtimemjd",&eventtimemjd);
333   fTree->SetBranchAddress("eventtimemjds",&eventtimemjds);
334   fTree->SetBranchAddress("eventtimemjdns",&eventtimemjdns);
335   fTree->SetBranchAddress("eventid",&eventid);
336   fTree->SetBranchAddress("runid",&runid);
337   fTree->SetBranchAddress("String",&String);
338   fTree->SetBranchAddress("OM",&OM);
339   fTree->SetBranchAddress("baseline",&baseline);
340   fTree->SetBranchAddress("type",&type);
341   fTree->SetBranchAddress("numberbins",&numberbins);
342   fTree->SetBranchAddress("starttime",&starttime);
343   fTree->SetBranchAddress("ntracks",&ntracks);
344
345   Int_t nmaxbins=(Int_t)fTree->GetLeaf("numberbins")->GetMaximum();
346   Float_t* wvform=new Float_t[nmaxbins+1];
347   fTree->SetBranchAddress("wvform",wvform);
348
349   Int_t nmaxtracks=(Int_t)fTree->GetLeaf("ntracks")->GetMaximum();
350   Float_t* trackx=new Float_t[nmaxtracks];
351   Float_t* tracky=new Float_t[nmaxtracks];
352   Float_t* trackz=new Float_t[nmaxtracks];
353   Float_t* trackzenith=new Float_t[nmaxtracks];
354   Float_t* trackazimuth=new Float_t[nmaxtracks];
355   Int_t* tracktype=new Int_t[nmaxtracks];
356   fTree->SetBranchAddress("trackx",trackx);
357   fTree->SetBranchAddress("tracky",tracky);
358   fTree->SetBranchAddress("trackz",trackz);
359   fTree->SetBranchAddress("trackzenith",trackzenith);
360   fTree->SetBranchAddress("trackazimuth",trackazimuth);
361   fTree->SetBranchAddress("tracktype",tracktype);
362
363   // Prepare for loop over entries
364   lastevent=0;
365
366   // Loop over waveforms in tree
367   for(Int_t ientry=0; ientry<fTree->GetEntries(); ientry++)
368   {
369    fTree->GetEntry(ientry);
370    
371    // If new event
372    if(eventid!=lastevent)
373    {
374     // Write old event to tree (if it contains data)
375     if(evt->GetDevices("IceAOM"))
376     {
377      // Invoke all available sub-tasks (if any) and write event to tree
378      CleanTasks();
379      ExecuteTasks(opt);
380      if(otree) otree->Fill();
381      if (fPrintfreq) { if (!(nevt%fPrintfreq)) evt->HeaderData(); }
382      // Update event counter
383      nevt++;
384      if (fMaxevt>-1 && nevt>=fMaxevt) 
385      {
386       evt->Reset();
387       break;
388      }
389     }
390
391     // Start new event
392     evt->Reset();
393     evt->SetRunNumber(runid);
394     evt->SetEventNumber(eventid);
395     evt->SetMJD(eventtimemjd,eventtimemjds,eventtimemjdns,0);
396     evt->AddDevice(daq);
397
398     // Loop over all the tracks and add them to the current event
399     for (Int_t itrack=0; itrack<ntracks; itrack++)
400     {
401      t.Reset();
402
403      // Beginpoint of the track
404      vec[0]=trackx[itrack];
405      vec[1]=tracky[itrack];
406      vec[2]=trackz[itrack];
407      r.SetPosition(vec,"car");
408      t.SetBeginPoint(r);
409
410      // Momentum in GeV/c
411      vec[0]=1; // Energy is by convention set to 1 
412      vec[1]=pi-trackzenith[itrack];
413      vec[2]=trackazimuth[itrack]+pi;
414      if(vec[2]>=2*pi) vec[2]-=2*pi;
415      p.SetVector(vec,"sph");
416      t.Set3Momentum(p);
417
418      // Check for unreconstructed tracks (NaN values) and reset track if necessary
419      if(trackx[itrack]!=trackx[itrack] || tracky[itrack]!=tracky[itrack] || trackz[itrack]!=trackz[itrack] ||
420         trackzenith[itrack]!=trackzenith[itrack] || trackazimuth[itrack]!=trackazimuth[itrack]){
421       t.Reset();
422      }
423
424      // Track ID and name
425      t.SetId(itrack+1);
426      if(tracktype[itrack]==1) { t.SetName("I3CFIRST"); t.SetTitle("IceTray CFIRST track"); }
427      else if(tracktype[itrack]==2) { t.SetName("I3direct"); t.SetTitle("IceTray direct walk track"); }
428      else if(tracktype[itrack]==3) { t.SetName("I3Jams_cluster"); t.SetTitle("IceTray Jams cluster"); }
429      else if(tracktype[itrack]==4) { t.SetName("I3Jams_qual"); t.SetTitle("IceTray Jams quality"); }
430      else if(tracktype[itrack]==5) { t.SetName("I3TOIFit"); t.SetTitle("IceTray TOI fit"); }
431      else if(tracktype[itrack]==6) { t.SetName("I3line-direct"); t.SetTitle("IceTray line direct walk"); }
432      else t.SetName("Unknown");
433
434      // Add track to event
435      evt->AddTrack(t);
436     }
437
438     // Remember event nr
439     lastevent=eventid;
440    }
441
442    omid=firstomonstring[-String]+OM-1;
443    // Get corresponding device from the current event structure  
444    omx=(IceAOM*)evt->GetIdDevice(omid);
445    if (!omx)
446    {
447     om.Reset(1);
448     om.SetUniqueID(omid);
449     evt->AddDevice(om);
450     omx=(IceAOM*)evt->GetIdDevice(omid);
451    }
452    if (!omx) continue;
453
454    // Store baseline info
455    hname="BASELINE-WF";
456    hname+=omx->GetNwaveforms()+1;
457    omx->AddNamedSlot(hname);
458    omx->SetSignal(baseline,hname);
459
460    // Fill the waveform histogram with this fragment
461    hname="OM";
462    hname+=omid;
463    hname+="-WF";
464    hname+=omx->GetNwaveforms()+1;
465
466    histo.Reset();
467    histo.SetName(hname.Data());
468    histo.SetBins(numberbins,starttime,starttime+numberbins*binwidth);
469
470    for (Int_t jbin=1; jbin<=numberbins; jbin++)
471    {
472     histo.SetBinContent(jbin,baseline-wvform[jbin-1]);
473    }
474
475    omx->SetWaveform(&histo,omx->GetNwaveforms()+1);
476
477   }
478
479   // Write last event to tree
480   if(evt->GetDevices("IceAOM"))
481   {
482    CleanTasks();
483    ExecuteTasks(opt);
484    if (otree) otree->Fill();
485    if (fPrintfreq) { if (!(nevt%fPrintfreq)) evt->HeaderData(); }
486
487    // Update event counter
488    nevt++;
489    if (fMaxevt>-1 && nevt>=fMaxevt) break;
490   }
491
492   // Clean up
493   delete[] wvform;
494   delete[] trackx;
495   delete[] tracky;
496   delete[] trackz;
497   delete[] trackzenith;
498   delete[] trackazimuth;
499   delete[] tracktype;
500  } // End of input file loop
501
502  // Flush possible memory resident data to the output file
503  if (fOutfile) fOutfile->Write();
504
505  // Remove the IceEvent object from the environment
506  // and delete it as well
507  if (evt)
508  {
509   RemoveObject(evt);
510   delete evt;
511  }
512 }
513 ///////////////////////////////////////////////////////////////////////////