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
  60 // q.SetInputFile("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  fInfile="";
 127  fOutfile=0;
 128
 129  fPdg=0;
 130  fOmdb=0;
 131  fFitdefs=0;
 132 }
 133 ///////////////////////////////////////////////////////////////////////////
 134 IceF2k::~IceF2k()
 135 {
 136 // Default destructor.
 137
 138  if (fPdg)
 139  {
 140   delete fPdg;
 141   fPdg=0;
 142  }
 143
 144  if (fOmdb)
 145  {
 146   delete fOmdb;
 147   fOmdb=0;
 148  }
 149
 150  if (fFitdefs)
 151  {
 152   delete fFitdefs;
 153   fFitdefs=0;
 154  }
 155 }
 156 ///////////////////////////////////////////////////////////////////////////
 157 void IceF2k::SetMaxEvents(Int_t n)
 158 {
 159 // Set the maximum number of events to be processed.
 160 // n=-1 implies processing of the complete input file, which is the default
 161 // initialisation in the constructor.
 162  fMaxevt=n;
 163 }
 164 ///////////////////////////////////////////////////////////////////////////
 165 void IceF2k::SetPrintFreq(Int_t f)
 166 {
 167 // Set the printfrequency to produce info every f events.
 168 // f=1 is the default initialisation in the constructor.
 169  if (f>=0) fPrintfreq=f;
 170 }
 171 ///////////////////////////////////////////////////////////////////////////
 172 void IceF2k::SetSplitLevel(Int_t split)
 173 {
 174 // Set the split level for the ROOT data file.
 175 // split=99 is the default initialisation in the constructor.
 176  if (split>=0) fSplit=split;
 177 }
 178 ///////////////////////////////////////////////////////////////////////////
 179 void IceF2k::SetBufferSize(Int_t bsize)
 180 {
 181 // Set the buffer size for the ROOT data file.
 182 // bsize=32000 is the default initialisation in the constructor.
 183  if (bsize>=0) fBsize=bsize;
 184 }
 185 ///////////////////////////////////////////////////////////////////////////
 186 void IceF2k::SetInputFile(TString name)
 187 {
 188 // Set the name of the F2K input file.
 189  fInfile=name;
 190 }
 191 ///////////////////////////////////////////////////////////////////////////
 192 void IceF2k::SetOutputFile(TFile* ofile)
 193 {
 194 // Set the output file for the ROOT data.
 195  if (fOutfile) delete fOutfile;
 196  fOutfile=ofile;
 197 }
 198 ///////////////////////////////////////////////////////////////////////////
 199 void IceF2k::SetOutputFile(TString name)
 200 {
 201 // Create the output file for the ROOT data.
 202  if (fOutfile) delete fOutfile;
 203  fOutfile=new TFile(name.Data(),"RECREATE","F2K data in IceEvent structure");
 204 }
 205 ///////////////////////////////////////////////////////////////////////////
 206 TFile* IceF2k::GetOutputFile()
 207 {
 208 // Provide pointer to the ROOT output file.
 209  return fOutfile;
 210 }
 211 ///////////////////////////////////////////////////////////////////////////
 212 TDatabasePDG* IceF2k::GetPDG()
 213 {
 214 // Provide pointer to the PDG database
 215  return fPdg;
 216 }
 217 ///////////////////////////////////////////////////////////////////////////
 218 AliObjMatrix* IceF2k::GetOMdbase()
 219 {
 220 // Provide pointer to the OM geometry, calib. etc... database
 221  return fOmdb;
 222 }
 223 ///////////////////////////////////////////////////////////////////////////
 224 AliDevice* IceF2k::GetFitdefs()
 225 {
 226 // Provide pointer to the fit definitions
 227  return fFitdefs;
 228 }
 229 ///////////////////////////////////////////////////////////////////////////
 230 void IceF2k::Exec(Option_t* opt)
 231 {
 232 // Job to loop over the specified number of events and convert the
 233 // F2K data into the IceEvent structure.
 234 // If maxevents<0 (default) all the entries of the input file
 235 // will be processed.
 236 // Every "printfreq" events a short event summary will be printed.
 237 // The default value is printfreq=1.
 238 // The output will be written on a standard output tree named "T".
 239 //
 240 // Notes :
 241 // -------
 242 // 1) This class is derived from AliJob, allowing a task based processing.
 243 //    After the conversion of an F2K event into an IceEvent structure,
 244 //    the processing of all available sub-tasks (if any) is invoked.
 245 //    This provides an event-by-event (sub)task processing before the
 246 //    final data structures are written out.
 247 // 2) The main object in this job environment is an IceEvent* pointer.
 248
 249  if (fInfile=="")
 250  {
 251   cout << " *IceF2k Exec* No data input file specified." << endl;
 252   return;
 253  }
 254
 255  // Open the input file in the default ascii format (autodetection) for reading
 256  fInput=rdmc_mcopen(fInfile.Data(),"r",RDMC_DEFAULT_ASCII_F);
 257
 258  if (!fInput)
 259  {
 260   cout << " *IceF2k Exec* No input file found with name : " << fInfile.Data() << endl;
 261   return;
 262  }
 263
 264  // Initialise the event structure
 265  rdmc_init_mevt(&fEvent);
 266
 267  // Read the file header information
 268  rdmc_rarr(fInput,&fHeader);
 269
 270  TTree* otree=0;
 271  if (fOutfile)
 272  {
 273   otree=new TTree("T","F2K Data converted to IceEvent structures");
 274   otree->SetDirectory(fOutfile);
 275  }
 276
 277  IceEvent* evt=new IceEvent();
 278  evt->SetTrackCopy(1);
 279  evt->SetDevCopy(1);
 280
 281  // Branch in the tree for the event structure
 282  if (otree) otree->Branch("IceEvent","IceEvent",&evt,fBsize,fSplit);
 283
 284  // Create the particle database and extend it with some F2000 specific definitions
 285  if (!fPdg) fPdg=new TDatabasePDG();
 286  Double_t me=fPdg->GetParticle(11)->Mass();
 287  fPdg->AddParticle("brems"   ,"brems"   ,0,1,0,0,"none",10001001,0,0);
 288  fPdg->AddParticle("deltae"  ,"deltae"  ,me,1,0,-3,"Lepton",10001002,0,0);
 289  fPdg->AddParticle("pairprod","pairprod",0,1,0,0,"none",10001003,0,0);
 290  fPdg->AddParticle("nucl_int","nucl_Int",0,1,0,0,"none",10001004,0,0);
 291  fPdg->AddParticle("mu_pair" ,"mu_pair" ,0,1,0,0,"none",10001005,0,0);
 292  fPdg->AddParticle("hadrons" ,"hadrons" ,0,1,0,0,"none",10001006,0,0);
 293  fPdg->AddParticle("fiberlaser","fiberlaser",0,1,0,0,"none",10002100,0,0);
 294  fPdg->AddParticle("n2laser"   ,"n2laser"   ,0,1,0,0,"none",10002101,0,0);
 295  fPdg->AddParticle("yaglaser"  ,"yaglaser"  ,0,1,0,0,"none",10002201,0,0);
 296  fPdg->AddParticle("z_primary","z_primary",0,1,0,0,"none",10003000,0,0);
 297  fPdg->AddParticle("a_primary","a_primary",0,1,0,0,"none",10003500,0,0);
 298
 299  // Fill the database with geometry, calib. etc... parameters
 300  // for all the devices
 301  FillOMdbase();
 302
 303  // Set the fit definitions according to the F2000 header info
 304  SetFitdefs();
 305
 306  // Initialise the job working environment
 307  SetMainObject(evt);
 308  if (fOutfile)
 309  {
 310   AddObject(fOutfile);
 311   AddObject(otree);
 312  }
 313
 314  cout << " ***" << endl;
 315  cout << " *** Start processing of job " << GetName() << " ***" << endl;
 316  cout << " ***" << endl;
 317  cout << " F2K input file : " << fInfile.Data() << endl;
 318  cout << " Maximum number of events to be processed : " << fMaxevt << endl;
 319  cout << " Print frequency : " << fPrintfreq << endl;
 320  if (fOutfile)
 321  {
 322   cout << " ROOT output file : " << fOutfile->GetName() << endl;
 323   cout << " Output characteristics : splitlevel = " << fSplit << " buffersize = " << fBsize << endl;
 324  }
 325
 326  ListEnvironment();
 327
 328  Int_t nevt=0;
 329  while (!rdmc_revt(fInput,&fHeader,&fEvent))
 330  {
 331   if (fMaxevt>-1 && nevt>=fMaxevt) break;
 332
 333   // Reset the complete Event structure
 334   evt->Reset();
 335
 336   evt->SetRunNumber(fEvent.nrun);
 337   evt->SetEventNumber(fEvent.enr);
 338   evt->SetMJD(fEvent.mjd,fEvent.secs,fEvent.nsecs);
 339
 340   PutMcTracks();
 341
 342   PutRecoTracks();
 343
 344   PutHits();
 345
 346   // Invoke all available sub-tasks (if any)
 347   CleanTasks();
 348   ExecuteTasks(opt);
 349
 350   if (fPrintfreq)
 351   {
 352    if (!(nevt%fPrintfreq)) evt->HeaderData();
 353   }
 354
 355   // Write the complete structure to the output Tree
 356   if (otree) otree->Fill();
 357
 358   // Update event counter
 359   nevt++;
 360  }
 361
 362  // Flush possible memory resident data to the output file
 363  if (fOutfile) fOutfile->Write();
 364
 365  // Remove the IceEvent object from the environment
 366  // and delete it as well
 367  if (evt)
 368  {
 369   RemoveObject(evt);
 370   delete evt;
 371  }
 372 }
 373 ///////////////////////////////////////////////////////////////////////////
 374 void IceF2k::FillOMdbase()
 375 {
 376 // Fill the database with geometry, calib. etc... parameters
 377 // for all the devices.
 378
 379  if (fHeader.nch<=0) return;
 380
 381  Int_t adccal=fHeader.is_calib.adc;
 382  Int_t tdccal=fHeader.is_calib.tdc;
 383  Int_t totcal=fHeader.is_calib.tot;
 384
 385  TF1 fadccal("fadccal","(x-[1])*[0]");
 386  TF1 fadcdecal("fadcdecal","(x/[0])+[1]");
 387  fadccal.SetParName(0,"BETA-ADC");
 388  fadccal.SetParName(1,"PED-ADC");
 389  fadcdecal.SetParName(0,"BETA-ADC");
 390  fadcdecal.SetParName(1,"PED-ADC");
 391
 392  TF1 ftdccal("ftdccal","(x*[0])-[1]-([0]-1.)*32767.-[2]/sqrt([3])");
 393  TF1 ftdcdecal("ftdcdecal","(x+([0]-1.)*32767.+[1]+[2]/sqrt([3]))/[0]");
 394  ftdccal.SetParName(0,"BETA-TDC");
 395  ftdccal.SetParName(1,"T0");
 396  ftdccal.SetParName(2,"ALPHA-TDC");
 397  ftdccal.SetParName(3,"ADC-SLEW");
 398  ftdcdecal.SetParName(0,"BETA-TDC");
 399  ftdcdecal.SetParName(1,"T0");
 400  ftdcdecal.SetParName(2,"ALPHA-TDC");
 401  ftdcdecal.SetParName(3,"ADC-SLEW");
 402
 403  TF1 ftotcal("ftotcal","x*[0]");
 404  TF1 ftotdecal("ftotdecal","x/[0]");
 405  ftotcal.SetParName(0,"BETA-TOT");
 406  ftotdecal.SetParName(0,"BETA-TOT");
 407
 408  if (fOmdb)
 409  {
 410   fOmdb->Reset();
 411  }
 412  else
 413  {
 414   fOmdb=new AliObjMatrix();
 415   fOmdb->SetNameTitle("OMDBASE","The OM geometry, calib. etc... database");
 416   fOmdb->SetOwner();
 417  }
 418
 419  IceAOM* dev=0;
 420  Double_t pos[3]={0,0,0};
 421  for (Int_t i=0; i<fHeader.nch; i++)
 422  {
 423   dev=new IceAOM();
 424   dev->SetUniqueID(i+1);
 425
 426   dev->SetSlotName("ADC",1);
 427   dev->SetSlotName("LE",2);
 428   dev->SetSlotName("TOT",3);
 429
 430   dev->SetSlotName("TYPE",4);
 431   dev->SetSlotName("ORIENT",5);
 432   dev->SetSlotName("THRESH",6);
 433   dev->SetSlotName("SENSIT",7);
 434   dev->SetSlotName("BETA-TDC",8);
 435   dev->SetSlotName("T0",9);
 436   dev->SetSlotName("ALPHA-TDC",10);
 437   dev->SetSlotName("PED-ADC",11);
 438   dev->SetSlotName("BETA-ADC",12);
 439   dev->SetSlotName("KAPPA-ADC",13);
 440   dev->SetSlotName("PED-TOT",14);
 441   dev->SetSlotName("BETA-TOT",15);
 442   dev->SetSlotName("KAPPA-TOT",16);
 443
 444   pos[0]=fHeader.x[i];
 445   pos[1]=fHeader.y[i];
 446   pos[2]=fHeader.z[i];
 447   dev->SetPosition(pos,"car");
 448
 449   fadccal.SetParameter(0,fHeader.cal[i].beta_a);
 450   fadccal.SetParameter(1,fHeader.cal[i].ped);
 451   fadcdecal.SetParameter(0,fHeader.cal[i].beta_a);
 452   if (!fHeader.cal[i].beta_a) fadcdecal.SetParameter(0,1);
 453   fadcdecal.SetParameter(1,fHeader.cal[i].ped);
 454
 455   ftdccal.SetParameter(0,fHeader.cal[i].beta_t);
 456   ftdccal.SetParameter(1,fHeader.cal[i].t_0);
 457   ftdccal.SetParameter(2,fHeader.cal[i].alpha_t);
 458   ftdccal.SetParameter(3,1.e20);
 459   ftdcdecal.SetParameter(0,fHeader.cal[i].beta_t);
 460   if (!fHeader.cal[i].beta_t) ftdcdecal.SetParameter(0,1);
 461   ftdcdecal.SetParameter(1,fHeader.cal[i].t_0);
 462   ftdcdecal.SetParameter(2,fHeader.cal[i].alpha_t);
 463   ftdcdecal.SetParameter(3,1.e20);
 464
 465   ftotcal.SetParameter(0,fHeader.cal[i].beta_tot);
 466   ftotdecal.SetParameter(0,fHeader.cal[i].beta_tot);
 467   if (!fHeader.cal[i].beta_tot) ftotdecal.SetParameter(0,1);
 468
 469   if (adccal)
 470   {
 471    dev->SetDecalFunction(&fadcdecal,1);
 472   }
 473   else
 474   {
 475    dev->SetCalFunction(&fadccal,1);
 476   }
 477
 478   if (tdccal)
 479   {
 480    dev->SetDecalFunction(&ftdcdecal,2);
 481   }
 482   else
 483   {
 484    dev->SetCalFunction(&ftdccal,2);
 485   }
 486
 487   if (totcal)
 488   {
 489    dev->SetDecalFunction(&ftotdecal,3);
 490   }
 491   else
 492   {
 493    dev->SetCalFunction(&ftotcal,3);
 494   }
 495
 496   dev->SetSignal(fHeader.type[i],4);
 497   dev->SetSignal((Float_t)fHeader.costh[i],5);
 498   dev->SetSignal(fHeader.thresh[i],6);
 499   dev->SetSignal(fHeader.sensit[i],7);
 500   dev->SetSignal(fHeader.cal[i].beta_t,8);
 501   dev->SetSignal(fHeader.cal[i].t_0,9);
 502   dev->SetSignal(fHeader.cal[i].alpha_t,10);
 503   dev->SetSignal(fHeader.cal[i].ped,11);
 504   dev->SetSignal(fHeader.cal[i].beta_a,12);
 505   dev->SetSignal(fHeader.cal[i].kappa,13);
 506   dev->SetSignal(fHeader.cal[i].ped_tot,14);
 507   dev->SetSignal(fHeader.cal[i].beta_tot,15);
 508   dev->SetSignal(fHeader.cal[i].kappa_tot,16);
 509
 510   fOmdb->EnterObject(i+1,1,dev);
 511  }
 512 }
 513 ///////////////////////////////////////////////////////////////////////////
 514 void IceF2k::SetFitdefs()
 515 {
 516 // Obtain the names of the variables for each fit procedure from the
 517 // f2000 header. Each different fit procedure is then stored as a separate
 518 // hit of an AliDevice object and the various fit variables are stored
 519 // as separate signal slots of the corresponding hit.
 520 // Via the GetFitdefs() memberfunction this AliDevice object can be
 521 // retrieved and stored in the ROOT output file if wanted.
 522 // The name of the object is FitDefinitions and the stored data can be
 523 // inspected via the AliDevice::Data() memberfunction and looks as follows :
 524 //
 525 //  *AliDevice::Data* Id :0 Name : FitDefinitions
 526 //    Position Vector in car coordinates : 0 0 0
 527 //    Err. in car coordinates : 0 0 0
 528 //  The following 8 hits are registered :
 529 //  *AliSignal::Data* Id :0
 530 //    Position Vector in car coordinates : 0 0 0
 531 //    Err. in car coordinates : 0 0 0
 532 //    Owned by device : AliDevice Name : FitDefinitions
 533 //    Slot : 1 Signal value : 1 name : id
 534 //    Slot : 2 Signal value : 2 name : rchi2
 535 //    Slot : 3 Signal value : 3 name : prob
 536 //    Slot : 4 Signal value : 4 name : sigth
 537 //    Slot : 5 Signal value : 5 name : covmin
 538 //    Slot : 6 Signal value : 6 name : covmax
 539 //    Slot : 7 Signal value : 7 name : cutflag
 540 //    Slot : 8 Signal value : 8 name : chi2
 541 //  *AliSignal::Data* Id :1
 542 //    Position Vector in car coordinates : 0 0 0
 543 //    Err. in car coordinates : 0 0 0
 544 //    Owned by device : AliDevice Name : FitDefinitions
 545 //    Slot : 1 Signal value : 1 name : id
 546 //    Slot : 2 Signal value : 2 name : rchi2
 547 //    Slot : 3 Signal value : 3 name : prob
 548 // etc....
 549 //
 550 // This memberfunction is based on the original idea/code by Adam Bouchta.
 551
 552  if (fHeader.n_fit<=0) return;
 553
 554  if (fFitdefs)
 555  {
 556   fFitdefs->Reset(1);
 557  }
 558  else
 559  {
 560   fFitdefs=new AliDevice();
 561  }
 562
 563  fFitdefs->SetName("FitDefinitions");
 564  fFitdefs->SetHitCopy (1);
 565
 566  AliSignal s;
 567  s.Reset();
 568
 569  for (Int_t i=0; i<fHeader.n_fit; i++)
 570  {
 571   s.SetUniqueID(fHeader.def_fit[i].id);
 572
 573   for (Int_t j=0; j<fHeader.def_fit[i].nwords; j++)
 574   {
 575    s.SetSlotName(TString(fHeader.def_fit[i].words[j]),j+1);
 576    s.SetSignal(j+1,j+1);
 577   }
 578
 579   fFitdefs->AddHit(s);
 580   s.Reset(1);
 581  }
 582 }
 583 ///////////////////////////////////////////////////////////////////////////
 584 void IceF2k::PutMcTracks()
 585 {
 586 // Get the MC tracks from the F2000 file into the IcePack structure.
 587 // Note : MC tracks are given negative track id's in the event structure.
 588 // This memberfunction is based on the original code by Adam Bouchta.
 589
 590  IceEvent* evt=(IceEvent*)GetMainObject();
 591  if (!evt || fEvent.ntrack<=0) return;
 592
 593  // Loop over all the tracks and add them to the current event
 594  AliTrack t;
 595  Double_t vec[3];
 596  AliPosition r;
 597  Ali3Vector p;
 598  Int_t tid=0;
 599  Int_t idpdg=0;
 600  Int_t idf2k=0;
 601  for (Int_t i=0; i<fEvent.ntrack; i++)
 602  {
 603   t.Reset ();
 604
 605   // Beginpoint of the track
 606   vec[0]=fEvent.gen[i].x;
 607   vec[1]=fEvent.gen[i].y;
 608   vec[2]=fEvent.gen[i].z;
 609   r.SetPosition(vec,"car");
 610   t.SetBeginPoint(r);
 611
 612   // Endpoint of the track
 613   vec[0]+=fEvent.gen[i].length*fEvent.gen[i].px;
 614   vec[1]+=fEvent.gen[i].length*fEvent.gen[i].py;
 615   vec[2]+=fEvent.gen[i].length*fEvent.gen[i].pz;
 616   r.SetPosition(vec,"car");
 617   t.SetEndPoint(r);
 618
 619   // Momentum in GeV/c
 620   vec[0]=fEvent.gen[i].e*fEvent.gen[i].px*1e-3;
 621   vec[1]=fEvent.gen[i].e*fEvent.gen[i].py*1e-3;
 622   vec[2]=fEvent.gen[i].e*fEvent.gen[i].pz*1e-3;
 623   p.SetVector (vec,"car");
 624   t.Set3Momentum(p);
 625
 626   // MC tracks are indicated by negative track id's
 627   tid=fEvent.gen[i].tag;
 628   t.SetId(-abs(tid));
 629
 630   idf2k=fEvent.gen[i].id;
 631   idpdg=0;
 632   if (idf2k>1000)
 633   {
 634    idpdg=idf2k+10000000;
 635   }
 636   else if (idf2k <= 48)
 637   {
 638    idpdg=fPdg->ConvertGeant3ToPdg(idf2k);
 639   }
 640   else
 641   {
 642    if (idf2k==201) idpdg=12;
 643    if (idf2k==202) idpdg=14;
 644    if (idf2k==203) idpdg=16;
 645    if (idf2k==204) idpdg=-12;
 646    if (idf2k==205) idpdg=-14;
 647    if (idf2k==206) idpdg=-16;
 648   }
 649
 650   t.SetParticleCode(idpdg);
 651   t.SetName(fPdg->GetParticle(idpdg)->GetName());
 652   t.SetTitle("MC track");
 653   t.SetMass(fPdg->GetParticle(idpdg)->Mass());
 654   t.SetCharge(fPdg->GetParticle(idpdg)->Charge()/3.);
 655
 656   evt->AddTrack(t);
 657  }
 658
 659  // Create the pointers to each particle's parent particle.
 660  Int_t txid=0;
 661  Int_t parid=0;
 662  for (Int_t itk=1; itk<=evt->GetNtracks (); itk++)
 663  {
 664   AliTrack* tx=evt->GetTrack(itk);
 665
 666   if (!tx) continue;
 667
 668   txid=tx->GetId();
 669
 670   parid=-1;
 671   for (Int_t j=0; j<fEvent.ntrack; j++)
 672   {
 673    tid=fEvent.gen[j].tag;
 674    if (-abs(tid) == txid) parid=fEvent.gen[j].parent;
 675   }
 676
 677   if (parid<0) continue;
 678
 679   AliTrack* tpar=evt->GetIdTrack(-abs(parid));
 680
 681   if (!tpar) continue;
 682
 683   tx->SetParentTrack(tpar);
 684  }
 685 }
 686 ///////////////////////////////////////////////////////////////////////////
 687 void IceF2k::PutRecoTracks()
 688 {
 689 // Get the reconstructed tracks from the F2000 file into the IcePack structure.
 690 // Note : Reco tracks are given positive track id's in the event structure.
 691 // This memberfunction is based on the original code by Adam Bouchta.
 692
 693  IceEvent* evt=(IceEvent*)GetMainObject();
 694  if (!evt || fEvent.nfit<=0) return;
 695
 696  // Loop over all the tracks and add them to the current event
 697  AliTrack t;
 698  Double_t vec[3];
 699  AliPosition r;
 700  Ali3Vector p;
 701  Int_t tid=0;
 702  Int_t idpdg=0;
 703  Int_t idf2k=0;
 704  for (Int_t i=0; i<fEvent.nfit; i++)
 705  {
 706   t.Reset ();
 707
 708   // Beginpoint of the track
 709   vec[0]=fEvent.rec[i].x;
 710   vec[1]=fEvent.rec[i].y;
 711   vec[2]=fEvent.rec[i].z;
 712   r.SetPosition(vec,"car");
 713   t.SetBeginPoint(r);
 714
 715   // Endpoint of the track
 716   vec[0]+=fEvent.rec[i].length*fEvent.rec[i].px;
 717   vec[1]+=fEvent.rec[i].length*fEvent.rec[i].py;
 718   vec[2]+=fEvent.rec[i].length*fEvent.rec[i].pz;
 719   r.SetPosition(vec,"car");
 720   t.SetEndPoint(r);
 721
 722   // Momentum in GeV/c
 723   if (fEvent.rec[i].e > 0)
 724   {
 725    vec[0]=fEvent.rec[i].e*fEvent.rec[i].px*1e-3;
 726    vec[1]=fEvent.rec[i].e*fEvent.rec[i].py*1e-3;
 727    vec[2]=fEvent.rec[i].e*fEvent.rec[i].pz*1e-3;
 728   }
 729   else // Give the track a nominal momentum of 1 GeV/c
 730   {
 731    vec[0]=fEvent.rec[i].px;
 732    vec[1]=fEvent.rec[i].py;
 733    vec[2]=fEvent.rec[i].pz;
 734   }
 735   p.SetVector (vec,"car");
 736   t.Set3Momentum(p);
 737
 738   // Use the fit number as track id
 739   tid=fEvent.rec[i].tag;
 740   t.SetId(abs(tid));
 741
 742   idf2k=fEvent.rec[i].id;
 743   idpdg=0;
 744   if (idf2k>1000)
 745   {
 746    idpdg=idf2k+10000000;
 747   }
 748   else if (idf2k <= 48)
 749   {
 750    idpdg=fPdg->ConvertGeant3ToPdg(idf2k);
 751   }
 752   else
 753   {
 754    if (idf2k==201) idpdg=12;
 755    if (idf2k==202) idpdg=14;
 756    if (idf2k==203) idpdg=16;
 757    if (idf2k==204) idpdg=-12;
 758    if (idf2k==205) idpdg=-14;
 759    if (idf2k==206) idpdg=-16;
 760   }
 761
 762   t.SetParticleCode(idpdg);
 763   t.SetName(fPdg->GetParticle(idpdg)->GetName());
 764   t.SetTitle("RECO track");
 765   t.SetMass(fPdg->GetParticle(idpdg)->Mass());
 766   t.SetCharge(fPdg->GetParticle(idpdg)->Charge()/3.);
 767
 768   // Retrieve the various fit parameters for this track
 769   AliSignal* fitdata=fFitdefs->GetIdHit(i);
 770   for (Int_t jval=0; jval<fEvent.fresult[i].nval; jval++)
 771   {
 772    fitdata->SetSignal(fEvent.fresult[i].val[jval],jval+1);
 773   }
 774
 775   // Store the various fit parameters for this track
 776   t.SetFitDetails(fitdata);
 777
 778   // Store the various reco tracks as track hypotheses.
 779   // A copy of the first reco track is entered as a new track instance
 780   // into the event and all reco tracks (incl. the first one) are
 781   // stored as hypotheses linked to this new reco track.
 782   if (i==0)
 783   {
 784    evt->AddTrack(t);
 785    AliTrack* tx=evt->GetTrack(evt->GetNtracks());
 786    Int_t nrec=evt->GetNtracks(1);
 787    tx->SetId(nrec+1);
 788   }
 789   AliTrack* tx=evt->GetTrack(evt->GetNtracks());
 790   if (tx) tx->AddTrackHypothesis(t);
 791  }
 792 }
 793 ///////////////////////////////////////////////////////////////////////////
 794 void IceF2k::PutHits()
 795 {
 796 // Get the hit and waveform info from the F2000 file into the IcePack structure.
 797 // This memberfunction is based on the original code by Adam Bouchta.
 798
 799  IceEvent* evt=(IceEvent*)GetMainObject();
 800  if (!evt) return;
 801
 802  // Loop over all the hits and add them to the current event
 803  IceAOM om;
 804  AliSignal s;
 805  s.SetSlotName("ADC",1);
 806  s.SetSlotName("LE",2);
 807  s.SetSlotName("TOT",3);
 808  Int_t chan=0;
 809  Int_t maxchan=800;
 810  if (fOmdb) maxchan=fHeader.nch;
 811  IceAOM* omx=0;
 812  AliSignal* sx=0;
 813  Int_t tid=0;
 814  AliTrack* tx=0;
 815  Float_t adc=0;
 816  for (Int_t i=0; i<fEvent.nhits; i++)
 817  {
 818   chan=fEvent.h[i].ch+1;
 819   if (chan>maxchan) continue; // Channels 9001, 9002 etc are trigger channels
 820
 821   // Get corresponding device from the current event structure
 822   omx=(IceAOM*)evt->GetIdDevice(chan);
 823   if (!omx)
 824   {
 825    if (fOmdb)
 826    {
 827     omx=(IceAOM*)fOmdb->GetObject(chan,1);
 828     evt->AddDevice(omx);
 829    }
 830    else
 831    {
 832     om.Reset(1);
 833     om.SetUniqueID(chan);
 834     evt->AddDevice(om);
 835    }
 836    omx=(IceAOM*)evt->GetIdDevice(chan);
 837   }
 838
 839   if (!omx) continue;
 840
 841   s.Reset();
 842   s.SetUniqueID(fEvent.h[i].id);
 843   s.SetSignal(fEvent.h[i].amp,1);
 844   s.SetSignal(fEvent.h[i].t,2);
 845   s.SetSignal(fEvent.h[i].tot,3);
 846
 847   omx->AddHit(s);
 848
 849   sx=omx->GetHit(omx->GetNhits());
 850   if (!sx) continue;
 851
 852   // ADC dependent TDC (de)calibration function for this hit
 853   TF1* fcal=omx->GetCalFunction("LE");
 854   TF1* fdecal=omx->GetDecalFunction("LE");
 855   if (fcal) sx->SetCalFunction(fcal,2);
 856   if (fdecal) sx->SetDecalFunction(fdecal,2);
 857   fcal=sx->GetCalFunction(2);
 858   fdecal=sx->GetDecalFunction(2);
 859   adc=sx->GetSignal(1,-4);
 860   if (adc>0)
 861   {
 862    if (fcal) fcal->SetParameter(3,adc);
 863    if (fdecal) fdecal->SetParameter(3,adc);
 864   }
 865   else
 866   {
 867    if (fcal) fcal->SetParameter(3,1.e20);
 868    if (fdecal) fdecal->SetParameter(3,1.e20);
 869   }
 870
 871   // Bi-directional link between this hit and the track that caused the ADC value.
 872   // This F2K info is probably only present for MC tracks.
 873   tid=fEvent.h[i].ma;
 874   if (tid > 0)
 875   {
 876    tx=evt->GetIdTrack(tid); // Reco tracks
 877    if (!tx) tx=evt->GetIdTrack(-tid); // MC tracks
 878    if (tx) sx->AddLink(tx);
 879   }
 880   else
 881   {
 882    if (tid == -2) s.SetNameTitle("N","Noise");
 883    if (tid == -3) s.SetNameTitle("A","Afterpulse");
 884   }
 885  }
 886
 887  // Loop over all the waveforms and add the histo(s) to the corresponding OM's
 888  TH1F histo;
 889  Int_t nbins=0;
 890  Float_t xlow=0;
 891  Float_t xup=0;
 892  TString hname;
 893  for (Int_t iwf=0; iwf<fEvent.nwf; iwf++)
 894  {
 895   chan=fEvent.wf[iwf].om;
 896   if (chan<=0 || chan>maxchan) continue; // Skip trigger channels
 897
 898   // Get corresponding device from the current event structure
 899   omx=(IceAOM*)evt->GetIdDevice(chan);
 900   if (!omx)
 901   {
 902    if (fOmdb)
 903    {
 904     omx=(IceAOM*)fOmdb->GetObject(chan,1);
 905     evt->AddDevice(omx);
 906    }
 907    else
 908    {
 909     om.Reset(1);
 910     om.SetUniqueID(chan);
 911     evt->AddDevice(om);
 912    }
 913    omx=(IceAOM*)evt->GetIdDevice(chan);
 914   }
 915
 916   if (!omx) continue;
 917
 918   omx->SetSlotName("BASELINE",omx->GetNnames()+1);
 919   omx->SetSignal(fEvent.wf[iwf].baseline,"BASELINE");
 920
 921   // Fill the waveform histogram
 922   hname="OM";
 923   hname+=chan;
 924   hname+="-WF";
 925   hname+=omx->GetNwaveforms()+1;
 926
 927   histo.Reset();
 928   histo.SetName(hname.Data());
 929   nbins=fEvent.wf[iwf].ndigi;
 930   xlow=fEvent.wf[iwf].t_start;
 931   xup=xlow+float(nbins)*fEvent.wf[iwf].t_bin;
 932   histo.SetBins(nbins,xlow,xup);
 933
 934   for (Int_t jbin=1; jbin<=fEvent.wf[iwf].ndigi; jbin++)
 935   {
 936    histo.SetBinContent(jbin,fEvent.wf[iwf].digi[jbin-1]);
 937   }
 938
 939   omx->SetWaveform(&histo,omx->GetNwaveforms()+1);
 940  }
 941
 942  // Set bi-directional links between hits and reco track hypotheses.
 943  // Note : Reco tracks are recognised by their positive id.
 944  Int_t hid=0;
 945  TObjArray* rectracks=evt->GetTracks(1);
 946  for (Int_t jtk=0; jtk<rectracks->GetEntries(); jtk++)
 947  {
 948   tx=(AliTrack*)rectracks->At(jtk);
 949   if (!tx) continue;
 950
 951   for (Int_t jhyp=1; jhyp<=tx->GetNhypotheses(); jhyp++)
 952   {
 953    AliTrack* hypx=tx->GetTrackHypothesis(jhyp);
 954    if (!hypx) continue;
 955
 956    // Loop over all combinations of F2K fits and used OM hits
 957    for (Int_t k=0; k<fEvent.nfit_uses; k++)
 958    {
 959     if (fEvent.fit_uses[k].useid != hypx->GetId()) continue;
 960     hid=fEvent.fit_uses[k].hitid;
 961     sx=evt->GetIdHit(hid,"IceAOM");
 962     if (sx) sx->AddLink(hypx);
 963    }
 964   }
 965  }
 966 }
 967 ///////////////////////////////////////////////////////////////////////////