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 //
  28 // Usage example :
  29 // ---------------
  30 //
  31 // gSystem->Load("ralice");
  32 // gSystem->Load("icepack");
  33 // gSystem->Load("iceconvert");
  34 //
  35 // IceF2k q("IceF2k","F2K to IcePack data structure conversion");
  36 //
  37 // // Limit the number of entries for testing
  38 // q.SetMaxEvents(10);
  39 //
  40 // // Print frequency to produce a short summary print every printfreq events
  41 // q.SetPrintFreq(1);
  42 //
  43 // // Split level for the output structures
  44 // q.SetSplitLevel(2);
  45 //
  46 // // Buffer size for the output structures
  47 // q.SetBufferSize(32000);
  48 //
  49 // // The F2K input filename
  50 // q.SetInputFile("run7825.f2k");
  51 //
  52 // // Output file for the event structures
  53 // TFile* ofile=new TFile("events.root","RECREATE","F2K data in IceEvent structure");
  54 // q.SetOutputFile(ofile);
  55 //
  56 // ///////////////////////////////////////////////////////////////////
  57 // // Here the user can specify his/her sub-tasks to be executed
  58 // // on an event-by-event basis after the IceEvent structure
  59 // // has been filled and before the data is written out.
  60 // // Sub-tasks (i.e. a user classes derived from TTask) are entered
  61 // // as follows :
  62 // //
  63 // //    MyXtalk task1("task1","Cross talk correction");
  64 // //    MyClean task2("task2","Hit cleaning");
  65 // //    q.Add(&task1);
  66 // //    q.Add(&task2);
  67 // //
  68 // // The sub-tasks will be executed in the order as they are entered.
  69 // ///////////////////////////////////////////////////////////////////
  70 //
  71 // // Perform the conversion and execute subtasks (if any)
  72 // // on an event-by-event basis
  73 // q.ExecuteJob();
  74 //
  75 // // Select various objects to be added to the output file
  76 //
  77 // AliObjMatrix* omdb=q.GetOMdbase();
  78 // if (omdb) omdb->Write();
  79 //
  80 // AliDevice* fitdefs=q.GetFitdefs();
  81 // if (fitdefs) fitdefs->Write();
  82 //
  83 // TDatabasePDG* pdg=q.GetPDG();
  84 // if (pdg) pdg->Write();
  85 //
  86 // // Close output file
  87 // ofile->Write();
  88 // ofile->Close();
  89 //
  90 //--- Author: Nick van Eijndhoven 11-mar-2005 Utrecht University
  91 //- Modified: NvE $Date$ Utrecht University
  92 ///////////////////////////////////////////////////////////////////////////
  93
  94 #include "IceF2k.h"
  95 #include "Riostream.h"
  96
  97 ClassImp(IceF2k) // Class implementation to enable ROOT I/O
  98
  99 IceF2k::IceF2k(const char* name,const char* title) : AliJob(name,title)
 100 {
 101 // Default constructor.
 102 // By default maxevent=-1, split=99, bsize=32000, printfreq=1.
 103
 104  fSplit=99;
 105  fBsize=32000;
 106  fMaxevt=-1;
 107  fPrintfreq=1;
 108  fInfile="";
 109  fOutfile=0;
 110
 111  fPdg=0;
 112  fOmdb=0;
 113  fFitdefs=0;
 114 }
 115 ///////////////////////////////////////////////////////////////////////////
 116 IceF2k::~IceF2k()
 117 {
 118 // Default destructor.
 119  IceEvent* evt=(IceEvent*)GetMainObject();
 120  if (evt) delete evt;
 121
 122  if (fPdg)
 123  {
 124   delete fPdg;
 125   fPdg=0;
 126  }
 127
 128  if (fOmdb)
 129  {
 130   delete fOmdb;
 131   fOmdb=0;
 132  }
 133
 134  if (fFitdefs)
 135  {
 136   delete fFitdefs;
 137   fFitdefs=0;
 138  }
 139 }
 140 ///////////////////////////////////////////////////////////////////////////
 141 void IceF2k::SetMaxEvents(Int_t n)
 142 {
 143 // Set the maximum number of events to be processed.
 144 // n=-1 implies processing of the complete input file, which is the default
 145 // initialisation in the constructor.
 146  fMaxevt=n;
 147 }
 148 ///////////////////////////////////////////////////////////////////////////
 149 void IceF2k::SetPrintFreq(Int_t f)
 150 {
 151 // Set the printfrequency to produce info every f events.
 152 // f=1 is the default initialisation in the constructor.
 153  if (f>0) fPrintfreq=f;
 154 }
 155 ///////////////////////////////////////////////////////////////////////////
 156 void IceF2k::SetSplitLevel(Int_t split)
 157 {
 158 // Set the split level for the ROOT data file.
 159 // split=99 is the default initialisation in the constructor.
 160  if (split>=0) fSplit=split;
 161 }
 162 ///////////////////////////////////////////////////////////////////////////
 163 void IceF2k::SetBufferSize(Int_t bsize)
 164 {
 165 // Set the buffer size for the ROOT data file.
 166 // bsize=32000 is the default initialisation in the constructor.
 167  if (bsize>=0) fBsize=bsize;
 168 }
 169 ///////////////////////////////////////////////////////////////////////////
 170 void IceF2k::SetInputFile(TString name)
 171 {
 172 // Set the name of the F2K input file.
 173  fInfile=name;
 174 }
 175 ///////////////////////////////////////////////////////////////////////////
 176 void IceF2k::SetOutputFile(TFile* ofile)
 177 {
 178 // Set the output file for the ROOT data.
 179  fOutfile=ofile;
 180 }
 181 ///////////////////////////////////////////////////////////////////////////
 182 TDatabasePDG* IceF2k::GetPDG()
 183 {
 184 // Provide pointer to the PDG database
 185  return fPdg;
 186 }
 187 ///////////////////////////////////////////////////////////////////////////
 188 AliObjMatrix* IceF2k::GetOMdbase()
 189 {
 190 // Provide pointer to the OM geometry, calib. etc... database
 191  return fOmdb;
 192 }
 193 ///////////////////////////////////////////////////////////////////////////
 194 AliDevice* IceF2k::GetFitdefs()
 195 {
 196 // Provide pointer to the fit definitions
 197  return fFitdefs;
 198 }
 199 ///////////////////////////////////////////////////////////////////////////
 200 void IceF2k::Exec(Option_t* opt)
 201 {
 202 // Job to loop over the specified number of events and convert the
 203 // F2K data into the IceEvent structure.
 204 // If maxevents<0 (default) all the entries of the input file
 205 // will be processed.
 206 // Every "printfreq" events a short event summary will be printed.
 207 // The default value is printfreq=1.
 208 // The output will be written on a standard output tree named "T".
 209 //
 210 // Notes :
 211 // -------
 212 // 1) This class is derived from AliJob, allowing a task based processing.
 213 //    After the conversion of an F2K event into an IceEvent structure,
 214 //    the processing of all available sub-tasks (if any) is invoked.
 215 //    This provides an event-by-event (sub)task processing before the
 216 //    final data structures are written out.
 217 // 2) The main object in this job environment is an IceEvent* pointer.
 218
 219  if (fInfile=="")
 220  {
 221   cout << " *IceF2k Exec* No data input file specified." << endl;
 222   return;
 223  }
 224
 225  // Open the input file in the default ascii format (autodetection) for reading
 226  fInput=rdmc_mcopen(fInfile.Data(),"r",RDMC_DEFAULT_ASCII_F);
 227
 228  if (!fInput)
 229  {
 230   cout << " *IceF2k Exec* No input file found with name : " << fInfile.Data() << endl;
 231   return;
 232  }
 233
 234  // Initialise the event structure
 235  rdmc_init_mevt(&fEvent);
 236
 237  // Read the file header information
 238  rdmc_rarr(fInput,&fHeader);
 239
 240  if (!fOutfile)
 241  {
 242   cout << " *IceF2k Exec* No ROOT output file specified." << endl;
 243   return;
 244  }
 245
 246  TTree* otree=new TTree("T","F2K Data converted to IceEvent structures");
 247
 248  IceEvent* evt=new IceEvent();
 249  evt->SetTrackCopy(1);
 250  evt->SetDevCopy(1);
 251
 252  // Branch in the tree for the event structure
 253  otree->Branch("IceEvent","IceEvent",&evt,fBsize,fSplit);
 254
 255  // Create the particle database and extend it with some F2000 specific definitions
 256  if (!fPdg) fPdg=new TDatabasePDG();
 257  Double_t me=fPdg->GetParticle(11)->Mass();
 258  fPdg->AddParticle("brems"   ,"brems"   ,0,1,0,0,"none",10001001,0,0);
 259  fPdg->AddParticle("deltae"  ,"deltae"  ,me,1,0,-3,"Lepton",10001002,0,0);
 260  fPdg->AddParticle("pairprod","pairprod",0,1,0,0,"none",10001003,0,0);
 261  fPdg->AddParticle("nucl_int","nucl_Int",0,1,0,0,"none",10001004,0,0);
 262  fPdg->AddParticle("mu_pair" ,"mu_pair" ,0,1,0,0,"none",10001005,0,0);
 263  fPdg->AddParticle("hadrons" ,"hadrons" ,0,1,0,0,"none",10001006,0,0);
 264  fPdg->AddParticle("fiberlaser","fiberlaser",0,1,0,0,"none",10002100,0,0);
 265  fPdg->AddParticle("n2laser"   ,"n2laser"   ,0,1,0,0,"none",10002101,0,0);
 266  fPdg->AddParticle("yaglaser"  ,"yaglaser"  ,0,1,0,0,"none",10002201,0,0);
 267  fPdg->AddParticle("z_primary","z_primary",0,1,0,0,"none",10003000,0,0);
 268  fPdg->AddParticle("a_primary","a_primary",0,1,0,0,"none",10003500,0,0);
 269
 270  // Fill the database with geometry, calib. etc... parameters
 271  // for all the devices
 272  FillOMdbase();
 273
 274  // Set the fit definitions according to the F2000 header info
 275  SetFitdefs();
 276
 277  // Initialise the job working environment
 278  SetMainObject(evt);
 279  AddObject(fOutfile);
 280  AddObject(otree);
 281
 282  cout << " ***" << endl;
 283  cout << " *** Start processing of job " << GetName() << " ***" << endl;
 284  cout << " ***" << endl;
 285  cout << " F2K input file : " << fInfile.Data() << endl;
 286  cout << " Maximum number of events to be processed : " << fMaxevt << endl;
 287  cout << " Print frequency : " << fPrintfreq << endl;
 288  cout << " ROOT output file : " << fOutfile->GetName() << endl;
 289  cout << " Output characteristics : splitlevel = " << fSplit << " buffersize = " << fBsize << endl;
 290
 291  ListEnvironment();
 292
 293  Int_t nevt=0;
 294  while (!rdmc_revt(fInput,&fHeader,&fEvent))
 295  {
 296   if (fMaxevt>-1 && nevt>=fMaxevt) break;
 297
 298   // Reset the complete Event structure
 299   evt->Reset();
 300
 301   evt->SetRunNumber(fEvent.nrun);
 302   evt->SetEventNumber(fEvent.enr);
 303   evt->SetMJD(fEvent.mjd,fEvent.secs,fEvent.nsecs);
 304
 305   PutMcTracks();
 306
 307   PutRecoTracks();
 308
 309   PutHits();
 310
 311   // Invoke all available sub-tasks (if any)
 312   ExecuteTasks(opt);
 313
 314   if (!(nevt%fPrintfreq)) evt->HeaderData();
 315
 316   // Write the complete structure to the output Tree
 317   otree->Fill();
 318
 319   // Update event counter
 320   nevt++;
 321  }
 322 }
 323 ///////////////////////////////////////////////////////////////////////////
 324 void IceF2k::FillOMdbase()
 325 {
 326 // Fill the database with geometry, calib. etc... parameters
 327 // for all the devices.
 328
 329  if (fHeader.nch<=0) return;
 330
 331  if (fOmdb)
 332  {
 333   fOmdb->Reset();
 334  }
 335  else
 336  {
 337   fOmdb=new AliObjMatrix();
 338   fOmdb->SetNameTitle("OMDBASE","The OM geometry, calib. etc... database");
 339   fOmdb->SetOwner();
 340  }
 341
 342  IceAOM* dev=0;
 343  Double_t pos[3]={0,0,0};
 344  for (Int_t i=0; i<fHeader.nch; i++)
 345  {
 346   dev=new IceAOM();
 347   dev->SetUniqueID(i+1);
 348   dev->SetSlotName("TYPE",1);
 349   dev->SetSlotName("ORIENT",2);
 350   dev->SetSlotName("T0",3);
 351   dev->SetSlotName("ALPHA",4);
 352   dev->SetSlotName("KADC",5);
 353   dev->SetSlotName("KTOT",6);
 354   dev->SetSlotName("KTDC",7);
 355
 356   pos[0]=fHeader.x[i];
 357   pos[1]=fHeader.y[i];
 358   pos[2]=fHeader.z[i];
 359   dev->SetPosition(pos,"car");
 360   dev->SetSignal(fHeader.type[i],1);
 361   dev->SetSignal((Float_t)fHeader.costh[i],2);
 362   dev->SetSignal(fHeader.cal[i].t_0,3);
 363   dev->SetSignal(fHeader.cal[i].alpha_t,4);
 364   dev->SetSignal(fHeader.cal[i].beta_a,5);
 365   dev->SetSignal(fHeader.cal[i].beta_tot,6);
 366   dev->SetSignal(fHeader.cal[i].beta_t,7);
 367   fOmdb->EnterObject(i+1,1,dev);
 368  }
 369 }
 370 ///////////////////////////////////////////////////////////////////////////
 371 void IceF2k::SetFitdefs()
 372 {
 373 // Obtain the names of the variables for each fit procedure from the
 374 // f2000 header. Each different fit procedure is then stored as a separate
 375 // hit of an AliDevice object and the various fit variables are stored
 376 // as separate signal slots of the corresponding hit.
 377 // Via the GetFitdefs() memberfunction this AliDevice object can be
 378 // retrieved and stored in the ROOT output file if wanted.
 379 // The name of the object is FitDefinitions and the stored data can be
 380 // inspected via the AliDevice::Data() memberfunction and looks as follows :
 381 //
 382 //  *AliDevice::Data* Id :0 Name : FitDefinitions
 383 //    Position Vector in car coordinates : 0 0 0
 384 //    Err. in car coordinates : 0 0 0
 385 //  The following 8 hits are registered :
 386 //  *AliSignal::Data* Id :0
 387 //    Position Vector in car coordinates : 0 0 0
 388 //    Err. in car coordinates : 0 0 0
 389 //    Owned by device : AliDevice Name : FitDefinitions
 390 //    Slot : 1 Signal value : 1 name : id
 391 //    Slot : 2 Signal value : 2 name : rchi2
 392 //    Slot : 3 Signal value : 3 name : prob
 393 //    Slot : 4 Signal value : 4 name : sigth
 394 //    Slot : 5 Signal value : 5 name : covmin
 395 //    Slot : 6 Signal value : 6 name : covmax
 396 //    Slot : 7 Signal value : 7 name : cutflag
 397 //    Slot : 8 Signal value : 8 name : chi2
 398 //  *AliSignal::Data* Id :1
 399 //    Position Vector in car coordinates : 0 0 0
 400 //    Err. in car coordinates : 0 0 0
 401 //    Owned by device : AliDevice Name : FitDefinitions
 402 //    Slot : 1 Signal value : 1 name : id
 403 //    Slot : 2 Signal value : 2 name : rchi2
 404 //    Slot : 3 Signal value : 3 name : prob
 405 // etc....
 406 //
 407 // This memberfunction is based on the original idea/code by Adam Bouchta.
 408
 409  if (fHeader.n_fit<=0) return;
 410
 411  if (fFitdefs)
 412  {
 413   fFitdefs->Reset(1);
 414  }
 415  else
 416  {
 417   fFitdefs=new AliDevice();
 418  }
 419
 420  fFitdefs->SetName("FitDefinitions");
 421  fFitdefs->SetHitCopy (1);
 422
 423  AliSignal s;
 424  s.Reset();
 425
 426  for (Int_t i=0; i<fHeader.n_fit; i++)
 427  {
 428   s.SetUniqueID(fHeader.def_fit[i].id);
 429
 430   for (Int_t j=0; j<fHeader.def_fit[i].nwords; j++)
 431   {
 432    s.SetSlotName(TString(fHeader.def_fit[i].words[j]),j+1);
 433    s.SetSignal(j+1,j+1);
 434   }
 435
 436   fFitdefs->AddHit(s);
 437   s.Reset(1);
 438  }
 439 }
 440 ///////////////////////////////////////////////////////////////////////////
 441 void IceF2k::PutMcTracks()
 442 {
 443 // Get the MC tracks from the F2000 file into the IcePack structure.
 444 // Note : MC tracks are given negative track id's in the event structure.
 445 // This memberfunction is based on the original code by Adam Bouchta.
 446
 447  IceEvent* evt=(IceEvent*)GetMainObject();
 448  if (!evt || fEvent.ntrack<=0) return;
 449
 450  // Loop over all the tracks and add them to the current event
 451  AliTrack t;
 452  Double_t vec[3];
 453  AliPosition r;
 454  Ali3Vector p;
 455  Int_t tid=0;
 456  Int_t idpdg=0;
 457  Int_t idf2k=0;
 458  for (Int_t i=0; i<fEvent.ntrack; i++)
 459  {
 460   t.Reset ();
 461
 462   // Beginpoint of the track
 463   vec[0]=fEvent.gen[i].x;
 464   vec[1]=fEvent.gen[i].y;
 465   vec[2]=fEvent.gen[i].z;
 466   r.SetPosition(vec,"car");
 467   t.SetBeginPoint(r);
 468
 469   // Endpoint of the track
 470   vec[0]+=fEvent.gen[i].length*fEvent.gen[i].px;
 471   vec[1]+=fEvent.gen[i].length*fEvent.gen[i].py;
 472   vec[2]+=fEvent.gen[i].length*fEvent.gen[i].pz;
 473   r.SetPosition(vec,"car");
 474   t.SetEndPoint(r);
 475
 476   // Momentum in GeV/c
 477   vec[0]=fEvent.gen[i].e*fEvent.gen[i].px*1e-3;
 478   vec[1]=fEvent.gen[i].e*fEvent.gen[i].py*1e-3;
 479   vec[2]=fEvent.gen[i].e*fEvent.gen[i].pz*1e-3;
 480   p.SetVector (vec,"car");
 481   t.Set3Momentum(p);
 482
 483   // MC tracks are indicated by negative track id's
 484   tid=fEvent.gen[i].tag;
 485   t.SetId(-abs(tid));
 486
 487   idf2k=fEvent.gen[i].id;
 488   idpdg=0;
 489   if (idf2k>1000)
 490   {
 491    idpdg=idf2k+10000000;
 492   }
 493   else if (idf2k <= 48)
 494   {
 495    idpdg=fPdg->ConvertGeant3ToPdg(idf2k);
 496   }
 497   else
 498   {
 499    if (idf2k==201) idpdg=12;
 500    if (idf2k==202) idpdg=14;
 501    if (idf2k==203) idpdg=16;
 502    if (idf2k==204) idpdg=-12;
 503    if (idf2k==205) idpdg=-14;
 504    if (idf2k==206) idpdg=-16;
 505   }
 506
 507   t.SetParticleCode(idpdg);
 508   t.SetName(fPdg->GetParticle(idpdg)->GetName());
 509   t.SetTitle("MC track");
 510   t.SetMass(fPdg->GetParticle(idpdg)->Mass());
 511   t.SetCharge(fPdg->GetParticle(idpdg)->Charge()/3.);
 512
 513   evt->AddTrack(t);
 514  }
 515
 516  // Create the pointers to each particle's parent particle.
 517  Int_t txid=0;
 518  Int_t parid=0;
 519  for (Int_t itk=1; itk<=evt->GetNtracks (); itk++)
 520  {
 521   AliTrack* tx=evt->GetTrack(itk);
 522
 523   if (!tx) continue;
 524
 525   txid=tx->GetId();
 526
 527   parid=-1;
 528   for (Int_t j=0; j<fEvent.ntrack; j++)
 529   {
 530    tid=fEvent.gen[j].tag;
 531    if (-abs(tid) == txid) parid=fEvent.gen[j].parent;
 532   }
 533
 534   if (parid<0) continue;
 535
 536   AliTrack* tpar=evt->GetIdTrack(-abs(parid));
 537
 538   if (!tpar) continue;
 539
 540   tx->SetParentTrack(tpar);
 541  }
 542 }
 543 ///////////////////////////////////////////////////////////////////////////
 544 void IceF2k::PutRecoTracks()
 545 {
 546 // Get the reconstructed tracks from the F2000 file into the IcePack structure.
 547 // Note : Reco tracks are given positive track id's in the event structure.
 548 // This memberfunction is based on the original code by Adam Bouchta.
 549
 550  IceEvent* evt=(IceEvent*)GetMainObject();
 551  if (!evt || fEvent.nfit<=0) return;
 552
 553  // Loop over all the tracks and add them to the current event
 554  AliTrack t;
 555  Double_t vec[3];
 556  AliPosition r;
 557  Ali3Vector p;
 558  Int_t tid=0;
 559  Int_t idpdg=0;
 560  Int_t idf2k=0;
 561  for (Int_t i=0; i<fEvent.nfit; i++)
 562  {
 563   t.Reset ();
 564
 565   // Beginpoint of the track
 566   vec[0]=fEvent.rec[i].x;
 567   vec[1]=fEvent.rec[i].y;
 568   vec[2]=fEvent.rec[i].z;
 569   r.SetPosition(vec,"car");
 570   t.SetBeginPoint(r);
 571
 572   // Endpoint of the track
 573   vec[0]+=fEvent.rec[i].length*fEvent.rec[i].px;
 574   vec[1]+=fEvent.rec[i].length*fEvent.rec[i].py;
 575   vec[2]+=fEvent.rec[i].length*fEvent.rec[i].pz;
 576   r.SetPosition(vec,"car");
 577   t.SetEndPoint(r);
 578
 579   // Momentum in GeV/c
 580   if (fEvent.rec[i].e > 0)
 581   {
 582    vec[0]=fEvent.rec[i].e*fEvent.rec[i].px*1e-3;
 583    vec[1]=fEvent.rec[i].e*fEvent.rec[i].py*1e-3;
 584    vec[2]=fEvent.rec[i].e*fEvent.rec[i].pz*1e-3;
 585   }
 586   else // Give the track a nominal momentum of 1 GeV/c
 587   {
 588    vec[0]=fEvent.rec[i].px;
 589    vec[1]=fEvent.rec[i].py;
 590    vec[2]=fEvent.rec[i].pz;
 591   }
 592   p.SetVector (vec,"car");
 593   t.Set3Momentum(p);
 594
 595   // Use the fit number as track id
 596   tid=fEvent.rec[i].tag;
 597   t.SetId(abs(tid));
 598
 599   idf2k=fEvent.rec[i].id;
 600   idpdg=0;
 601   if (idf2k>1000)
 602   {
 603    idpdg=idf2k+10000000;
 604   }
 605   else if (idf2k <= 48)
 606   {
 607    idpdg=fPdg->ConvertGeant3ToPdg(idf2k);
 608   }
 609   else
 610   {
 611    if (idf2k==201) idpdg=12;
 612    if (idf2k==202) idpdg=14;
 613    if (idf2k==203) idpdg=16;
 614    if (idf2k==204) idpdg=-12;
 615    if (idf2k==205) idpdg=-14;
 616    if (idf2k==206) idpdg=-16;
 617   }
 618
 619   t.SetParticleCode(idpdg);
 620   t.SetName(fPdg->GetParticle(idpdg)->GetName());
 621   t.SetTitle("RECO track");
 622   t.SetMass(fPdg->GetParticle(idpdg)->Mass());
 623   t.SetCharge(fPdg->GetParticle(idpdg)->Charge()/3.);
 624
 625   // Retrieve the various fit parameters for this track
 626   AliSignal* fitdata=fFitdefs->GetIdHit(i);
 627   for (Int_t jval=0; jval<fEvent.fresult[i].nval; jval++)
 628   {
 629    fitdata->SetSignal(fEvent.fresult[i].val[jval],jval+1);
 630   }
 631
 632   // Store the various fit parameters for this track
 633   t.SetFitDetails(fitdata);
 634
 635   // Store the various reco tracks as track hypotheses.
 636   // A copy of the first reco track is entered as a new track instance
 637   // into the event and all reco tracks (incl. the first one) are
 638   // stored as hypotheses linked to this new reco track.
 639   if (i==0)
 640   {
 641    evt->AddTrack(t);
 642    AliTrack* tx=evt->GetTrack(evt->GetNtracks());
 643    Int_t nrec=evt->GetNtracks(1);
 644    tx->SetId(nrec+1);
 645   }
 646   AliTrack* tx=evt->GetTrack(evt->GetNtracks());
 647   if (tx) tx->AddTrackHypothesis(t);
 648  }
 649 }
 650 ///////////////////////////////////////////////////////////////////////////
 651 void IceF2k::PutHits()
 652 {
 653 // Get the hit and waveform info from the F2000 file into the IcePack structure.
 654 // This memberfunction is based on the original code by Adam Bouchta.
 655
 656  IceEvent* evt=(IceEvent*)GetMainObject();
 657  if (!evt) return;
 658
 659  // Loop over all the hits and add them to the current event
 660  IceAOM om;
 661  AliSignal s;
 662  s.SetSlotName("ADC",1);
 663  s.SetSlotName("LE",2);
 664  s.SetSlotName("TOT",3);
 665  Int_t chan=0;
 666  Int_t maxchan=800;
 667  if (fOmdb) maxchan=fHeader.nch;
 668  IceAOM* omx=0;
 669  AliSignal* sx=0;
 670  Int_t tid=0;
 671  AliTrack* tx=0;
 672  for (Int_t i=0; i<fEvent.nhits; i++)
 673  {
 674   chan=fEvent.h[i].ch+1;
 675   if (chan>maxchan) continue; // Channels 9001, 9002 etc are trigger channels
 676
 677   // Get corresponding device from the current event structure
 678   omx=(IceAOM*)evt->GetIdDevice(chan);
 679   if (!omx)
 680   {
 681    if (fOmdb)
 682    {
 683     omx=(IceAOM*)fOmdb->GetObject(chan,1);
 684     evt->AddDevice(omx);
 685    }
 686    else
 687    {
 688     om.Reset(1);
 689     om.SetUniqueID(chan);
 690     evt->AddDevice(om);
 691    }
 692    omx=(IceAOM*)evt->GetIdDevice(chan);
 693   }
 694
 695   if (!omx) continue;
 696
 697   s.Reset();
 698   s.SetUniqueID(fEvent.h[i].id);
 699   s.SetSignal(fEvent.h[i].amp,1);
 700   s.SetSignal(fEvent.h[i].t,2);
 701   s.SetSignal(fEvent.h[i].tot,3);
 702
 703   omx->AddHit(s);
 704
 705   sx=omx->GetHit(omx->GetNhits());
 706   if (!sx) continue;
 707
 708   // Bi-directional link between this hit and the track that caused the ADC value.
 709   // This F2K info is probably only present for MC tracks.
 710   tid=fEvent.h[i].ma;
 711   if (tid > 0)
 712   {
 713    tx=evt->GetIdTrack(tid); // Reco tracks
 714    if (!tx) tx=evt->GetIdTrack(-tid); // MC tracks
 715    if (tx) sx->AddLink(tx);
 716   }
 717   else
 718   {
 719    if (tid == -2) s.SetNameTitle("N","Noise");
 720    if (tid == -3) s.SetNameTitle("A","Afterpulse");
 721   }
 722  }
 723
 724  // Loop over all the waveforms and add the histo(s) to the corresponding OM's
 725  TH1F histo;
 726  Int_t nbins=0;
 727  Float_t xlow=0;
 728  Float_t xup=0;
 729  TString hname;
 730  for (Int_t iwf=0; iwf<fEvent.nwf; iwf++)
 731  {
 732   chan=fEvent.wf[iwf].om;
 733   if (chan<=0 || chan>maxchan) continue; // Skip trigger channels
 734
 735   // Get corresponding device from the current event structure
 736   omx=(IceAOM*)evt->GetIdDevice(chan);
 737   if (!omx)
 738   {
 739    if (fOmdb)
 740    {
 741     omx=(IceAOM*)fOmdb->GetObject(chan,1);
 742     evt->AddDevice(omx);
 743    }
 744    else
 745    {
 746     om.Reset(1);
 747     om.SetUniqueID(chan);
 748     evt->AddDevice(om);
 749    }
 750    omx=(IceAOM*)evt->GetIdDevice(chan);
 751   }
 752
 753   if (!omx) continue;
 754
 755   omx->SetSlotName("BASELINE",omx->GetNnames()+1);
 756   omx->SetSignal(fEvent.wf[iwf].baseline,"BASELINE");
 757
 758   // Fill the waveform histogram
 759   hname="OM";
 760   hname+=chan;
 761   hname+="-WF";
 762   hname+=omx->GetNwaveforms()+1;
 763
 764   histo.Reset();
 765   histo.SetName(hname.Data());
 766   nbins=fEvent.wf[iwf].ndigi;
 767   xlow=fEvent.wf[iwf].t_start;
 768   xup=xlow+float(nbins)*fEvent.wf[iwf].t_bin;
 769   histo.SetBins(nbins,xlow,xup);
 770
 771   for (Int_t jbin=1; jbin<=fEvent.wf[iwf].ndigi; jbin++)
 772   {
 773    histo.SetBinContent(jbin,fEvent.wf[iwf].digi[jbin]);
 774   }
 775
 776   omx->SetWaveform(&histo,omx->GetNwaveforms()+1);
 777  }
 778
 779  // Set bi-directional links between hits and reco track hypotheses.
 780  // Note : Reco tracks are recognised by their positive id.
 781  Int_t hid=0;
 782  TObjArray* rectracks=evt->GetTracks(1);
 783  for (Int_t jtk=0; jtk<rectracks->GetEntries(); jtk++)
 784  {
 785   tx=(AliTrack*)rectracks->At(jtk);
 786   if (!tx) continue;
 787
 788   for (Int_t jhyp=1; jhyp<=tx->GetNhypotheses(); jhyp++)
 789   {
 790    AliTrack* hypx=tx->GetTrackHypothesis(jhyp);
 791    if (!hypx) continue;
 792
 793    // Loop over all combinations of F2K fits and used OM hits
 794    for (Int_t k=0; k<fEvent.nfit_uses; k++)
 795    {
 796     if (fEvent.fit_uses[k].useid != hypx->GetId()) continue;
 797     hid=fEvent.fit_uses[k].hitid;
 798     sx=evt->GetIdHit(hid,"IceAOM");
 799     if (sx) sx->AddLink(hypx);
 800    }
 801   }
 802  }
 803 }
 804 ///////////////////////////////////////////////////////////////////////////