/////////////////////////////////////////////////////////////////////////////// // fill RPC branches of GEANT ROOT tree // // Modified on 05/12/2008 by I. Koenig. // Modified on 23/12/2009 by D. Gonzalez-Diaz. // Implemented a fast algorithm for filtering // tracks with only valid hits at the RPC box. // Modified on 19/07/2012 by Alessio Mangiarotti. // The RPC geometry file in Geant has the lateral // columns exchanged with respect to the accepted // convention and Hydra. As suggested by Jochen // Markert, the columns are swapped here and no // further action is necessary in the digitiser or // in the hit finder. // Modified on 20/09/2012 by Alessio Mangiarotti. // Removed bug that lead to a loss of good hits // (i.e. hits in RPC gaps) by secondary particles // produced in materials located inside the RPC box. // Modified on 14/03/2013 by Alessio Mangiarotti. // Fillrpc adapted to the new HGeantRpc structure // which is organised in a cellwise fashion. The // combinations of all four gaps into a cell is now // carried on here (before it was performed in the // digitiser). The code is very similar to the one // in the digitiser: two new internal structs // celltrack and celldat have been added to handle // the reorganisation of the gaps into cells. // The use of the old HGeantRpc structure (but now // filled with average per cell information) can still // be restored by commenting the "NEW after Jan 2013" // and uncommenting the "OLD before Jan 2013" parts // of the code below. A further reduction in output // size can be obtained by saving only the average gap // information (see full document for performance // details). To obtain this configuration comment // the "ALL information" and uncomment // "ONLY average information". /////////////////////////////////////////////////////////////////////////////// #include "fillrpc.h" #include "cfortran.h" #include "hgeantrpc.h" #include "hades.h" #include "hlocation.h" #include "hcategory.h" #include "hrecevent.h" #include "hpartialevent.h" #include "hgeantkine.h" #include "hgeantmaxtrk.h" struct RPC_DEF { // structure mapped to the common block RPCTUPLE Int_t ntrk; Int_t rpctrk[MAXTRKRPC]; Int_t rpcdet[MAXTRKRPC]; Float_t rpce[MAXTRKRPC]; Float_t rpcx[MAXTRKRPC]; Float_t rpcy[MAXTRKRPC]; Float_t rpcz[MAXTRKRPC]; Float_t rpctheta[MAXTRKRPC]; Float_t rpcphi[MAXTRKRPC]; Float_t rpctof[MAXTRKRPC]; Float_t rpcmom[MAXTRKRPC]; Float_t rpclen[MAXTRKRPC]; Float_t rpcloclen[MAXTRKRPC]; }; #define RPCTUPLE COMMON_BLOCK(RPCTUPLE,rpctuple) COMMON_BLOCK_DEF(RPC_DEF,RPCTUPLE); typedef struct celltrack { // stores all needed info per gap Float_t gape[4]; // GEANT energy loss Float_t gapx[4]; // GEANT x position Float_t gapy[4]; // GEANT y position Float_t gapz[4]; // GEANT z position Float_t gaptheta[4]; // GEANT theta Float_t gapphi[4]; // GEANT phi Float_t gaptof[4]; // GEANT time Float_t gapmom[4]; // GEANT momentum Float_t gaplen[4]; // GEANT track length Float_t gaploclen[4]; // GEANT local track length Float_t ene; // GEANT average energy loss Float_t posx; // GEANT average x position Float_t posy; // GEANT average y position Float_t posz; // GEANT average z position Float_t theta; // GEANT average theta Float_t phi; // GEANT average phi Float_t tof; // GEANT average time Float_t mom; // GEANT average momentum Float_t len; // GEANT average track length Float_t loclen; // GEANT average local track length Int_t det; // GEANT detector index Int_t track; // GEANT track number at the box // Reset initial values. void reset(){ for(UInt_t i=0;i<4;i++){ gape[i] = 0.; gapx[i] = 0.; gapy[i] = 0.; gapz[i] = 0.; gaptheta[i] = 0.; gapphi[i] = 0.; gaptof[i] = 0.; gapmom[i] = 0.; gaplen[i] = 0.; gaploclen[i] = 0.; } track = -1; ene = 0.; posx = 0.; posy = 0.; posz = 0.; theta = 0.; phi = 0.; tof = 0.; mom = 0.; len = 0.; loclen = 0.; det = -999; } // Calculate the mean of all crossed gaps void calcMeans(){ UInt_t ii=0; for(UInt_t i=0;i<4;i++){ if(gaptof[i]>0.) { ene += gape[i]; posx += gapx[i]; posy += gapy[i]; posz += gapz[i]; theta += gaptheta[i]; phi += gapphi[i]; tof += gaptof[i]; mom += gapmom[i]; len += gaplen[i]; loclen += gaploclen[i]; ii++; } } if(ii>0) { Float_t aii=(Float_t)ii; ene /= aii; posx /= aii; posy /= aii; posz /= aii; theta /= aii; phi /= aii; tof /= aii; mom /= aii; len /= aii; loclen /= aii; } } } celltrack; typedef struct { // collection of all celltracks per cell //Short_t linIndex; // linear index in working array // (= column * maxCell + cell ) Short_t col; // column Short_t cell; // cell vector celltr;// gaptracks stored according to left side of the cell void reset(){ // Reset initial values. deletes right/left vector //linIndex= -1; col = -999; cell = -999; for(UInt_t i=0;igetCurrentEvent()); HPartialEvent* pSimul = pHGeantEvent->getPartialEvent(catSimul); HCategory* pRpcCat = pSimul->getCategory(catRpcGeantRaw); HCategory* pKin = pSimul->getCategory(catGeantKine); if(*sector == 1) pRpcCat->Clear(); //--------------------------------------------- // create working arrays (vector of hyp objects hold dynamic lists on hyp) // (to be cleared at end of function, otherwise meomory leak!) vector cellobjects; //! temporary working array for regrouping gaps into cells cellobjects.resize(maxTot+1, 0 ); // size is constant, but objects not yet created! //--------------------------------------------- Int_t index=-1, column=-999, cell=-999, gap=-999, ind=0; Int_t track=0, track_sec=0, track_sec2=0, track_main=0, trk; Int_t Id, track_tmp; Bool_t hasNoDaughterInCell[MAXTRKRPC], hasNot=kFALSE, WasAlreadyTagged=kFALSE; Float_t Time_prev=0, Time=0; //Beginning of loop for tagging tracks at the entrance of the RPC box //with daughters having a Geant hit in an RPC cell. //The main purpose of the tagging is to REMOVE tracks that cross the //RPC box but do not leave any hit to save disk space. for(Int_t i=0; igetObject(track-1))->getParentTrack(); j++; } if(j==nbranches) cout << "fillrpc: exhausted depth of search " << endl; if(!hasNot) { // Bug fix AM 20-09-2012 ! if (track_sec!=track_main) { for(Int_t trk_sec2=0; trk_sec2=0) { if (hasNoDaughterInCell[trk_sec2]) { cout << "fillrpc: something strange with flagging for track " << RPCTUPLE.rpctrk[trk_sec2] << endl; // Uncomment this to print out more information! //for(Int_t trk_sec3=0; trk_sec3getObject(track-1))->getParentTrack(); // cout << RPCTUPLE.rpctrk[trk_sec3] // << " " << RPCTUPLE.rpcdet[trk_sec3] // << " " << hasNoDaughterInCell[trk_sec3] // << " " << RPCTUPLE.rpctof[trk_sec3] // << " " << track << endl; //} } } } // Keep everything for(Int_t trk_sec2=0; trk_sec2=0) { // RPC CELL column = RPCTUPLE.rpcdet[trk]/(35*4); cell = (RPCTUPLE.rpcdet[trk] - column*(35*4))/4; gap = RPCTUPLE.rpcdet[trk] - column*(35*4) - cell*4; // // ATTENTION ATTENTION ATTENTION ATTENTION ATTENTION ATTENTION // Correction for the mistake in the Geant RPC geometry file // ATTENTION ATTENTION ATTENTION ATTENTION ATTENTION ATTENTION // if (column==0) { column=4; } else if (column==1) { column=5; } else if (column==4) { column=0; } else if (column==5) { column=1; } // // ATTENTION ATTENTION ATTENTION ATTENTION ATTENTION ATTENTION // // if the hit is in an roc cell, objects are stored linear ind = column * maxCell + cell; } else { // EBOX column=RPCTUPLE.rpcdet[trk]; cell=RPCTUPLE.rpcdet[trk]; gap=0; // if hit is at ebox, use special position at maxCol*maxCell+1 ind = maxTot; } if(cellobjects[ind] == 0) { // first time we have to create the object cellobjects[ind] = new celldat; //rpcobjects[ind]->linIndex= ind; cellobjects[ind]->col = column; cellobjects[ind]->cell = cell; } celldat* cdat = cellobjects[ind]; //--------------------------------------------------------------- // create new objects for the track crossing the cell // and add them to the list of tracks for this cell celltrack* celltr; // CASE I: no tacks in this cell yet or track crossing ebox if((cdat->celltr.size()==0)||(column<0)) { celltr = new celltrack; celltr->reset(); cdat->celltr.push_back(celltr); // CASE II: this cell has already tracks associated // with it } else { // Look if current track is already present in the list for(UInt_t i=0;icelltr.size();i++){ celltr = cdat->celltr[i]; // CASE II a: the current track is already associated, // so select it if(celltr->track == track) break; } // CASE II b: the current track is not associated, // so add it if(celltr->track != track) { celltr = new celltrack; celltr->reset(); cdat->celltr.push_back(celltr); } } //--------------------------------------------------------------- // The condition on the time takes care of possible multiple // re-entries of the same track. if(celltr->gaptof[gap]==0.) { celltr->gape[gap] = RPCTUPLE.rpce[trk]; celltr->gapx[gap] = RPCTUPLE.rpcx[trk]; celltr->gapy[gap] = RPCTUPLE.rpcy[trk]; celltr->gapz[gap] = RPCTUPLE.rpcz[trk]; celltr->gaptheta[gap] = RPCTUPLE.rpctheta[trk]; celltr->gapphi[gap] = RPCTUPLE.rpcphi[trk]; celltr->gaptof[gap] = RPCTUPLE.rpctof[trk]; celltr->gapmom[gap] = RPCTUPLE.rpcmom[trk]; celltr->gaplen[gap] = RPCTUPLE.rpclen[trk]; celltr->gaploclen[gap] = RPCTUPLE.rpcloclen[trk]; celltr->track = track; celltr->det = RPCTUPLE.rpcdet[trk]; } else if(celltr->gaptof[gap]gape[gap] = RPCTUPLE.rpce[trk]; celltr->gapx[gap] = RPCTUPLE.rpcx[trk]; celltr->gapy[gap] = RPCTUPLE.rpcy[trk]; celltr->gapz[gap] = RPCTUPLE.rpcz[trk]; celltr->gaptheta[gap] = RPCTUPLE.rpctheta[trk]; celltr->gapphi[gap] = RPCTUPLE.rpcphi[trk]; celltr->gaptof[gap] = RPCTUPLE.rpctof[trk]; celltr->gaplen[gap] = RPCTUPLE.rpclen[trk]; celltr->gaploclen[gap] = RPCTUPLE.rpcloclen[trk]; celltr->track = track; celltr->det = RPCTUPLE.rpcdet[trk]; //printf("Reentry of track=%i in RPC cell Sector=%i Module=%i Cell=%i\n", // track,*sector-1,column,cell); } } // End of loop for filling celldat // This is the second loop to save all cellwise information // in HGeantRpc. // Counter of stored cells per sector trk=0; for(UInt_t i = 0; i < cellobjects.size(); i ++) { celldat* cdat = cellobjects[i]; if(cdat) { column = cdat->col; cell = cdat->cell; for(UInt_t j = 0; j < cdat->celltr.size(); j ++) { celltrack* celltr = cdat->celltr[j]; celltr->calcMeans(); track = celltr->track; loc.set(2,(*sector-1),trk); HGeantRpc* pRpc = (HGeantRpc*)pRpcCat->getSlot(loc,&index); if(pRpc==NULL) { printf("\n*** FILLRPC: HGeantRpc location not found ! ***\n"); printf("sector=%d track=%d\n",(*sector-1),trk); return; } pRpc = new(pRpc) HGeantRpc; pRpc->setTrack(track); // // Old before Jan 2013 // //pRpc->setHit(celltr->posx,celltr->posy,celltr->posz, // celltr->tof,celltr->mom,celltr->ene); //pRpc->setTLength(celltr->len,celltr->loclen); // // New after Jan 2013 // // Use this to save ONLY average information //pRpc->setHit(celltr->posx,celltr->posy,celltr->posz,celltr->tof, // celltr->mom,celltr->ene,celltr->loclen); //// Use this to save ONLY average information ////for(UInt_t k = 0; k < 4; k ++) { //// if(k==0) { //// pRpc->setGap(k,celltr->posx,celltr->posy,celltr->mom, //// celltr->ene,celltr->loclen); //// } else { //// pRpc->setGap(k,0.,0.,0.,0.,0.); //// } ////} // Use this to save ALL information for(UInt_t k = 0; k < 4; k ++) { pRpc->setGap(k,celltr->gapx[k],celltr->gapy[k],celltr->gapmom[k], celltr->gape[k],celltr->gaploclen[k]); } // Common to all NEW after Jan 2013 pRpc->setZHit(celltr->posz); pRpc->setTofHit(celltr->tof); pRpc->setTLengthHit(celltr->len); pRpc->setVersion(5); // // Common to OLD before Jan 2013 and NEW after Jan 2013 // pRpc->setIncidence(celltr->theta,celltr->phi); if(celltr->det<0) pRpc->setDetectorID(celltr->det); // Gap is always zero for hit in EBOX else pRpc->setAddress(*sector-1, column, cell, 0); // Add TClonesArray index to access from HGeantKine if(pKin != NULL && track > 0) { ((HGeantKine*)pKin->getObject(track-1))->setRpcHitIndex(index); pRpc->sortVariable = celltr->tof; } // Increment counter of stored cells per sector trk++; } } } //---------------------------------------------- // clean the working array for(UInt_t i = 0; i < cellobjects.size(); i ++) { celldat* cdat = cellobjects[i]; if(cdat) { cdat->reset(); } delete cdat; } cellobjects.clear(); //---------------------------------------------- // End of loop for filling HGeantRpc // Sort linked lists by sortVariable = celltr->tof if(pKin != NULL && *sector == 6) { Int_t nTracks = pKin->getEntries(); for(Int_t i=0;igetObject(i))->sortRpcHits(); } } } // wrapping of the c++ routine as Fortran-subroutine using cfortran.h // of Cern-library FCALLSCSUB1(fillrpc,FILLRPC,fillrpc,PINT)