PHOS/AliPHOS.cxx

   1 ////////////////////////////////////////////////
   2 //  Manager and hits classes for set:PHOS     //
   3 ////////////////////////////////////////////////
   4
   5 // --- ROOT system ---
   6 #include "TH1.h"
   7 #include "TRandom.h"
   8 #include "TFile.h"
   9 #include "TTree.h"
  10 #include "TBRIK.h"
  11 #include "TNode.h"
  12 #include "TMath.h"
  13
  14 // --- Standard library ---
  15 #include <stdio.h>
  16 #include <string.h>
  17 #include <stdlib.h>
  18
  19 // --- galice header files ---
  20 #include "AliPHOS.h"
  21 #include "AliRun.h"
  22
  23 //______________________________________________________________________________
  24
  25
  26 ClassImp(AliPHOS)
  27
  28 //______________________________________________________________________________
  29
  30 AliPHOS::~AliPHOS(void)
  31 {
  32   delete fHits;                 // 28.12.1998
  33   delete fTreePHOS;             // 28.12.1998
  34   fCradles->Delete();
  35   delete fCradles;
  36 }
  37
  38 //______________________________________________________________________________
  39
  40 AliPHOS::AliPHOS() :
  41          fDebugLevel            (0),
  42          fTreePHOS              (NULL),
  43          fBranchNameOfCradles   ("AliPHOSCradles"),
  44          fTreeName              ("PHOS")
  45 {
  46    fIshunt   = 0;
  47
  48   if( NULL==(fCradles=new TObjArray) )
  49   {
  50     Error("AliPHOS","Can not create fCradles");
  51     exit(1);
  52   }
  53   DefPars();
  54 }
  55
  56 //______________________________________________________________________________
  57
  58 AliPHOS::AliPHOS(const char *name, const char *title)
  59        : AliDetector            (name,title),
  60          fDebugLevel            (0),
  61          fTreePHOS              (NULL),
  62          fBranchNameOfCradles   ("AliPHOSCradles"),
  63          fTreeName              ("PHOS")
  64 {
  65 //Begin_Html
  66 /*
  67 <img src="picts/aliphos.gif">
  68 */
  69 //End_Html
  70
  71    fHits   = new TClonesArray("AliPHOShit",  405);
  72
  73    fIshunt     =  0;
  74
  75    SetMarkerColor(kGreen);
  76    SetMarkerStyle(2);
  77    SetMarkerSize(0.4);
  78
  79   if( NULL==(fCradles=new TObjArray) ) {
  80      Error("AliPHOS","Can not create fCradles");
  81      exit(1);
  82   }
  83   DefPars();
  84 }
  85
  86 //______________________________________________________________________________
  87
  88 void AliPHOS::DefPars()
  89 {
  90       PHOSflags[0]=0;
  91       PHOSflags[1]=1;
  92       PHOSflags[2]=0;
  93       PHOSflags[3]=0;
  94       PHOSflags[4]=0;
  95       PHOSflags[5]=0;
  96       PHOSflags[6]=0;
  97       PHOSflags[7]=0;
  98       PHOSflags[8]=0;
  99       PHOScell[0]=2.2;
 100       PHOScell[1]=18.;
 101       PHOScell[2]=0.01;
 102       PHOScell[3]=0.01;
 103       PHOScell[4]=1.0;
 104       PHOScell[5]=0.1;
 105       PHOScell[6]=0.;
 106       PHOScell[7]=0.;
 107       PHOScell[8]=0.;
 108       PHOSradius=460.;
 109       PHOSsize[0]=104;
 110       PHOSsize[1]=88;
 111       PHOSsize[2]=4;
 112       PHOScradlesA=0.;
 113       PHOSextra[0]=0.001;
 114       PHOSextra[1]=6.95;
 115       PHOSextra[2]=4.;
 116       PHOSextra[3]=5.;
 117       PHOSextra[4]=2.;
 118       PHOSextra[5]=0.06;
 119       PHOSextra[6]=10.;
 120       PHOSextra[7]=3.;
 121       PHOSextra[8]=1.;
 122       PHOSTXW[0]=209.;
 123       PHOSTXW[1]=71.;
 124       PHOSTXW[2]=250.;
 125       PHOSAIR[0]=206.;
 126       PHOSAIR[1]=66.;
 127       PHOSAIR[2]=244.;
 128       PHOSFTI[0]=214.6;
 129       PHOSFTI[1]=80.;
 130       PHOSFTI[2]=260.;
 131       PHOSFTI[3]=467.;
 132 }
 133 //______________________________________________________________________________
 134
 135 void AliPHOS::AddHit(Int_t track, Int_t *vol, Float_t *hits)
 136 {
 137   TClonesArray &lhits = *fHits;
 138   new(lhits[fNhits++]) AliPHOShit(fIshunt,track,vol,hits);
 139 }
 140
 141 //___________________________________________
 142 void AliPHOS::BuildGeometry()
 143 {
 144
 145   TNode *Node, *Top;
 146
 147   const int kColorPHOS = kRed;
 148   //
 149   Top=gAlice->GetGeometry()->GetNode("alice");
 150
 151
 152   // PHOS
 153   Float_t pphi=12.9399462;
 154   new TRotMatrix("rot988","rot988",90,-3*pphi,90,90-3*pphi,0,0);
 155   new TRotMatrix("rot989","rot989",90,-  pphi,90,90-  pphi,0,0);
 156   new TRotMatrix("rot990","rot990",90,   pphi,90,90+  pphi,0,0);
 157   new TRotMatrix("rot991","rot991",90, 3*pphi,90,90+3*pphi,0,0);
 158   new TBRIK("S_PHOS","PHOS box","void",107.3,40,130);
 159   Top->cd();
 160   Node = new TNode("PHOS1","PHOS1","S_PHOS",-317.824921,-395.014343,0,"rot988");
 161   Node->SetLineColor(kColorPHOS);
 162   fNodes->Add(Node);
 163   Top->cd();
 164   Node = new TNode("PHOS2","PHOS2","S_PHOS",-113.532333,-494.124908,0,"rot989");
 165   fNodes->Add(Node);
 166   Node->SetLineColor(kColorPHOS);
 167   Top->cd();
 168   Node = new TNode("PHOS3","PHOS3","S_PHOS", 113.532333,-494.124908,0,"rot990");
 169   Node->SetLineColor(kColorPHOS);
 170   fNodes->Add(Node);
 171   Top->cd();
 172   Node = new TNode("PHOS4","PHOS4","S_PHOS", 317.824921,-395.014343,0,"rot991");
 173   Node->SetLineColor(kColorPHOS);
 174   fNodes->Add(Node);
 175 }
 176
 177 //___________________________________________
 178 void AliPHOS::CreateMaterials()
 179 {
 180 // *** DEFINITION OF AVAILABLE PHOS MATERIALS ***
 181
 182 // CALLED BY : PHOS_MEDIA
 183 // ORIGIN    : NICK VAN EIJNDHOVEN
 184
 185
 186
 187     Int_t   ISXFLD = gAlice->Field()->Integ();
 188     Float_t SXMGMX = gAlice->Field()->Max();
 189
 190 // --- The PbWO4 crystals ---
 191     Float_t ax[3] = { 207.19,183.85,16. };
 192     Float_t zx[3] = { 82.,74.,8. };
 193     Float_t wx[3] = { 1.,1.,4. };
 194     Float_t dx    = 8.28;
 195 // --- Stainless Steel ---
 196     Float_t as[5] = { 55.847,12.011,51.9961,58.69,28.0855 };
 197     Float_t zs[5] = { 26.,6.,24.,28.,14. };
 198     Float_t ws[5] = { .6392,8e-4,.2,.14,.02 };
 199     Float_t ds    = 8.;
 200 // --- The polysterene scintillator (CH) ---
 201     Float_t ap[2] = { 12.011,1.00794 };
 202     Float_t zp[2] = { 6.,1. };
 203     Float_t wp[2] = { 1.,1. };
 204     Float_t dp    = 1.032;
 205 // --- Tyvek (CnH2n)
 206     Float_t at[2] = { 12.011,1.00794 };
 207     Float_t zt[2] = { 6.,1. };
 208     Float_t wt[2] = { 1.,2. };
 209     Float_t dt    = .331;
 210 // --- Polystyrene foam ---
 211     Float_t af[2] = { 12.011,1.00794 };
 212     Float_t zf[2] = { 6.,1. };
 213     Float_t wf[2] = { 1.,1. };
 214     Float_t df    = .12;
 215 //--- Foam thermo insulation (actual chemical composition unknown yet!) ---
 216     Float_t ati[2] = { 12.011,1.00794 };
 217     Float_t zti[2] = { 6.,1. };
 218     Float_t wti[2] = { 1.,1. };
 219     Float_t dti    = .1;
 220 // --- Textolit (actual chemical composition unknown yet!) ---
 221     Float_t atx[2] = { 12.011,1.00794 };
 222     Float_t ztx[2] = { 6.,1. };
 223     Float_t wtx[2] = { 1.,1. };
 224     Float_t dtx    = 1.83;
 225
 226     Int_t *idtmed = fIdtmed->GetArray()-699;
 227
 228     AliMixture(  0, "PbWO4$",          ax, zx, dx, -3, wx);
 229     AliMixture(  1, "Polystyrene$",    ap, zp, dp, -2, wp);
 230     AliMaterial( 2, "Al$",             26.98, 13., 2.7, 8.9, 999);
 231 // ---                                Absorption length^ is ignored ---
 232     AliMixture(  3, "Tyvek$",           at, zt, dt, -2, wt);
 233     AliMixture(  4, "Foam$",            af, zf, df, -2, wf);
 234     AliMixture(  5, "Stainless Steel$", as, zs, ds, 5, ws);
 235     AliMaterial( 6, "Si$",              28.09, 14., 2.33, 9.36, 42.3);
 236     AliMixture(  7, "Thermo Insul.$",   ati, zti, dti, -2, wti);
 237     AliMixture(  8, "Textolit$",        atx, ztx, dtx, -2, wtx);
 238     AliMaterial(99, "Air$",             14.61, 7.3, .001205, 30420., 67500);
 239
 240     AliMedium(0, "PHOS Xtal    $", 0, 1, ISXFLD, SXMGMX, 10., .1, .1, .1, .1);
 241     AliMedium(2, "Al parts     $", 2, 0, ISXFLD, SXMGMX, 10., .1, .1, .001, .001);
 242     AliMedium(3, "Tyvek wrapper$", 3, 0, ISXFLD, SXMGMX, 10., .1, .1, .001, .001);
 243     AliMedium(4, "Polyst. foam $", 4, 0, ISXFLD, SXMGMX, 10., .1, .1, .1, .1);
 244     AliMedium(5, "Steel cover  $", 5, 0, ISXFLD, SXMGMX, 10., .1, .1, 1e-4, 1e-4);
 245     AliMedium(6, "Si PIN       $", 6, 0, ISXFLD, SXMGMX, 10., .1, .1, .01, .01);
 246     AliMedium(7, "Thermo Insul.$", 7, 0, ISXFLD, SXMGMX, 10., .1, .1, .1, .1);
 247     AliMedium(8, "Textolit     $", 8, 0, ISXFLD, SXMGMX, 10., .1, .1, .1, .1);
 248     AliMedium(99, "Air          $",99, 0, ISXFLD, SXMGMX, 10., 1., .1, .1, 10);
 249
 250 // --- Generate explicitly delta rays in the steel cover ---
 251     gMC->Gstpar(idtmed[704], "LOSS", 3.);
 252     gMC->Gstpar(idtmed[704], "DRAY", 1.);
 253 // --- and in aluminium parts ---
 254     gMC->Gstpar(idtmed[701], "LOSS", 3.);
 255     gMC->Gstpar(idtmed[701], "DRAY", 1.);
 256 }
 257
 258 //______________________________________________________________________________
 259
 260 void AliPHOS::AddPHOSCradles()
 261 {
 262   Int_t i;
 263   for(i=0;i<GetCradlesAmount();i++) {
 264
 265     int n = fCradles->GetEntries();
 266     fCradles->Add(new AliPHOSCradle( IsVersion(),            // geometry.
 267                                      GetCrystalSideSize    (),
 268                                      GetCrystalLength      (),
 269                                      GetWrapThickness      (),
 270                                      GetAirThickness       (),
 271                                      GetPIN_SideSize       (),
 272                                      GetPIN_Length         (),
 273                                      GetRadius             (),
 274                                      GetNz                 (),
 275                                      GetNphi               (),
 276                                      GetCradleAngle        (i)));
 277
 278     if( n+1 != fCradles->GetEntries() || NULL == fCradles->At(n) )
 279       {
 280         cout << "  Can not create or add AliPHOSCradle.\n";
 281         exit(1);
 282       }
 283   }
 284 }
 285
 286 //______________________________________________________________________________
 287
 288 Int_t AliPHOS::DistancetoPrimitive(Int_t , Int_t )
 289 {
 290    return 9999;
 291 }
 292
 293 //___________________________________________
 294 void AliPHOS::Init()
 295 {
 296   Int_t i;
 297   //
 298   printf("\n");
 299   for(i=0;i<35;i++) printf("*");
 300   printf(" PHOS_INIT ");
 301   for(i=0;i<35;i++) printf("*");
 302   printf("\n");
 303   //
 304   // Here the ABSO initialisation code (if any!)
 305   for(i=0;i<80;i++) printf("*");
 306   printf("\n");
 307 }
 308
 309 //______________________________________________________________________________
 310
 311 void AliPHOS::MakeBranch(Option_t *)
 312 {
 313 // ROOT output initialization to ROOT file.
 314 //
 315 // AliDetector::MakeBranch()  is always called.
 316 //
 317 // There will be also special tree "PHOS" with one branch "AliPHOSCradles"
 318 // if it was set next flag in the galice card file:
 319 //  * PHOSflags:    YES: X<>0   NO: X=0
 320 //  * PHOSflags(1) : -----X.  Create branch for TObjArray of AliPHOSCradle
 321 //     Examples:
 322 //     PHOSflags      1.
 323 //     PHOSflags 636301.
 324 // In that case special bit CradlesBranch_Bit will be set for AliPHOS
 325
 326   AliDetector::MakeBranch();
 327
 328   int i;
 329   float t = GetPHOS_flag(0)/10;
 330   i = (int) t;
 331   i = (int) ((t-i)*10);
 332   if( !i )
 333     return;
 334
 335   SetBit(CradlesBranch_Bit);
 336
 337   if( NULL==(fTreePHOS=new TTree(fTreeName.Data(),"PHOS events tree")) )
 338   {
 339     Error("MakeBranch","Can not create TTree");
 340     exit(1);
 341   }
 342
 343   if( NULL==fTreePHOS->GetCurrentFile() )
 344   {
 345     Error("MakeBranch","There is no opened ROOT file");
 346     exit(1);
 347   }
 348
 349   // Create a new branch in the current Root Tree.
 350
 351   if( NULL==fTreePHOS->Branch(fBranchNameOfCradles.Data(),"TObjArray",&fCradles,4000,0) )
 352   {
 353     Error("MakeBranch","Can not create branch");
 354     exit(1);
 355   }
 356
 357   printf("The branch %s has been created\n",fBranchNameOfCradles.Data());
 358 }
 359
 360 //______________________________________________________________________________
 361
 362 void AliPHOS::SetTreeAddress(void)
 363 {
 364 // ROOT input initialization.
 365 //
 366 // AliDetector::SetTreeAddress()  is always called.
 367 //
 368 // If CradlesBranch_Bit is set (see AliPHOS::MakeBranch) than fTreePHOS is
 369 // initilized.
 370
 371   AliDetector::SetTreeAddress();
 372
 373   if( !TestBit(CradlesBranch_Bit) )
 374     return;
 375
 376   if( NULL==(fTreePHOS=(TTree*)gDirectory->Get((char*)(fTreeName.Data()))  ) )
 377   {
 378     Error("SetTreeAddress","Can not find Tree \"%s\"\n",fTreeName.Data());
 379     exit(1);
 380   }
 381
 382   TBranch *branch = fTreePHOS->GetBranch(fBranchNameOfCradles.Data());
 383   if( NULL==branch )
 384   {
 385     Error("SetTreeAddress","Can not find branch %s in TTree:%s",fBranchNameOfCradles.Data(),fTreeName.Data());
 386     exit(1);
 387   }
 388
 389   branch->SetAddress(&fCradles);
 390 }
 391
 392 //______________________________________________________________________________
 393
 394 AliPHOSCradle *AliPHOS::GetCradleOfTheParticle(const TVector3 &p,const TVector3 &v) const
 395 {
 396 // For a given direction 'p' and source point 'v' returns pointer to AliPHOSCradle
 397 // in that direction or NULL if AliPHOSCradle was not found.
 398
 399   for( int m=0; m<fCradles->GetEntries(); m++ )
 400   {
 401     AliPHOS *PHOS = (AliPHOS *)this;     // Removing 'const'...
 402     AliPHOSCradle *cradle = (AliPHOSCradle *)PHOS->fCradles->operator[](m);
 403
 404     float x,y,l;
 405     const float d = cradle->GetRadius();
 406     cradle->GetXY(p,v,d,x,y,l);
 407
 408     if( l>0 && TMath::Abs(x)<cradle->GetNz  ()*cradle->GetCellSideSize()/2
 409             && TMath::Abs(y)<cradle->GetNphi()*cradle->GetCellSideSize()/2 )
 410       return cradle;
 411   }
 412
 413   return NULL;
 414 }
 415
 416 //______________________________________________________________________________
 417
 418 void AliPHOS::Reconstruction(Float_t signal_step, UInt_t min_signal_reject)
 419 {
 420 // Call AliPHOSCradle::Reconstruction(Float_t signal_step, UInt_t min_signal_reject)
 421 // for all AliPHOSCradles.
 422
 423   for( int i=0; i<fCradles->GetEntries(); i++ )
 424     GetCradle(i).Reconstruction(signal_step,min_signal_reject);
 425 }
 426
 427 //______________________________________________________________________________
 428
 429 void AliPHOS::ResetDigits(void)
 430 {
 431   AliDetector::ResetDigits();
 432
 433   for( int i=0; i<fCradles->GetEntries(); i++ )
 434     ((AliPHOSCradle*)(*fCradles)[i]) -> Clear();
 435 }
 436
 437 //______________________________________________________________________________
 438
 439 void AliPHOS::FinishEvent(void)
 440 {
 441 // Called at the end of each 'galice' event.
 442
 443   if( NULL!=fTreePHOS )
 444     fTreePHOS->Fill();
 445 }
 446
 447 //______________________________________________________________________________
 448
 449 void AliPHOS::FinishRun(void)
 450 {
 451 }
 452
 453 //______________________________________________________________________________
 454
 455 void AliPHOS::Print(Option_t *opt)
 456 {
 457 // Print PHOS information.
 458 // For each AliPHOSCradle the function AliPHOSCradle::Print(opt) is called.
 459
 460   AliPHOS &PHOS = *(AliPHOS *)this;     // Removing 'const'...
 461
 462   for( int i=0; i<fCradles->GetEntries(); i++ )
 463   {
 464     printf("PHOS cradle %d from %d\n",i+1, fCradles->GetEntries());
 465     PHOS.GetCradle(i).Print(opt);
 466     printf( "---------------------------------------------------\n");
 467   }
 468 }
 469
 470 //______________________________________________________________________________
 471 void AliPHOS::SetFlags(Float_t p1,Float_t p2,Float_t p3,Float_t p4,
 472                        Float_t p5,Float_t p6,Float_t p7,Float_t p8,Float_t p9)
 473 {
 474   PHOSflags[0]=p1;
 475   PHOSflags[1]=p2;
 476   PHOSflags[2]=p3;
 477   PHOSflags[3]=p4;
 478   PHOSflags[4]=p5;
 479   PHOSflags[5]=p6;
 480   PHOSflags[6]=p7;
 481   PHOSflags[7]=p8;
 482   PHOSflags[8]=p9;
 483 }
 484
 485 //______________________________________________________________________________
 486 void AliPHOS::SetCell(Float_t p1,Float_t p2,Float_t p3,Float_t p4,
 487                        Float_t p5,Float_t p6,Float_t p7,Float_t p8,Float_t p9)
 488 {
 489   PHOScell[0]=p1;
 490   PHOScell[1]=p2;
 491   PHOScell[2]=p3;
 492   PHOScell[3]=p4;
 493   PHOScell[4]=p5;
 494   PHOScell[5]=p6;
 495   PHOScell[6]=p7;
 496   PHOScell[7]=p8;
 497   PHOScell[8]=p9;
 498 }
 499
 500 //______________________________________________________________________________
 501 void AliPHOS::SetRadius(Float_t radius)
 502 {
 503    PHOSradius=radius;
 504 }
 505
 506 //______________________________________________________________________________
 507 void AliPHOS::SetCradleSize(Int_t nz, Int_t nphi, Int_t ncradles)
 508 {
 509    PHOSsize[0]=nz;
 510    PHOSsize[1]=nphi;
 511    PHOSsize[2]=ncradles;
 512 }
 513
 514 //______________________________________________________________________________
 515 void AliPHOS::SetCradleA(Float_t angle)
 516 {
 517    PHOScradlesA=angle;
 518 }
 519
 520 //______________________________________________________________________________
 521 void AliPHOS::SetExtra(Float_t p1,Float_t p2,Float_t p3,Float_t p4,
 522                        Float_t p5,Float_t p6,Float_t p7,Float_t p8,Float_t p9)
 523 {
 524    PHOSextra[0] = p1;
 525    PHOSextra[1] = p2;
 526    PHOSextra[2] = p3;
 527    PHOSextra[3] = p4;
 528    PHOSextra[4] = p5;
 529    PHOSextra[5] = p6;
 530    PHOSextra[6] = p7;
 531    PHOSextra[7] = p8;
 532    PHOSextra[8] = p9;
 533 }
 534
 535 //______________________________________________________________________________
 536 void AliPHOS::SetTextolitWall(Float_t dx, Float_t dy, Float_t dz)
 537 {
 538    PHOSTXW[0] = dx;
 539    PHOSTXW[1] = dy;
 540    PHOSTXW[2] = dz;
 541 }
 542
 543 //______________________________________________________________________________
 544 void AliPHOS::SetInnerAir(Float_t dx, Float_t dy, Float_t dz)
 545 {
 546    PHOSAIR[0] = dx;
 547    PHOSAIR[1] = dy;
 548    PHOSAIR[2] = dz;
 549 }
 550
 551 //______________________________________________________________________________
 552 void AliPHOS::SetFoam(Float_t dx, Float_t dy, Float_t dz, Float_t dr)
 553 {
 554    PHOSFTI[0] = dx;
 555    PHOSFTI[1] = dy;
 556    PHOSFTI[2] = dz;
 557    PHOSFTI[3] = dr;
 558 }
 559
 560 ClassImp(AliPHOSCradle)
 561
 562 //______________________________________________________________________________
 563
 564 AliPHOSCradle::AliPHOSCradle(void) {}
 565
 566 //______________________________________________________________________________
 567
 568 AliPHOSCradle::AliPHOSCradle( int   Geometry           ,
 569                               float CrystalSideSize    ,
 570                               float CrystalLength      ,
 571                               float WrapThickness      ,
 572                               float AirThickness       ,
 573                               float PIN_SideSize       ,
 574                               float PIN_Length         ,
 575                               float Radius             ,
 576                               int   Nz                 ,
 577                               int   Nphi               ,
 578                               float Angle              ) :
 579     fGeometry                   (Geometry),
 580 //  fCellEnergy                 (),
 581 //  fChargedTracksInPIN         (),
 582     fCrystalSideSize            (CrystalSideSize),
 583     fCrystalLength              (CrystalLength),
 584     fWrapThickness              (WrapThickness),
 585     fAirThickness               (AirThickness),
 586     fPIN_SideSize               (PIN_SideSize),
 587     fPIN_Length                 (PIN_Length),
 588     fRadius                     (Radius),
 589     fNz                         (Nz),
 590     fNphi                       (Nphi),
 591     fPhi                        (Angle)
 592 {
 593         fCellEnergy         = TH2F("CellE","Energy deposition in a cells",fNz,0,fNz,fNphi,0,fNphi);
 594         fCellEnergy           .SetDirectory(0);
 595         fChargedTracksInPIN = TH2S("PINCtracks","Amount of charged tracks in PIN",fNz,0,fNz,fNphi,0,fNphi);
 596         fChargedTracksInPIN   .SetDirectory(0);
 597 }
 598
 599 //______________________________________________________________________________
 600
 601 AliPHOSCradle::~AliPHOSCradle(void)        // 28.12.1998
 602 {
 603   fGammasReconstructed.Delete();
 604   fParticles          .Delete();
 605 }
 606
 607 //______________________________________________________________________________
 608
 609 void AliPHOSCradle::Clear(Option_t *)
 610 {
 611 // Clear digit. information.
 612
 613   fCellEnergy              .Reset();
 614   fChargedTracksInPIN      .Reset();
 615   GetParticles()           .Delete();
 616   GetParticles()           .Compress();
 617   GetGammasReconstructed() .Delete();
 618   GetGammasReconstructed() .Compress();
 619
 620 }
 621
 622 //______________________________________________________________________________
 623
 624 void AliPHOSCradle::GetXY(const TVector3 &p,const TVector3 &v,float R,float &x,float &y,float &l) const
 625 {
 626 // This function calculates hit position (x,y) in the CRADLE cells plain from particle in
 627 // the direction given by 'p' (not required to be normalized) and start point
 628 // given by 3-vector 'v'. So the particle trajectory is   t(l) = v + p*l
 629 // were 'l' is a number (distance from 'v' to CRADLE cells plain) and 't' is resulting
 630 // three-vector of trajectory point.
 631 //
 632 // After the call to this function user should test that l>=0 (the particle HITED the
 633 // plain) and (x,y) are in the region of CRADLE:
 634 //
 635 // Example:
 636 //   AliPHOSCradle cradle(......);
 637 //   TVector3 p(....), v(....);
 638 //   Float_t x,y,l;
 639 //   cradle.GetXY(p,v,x,y,l);
 640 //   if( l<0 || TMath::Abs(x)>cradle.GetNz()  *cradle.GetCellSideSize()/2
 641 //           || TMath::Abs(y)>cradle.GetNphi()*cradle.GetCellSideSize()/2 )
 642 //     cout << "Outside the CRADLE.\n";
 643
 644   // We have to create three vectors:
 645   //    s  - central point on the PHOS surface
 646   //    n1 - first vector in CRADLE plain
 647   //    n2 - second vector in CRADLE plain
 648   // This three vectors are orthonormalized.
 649
 650   double phi = fPhi/180*TMath::Pi();
 651   TVector3        n1(   0.0      ,   0.0      , 1.0 ),   // Z direction (X)
 652                   n2(  -sin(phi) ,   cos(phi) , 0 ),   // around beam (Y)
 653                   s ( R*cos(phi) , R*sin(phi) , 0 );   // central point
 654
 655   const double l1_min = 1e-2;
 656   double l1,
 657          p_n1 = p*n1,        // * - scalar product.
 658          p_n2 = p*n2,
 659          v_n1 = v*n1,
 660          v_n2 = v*n2,
 661          s_n1 = s*n1, // 0
 662          s_n2 = s*n2; // 0
 663
 664   if      ( TMath::Abs(l1=p.X()-n1.X()*p_n1-n2.X()*p_n2)>l1_min )
 665     { l = (-v.X()+s.X()+n1.X()*(v_n1-s_n1)+n2.X()*(v_n2-s_n2))/l1; }
 666   else if ( TMath::Abs(l1=p.Y()-n1.Y()*p_n1-n2.Y()*p_n2)>l1_min )
 667     { l = (-v.Y()+s.Y()+n1.Y()*(v_n1-s_n1)+n2.Y()*(v_n2-s_n2))/l1; }
 668   else if ( TMath::Abs(l1=p.Z()-n1.Z()*p_n1-n2.Z()*p_n2)>l1_min )
 669     { l = (-v.Z()+s.Z()+n1.Z()*(v_n1-s_n1)+n2.Z()*(v_n2-s_n2))/l1; }
 670
 671 //         double lx = (-v.X()+s.X()+n1.X()*(v.dot(n1)-s.dot(n1))+n2.X()*(v.dot(n2)-s.dot(n2)))/
 672 //                     (p.X()-n1.X()*p.dot(n1)-n2.X()*p.dot(n2)),
 673 //                ly = (-v.Y()+s.Y()+n1.Y()*(v.dot(n1)-s.dot(n1))+n2.Y()*(v.dot(n2)-s.dot(n2)))/
 674 //                     (p.Y()-n1.Y()*p.dot(n1)-n2.Y()*p.dot(n2)),
 675 //                lz = (-v.Z()+s.Z()+n1.Z()*(v.dot(n1)-s.dot(n1))+n2.Z()*(v.dot(n2)-s.dot(n2)))/
 676 //                     (p.Z()-n1.Z()*p.dot(n1)-n2.Z()*p.dot(n2));
 677 //         cout.form("x: %g %g %g %g\n",lx,-v.X()+s.X()+n1.X()*(v.dot(n1)-s.dot(n1))+n2.X()*(v.dot(n2)-s.dot(n2)),p.X()-n1.X()*p.dot(n1)-n2.X()*p.dot(n2));
 678 //         cout.form("y: %g %g %g %g\n",lx,-v.Y()+s.Y()+n1.Y()*(v.dot(n1)-s.dot(n1))+n2.Y()*(v.dot(n2)-s.dot(n2)),p.Y()-n1.Y()*p.dot(n1)-n2.Y()*p.dot(n2));
 679 //         cout.form("z: %g %g %g %g\n",lx,-v.Z()+s.Z()+n1.Z()*(v.dot(n1)-s.dot(n1))+n2.Z()*(v.dot(n2)-s.dot(n2)),p.Z()-n1.Z()*p.dot(n1)-n2.Z()*p.dot(n2));
 680 //         cout.form("lx,ly,lz =   %g,%g,%g\n",lx,ly,lz);
 681
 682   x = p_n1*l + v_n1 - s_n1;
 683   y = p_n2*l + v_n2 - s_n2;
 684 }
 685
 686 //______________________________________________________________________________
 687
 688 void AliPHOSCradle::Print(Option_t *opt)
 689 {
 690 // Print AliPHOSCradle information.
 691 //
 692 // options:  'd' - print energy deposition for EVERY cell
 693 //           'p' - print particles list that hit the cradle
 694 //           'r' - print list of reconstructed particles
 695
 696   AliPHOSCradle *cr = (AliPHOSCradle *)this;     // Removing 'const'...
 697
 698   printf("AliPHOSCradle:  Nz=%d  Nphi=%d, fPhi=%f, E=%g\n",fNz,fNphi,fPhi,
 699        cr->fCellEnergy.GetSumOfWeights());
 700
 701   if( NULL!=strchr(opt,'d') )
 702   {
 703     printf("\n\nCells Energy (in MeV):\n\n   |");
 704     for( int x=0; x<fNz; x++ )
 705       printf(" %4d|",x+1);
 706     printf("\n");
 707
 708     for( int y=fNphi-1; y>=0; y-- )
 709     {
 710       printf("%3d|",y+1);
 711       for( int x=0; x<fNz; x++ )
 712         printf("%6d",(int)(cr->fCellEnergy.GetBinContent(cr->fCellEnergy.GetBin(x,y))*1000));
 713       printf("\n");
 714     }
 715     printf("\n");
 716   }
 717
 718   if( NULL!=strchr(opt,'p') )
 719   {
 720     printf("This cradle was hit by %d particles\n",
 721          ((AliPHOSCradle*)this)->GetParticles().GetEntries());
 722     TObjArray &p=((AliPHOSCradle*)this)->GetParticles();
 723     for( int i=0; i<p.GetEntries(); i++ )
 724       ((AliPHOSgamma*)(p[i]))->Print();
 725   }
 726
 727   if( NULL!=strchr(opt,'p') )
 728   {
 729     printf("Amount of reconstructed gammas is %d\n",
 730          ((AliPHOSCradle*)this)->GetGammasReconstructed().GetEntries());
 731
 732     TObjArray &p=((AliPHOSCradle*)this)->GetGammasReconstructed();
 733     for( int i=0; i<p.GetEntries(); i++ )
 734       ((AliPHOSgamma*)(p[i]))->Print();
 735   }
 736 }
 737
 738 //______________________________________________________________________________
 739
 740 void AliPHOSCradle::Distortion(const TH2F *Noise, const TH2F *Stochastic, const TH2F *Calibration)
 741 {
 742 // This function changes histogram of cell energies fCellEnergy on the base of input
 743 // histograms Noise, Stochastic, Calibration. The histograms must have
 744 // size Nz x Nphi.
 745
 746   //////////////////////////////////
 747   // Testing the histograms size. //
 748   //////////////////////////////////
 749
 750   if( fNz!=fCellEnergy.GetNbinsX() || fNphi!=fCellEnergy.GetNbinsY() )
 751   {
 752     printf      ("Bad size of CellEnergy!   Must be:   Nz x Nphi = %d x %d\n"
 753                  "but size of CellEnergy is:  %d x %d\n",
 754                  fNz,fNphi,fCellEnergy.GetNbinsX(),fCellEnergy.GetNbinsY());
 755     exit(1);
 756   }
 757
 758   if( fNz!=fChargedTracksInPIN.GetNbinsX() || fNphi!=fChargedTracksInPIN.GetNbinsY() )
 759   {
 760     printf      ("Bad size of ChargedTracksInPIN!   Must be:   Nz x Nphi = %d x %d\n"
 761                  "but size of ChargedTracksInPIN is:  %d x %d\n",
 762                  fNz,fNphi,fChargedTracksInPIN.GetNbinsX(),fChargedTracksInPIN.GetNbinsY());
 763     exit(1);
 764   }
 765
 766   if( NULL!=Noise && (fNz!=Noise->GetNbinsX() || fNphi!=Noise->GetNbinsX()) )
 767   {
 768     printf      ("Bad size of Noise!   Must be:   Nz x Nphi = %d x %d\n"
 769                  "but size of Noise is:  %d x %d\n",
 770                  fNz,fNphi,fChargedTracksInPIN.GetNbinsX(),fChargedTracksInPIN.GetNbinsY());
 771     exit(1);
 772   }
 773
 774   if( NULL!=Stochastic && (fNz!=Stochastic->GetNbinsX() || fNphi!=Stochastic->GetNbinsX()) )
 775   {
 776     printf      ("Bad size of Stochastic!   Must be:   Nz x Nphi = %d x %d\n"
 777                  "but size of Stochastic is:  %d x %d\n",
 778                  fNz,fNphi,fChargedTracksInPIN.GetNbinsX(),fChargedTracksInPIN.GetNbinsY());
 779     exit(1);
 780   }
 781
 782   if( NULL!=Calibration && (fNz!=Calibration->GetNbinsX() || fNphi!=Calibration->GetNbinsX()) )
 783   {
 784     printf      ("Bad size of Calibration!   Must be:   Nz x Nphi = %d x %d\n"
 785                  "but size of Calibration is:  %d x %d\n",
 786                  fNz,fNphi,fChargedTracksInPIN.GetNbinsX(),fChargedTracksInPIN.GetNbinsY());
 787     exit(1);
 788   }
 789
 790   ////////////////////
 791   // Do distortion! //
 792   ////////////////////
 793
 794   for( int y=0; y<fNphi; y++ )
 795     for( int x=0; x<fNz; x++ )
 796     {
 797       const int n = fCellEnergy.GetBin(x,y);   // Bin number
 798       static TRandom r;
 799
 800       Float_t   E_old=fCellEnergy.GetBinContent(n),   E_new=E_old;
 801
 802       if( NULL!=Stochastic )
 803         E_new   = r.Gaus(E_old,sqrt(E_old)*GetDistortedValue(Stochastic,n));
 804
 805       if( NULL!=Calibration )
 806         E_new  *=  GetDistortedValue(Calibration,n);
 807
 808       if( NULL!=Noise )
 809         E_new  +=  GetDistortedValue(Noise,n);
 810
 811       fCellEnergy.SetBinContent(n,E_new);
 812     }
 813 }
 814
 815 ////////////////////////////////////////////////////////////////////////////////
 816
 817 TH2F* AliPHOSCradle::CreateHistForDistortion(const char *name, const char *title,
 818                                              Int_t Nx, Int_t Ny,
 819                                              Float_t MU_mu,    Float_t MU_sigma,
 820                                              Float_t SIGMA_mu, Float_t SIGMA_sigma)
 821 {
 822 // Create (new TH2F(...)) histogram with information (for every bin) that will
 823 // be used for VALUE creation.
 824 // Two values will be created for each bin:
 825 // MU    = TRandom::Gaus(MU_mu,MU_sigma)
 826 // and
 827 // SIGMA = TRandom::Gaus(SIGMA_mu,SIGMA_sigma)
 828 // The VALUE in a particluar bin will be equal
 829 // VALUE = TRandom::Gaus(MU,SIGMA)
 830 //
 831 // Do not forget to delete the histogram at the end of the work.
 832
 833   TH2F *h = new TH2F( name,title, Nx,1,Nx, Ny,1,Ny );
 834   if( h==NULL )
 835   {
 836     Error("CreateHistForDistortion","Can not create the histogram");
 837     exit(1);
 838   }
 839   h->SetDirectory(0);
 840
 841   for( int y=0; y<Ny; y++ )
 842     for( int x=0; x<Nx; x++ )
 843     {
 844       const int n = h->GetBin(x,y);
 845       h->SetBinContent(n,r.Gaus(   MU_mu,   MU_sigma));
 846       h->SetBinError  (n,r.Gaus(SIGMA_mu,SIGMA_sigma));
 847     }
 848
 849   return h;
 850 }
 851
 852 ////////////////////////////////////////////////////////////////////////////////
 853
 854 Float_t AliPHOSCradle::GetDistortedValue(const TH2F *h, UInt_t n)
 855 {
 856   return r.Gaus(((TH2F*)h)->GetBinContent(n),n);
 857 }
 858
 859 ////////////////////////////////////////////////////////////////////////////////
 860 //______________________________________________________________________________
 861
 862 #ifdef WIN32
 863   #define common_for_event_storing COMMON_FOR_EVENT_STORING
 864 #else
 865   #define common_for_event_storing common_for_event_storing_
 866 #endif
 867
 868 extern "C" struct
 869 {
 870   enum { crystals_matrix_amount_max=4, crystals_in_matrix_amount_max=40000 };
 871
 872   // Event-independent information
 873   UShort_t      crystals_matrix_amount_PHOS,
 874                 crystal_matrix_type,
 875                 amount_of_crystals_on_Z,
 876                 amount_of_crystals_on_PHI;
 877   Float_t       radius,
 878                 crystal_size,
 879                 crystal_length,
 880                 matrix_coordinate_Z             [crystals_matrix_amount_max],
 881                 matrix_coordinate_PHI           [crystals_matrix_amount_max];
 882   UInt_t        event_number;
 883   UShort_t      crystals_amount_with_amplitudes [crystals_matrix_amount_max],
 884                 crystals_amplitudes_Iad         [crystals_matrix_amount_max]
 885                                                 [crystals_in_matrix_amount_max][2];
 886 } common_for_event_storing;
 887
 888 //       integer*4 crystals_amount_max,crystals_in_matrix_amount_max,
 889 //      +          crystals_matrix_amount_max
 890 //       parameter (crystals_matrix_amount_max=4)
 891 //       parameter (crystals_in_matrix_amount_max=40000)
 892 //       parameter (crystals_amount_max =crystals_matrix_amount_max*
 893 //      +                                crystals_in_matrix_amount_max)
 894 //
 895 // * All units are in GeV, cm, radian
 896 //       real       crystal_amplitudes_unit, radius_unit,
 897 //      +           crystal_size_unit, crystal_length_unit,
 898 //      +           matrix_coordinate_Z_unit, matrix_coordinate_PHI_unit
 899 //       integer    crystal_amplitudes_in_units_min
 900 //       parameter (crystal_amplitudes_in_units_min        = 1)
 901 //       parameter (crystal_amplitudes_unit                = 0.001 ) ! 1.0  MeV
 902 //       parameter (radius_unit                            = 0.1   ) ! 0.1  cm
 903 //       parameter (crystal_size_unit                      = 0.01  ) ! 0.01 cm
 904 //       parameter (crystal_length_unit                    = 0.01  ) ! 0.01 cm
 905 //       parameter (matrix_coordinate_Z_unit               = 0.1   ) ! 0.1  cm
 906 //       parameter (matrix_coordinate_PHI_unit             = 1e-4  ) ! 1e-4 radian
 907 //
 908 //       integer*2 crystals_matrix_amount_PHOS, crystal_matrix_type,
 909 //      +          amount_of_crystals_on_Z, amount_of_crystals_on_PHI,
 910 //      +          crystals_amount_with_amplitudes, crystals_amplitudes_Iad
 911 //       integer*4 event_number
 912 //
 913 //       real      radius, crystal_size, crystal_length,
 914 //      +          matrix_coordinate_Z, matrix_coordinate_PHI
 915 //
 916 //       real      crystals_amplitudes, crystals_energy_total
 917 //       integer   event_file_unit_number
 918 //
 919 //       common /common_for_event_storing/
 920 //      + ! Event-independent information
 921 //      +        crystals_matrix_amount_PHOS,
 922 //      +        crystal_matrix_type,
 923 //      +        amount_of_crystals_on_Z,
 924 //      +        amount_of_crystals_on_PHI,
 925 //      +        radius,
 926 //      +        crystal_size,
 927 //      +        crystal_length,
 928 //      +        matrix_coordinate_Z     (crystals_matrix_amount_max),
 929 //      +        matrix_coordinate_PHI   (crystals_matrix_amount_max),
 930 //      +
 931 //      + ! Event-dependent information
 932 //      +        event_number,
 933 //      +        crystals_amount_with_amplitudes
 934 //      +                                (crystals_matrix_amount_max),
 935 //      +        crystals_amplitudes_Iad (2,crystals_in_matrix_amount_max,
 936 //      +                                 crystals_matrix_amount_max),
 937 //      +
 938 //      + ! These information don't store in data file
 939 //      +        crystals_amplitudes     (crystals_amount_max),
 940 //      +        crystals_energy_total,
 941 //      +        event_file_unit_number
 942
 943
 944 //      parameter (NGp=1000,nsps=10,nvertmax=1000)
 945 //         COMMON /GAMMA/KG,MW(ngp),ID(ngp),JD(ngp),E(ngp),E4(ngp),
 946 //      ,  XW(ngp),YW(ngp),ES(nsps,ngp),ET(nsps,ngp),ISsd(ngp),
 947 //      ,  IGDEV(ngp),ZGDEV(ngp),sigexy(3,ngp),Emimx(2,nsps,ngp),
 948 //      ,  kgfix,igfix(ngp),cgfix(3,ngp),sgfix(3,ngp),hiw(ngp),
 949 //      ,  wsw(nsps,ngp),h1w(ngp),h0w(ngp),raxay(5,ngp),
 950 //      ,  sigmaes0(nsps,ngp),dispeces(nsps,ngp),
 951 //      ,  igamvert(ngp)
 952
 953
 954 #ifdef WIN32
 955 #define rcgamma RCGAMMA
 956 #else
 957 #define rcgamma rcgamma_
 958 #endif
 959
 960 extern "C" struct
 961 {
 962   enum {NGP=1000, nsps=10, nvertmax=1000};
 963   int   recons_gammas_amount, mw[NGP],ID[NGP],JD[NGP];
 964   float E[NGP], E4[NGP], XW[NGP], YW[NGP], ES[NGP][nsps],ET[NGP][nsps],ISsd[NGP],
 965         igdev[NGP],Zgdev[NGP];
 966 //      sigexy(3,ngp),Emimx(2,nsps,ngp),
 967 //   ,  kgfix,igfix(ngp),cgfix(3,ngp),sgfix(3,ngp),hiw(ngp),
 968 //   ,  wsw(nsps,ngp),h1w(ngp),h0w(ngp),raxay(5,ngp),
 969 //   ,  sigmaes0(nsps,ngp),dispeces(nsps,ngp),
 970 //   ,  igamvert(ngp)
 971 } rcgamma;
 972
 973 #ifdef WIN32
 974 #define reconsfirst RECONSFIRST
 975 #define type_of_call _stdcall
 976 #else
 977 #define reconsfirst reconsfirst_
 978 #define type_of_call
 979 #endif
 980
 981 extern "C" void type_of_call reconsfirst(const float &,const float &);
 982
 983 void AliPHOSCradle::Reconstruction(Float_t signal_step, UInt_t min_signal_reject)
 984 {
 985 // Call of PHOS reconstruction program.
 986 // signal_step=0.001  GeV (1MeV)
 987 // min_signal_reject = 15 or 30 MeV
 988
 989
 990   common_for_event_storing.event_number                       = 0;  // We do not know event number?
 991   common_for_event_storing.crystals_matrix_amount_PHOS        = 1;
 992   common_for_event_storing.crystal_matrix_type                = 1; // 1 - rectangular
 993   common_for_event_storing.amount_of_crystals_on_Z            = fNz;
 994   common_for_event_storing.amount_of_crystals_on_PHI          = fNphi;
 995
 996   common_for_event_storing.radius                             = fRadius;
 997   common_for_event_storing.crystal_size                       = GetCellSideSize();
 998   common_for_event_storing.crystal_length                     = fCrystalLength;
 999
1000   common_for_event_storing.matrix_coordinate_Z            [0] = 0;
1001   common_for_event_storing.matrix_coordinate_PHI          [0] = fPhi;
1002
1003   #define  k    common_for_event_storing.crystals_amount_with_amplitudes[0]
1004   k=0;
1005
1006   for( int y=0; y<fNphi; y++ )
1007     for( int x=0; x<fNz; x++ )
1008     {
1009       UInt_t    n       = fCellEnergy.GetBin(x,y);
1010       UInt_t    signal  = (int) (fCellEnergy.GetBinContent(n)/signal_step);
1011       if( signal>=min_signal_reject )
1012       {
1013         common_for_event_storing.crystals_amplitudes_Iad[0][k][0] = signal;
1014         common_for_event_storing.crystals_amplitudes_Iad[0][k][1] = x + y*fNz;
1015         k++;
1016       }
1017     }
1018   #undef  k
1019
1020   GetGammasReconstructed().Delete();
1021   GetGammasReconstructed().Compress();
1022
1023   const float   stochastic_term   = 0.03,        // per cents over sqrt(E);  E in GeV
1024                 electronic_noise  = 0.01;        // GeV
1025   reconsfirst(stochastic_term,electronic_noise); // Call of reconstruction program.
1026
1027   for( int i=0; i<rcgamma.recons_gammas_amount; i++ )
1028   {
1029 //     new (GetGammasReconstructed().UncheckedAt(i) ) AliPHOSgamma;
1030 //     AliPHOSgamma &g = *(AliPHOSgamma*)(GetGammasReconstructed().UncheckedAt(i));
1031
1032     AliPHOSgamma *gggg = new AliPHOSgamma;
1033     if( NULL==gggg )
1034     {
1035       Error("Reconstruction","Can not create AliPHOSgamma");
1036       exit(1);
1037     }
1038
1039     GetGammasReconstructed().Add(gggg);
1040     AliPHOSgamma &g=*gggg;
1041
1042     Float_t thetta, alpha, betta, R=fRadius+rcgamma.Zgdev[i]/10;
1043
1044     g.fX      = rcgamma.YW[i]/10;
1045     g.fY      = rcgamma.XW[i]/10;
1046     g.fE      = rcgamma.E [i];
1047
1048     thetta      = atan(g.fX/R);
1049
1050     alpha = atan(g.fY/R);
1051     betta = fPhi/180*TMath::Pi() + alpha;
1052
1053     g.fPx = g.fE * cos(thetta) * cos(betta);
1054     g.fPy = g.fE * cos(thetta) * sin(betta);
1055     g.fPz = g.fE * sin(thetta);
1056   }
1057 }
1058
1059 //______________________________________________________________________________
1060 //______________________________________________________________________________
1061 //______________________________________________________________________________
1062 //______________________________________________________________________________
1063 //______________________________________________________________________________
1064
1065 ClassImp(AliPHOSgamma)
1066
1067 //______________________________________________________________________________
1068
1069 void AliPHOSgamma::Print(Option_t *)
1070 {
1071   float mass = fE*fE - fPx*fPx - fPy*fPy - fPz*fPz;
1072
1073   if( mass>=0 )
1074     mass =  sqrt( mass);
1075   else
1076     mass = -sqrt(-mass);
1077
1078   printf("XY=(%+7.2f,%+7.2f)  (%+7.2f,%+7.2f,%+7.2f;%7.2f)  mass=%8.4f  Ipart=%2d\n",
1079           fX,fY,fPx,fPy,fPz,fE,mass,fIpart);
1080 }
1081
1082 //______________________________________________________________________________
1083
1084 AliPHOSgamma &AliPHOSgamma::operator=(const AliPHOSgamma &g)
1085 {
1086   fX           = g.fX;
1087   fY           = g.fY;
1088   fE           = g.fE;
1089   fPx          = g.fPx;
1090   fPy          = g.fPy;
1091   fPz          = g.fPz;
1092   fIpart       = g.fIpart;
1093
1094   return *this;
1095 }
1096
1097 //______________________________________________________________________________
1098 //______________________________________________________________________________
1099 //______________________________________________________________________________
1100 //______________________________________________________________________________
1101 //______________________________________________________________________________
1102
1103 ClassImp(AliPHOShit)
1104
1105 //______________________________________________________________________________
1106
1107 AliPHOShit::AliPHOShit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits):
1108 AliHit(shunt, track)
1109 {
1110    Int_t i;
1111    for (i=0;i<5;i++) fVolume[i] = vol[i];
1112    fX       = hits[0];
1113    fY       = hits[1];
1114    fZ       = hits[2];
1115    fELOS    = hits[3];
1116 }
1117
1118 //______________________________________________________________________________