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