/////////////////////////////////////////////////////////////////////////////// // fill Rich branches of GEANT ROOT tree // // last modified on 05/12/2008 by I. Koenig /////////////////////////////////////////////////////////////////////////////// #include "fillrich.h" #include "cfortran.h" #include "hgeantrich.h" #include "hgeantmaxtrk.h" #include "hades.h" #include "hlocation.h" #include "hcategory.h" #include "hrecevent.h" #include "hpartialevent.h" #include "hgeantkine.h" struct RICH_DEF { // structure mapped to RICHTUPS Fortran COMMON block Int_t ncer; Float_t xcer[MAXCKOVRICH]; Float_t ycer[MAXCKOVRICH]; Float_t ecer[MAXCKOVRICH]; Int_t parentcer[MAXCKOVRICH]; Int_t nhit; Float_t xhit[MAXPARTRICH]; Float_t yhit[MAXPARTRICH]; Float_t zhit[MAXPARTRICH]; Int_t parthit[MAXPARTRICH]; Float_t momhit[MAXPARTRICH]; Float_t thetahit[MAXPARTRICH]; Float_t phihit[MAXPARTRICH]; Float_t elosshit[MAXPARTRICH]; Float_t tlenhit[MAXPARTRICH]; Int_t trkhit[MAXPARTRICH]; Int_t nmir; Int_t numphot[MAXMIRRRICH]; Int_t leptrack[MAXMIRRRICH]; Float_t xring[MAXMIRRRICH]; Float_t yring[MAXMIRRRICH]; Float_t xlep[MAXMIRRRICH]; Float_t ylep[MAXMIRRRICH]; Float_t zlep[MAXMIRRRICH]; Int_t lepid[MAXMIRRRICH]; }; struct RICH_FLAG { // structure mapped to RICHFLAG common Int_t isPMT; }; #define RICHTUPLE COMMON_BLOCK(RICHTUPLE,richtups) COMMON_BLOCK_DEF(RICH_DEF,RICHTUPLE); #define RICHFLAG COMMON_BLOCK(RICHFLAG,richflag) COMMON_BLOCK_DEF(RICH_FLAG,RICHFLAG); void fillrich(int* sector) { HLocation loc; HRecEvent* pHGeantEvent = (HRecEvent*)(gHades->getCurrentEvent()); // get pointer to event HPartialEvent* pSimul = pHGeantEvent->getPartialEvent(catSimul); HCategory* pRichCat1 = pSimul->getCategory(catRichGeantRaw); // get cat 1 if(*sector == 1) pRichCat1->Clear(); // clear category if 1st sector HCategory* pKin = pSimul->getCategory(catGeantKine); // get KINE cat //---------------------------------------------------------------------------------- if(RICHFLAG.isPMT==1){ // NEW RICH Int_t pmtId = *sector; Int_t sectorId = 0; // TODO set proper sector Id (*sector - 1); // Fortran index starts at 1! Int_t index; for(int trk=0; trkgetNewSlot(loc,&index); // get slot and its index if(pRichPhoton==NULL) { return; } //---------------------------------------------------------------------------------- // In PMT geometry there are no sectors any more ... lets take it for phi Float_t phi = TMath::ACos(RICHTUPLE.xcer[trk]/ sqrt(RICHTUPLE.xcer[trk]*RICHTUPLE.xcer[trk]+ RICHTUPLE.ycer[trk]*RICHTUPLE.ycer[trk])); if (RICHTUPLE.ycer[trk]<0) phi = 2.*TMath::Pi()-phi; Int_t sectorPhi = (Int_t)(phi/1.0471975)-1; // get sector from phi angle if (sectorPhi == -1) sectorPhi = 5; // 1st sector is along y axis! //---------------------------------------------------------------------------------- pRichPhoton = new(pRichPhoton) HGeantRichPhoton; // RICH photon data Int_t track = RICHTUPLE.parentcer[trk]; pRichPhoton->setTrack(track); pRichPhoton->setHit(RICHTUPLE.xcer[trk],RICHTUPLE.ycer[trk],RICHTUPLE.ecer[trk]); pRichPhoton->setAddress((Char_t)sectorPhi); // TODO: set sector number pRichPhoton->setPmtId(pmtId); if(pKin != NULL && track > 0) { // add hit to linked list of track ((HGeantKine*)pKin->getObject(track-1))->setRichHitIndex(index); pRichPhoton->sortVariable = RICHTUPLE.ecer[trk]; } } HCategory* pRichCat2 = pSimul->getCategory(catRichGeantRaw+1); // get cat 2 if(*sector == 1) pRichCat2->Clear(); // clear category if 1st sector for(int trk=0; trkgetNewSlot(loc); // get slot if(pRichDirect==NULL) { printf("\n*** FILLRICH: HGeantRichDirect location not found ! ***\n"); printf("sec=%d trk=%d\n",(*sector-1),trk); return; } //---------------------------------------------------------------------------------- // In PMT geometry there are no sectors any more ... lets take it for phi Float_t phi = TMath::ACos(RICHTUPLE.xhit[trk]/ sqrt(RICHTUPLE.xhit[trk]*RICHTUPLE.xhit[trk]+ RICHTUPLE.yhit[trk]*RICHTUPLE.yhit[trk])); if (RICHTUPLE.yhit[trk]<0) phi = 2.*TMath::Pi()-phi; Int_t sectorPhi = (Int_t)(phi/1.0471975)-1; // get sector from phi angle if (sectorPhi == -1) sectorPhi = 5; // 1st sector is along y axis! //---------------------------------------------------------------------------------- pRichDirect = new(pRichDirect) HGeantRichDirect; // RICH direct hit data pRichDirect->setTrack(RICHTUPLE.trkhit[trk],RICHTUPLE.parthit[trk]); pRichDirect->setHit(RICHTUPLE.xhit[trk],RICHTUPLE.yhit[trk],RICHTUPLE.zhit[trk]); pRichDirect->setMomentum(RICHTUPLE.momhit[trk],RICHTUPLE.thetahit[trk],RICHTUPLE.phihit[trk]); pRichDirect->setELoss(RICHTUPLE.elosshit[trk],RICHTUPLE.tlenhit[trk]); pRichDirect->setAddress((Char_t)sectorPhi); // TODO: set sector number pRichDirect->setPmtId(pmtId); } } else { // OLD RICH Int_t index; for(int trk=0; trkgetSlot(loc,&index); // get slot and its index if(pRichPhoton==NULL) { return; } pRichPhoton = new(pRichPhoton) HGeantRichPhoton; // RICH photon data Int_t track = RICHTUPLE.parentcer[trk]; pRichPhoton->setTrack(track); pRichPhoton->setHit(RICHTUPLE.xcer[trk],RICHTUPLE.ycer[trk], RICHTUPLE.ecer[trk]); pRichPhoton->setAddress((Char_t)(*sector-1)); if(pKin != NULL && track > 0) { // add hit to linked list of track ((HGeantKine*)pKin->getObject(track-1))->setRichHitIndex(index); pRichPhoton->sortVariable = RICHTUPLE.ecer[trk]; } } HCategory* pRichCat2 = pSimul->getCategory(catRichGeantRaw+1); // get cat 2 if(*sector == 1) pRichCat2->Clear(); // clear category if 1st sector for(int trk=0; trkgetSlot(loc); // get slot if(pRichDirect==NULL) { printf("\n*** FILLRICH: HGeantRichDirect location not found ! ***\n"); printf("sec=%d trk=%d\n",(*sector-1),trk); return; } pRichDirect = new(pRichDirect) HGeantRichDirect; // RICH direct hit data pRichDirect->setTrack(RICHTUPLE.trkhit[trk],RICHTUPLE.parthit[trk]); pRichDirect->setHit(RICHTUPLE.xhit[trk],RICHTUPLE.yhit[trk], RICHTUPLE.zhit[trk]); pRichDirect->setMomentum(RICHTUPLE.momhit[trk],RICHTUPLE.thetahit[trk], RICHTUPLE.phihit[trk]); pRichDirect->setELoss(RICHTUPLE.elosshit[trk],RICHTUPLE.tlenhit[trk]); pRichDirect->setAddress((Char_t)(*sector-1)); } } // end OLD RICH } void fillrichmirror(int* zombie) { HLocation loc; HRecEvent* pHGeantEvent = (HRecEvent*)(gHades->getCurrentEvent()); // get pointer to event HPartialEvent* pSimul = pHGeantEvent->getPartialEvent(catSimul); HCategory* pRichCat3 = pSimul->getCategory(catRichGeantRaw+2); // get cat 3 pRichCat3->Clear(); // clear category // HCategory* pKin = pSimul->getCategory(catGeantKine); // get KINE cat for(int trk=0; trk0) { Float_t phi = TMath::ACos(RICHTUPLE.xlep[trk]/ sqrt(RICHTUPLE.xlep[trk]*RICHTUPLE.xlep[trk]+ RICHTUPLE.ylep[trk]*RICHTUPLE.ylep[trk])); if (RICHTUPLE.ylep[trk]<0) phi = 2.*TMath::Pi()-phi; Int_t sector = (Int_t)(phi/1.0471975)-1; // get sector from phi angle if (sector == -1) sector = 5; // 1st sector is along y axis! loc.set(2,sector,trk); HGeantRichMirror* pRichMirror = (HGeantRichMirror*)pRichCat3->getSlot(loc); // get slot if (pRichMirror==NULL) { printf("\n*** FILLRICH: HGeantRichMirror location not found ! ***\n"); printf("sec=%d trk=%d\n",sector,trk); return; } pRichMirror = new(pRichMirror) HGeantRichMirror; // RICH Mirror hit data pRichMirror->setTrack(RICHTUPLE.leptrack[trk],RICHTUPLE.lepid[trk]); pRichMirror->setHit(RICHTUPLE.xlep[trk],RICHTUPLE.ylep[trk], RICHTUPLE.zlep[trk]); pRichMirror->setXYring(RICHTUPLE.xring[trk],RICHTUPLE.yring[trk]); pRichMirror->setAddress((Char_t)sector); pRichMirror->setNumPhot(RICHTUPLE.numphot[trk]); } } } // wrapping of the c++ routine as Fortran-subroutine using cfortran.h // of Cern-library FCALLSCSUB1(fillrich,FILLRICH,fillrich,PINT) FCALLSCSUB1(fillrichmirror,FILLRICHMIRROR,fillrichmirror,PINT)