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 ///////////////////////////////////////////////////////////////////////////