//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////
//                                                          //
//  This program has been generated automatically by the    //
//  Root utility gh2root from a Geant RZ geometry file.     //
//                                                          //
//  The class XRun is the main processor.                  //
//  It contains member functions to:                        //
//    - Initialize Root and all detector makers.            //
//    - Control the processing of one event.                //
//    - Simulate random tracks and hits.                    //
//    - Import a Zebra KINE,HITS,DIGI structures from FZ.   //
//    - Create and Fill the Root Trees with particles/hits. //
//    - Browse the generated Root file.                     //
//    - Control functions for graphics:                     //
//        . DistancetoPrimitive                             //
//        . Paint                                           //
//                                                          //
//  Some Root macros are also generated by gh2root:         //
//    - fake.C    : to generate random particles and hits.  //
//    - run.C     : to convert a Geant FZ file into Root.   //
//    - browse.C  : to browse a Root file with hits/digits. //
//    - analyze.C : example of analysis code to loop on     //
//                  generated events and make histograms.   //
//                                                          //
//         Author:  Rene Brun                               //
//////////////////////////////////////////////////////////////
 
#include <TRandom.h>
#include "GParticle.h"
#include "StarRun.h"
#include "StarMaker.h"
#include "StarSVTH.h"
#include "StarTPCH.h"
#include "StarFTPH.h"
#include "StarBTOH.h"
#include "StarCALH.h"
 
 
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
 
//--------------Declarations for ZEBRA---------------------
#ifndef WIN32
#  define quest quest_
#  define gclink gclink_
int quest[100];
int gclink[20];
#else
#  define quest  QUEST
#  define gclink GCLINK
extern "C" int quest[100];
extern "C" int gclink[20];
#endif
 
int *iq, *lq;
float *q;
int *iquest;
//  Define the names of the Fortran subroutine and functions for the different OSs
 
#ifndef WIN32
# define mzebra  mzebra_
# define mzstor  mzstor_
# define mzlink  mzlink_
# define mzwipe  mzwipe_
# define fzin    fzin_
# define fzfile  fzfile_
# define dzshow  dzshow_
# define cfopen  cfopen_
 
# define type_of_call
# define DEFCHAR  const char*
# define PASSCHAR(string) string
#else
# define mzebra  MZEBRA
# define mzstor  MZSTOR
# define mzlink  MZLINK
# define mzwipe  MZWIPE
# define fzin    FZIN
# define fzfile  FZFILE
# define dzshow  DZSHOW
# define cfopen  CFOPEN
#endif
 
extern "C" void  type_of_call mzebra(const int&);
extern "C" void  type_of_call mzstor(const int&,DEFCHAR,DEFCHAR, int*, int*, int*, int*, int*, int*
#ifndef WIN32
                        ,const int, const int
#endif
                        );
extern "C" void  type_of_call mzlink(const int&,DEFCHAR, int*, int*, int*
#ifndef WIN32
                        ,const int
#endif
                        );
extern "C" void  type_of_call mzwipe(const int&);
extern "C" void  type_of_call dzshow(DEFCHAR,const int&,const int&,DEFCHAR,const int&, const int&, const int&, const int&
#ifndef WIN32
           ,const int,const int
#endif
                        );
extern "C" void  type_of_call fzin(const int&,const int&, int*,const int&,
                                   DEFCHAR,const int&, int*
#ifndef WIN32
                        ,const int
#endif
                        );
extern "C" void  type_of_call cfopen(const int&,const int&,const int&,
                                   DEFCHAR,const int&,DEFCHAR,const int&
#ifndef WIN32
                        ,const int,const int
#endif
                        );
extern "C" void  type_of_call fzfile(const int&,const int&,DEFCHAR
#ifndef WIN32
                        ,const int
#endif
                        );
 
int lunfz;
//      COMMON/GCLINK/JDIGI ,JDRAW ,JHEAD ,JHITS ,JKINE ,JMATE ,JPART
//     +      ,JROTM ,JRUNG ,JSET  ,JSTAK ,JGSTAT,JTMED ,JTRACK,JVERTX
//     +      ,JVOLUM,JXYZ  ,JGPAR ,JGPAR2,JSKLT
//*
   int *ldigi, *lhead, *lhits, *lkine, *lvertx, *lpart, *lsklt;
//--------------End of Declarations for ZEBRA------------------
 
StarRun *Run;
 
ClassImp(StarRun)
 
//___________________________________________
 StarRun::StarRun()
{
   fMakers    = 0;
   fTreeK     = 0;
   fTreeH     = 0;
   fTreeD     = 0;
   fParticles = 0;
}
 
//___________________________________________
 StarRun::StarRun(const char *name, const char *title)
        : TNamed(name,title)
{
//  Constructor for the main processor.
//  Creates the list of Makers.
//  Creates the list of particles.
 
   Run    = this;
   fTreeK = 0;
   fTreeH = 0;
   fTreeD = 0;
 
   gROOT->GetListOfBrowsables()->Add(this,name);
 
// create the support list for the various Makers
   fMakers = new TList();
 
// create "standard" makers and add them to the list of makers (in XMaker constructor
// Note that the order in which makers are added to the list of makers is important
// makers will be processed in this order !!
 
   new StarSVTHSVTDMaker("SVTHSVTDMaker");
   new StarSVTHSFSDMaker("SVTHSFSDMaker");
   new StarTPCHTPADMaker("TPCHTPADMaker");
   new StarTPCHTPAIMaker("TPCHTPAIMaker");
   new StarTPCHTMSEMaker("TPCHTMSEMaker");
   new StarFTPHFSECMaker("FTPHFSECMaker");
   new StarBTOHBCSBMaker("BTOHBCSBMaker");
   new StarCALHCSUPMaker("CALHCSUPMaker");
   new StarCALHCSDAMaker("CALHCSDAMaker");
 
// create array for particles (conversion from JVERTX,JKINE to TClonesArray)
   fParticles = new TClonesArray("GParticle",100);
}
 
//___________________________________________
 void StarRun::Browse(TBrowser *b)
{
// Called when the item "Run" is clicked on the left pane
// of the Root browser.
// It displays the Root Trees and all detector makers.
   if (fTreeK) b->Add(fTreeK,fTreeK->GetName());
   if (fTreeH) b->Add(fTreeH,fTreeH->GetName());
   if (fTreeD) b->Add(fTreeD,fTreeD->GetName());
 
   TIter next(fMakers);
   StarMaker *maker;
   while((maker = (StarMaker*)next())) {
      b->Add(maker,maker->GetName());
   }
}
 
//___________________________________________
 void StarRun::Clear(Option_t *option)
{
//  Reset all Makers. Called at the end of each event
   TIter next(fMakers);
   StarMaker *maker;
   while((maker = (StarMaker*)next())) {
      maker->Clear(option);
   }
 
   fParticles->Clear();
}
 
//___________________________________________
 Int_t StarRun::DistancetoPrimitive(Int_t, Int_t)
{
   return 9999;
}
 
//___________________________________________
 void StarRun::FakeDigits(Int_t nd, Int_t snd)
{
// Generates random digits.
// Loop on all detector makers
 
   TIter next(fMakers);
   StarMaker *maker;
   while((maker = (StarMaker*)next())) {
      maker->FakeDigits(nd, snd);
   }
}
 
//___________________________________________
 void StarRun::FakeHits(Int_t nh, Int_t snh)
{
// Generates random hits.
// Loop on all detector makers
 
   TIter next(fMakers);
   StarMaker *maker;
   while((maker = (StarMaker*)next())) {
      maker->FakeHits(nh, snh);
   }
}
 
//___________________________________________
 void StarRun::FakeKine(Int_t nt, Int_t snt)
{
// Generates random particles following a gaussian
// with mean value nt and standard deviation nh.
// Generated particles are added to the list of particles.
   Int_t it,ipart,ks,kfcode,parent;
   Int_t firstchild,lastchild;
   Float_t xmass, tlife, time;
   Float_t E,px,py,pz,vx,vy,vz;
 
   TClonesArray &particles = *fParticles;
   Int_t ntrack = Int_t(gRandom->Gaus(nt, snt));
 
   for (it=1;it<=ntrack;it++) {
      ipart  = Int_t(gRandom->Rndm());
      xmass  = 10*gRandom->Rndm();
      tlife  = 1.e-9*gRandom->Rndm();
      px     = gRandom->Gaus(100,200);
      py     = gRandom->Gaus(100,200);
      pz     = gRandom->Gaus(500,1000);
      E      = TMath::Sqrt(xmass*xmass+px*px+py*py+pz*pz);
      ks     = 0;
      kfcode = ipart;
      parent = 0;
      time   = 0;
      vx     = gRandom->Gaus(0,2);
      vy     = gRandom->Gaus(0,2);
      vz     = gRandom->Gaus(0,100);
      firstchild = 0;
      lastchild  = 0;
 
      new(particles[it-1]) GParticle(ks,kfcode,parent,firstchild,lastchild,
         px,py,pz,E,xmass,vx,vy,vz,time,tlife);
   }
}
 
//___________________________________________
 void StarRun::FillTree()
{
   if (fTreeK) fTreeK->Fill();
   if (fTreeH) fTreeH->Fill();
   if (fTreeD) fTreeD->Fill();
}
 
//___________________________________________
 void StarRun::Finish()
{
// Called at the end of the run.
   TIter next(fMakers);
   StarMaker *maker;
   while((maker = (StarMaker*)next())) {
      maker->Finish();
   }
}
 
//___________________________________________
 Int_t StarRun::GetEvent(const Int_t event)
{
// Returns event number event from the Root Trees.
   fEvent = event;
   Int_t nbytes = 0;
   if (fTreeK) nbytes += fTreeK->GetEvent(event);
   if (fTreeH) nbytes += fTreeH->GetEvent(event);
   if (fTreeD) nbytes += fTreeD->GetEvent(event);
   return nbytes;
}
 
//___________________________________________
 Int_t StarRun::GetGeantEvent(Option_t *option)
{
// Read the next event from the Geant FZ file.
// Search for the next HEAD record.
// If found, read VERT,KINE,HITS and DIGI records
//    returns 0 if nothing read or EOF reached.
//    returns 1 otherwise.
   const int maxhead = 100;
   int iuhead[maxhead];
   char *H = strstr(option,"H");
   char *D = strstr(option,"D");
   int nuhead;
 
   if (fDebug > 0) printf("GetGeantEvent startingn");
   mzwipe(0);
 
//       Get next event header
   while(1) {
      nuhead = maxhead;
      fzin(lunfz,0,0,0,"S",nuhead,iuhead,1);
      if (fDebug > 2) printf("iquest(1)=%d,nuhead=%d,iuhead(1)=%d,iuhead(2)=%d,iuhead(3)=%4sn",iquest[0],nuhead,iuhead[0],iuhead[1],(char*)&iuhead[2]);
      if (iquest[0] >= 2) return 0;
      if (nuhead >= 3) {
         if (strncmp((char*)&iuhead[2],"HEAD",4) == 0) break;
      }
   }
 
//        header found, read all events data structures
   fzin(lunfz,0,lhead,1,"A",nuhead,iuhead,1);
   if (fDebug > 1) printf("JHEAD structure read OKn");
   while(1) {
      nuhead = maxhead;
      fzin(lunfz,0,0,0,"S",nuhead,iuhead,1);
      if (iquest[0] >= 2) return 0;
      if (fDebug > 2) printf("iquest(1)=%d,nuhead=%d,iuhead(1)=%d,iuhead(2)=%d,iuhead(3)=%4sn",iquest[0],nuhead,iuhead[0],iuhead[1],(char*)&iuhead[2]);
      if (strncmp((char*)&iuhead[2],"VERT",4) == 0) {
         fzin(lunfz,0,lvertx,1,"A",nuhead,iuhead,1);
         if (fDebug > 1) printf("JVERTX structure read OKn");
      } else if (strncmp((char*)&iuhead[2],"KINE",4) == 0) {
         fzin(lunfz,0,lkine,1,"A",nuhead,iuhead,1);
         if (fDebug > 1) printf("JKINE structure read OKn");
         if (fDebug > 4) dzshow("kine",0,*lkine,"BSLV",0,0,0,0,4,4);
         if (!H && !D) return 1;
      } else if (strncmp((char*)&iuhead[2],"HITS",4) == 0) {
         fzin(lunfz,0,lhits,1,"A",nuhead,iuhead,1);
         if (fDebug > 1) printf("JHITS structure read OKn");
         if (fDebug > 4) dzshow("hits",0,*lhits,"BSLV",0,0,0,0,4,4);
         if (!D) return 1;
      } else if (strncmp((char*)&iuhead[2],"DIGI",4) == 0) {
         fzin(lunfz,0,ldigi,1,"A",nuhead,iuhead,1);
         if (fDebug > 1) printf("JDIGI structure read OKn");
         if (fDebug > 4) dzshow("digi",0,*ldigi,"BSLV",0,0,0,0,4,4);
         return 1;
      }
   }
}
 
//___________________________________________
 Int_t StarRun::ImportGeantDigits()
{
   TIter next(fMakers);
   StarMaker *maker;
   Int_t totaldigits = 0;
   while((maker = (StarMaker*)next())) {
      totaldigits += maker->ImportGeantDigits();
   }
   if (fDebug > 1) printf(" Total number of imported digits=%dn",totaldigits);
   return totaldigits;
}
 
//___________________________________________
 Int_t  StarRun::ImportGeantDigits(StarMaker *maker, Int_t iset, Int_t ngdet, Int_t *gdet)
{
// Copy Geant3 digits for set iset and all ngdet detectors of array gdet in TClonesArray fDigits.
//
// The following is a description of the original Geant3
// JDIGI Zebra data structure.
//
//-----------------------------------------------
   Int_t i,jdigi,idet;
   Int_t jdi,jdid,ilast,nv,nd,nw,nwdi,ntrt,nwtr;
   Int_t itr,ntrm1,mk,nb,mask,nk,k,iv,ih;
   Int_t setdigits, detdigits;
   const Int_t kMAXVOL = 20;
   const Int_t kMAXTRK = 20;
   Int_t vol[kMAXVOL];
   Int_t digits[kMAXVOL];
   Int_t tracks[kMAXTRK];
   Int_t *nbitsv, *nbitsd;
   jdigi = *ldigi;
   if (jdigi <= 0) return 0;
 
   jdi  = lq[jdigi-iset];
   if (jdi <= 0) return 0;
 
   setdigits = 0;
   for (Int_t id=0;id<ngdet;id++) {
      detdigits = 0;
      idet = gdet[id];
      if (idet <= 0) continue;
 
      jdid = lq[jdi-idet];
      if (jdid <= 0) continue;
 
      ilast = iq[jdi+idet];
      if (ilast == 0) continue;
      nv     = maker->GetNV();
      nd     = maker->GetND();
      nw     = maker->GetNW();
      nbitsv = maker->GetNbitsV();
      nbitsd = maker->GetNbitsD();
      if (nv <= 0) continue;
 
             // Loop on all digits
 
      i    = 0;
      nwdi = 0;
      while (i < ilast) {
         i += nwdi;
         if (i >= ilast) break;
         nwdi  = iq[jdid+i+1];
         if (nwdi <= 0) break;
         nk    = 2;
         ntrm1 = iq[jdid+i+nk] & 65535;
         ntrt  = ntrm1 + 1;
         nwtr  = ntrt/2 + 1;
         nk   += nwtr;
         if (nv > 0) {
            k = 0;
            for (iv=0;iv<nv;iv++) {
               nb = nbitsv[iv];
               if (nb <= 0) {
                  if (k > 0) { k=0; nk++;}
                  vol[iv] = iq[jdid+i+nk];
                  if (iv != nv-1) nk++;
               } else {
                  mask = (1<<nb) -1;
                  if (k+nb > 32) { k=0; nk++;}
                  vol[iv] = (iq[jdid+i+nk]>>k) & mask;
                  k += nb;
               }
            }
            nk++;
         }
 
           // Get track numbers
 
         mk = nk;
         nk = 2;
         if (ntrt > 0) {
            if (ntrm1 >= 1) {
               for (itr=1;itr<=ntrm1;itr+=2) {
                  if (itr <= kMAXTRK) tracks[itr-1] = (iq[jdid+i+nk]>>16) & 65535;
                  nk++;
                  if (itr <  kMAXTRK) tracks[itr]   = iq[jdid+i+nk] & 65535;
               }
            }
            if (ntrt <= kMAXTRK && (ntrt%2 ==1)) {
               tracks[ntrt-1] = (iq[jdid+i+nk]>>16) & 65535;
            }
         }
         nk = mk;
 
           // Get unpacked digits
 
         k = 0;
         Int_t kdigit;
         for (ih=0;ih<nd;ih++) {
            nb = nbitsd[ih];
            if (nb <= 0) {
               if (k > 0) {k = 0; nk++;}
               kdigit = iq[jdid+i+nk];
               nk++;
            } else {
               mask = (1<<nb) -1;
               if (k+nb > 32) { k=0; nk++;}
               kdigit = (iq[jdid+i+nk]>>k) & mask;
               k += nb;
            }
            digits[ih] = kdigit;
         }
         maker->AddDigit(idet,ntrt,tracks,vol,digits);
         detdigits++;
      }
      setdigits += detdigits;
      if (fDebug > 3) printf(" detector digits=%dn",detdigits);
   }
   if (fDebug > 2) printf(" Total digits in %s =%dn",maker->GetName(),setdigits);
   return setdigits;
}
 
//___________________________________________
 Int_t StarRun::ImportGeantHits()
{
   TIter next(fMakers);
   StarMaker *maker;
   Int_t totalhits = 0;
   while((maker = (StarMaker*)next())) {
      totalhits += maker->ImportGeantHits();
   }
   if (fDebug > 1) printf(" Total number of imported hits=%dn",totalhits);
   return totalhits;
}
 
//___________________________________________
 Int_t  StarRun::ImportGeantHits(StarMaker *maker, Int_t iset, Int_t ngdet, Int_t *gdet)
{
// Copy Geant3 hits for set iset and all ngdet detectors of array gdet in TClonesArray fHits.
//
// The following is a description of the original Geant3
// JSET and JHITS Zebra data structures.
//
//--------------------JSET-----------------------
// IUSET     set identifier (4 characters), user defined
// IUDET     detector identifier  (4 characters),   name of  an existing volume
// NV        number of volume descriptors
// NAMESV    vector of NV volume descriptors (4 characters)
// NBITSV    vector of  NV bit numbers  for packing  the volume numbers
// ISET      position of set in bank JSET
// IDET      position of detector in bank JS=IB(JSET-ISET)
// Remarks:
// - The vector NAMESV (length NV)  contains the list of volume
//   names which  permit the  identification of  every physical
//   detector with detector name IUDET.
// - Each  element of  the vector  NBITSV (length  NV)  is  the
//   number  of  bits  used  for  packing  the  number  of  the
//   corresponding volume,  when building the packed identifier
//   of a given physical detector.
//       IQ(JSET+ISET) = IUSET
//       JS = LQ(JSET-ISET) + pointer to set parameters
//       IQ(JS+IDET)=IUDET
//       JD= LQ(JS-1)  = pointer to detector IDET
//       IQ(JD+1)=Total number of words to store packed volumes
//       IQ(JD+2)=NV
//       IQ(JD+3)=Number of words required per hit
//       IQ(JD+4)=Number of different hits types
//       IQ(JD+5)=Number of words per digit
//       IQ(JD+6)=Number of different digit types
//       IQ(JD+7)=number of words for allocation of HITS banks
//       IQ(JD+8)=number of words for allocation of DIGI banks
//       IQ(JD+9)=Number of paths through the JVOLUM tree
//       IQ(JD+10)= For an alias only, IDET of main detector
//       IQ(JD+2*I+9) = name of volume i = NAMESV(I)
//       IQ(JD+2*I+10)= number of bits/volume = NBITSV(I)
//
//            The Detector Set data structure JSET
//            ------------------------------------
//
//                                        | JSET
//    NSET            ISET                v         NSET
//     ................................................
//     |              | |               |  | Set names|
//     ................................................
//                     | JS
//                     |
//    NDET       IDET  v                    NDET
//     ........................................
//     |        |  |  | | Detector names      |
//     ........................................
//               | JD
//      -3 -2 -1 v
//     ................................................
//     |  |  |  |  |   Parameters of GSDET            |
//     ................................................
//      |  |  |
//      |  |  |  JDH
//      |  |  |
//      |  |  |           .............................
//      |  |  ............| Parameters of GSDETH      |
//      |  |              .............................
//      |  |
//      |  | JDD
//      |  |
//      |  |              .............................
//      |  ...............|  Parameters of GSDETD     |
//      |                 .............................
//      |
//      |  JDU
//      |                 .............................
//      ..................| Parameters of GSDETU      |
//                        .............................
//
//  JS = LQ(JSET-ISET) pointer to detector set number ISET
// The JSET data structure is filled by GSDET, GSDETH,
// GSDETD, GSDETU and eventually by GSDETA.
//
//--------------------JHITS----------------------
//       JH  = LQ(JHITS-ISET)
//       JHD = LQ(JH-IDET)
//       IQ(JH+IDET) = pointer to LAST USED word in JHD
//         Each hit is packed into JHD in the following format
//        --Track number ITRA not packed
//        --Volume numbers packed
//        --Hits transformed and packed
//
//                The Hit data structure JHITS
//                ----------------------------
//
//                                            | JHITS
//    NSET                  ISET              v
//     ..........................................
//     |                   |  |              |  |
//     ..........................................
//                           |
//                           | JH
//      NDET    IDET         v               NDET
//        .....................................
//        |      |  |      |  |               |
//        .....................................
//                |
//                | JHD
//                v
//               .........................................
//               |  | 1st hit        | 2nd hit, etc.     |
//               .........................................
//                         Bank layout
// JH            =  LQ(JHITS-ISET,)  pointer  to  hits for  set  number ISET
// JHD           = LQ(JH-IDET),   pointer to  hits of  detector IDET of set ISET
// IQ(JH+IDET)    number of words  used so far for  storing the hits of detector IDET
// IQ(JHD+1)      1st word of 1st hit
// IQ(JHD+NWH+1)  1st word of 2nd hit
//                      JS=LQ(JSET-ISET)
//                      JD=LQ(JS-IDET)
//                      NWH=IQ(JD+3)
//
// The Geant3 common PAWC is defined as:
//     common/pawc/paw(nwpaw)
//     INTEGER   IQ(2), LQ(8000)
//     REAL            Q(2)
//     EQUIVALENCE (LQ(1),paw(11)),(IQ(1),paw(19)),(Q(1),paw(19))
//
// The Geant3 common GCLINK is defined as:
//     COMMON/GCLINK/JDIGI ,JDRAW ,JHEAD ,JHITS ,JKINE ,JMATE ,JPART
//     ,JROTM ,JRUNG ,JSET  ,JSTAK ,JGSTAT,JTMED ,JTRACK,JVERTX
//     ,JVOLUM,JXYZ  ,JGPAR ,JGPAR2,JSKLT
//-----------------------------------------------
   Int_t i,track,jhits,idet;
   Int_t jhd,ilast,nv,nh,nw;
   Int_t nb,mask,nk,k,iv,ih;
   Int_t *nbitsv, *nbitsh;
   Float_t *factor, *origin;
   jhits = *lhits;
   if (jhits <= 0) return 0;
 
   Int_t jh  = lq[jhits-iset];
   if (jh <= 0) return 0;
 
   Int_t sethits = 0;
   for (Int_t id=0;id<ngdet;id++) {
      idet = gdet[id];
      if (idet <= 0) return sethits;
 
      jhd = lq[jh-idet];
      if (jhd <= 0) return sethits;
 
      Int_t dethits = 0;
      ilast = iq[jh+idet];
      if (ilast == 0) return sethits;
      nv     = maker->GetNV();
      nh     = maker->GetNH();
      nw     = maker->GetNW();
      nbitsv = maker->GetNbitsV();
      nbitsh = maker->GetNbitsH();
      factor = maker->GetFactor();
      origin = maker->GetOrigin();
      Int_t vol[20];
      Float_t hits[20];
 
             // Loop on all hits
 
      for (i=1;i<=ilast;i +=nw) {
         track = iq[jhd+i];
         nk = 1;
         if (nv > 0) {
            k = 0;
            for (iv=0;iv<nv;iv++) {
               nb = nbitsv[iv];
               if (nb <= 0) {
                  if (k > 0) { k=0; nk++;}
                  vol[iv] = iq[jhd+i+nk];
                  if (iv != nv-1) nk++;
               } else {
                  mask = (1<<nb) -1;
                  if (k+nb > 32) { k=0; nk++;}
                  vol[iv] = (iq[jhd+i+nk]>>k) & mask;
                  k += nb;
               }
            }
            nk++;
         }
           // Get unpacked hits
           // Hits origin is shifted . Division by scale factor
//
         k = 0;
         Int_t khit;
         for (ih=0;ih<nh;ih++) {
            nb = nbitsh[ih];
            if (nb <= 0) {
               if (k > 0) {k = 0; nk++;}
               khit = iq[jhd+i+nk];
               nk++;
            } else {
               mask = (1<<nb) -1;
               if (k+nb > 32) { k=0; nk++;}
               khit = (iq[jhd+i+nk]>>k) & mask;
               k += nb;
            }
            hits[ih] = Float_t(khit)/factor[ih] - origin[ih];
         }
         maker->AddHit(idet,track,vol,hits);
         dethits++;
      }
      sethits += dethits;
      if (fDebug > 3) printf(" detector hits=%dn",dethits);
   }
   if (fDebug > 2) printf(" Total hits in %s =%dn",maker->GetName(),sethits);
   return sethits;
}
 
//___________________________________________
 Int_t StarRun::ImportGeantKine()
{
//  Copy the Geant data structures JVERTX and JKINE into a TClonesArray
//  of GParticles
 
   Int_t it,jk,ipart,ks,kfcode,ivor,jvor,parent;
   Int_t nvsons,ivson1,ivsonl,jvson1,jvsonl;
   Int_t ntson1,ntsonl,firstchild,lastchild;
   Float_t xmass, tlife, time;
   Float_t E,px,py,pz,vx,vy,vz;
 
   Int_t jkine   = *lkine;
// Int_t jpart   = *lpart;
   Int_t jvertx  = *lvertx;
   if (jkine  <= 0) return 0;
   if (jvertx <= 0) return 0;
 
   TClonesArray &particles = *fParticles;
   Int_t ntrack = iq[jkine+1];
 
   for (it=1;it<=ntrack;it++) {
      jk     = lq[jkine-it]; if (jk <= 0) continue;
      ipart  = Int_t(q[jk+5]);   //particle nr in JPART
//    jp     = lq[jpart-ipart];
      tlife  = 0; //should be q[jp+9]
      px     = q[jk+1];
      py     = q[jk+2];
      pz     = q[jk+3];
      E      = q[jk+4];
      xmass  = TMath::Sqrt(E*E -px*px-py*py-pz*pz);
      ks     = 0;
      kfcode = ipart;               //this should be set to LUJETS convention
      ivor   = Int_t(q[jk+6]);  //vertex number origin of this track
      if (ivor) {
         jvor   = lq[jvertx-ivor];
         parent = Int_t(q[jvor+5]);
         time   = q[jvor+4];    //time of vertex generation
         vx     = q[jvor+1];
         vy     = q[jvor+2];
         vz     = q[jvor+3];
      } else {
         parent = 0;
         time   = 0;
         vx = vy = vz = 0;
      }
      nvsons = Int_t(q[jk+7]);
      if (nvsons) {
         ivson1     = Int_t(q[jk+8]);             //first generated vertex
         ivsonl     = Int_t(q[jk+nvsons+7]);      //last generated vertex
         jvson1     = lq[jvertx-ivson1];          //pointer to first gen vertex
         jvsonl     = lq[jvertx-ivsonl];          //pointer to last gen vertex
         ntson1     = Int_t(q[jvson1+7]);         //number of tracks of first gen vertex
         firstchild = Int_t(q[jvson1+8]);         //first track of first gen vertex
         ntsonl     = Int_t(q[jvsonl+7]);         //number of tracks of last gen vertex
         lastchild  = Int_t(q[jvsonl+ntsonl+7]);  //last track of last gen vertex
      } else {
         firstchild = 0;
         lastchild  = 0;
      }
      new(particles[it-1]) GParticle(ks,kfcode,parent,firstchild,lastchild,
         px,py,pz,E,xmass,vx,vy,vz,time,tlife);
   }
   if (fDebug > 1) printf("ImportKine found %d particlesn",ntrack);
   return ntrack;
}
 
//___________________________________________
 void StarRun::Init()
{
   TIter next(fMakers);
   TObject *objfirst, *objlast;
   StarMaker *maker;
   while((maker = (StarMaker*)next())) {
     // save last created histogram in current Root directory
      objlast = gDirectory->GetList()->Last();
 
     // Initialise maker
      maker->Init();
 
     // Add Maker histograms in Maker list of histograms
      if (objlast) objfirst = gDirectory->GetList()->After(objlast);
      else         objfirst = gDirectory->GetList()->First();
      while (objfirst) {
         maker->Histograms()->Add(objfirst);
         objfirst = gDirectory->GetList()->After(objfirst);
      }
   }
}
 
//___________________________________________
 void StarRun::InitGeant()
{
}
 
//___________________________________________
 void StarRun::InitZebra(const Int_t nzebra)
{
// Initialize the Zebra package.
// The Zebra space is allocated dynamically
// nzebra is the size of the Zebra store in words.
 
   Int_t *zebra = new Int_t[nzebra];
   lq        = &zebra[9];
   iq        = &zebra[17];
   void *qq  = iq;
   q         = (float*)qq;
   iquest    = &quest[0];
   int ixdiv = 0;
 
   mzebra(-1);
   Int_t ixstore = 0;
   mzstor(ixstore,"z"," ",&zebra[5],&zebra[10],&zebra[11],&zebra[11],&zebra[2000],&zebra[nzebra],1,1);
 
//           Create a permanent link area for master pointers
 
   ldigi  = &gclink[0];
   lhead  = &gclink[2];
   lhits  = &gclink[3];
   lkine  = &gclink[4];
   lvertx = &gclink[14];
   lpart  = &gclink[6];
   lsklt  = &gclink[19];
 
   mzlink(ixdiv,"/GCLINK/",ldigi,lsklt,ldigi,8);
}
 
//___________________________________________
 void StarRun::MakeTree(Option_t *option)
{
//  Create the ROOT trees
//  Loop on all makers to create the Root branch (if any)
 
   char *K = strstr(option,"K");
   char *H = strstr(option,"H");
   char *D = strstr(option,"D");
   if (K && !fTreeK) fTreeK = new TTree("TK","Kinematics");
   if (H && !fTreeH) fTreeH = new TTree("TH","Hits");
   if (D && !fTreeD) fTreeD = new TTree("TD","Digits");
 
   TIter next(fMakers);
   StarMaker *maker;
   while((maker = (StarMaker*)next())) {
      maker->MakeBranch();
   }
 
//  Create a branch for particles
   if (fTreeK) fTreeK->Branch("Particles",&fParticles,4000);
}
 
//___________________________________________
 void StarRun::OpenGeantFile(const char *fname, Int_t lrecl, const char *option)
{
// This function assumes a Geant file in Zebra/FZ format
// created with the C I/O option "L"
// If not, one must create a Fortran routine with:
//      open(unit=1,file=....)
//      call fzfile(1,0,"xi")
 
   printf("Openning Geant file:%sn",fname);
   int ier;
   int len = strlen(fname);
   cfopen(lunfz,0,0,"r ",0,fname,ier,2,len);
   iquest[0] = lunfz;
   Int_t lopt = strlen(option);
   fzfile(lunfz,lrecl,option,lopt);
}
 
//___________________________________________
 void StarRun::Paint(Option_t *option)
{
   TIter next(fMakers);
   StarMaker *maker;
   while((maker = (StarMaker*)next())) {
      maker->Paint(option);
   }
}
 
//___________________________________________
 void StarRun::RunGeant()
{
}
 
//___________________________________________
 void StarRun::Streamer(TBuffer &R__b)
{
   // Stream an object of class XRun.
 
   if (R__b.IsReading()) {
      Version_t R__v = R__b.ReadVersion(); if (R__v) { }
      TNamed::Streamer(R__b);
      if (!Run) Run = this;
      gROOT->GetListOfBrowsables()->Add(this,"Run");
      R__b >> fRun;
      R__b >> fEvent;
      R__b >> fNvertex;
      R__b >> fNtrack;
      fTreeK = (TTree*)gDirectory->Get("TK");
      fTreeH = (TTree*)gDirectory->Get("TH");
      fTreeD = (TTree*)gDirectory->Get("TD");
      R__b >> fMakers;
      R__b >> fParticles;
      if (fTreeK == 0) return;
      TBranch *branch = fTreeK->GetBranch("Particles");
      if (branch)  branch->SetAddress(&fParticles);
   } else {
      R__b.WriteVersion(StarRun::IsA());
      TNamed::Streamer(R__b);
      R__b << fRun;
      R__b << fEvent;
      R__b << fNvertex;
      R__b << fNtrack;
      fTreeK->Write(); // added
      fTreeH->Write(); // added
      fTreeD->Write(); // added
      R__b << fMakers;
      R__b << fParticles;
   }
}
 


ROOT page - Class index - Top of the page

This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.