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