///////////////////////////////////////////////////////////////////////////// // // fill FRpc branches of GEANT ROOT tree // // last modified on 30/11/2015 by I. Koenig ///////////////////////////////////////////////////////////////////////////// #include "fillfrpc.h" #include "cfortran.h" #include "hgeantfrpc.h" #include "hgeantmaxtrk.h" #include "hades.h" #include "hlocation.h" #include "hcategory.h" #include "hrecevent.h" #include "hpartialevent.h" #include "hgeantkine.h" #include struct FRPC_DEF { // structure mapped to FRPCTUPLE Fortran COMMON block Int_t ntrk; Int_t frpctrk[MAXTRKFRPC]; Int_t frpcpart[MAXTRKFRPC]; Int_t frpcmod[MAXTRKFRPC]; Int_t frpclayer[MAXTRKFRPC]; Int_t frpccell[MAXTRKFRPC]; Int_t frpcsub[MAXTRKFRPC]; Float_t frpcx[MAXTRKFRPC]; Float_t frpcy[MAXTRKFRPC]; Float_t frpcz[MAXTRKFRPC]; Float_t frpcpx[MAXTRKFRPC]; Float_t frpcpy[MAXTRKFRPC]; Float_t frpcpz[MAXTRKFRPC]; Float_t frpctof[MAXTRKFRPC]; Float_t frpclen[MAXTRKFRPC]; Float_t frpce[MAXTRKFRPC]; }; #define FRPCTUPLE COMMON_BLOCK(FRPCTUPLE,frpctuple) COMMON_BLOCK_DEF(FRPC_DEF,FRPCTUPLE); void fillfrpc(void) { HRecEvent* pHGeantEvent = (HRecEvent*)(gHades->getCurrentEvent()); // get pointer to event HPartialEvent* pSimul = pHGeantEvent->getPartialEvent(catSimul); HCategory* pFRpcCat = pSimul->getCategory(catFRpcGeantRaw); // get category pFRpcCat->Clear(); // clear category if 1st sector HCategory* pKin = pSimul->getCategory(catGeantKine); // get KINE cat Int_t gmodule, glayer, gcell, gsub; Int_t old_hitid = -1; Hit old_hit; old_hit.track = -1; old_hit.module = -1; old_hit.layer = -1; old_hit.cell = -1; old_hit.subcell = 0; old_hit.mult = 0; // loop over tracks (TClonesArray, expands if necessary) for(Int_t idx = 0; idx < FRPCTUPLE.ntrk; ++idx) { gmodule = FRPCTUPLE.frpcmod[idx]-1; //GEANT indexing is from 1 glayer = FRPCTUPLE.frpclayer[idx]-1; gcell = FRPCTUPLE.frpccell[idx]-1; gsub = FRPCTUPLE.frpcsub[idx]; Hit hit; hit.module = -1; hit.layer = -1; hit.cell = -1; hit.subcell = 0; hit.mult = 1; hit.index = idx; hit.track = FRPCTUPLE.frpctrk[idx]; hit.px = FRPCTUPLE.frpcpx[idx]; hit.py = FRPCTUPLE.frpcpy[idx]; hit.pz = FRPCTUPLE.frpcpz[idx]; hit.x = FRPCTUPLE.frpcx[idx]; hit.y = FRPCTUPLE.frpcy[idx]; hit.z = FRPCTUPLE.frpcz[idx]; hit.tof = FRPCTUPLE.frpctof[idx]; hit.trklen = FRPCTUPLE.frpclen[idx]; hit.eloss = FRPCTUPLE.frpce[idx]; char _gsub = gsub-1; Int_t hitid = gcell + // rpc pad number 1..32 / 0..31 100*glayer + // rpc layer 1..4 / 0..3 1000*gmodule + // rpc module, currently always 7 / 6 10000*hit.track; hit.module = gmodule; hit.layer = glayer; hit.cell = gcell; if (hitid != old_hitid) { if (old_hitid != -1) { old_hit.px /= old_hit.mult; old_hit.py /= old_hit.mult; old_hit.pz /= old_hit.mult; old_hit.x /= old_hit.mult; old_hit.y /= old_hit.mult; old_hit.z /= old_hit.mult; fillFRpcHit(pKin, pFRpcCat, old_hit); } old_hitid = hitid; old_hit = hit; // see comment after "else" old_hit.subcell = (old_hit.subcell & ~(0b11 << _gsub*2)) | (1 << _gsub*2); } else { ++old_hit.mult; // subcell is a 8 bit word, where each two bits: 0b aa bb cc dd // store number o subcels fired for this track, usually 0, 1, 2, 3 or more char sc = (old_hit.subcell >> _gsub*2) & 0b11; if (sc < 3) ++sc; old_hit.subcell = (old_hit.subcell & ~(0b11 << _gsub*2)) | (sc << _gsub*2); old_hit.px += hit.px; old_hit.py += hit.py; old_hit.pz += hit.pz; old_hit.x += hit.x; old_hit.y += hit.y; old_hit.z += hit.z; if (old_hit.tof > hit.tof) old_hit.tof = hit.tof; if (old_hit.trklen > hit.trklen) old_hit.trklen = hit.trklen; old_hit.eloss += hit.eloss; } } if (old_hitid != -1) { old_hit.px /= old_hit.mult; old_hit.py /= old_hit.mult; old_hit.pz /= old_hit.mult; old_hit.x /= old_hit.mult; old_hit.y /= old_hit.mult; old_hit.z /= old_hit.mult; fillFRpcHit(pKin, pFRpcCat, old_hit); } } Bool_t fillFRpcHit(HCategory* pKin, HCategory* pFRpcCat, const Hit & hit) { HLocation loc; loc.set(0); HGeantFRpc* pFRpc = (HGeantFRpc*)pFRpcCat->getNewSlot(loc); // get slot if (pFRpc == NULL) { printf("\n*** FILLFRPC: HGeantFRpc location not found ! ***\n"); printf("index=%d\n", hit.index); return kFALSE; } pFRpc = new(pFRpc) HGeantFRpc; // FRpc track data if (!pFRpc) return kFALSE; pFRpc->setTrack(hit.track); pFRpc->setAddress(hit.module, hit.layer, hit.cell, hit.subcell); pFRpc->setHit(hit.x, hit.y, hit.z, hit.px, hit.py, hit.pz, hit.tof, hit.trklen, hit.eloss); // add hit to linked list of track if (pKin != NULL && hit.track > 0) { ((HGeantKine*)pKin->getObject(hit.track-1))->setFRpcHitIndex(hit.index); } return kTRUE; } // wrapping of the c++ routine as Fortran-subroutine using cfortran.h of Cern-library FCALLSCSUB0(fillfrpc,FILLFRPC,fillfrpc)