BmnZdc.cxx

Go to the documentation of this file.

 
/*****************************************************************************
 *
 *         Class BmnZdc
 *
 *  Adopted for BMN by:   Elena Litvinenko (EL)
 *  e-mail:   litvin@nf.jinr.ru
 *  Version:  06-11-2015
 *  Last update:  22-Feb-2012 (EL)
 *  Last update:  05-May-2016 (MG)
 *
 *****************************************************************************/
 
#include "BmnZdc.h"
 
#include "BmnZdcGeo.h"
#include "BmnZdcPoint.h"
#include "CbmStack.h"
#include "FairGeoInterface.h"
#include "FairGeoLoader.h"
#include "FairGeoNode.h"
#include "FairGeoRootBuilder.h"
#include "FairRootManager.h"
#include "FairVolume.h"
#include "TClonesArray.h"
#include "TGeoArb8.h"
#include "TGeoMCGeometry.h"
#include "TGeoManager.h"
#include "TLorentzVector.h"
#include "TParticle.h"
#include "TVirtualMC.h"
 
#include <iostream>
// add on for debug
// #include "FairGeoG3Builder.h"
#include "FairRun.h"
#include "FairRuntimeDb.h"
#include "TObjArray.h"
#include "TParticlePDG.h"
 
// -----   Default constructor   -------------------------------------------
BmnZdc::BmnZdc()
{
    fZdcCollection = new TClonesArray("BmnZdcPoint");
    volDetector = 0;
    fPosIndex = 0;
    // fpreflag = 0;
    // fpostflag = 0;
    fVSCVolId = 0;
    fVSCVolId_after = 0;
    fEventID = -1;
    fVerboseLevel = 1;
}
// -------------------------------------------------------------------------
 
// -----   Standard constructor   ------------------------------------------
BmnZdc::BmnZdc(const char* name, Bool_t active)
    : FairDetector(name, active)
{
    fZdcCollection = new TClonesArray("BmnZdcPoint");
    fPosIndex = 0;
    volDetector = 0;
    // fpreflag = 0;
    // fpostflag = 0;
    fVSCVolId = 0;
    fVSCVolId_after = 0;
    fEventID = -1;
    fVerboseLevel = 1;
}
// -------------------------------------------------------------------------
 
// -----   Destructor   ----------------------------------------------------
BmnZdc::~BmnZdc()
{
    if (fZdcCollection) {
        fZdcCollection->Delete();
        delete fZdcCollection;
    }
}
// -------------------------------------------------------------------------
 
// -----   Public method Intialize   ---------------------------------------
void BmnZdc::Initialize()
{
    // Init function
    // cout <<"Initialize marina " <<endl;
    FairDetector::Initialize();
    // FairRun* sim = FairRun::Instance();
    // FairRuntimeDb* rtdb=sim->GetRuntimeDb();
 
    fVSCVolId = gMC->VolId("zdc01s");
    fVSCVolId_after = gMC->VolId("zdc01s_after");
    // fVSCVETOVolId = gMC->VolId("ndet01s_VETO");
}
// -------------------------------------------------------------------------
void BmnZdc::BeginEvent()
{
    // Begin of the event
    // cout <<"BeginEvent marina " <<endl;
}
 
// -----   Public method ProcessHits  --------------------------------------
Bool_t BmnZdc::ProcessHits(FairVolume* vol)
{
 
    Int_t copyNoVSC /*,copyNoVTYVEC,copyNoVMOD,copyNoVZDC*/;
    // Int_t copyNoVSCNICA,copyNoVTYVECNICA,copyNoVMODNICA,copyNoVZDCNICA;
    Int_t copyNoVSC_after /*,copyNoVTYVECNICA,copyNoVMODNICA,copyNoVZDCNICA*/;
    // Int_t copyNoVSCVETO,copyNoVTYVECVETO,copyNoVMODVETO,copyNoVZDCVETO;
    // Int_t copyNoVSCCom,copyNoVTYVECCom,copyNoVMODCom,copyNoVZDCCom;
 
    Int_t ivol;
    TLorentzVector tPos1, tMom1;
 
    // Int_t module, module_nica;
    // Int_t slice, slice_nica;
 
    TLorentzVector tPos, tMom;
    Int_t copyNo;
    Int_t copyNo_slice;
 
    Double_t time = 0;
    Double_t length = 0;
 
    // TParticle* part;
    // Double_t charge;
 
    // Double_t timeCut = 45;
    // Double_t timeCut = 55;
    // Double_t  QCF=1; //quenching for Birk
    // Double_t  BirkConst = 12.6; //0.126 mm/MeV for polystyrene
    // 0.126 *(0.1/0.001) = 12.6 cm/GeV
    //(0.126 mm/MeV - from Wikipedia, 0.07943mm/MeV є– in Geant4)
 
    if (gMC->CurrentVolID(copyNoVSC) != fVSCVolId && gMC->CurrentVolID(copyNoVSC_after) != fVSCVolId_after) {
        return kFALSE;
    }
 
    ivol = vol->getMCid();
    fEventID = gMC->CurrentEvent();
 
    if (gMC->CurrentVolID(copyNoVSC) == fVSCVolId || gMC->CurrentVolID(copyNoVSC_after) == fVSCVolId_after) {
        // gMC->CurrentVolOffID(0, slice); //no tyvec
        // gMC->CurrentVolOffID(1, module);
        gMC->CurrentVolOffID(0, copyNo);   // planeNb 1-6
        gMC->CurrentVolOffID(1, copyNo_slice);
        // copyNoVTYVECCom = slice; copyNoVMODCom = module;
    }
 
    if (gMC->IsTrackEntering()) {
 
        ResetParameters();
        fELoss = 0.;
        time = 0.;
        length = 0.;
        gMC->TrackPosition(tPos);
        gMC->TrackMomentum(tMom);
 
        fTrackID = gMC->GetStack()->GetCurrentTrackNumber();
        // part    = gMC->GetStack()->GetCurrentTrack();
 
    }   // if (gMC->IsTrackEntering())
 
    if (gMC->IsTrackExiting()) {
 
        gMC->TrackPosition(tPos);
        gMC->TrackMomentum(tMom);
        // part    = gMC->GetStack()->GetCurrentTrack();
        // fTrackID = gMC->GetStack()->GetCurrentTrackNumber();
 
        // Create BmnNdetPoint
        fTrackID = gMC->GetStack()->GetCurrentTrackNumber();
        time = gMC->TrackTime() * 1.0e09;
        // length  += gMC->TrackLength();
 
        // if(fELoss>0) {
 
        length = 1;
 
        // cout <<"Exit_WALL " <<length <<" " <<gMC->CurrentVolPath() <<" " <<copyNo_slice <<" " <<copyNo <<" "
        // <<tPos.X() <<" " <<tPos.Y() <<" " <<tPos.Z() <<" " <<tMom.Px() <<" " <<tMom.Py() <<" " <<tMom.Pz() <<" "
        // <<time <<" " <<fELoss  <<" " <<part->GetPDG()->PdgCode() <<" " <<part->GetMother(0) <<endl;
        AddHit(fTrackID, ivol, copyNo_slice, copyNo, TVector3(tPos.X(), tPos.Y(), tPos.Z()),
               TVector3(tMom.Px(), tMom.Py(), tMom.Pz()), time, length, fELoss);
        /*
        Int_t points = gMC->GetStack()->GetCurrentTrack()->GetMother(1);
  //       Int_t nZdcPoints = (points & (1<<30)) >> 30;
  //       nZdcPoints ++;
  //       if (nZdcPoints > 1) nZdcPoints = 1;
  //      points = ( points & ( ~ (1<<30) ) ) | (nZdcPoints << 30);
        points = ( points & ( ~ (1<<30) ) ) | (1 << 30);
        gMC->GetStack()->GetCurrentTrack()->SetMother(1,points);
        */
        ((CbmStack*)gMC->GetStack())->AddPoint(kZDC);
 
        //}//if(fELoss>0)
    }   // if ( gMC->IsTrackExiting()) {
    /*
       Int_t points = gMC->GetStack()->GetCurrentTrack()->GetMother(1);
//       Int_t nZdcPoints = (points & (1<<30)) >> 30;
//       nZdcPoints ++;
//       if (nZdcPoints > 1) nZdcPoints = 1;
//      points = ( points & ( ~ (1<<30) ) ) | (nZdcPoints << 30);
      points = ( points & ( ~ (1<<30) ) ) | (1 << 30);
      gMC->GetStack()->GetCurrentTrack()->SetMother(1,points);
 
      ((CbmStack*)gMC->GetStack())->AddPoint(kZDC);
    */
 
    //////////////////////////////////////////////////
    // cout <<"ProcessHits marina " <<endl;
    //  if (TMath::Abs(gMC->TrackCharge()) <= 0) return kFALSE;
 
    // Int_t      ivol    = vol->getMCid();
    // TLorentzVector tPos1, tMom1;
 
    // commented for vacuum FWall
    /*
    #ifndef EDEBUG
          if (fELoss == 0. ) return kFALSE;
    #else
          if ((fELoss == 0. ) &&
          (!((gMC->GetStack()->GetCurrentTrack()->GetPdgCode()==2112)&&(gMC->GetStack()->GetCurrentTrack()->GetMother(0)==-1)))
    ) return kFALSE;
    #endif
    */
    // END comment for vacuum FWall
    /*
        if (gMC->IsTrackEntering()) {
 
          ResetParameters();
          fELoss = 0.;
        }
         if ( gMC->IsTrackExiting()    ||
         gMC->IsTrackStop()       ||
         gMC->IsTrackDisappeared()   ) {
    */
 
    /*
          //old geom
          Int_t copyNo;
          Int_t ivol1 = gMC->CurrentVolID(copyNo);
          //      ivol1 = vol->getVolumeId();
          Int_t iCell, iCell2 ;
          gMC->CurrentVolOffID(1, iCell);
          gMC->CurrentVolOffID(2, iCell2);
     */
    /*
          //new geom
          Int_t copyNo;
          Int_t copyNo_slice;
          //Int_t ivol1 = gMC->CurrentVolID(copyNo_slice);
          //      ivol1 = vol->getVolumeId();
          Int_t iCell, iCell2 ;
          gMC->CurrentVolOffID(1, copyNo_slice);
          gMC->CurrentVolOffID(0, copyNo);
 
          TParticle* part    = gMC->GetStack()->GetCurrentTrack();
          Double_t charge = part->GetPDG()->Charge() / 3. ;
 
     // Create BmnZdcPoint
         fTrackID = gMC->GetStack()->GetCurrentTrackNumber();
          Double_t time    = gMC->TrackTime() * 1.0e09;
          Double_t length  = gMC->TrackLength();
          TLorentzVector tPos, tMom;
          //gMC->TrackPosition(tPos);
          //gMC->TrackMomentum(tMom);
          Double_t ZZZ;
        Double_t eLoss   = gMC->Edep();
        if (eLoss != 0.)
          fELoss += eLoss;
 
        //if (gMC->IsTrackEntering()) {
 
        //ResetParameters();
        //fELoss = 0.;
        length = 0; //entering
        fTrackID = gMC->GetStack()->GetCurrentTrackNumber();
        time    = gMC->TrackTime() * 1.0e09;
        gMC->TrackPosition(tPos);
        gMC->TrackMomentum(tMom);
 
          //cout <<"Enter " <<gMC->CurrentVolPath() <<" " <<ivol1 <<" " <<copyNo_slice <<" "
     <<gMC->CurrentVolID(copyNo_slice) <<" " <<gMC->CurrentVolOffID(0, iCell) <<" " <<iCell <<" "
     <<gMC->CurrentVolOffID(1, copyNo) <<" " <<copyNo <<endl;
          //if(copyNo==1 && gMC->IsTrackEntering()) cout <<"Enter_WALL " <<length <<" " <<gMC->CurrentVolPath() <<" "
     <<copyNo_slice <<" " <<copyNo <<" " <<tPos.X() <<" " <<tPos.Y() <<" " <<tPos.Z() <<" " <<tMom.Px() <<" "
     <<tMom.Py() <<" " <<tMom.Pz() <<" " <<time <<" " <<fELoss <<endl; cout <<"Exit_WALL " <<length <<" "
     <<gMC->CurrentVolPath() <<" " <<copyNo_slice <<" " <<copyNo <<" " <<tPos.X() <<" " <<tPos.Y() <<" " <<tPos.Z() <<"
     " <<tMom.Px() <<" " <<tMom.Py() <<" " <<tMom.Pz() <<" " <<time <<" " <<fELoss  <<endl; AddHit(fTrackID, ivol,
     copyNo_slice, copyNo, TVector3(tPos.X(), tPos.Y(), tPos.Z()), TVector3(tMom.Px(), tMom.Py(), tMom.Pz()), time,
     length, fELoss);
        }//if (gMC->IsTrackEntering())
    */
 
    /*
    if ( gMC->IsTrackExiting()    ||
     gMC->IsTrackStop()       ||
     gMC->IsTrackDisappeared()   ) {
    */
 
    /*
    if ( gMC->IsTrackExiting() ) {
 
      //cout <<"Exit " <<gMC->CurrentVolPath() <<" " <<ivol1 <<" " <<copyNo_slice <<" "
<<gMC->CurrentVolID(copyNo_slice) <<" " <<iCell <<" " <<gMC->CurrentVolOffID(0, iCell) <<" " <<copyNo <<" "
<<gMC->CurrentVolOffID(1, copyNo) <<endl;
 
    part    = gMC->GetStack()->GetCurrentTrack();
 
// Create BmnZdcPoint
    fTrackID = gMC->GetStack()->GetCurrentTrackNumber();
    time    = gMC->TrackTime() * 1.0e09;
    //length  += gMC->TrackLength();
    length = 1; //exiting
    */
 
    /*
#ifdef EDEBUG
    static Bool_t already=0;
    if (lEDEBUGcounter<200) {
  std::cout << "EDEBUG-- BmnZdc::ProcessHits: TrackID:" << fTrackID <<
    //    " ELoss: " << fELoss <<
    //        "   particle: " << (part->GetName()) <<
    "   " << gMC->CurrentVolPath() << " " << tPos.X()<< " " << tPos.Y() << " " << tPos.Z() << "  " <<
part->GetPDG()->PdgCode()<<
    //          "   " << (gMC->GetStack()->GetCurrentTrack()->GetMother(1)) <<
    //        "   "  << ivol << "=="<< gMC->CurrentVolID(copyNo) << ","<< copyNo <<
    "   "  << ivol << "  " << vol->getRealName() << "   " <<
    "   "  << gMC->CurrentVolOffID(1,iCell) << " " <<  gMC->CurrentVolOffName(1) <<
    "   " << iCell << " " << iCell2 << "  "<< gMC->CurrentVolOffName(2) <<
    "   " << vol->getModId()<< " " <<
    //    " " <<  gMC->CurrentVolOffName(2) << " " <<  gMC->CurrentVolOffName(1) << " " << gMC->CurrentVolOffName(0) <<
" "
        // "   " << vol->getRealName() << "  " << gMC->CurrentVolPath() <<
        // "   ivol,iCell,copyNo= " << ivol << ","<< iCell << ","<< copyNo <<
        // "   " << vol->getRealName() << "  "<< gMC->CurrentVolName() << "  "<< gMC->CurrentVolPath() <<
        // "   "  << ivol << ","<< vol->getVolumeId() << " : "<< gMC->CurrentVolID(copyNo) << ","<< copyNo <<
            //  "  "<< gMC->CurrentVolOffName(2) << "  " <<
    //             "  "<< gMC->CurrentVolOffName(2) << "  "<< gMC->CurrentVolOffName(3) <<
    std::endl;
  //    vol->getGeoNode()->getMotherNode()->getMotherNode()->getCenterPosition().print();
  lEDEBUGcounter++;
    }
    if ((iCell==2)&&(lEDEBUGcounter>=100)&&(!already)) {
  already=1;
  lEDEBUGcounter=0;
    }
//       if ((part->GetPdgCode())==321) {
//  std::cout << "EDEBUG-- BmnZdc::ProcessHits(..)  K+:  " << fTrackID << "   " << (  gMC->IsTrackExiting()) << "  " <<
//    (gMC->IsTrackStop()) << "  " << (gMC->IsTrackDisappeared()) << "   " << fELoss << "  " << time << std::endl;
//       }
//#endif
 
//       if(copyNo==1)
//  AddHit(fTrackID, ivol, copyNo, iCell, TVector3(tPos1.X(), tPos1.Y(), tPos1.Z()),
//         TVector3(tMom1.Px(), tMom1.Py(), tMom1.Pz()),
//         time, length, fELoss);
//       else
 
  AddHit(fTrackID, ivol, copyNo_slice, copyNo, TVector3(tPos.X(), tPos.Y(), tPos.Z()),
         TVector3(tMom.Px(), tMom.Py(), tMom.Pz()),
         time, length, fELoss);
#else
*/
    // if(copyNo<=2) cout <<"Exit_WALL " <<length <<" " <<gMC->CurrentVolPath() <<" " <<copyNo_slice <<" " <<copyNo <<"
    // " <<tPos.Z() <<endl;
 
    /*
      AddHit(fTrackID, ivol, copyNo_slice, copyNo, TVector3(tPos.X(), tPos.Y(), tPos.Z()),
         TVector3(tMom.Px(), tMom.Py(), tMom.Pz()),
         time, length, fELoss);
//#endif
  }//if ( gMC->IsTrackExiting()) {
    */
    /*
      Int_t points = gMC->GetStack()->GetCurrentTrack()->GetMother(1);
//       Int_t nZdcPoints = (points & (1<<30)) >> 30;
//       nZdcPoints ++;
//       if (nZdcPoints > 1) nZdcPoints = 1;
//      points = ( points & ( ~ (1<<30) ) ) | (nZdcPoints << 30);
      points = ( points & ( ~ (1<<30) ) ) | (1 << 30);
      gMC->GetStack()->GetCurrentTrack()->SetMother(1,points);
 
      ((CbmStack*)gMC->GetStack())->AddPoint(kZDC);
    */
    //}
 
    return kTRUE;
 
    //  }
    // #undef EDEBUG
}
 
// ----------------------------------------------------------------------------
 
// -----   Public method EndOfEvent   -----------------------------------------
void BmnZdc::EndOfEvent()
{
    if (fVerboseLevel)
        Print("");
    Reset();
}
 
// -----   Public method Register   -------------------------------------------
void BmnZdc::Register()
{
    FairRootManager::Instance()->Register("ZdcPoint", "Zdc", fZdcCollection, kTRUE);
}
 
// -----   Public method GetCollection   --------------------------------------
TClonesArray* BmnZdc::GetCollection(Int_t iColl) const
{
    if (iColl == 0)
        return fZdcCollection;
 
    return NULL;
}
 
// -----   Public method Print   ----------------------------------------------
void BmnZdc::Print(Option_t*) const
{
    Int_t nHits = fZdcCollection->GetEntriesFast();
    cout << "-I- BmnZdc: " << nHits << " points registered in this event." << endl;
 
    if (fVerboseLevel > 1)
        for (Int_t i = 0; i < nHits; i++)
            (*fZdcCollection)[i]->Print();
}
 
// -----   Public method Reset   ----------------------------------------------
void BmnZdc::Reset()
{
    fZdcCollection->Delete();
 
    fPosIndex = 0;
}
 
// guarda in FairRootManager::CopyClones
// -----   Public method CopyClones   -----------------------------------------
void BmnZdc::CopyClones(TClonesArray* cl1, TClonesArray* cl2, Int_t offset)
{
    Int_t nEntries = cl1->GetEntriesFast();
    // cout << "-I- BmnZdc: " << nEntries << " entries to add." << endl;
    TClonesArray& clref = *cl2;
    BmnZdcPoint* oldpoint = NULL;
    for (Int_t i = 0; i < nEntries; i++) {
        oldpoint = (BmnZdcPoint*)cl1->At(i);
        Int_t index = oldpoint->GetTrackID() + offset;
        oldpoint->SetTrackID(index);
        new (clref[fPosIndex]) BmnZdcPoint(*oldpoint);
        fPosIndex++;
    }
    cout << " -I- BmnZdc: " << cl2->GetEntriesFast() << " merged entries." << endl;
}
 
// -----   Public method ConstructGeometry   ----------------------------------
 
void BmnZdc::ConstructGeometry()
{
 
    TString fileName = GetGeometryFileName();
    if (fileName.EndsWith(".root")) {
        // FairLogger::GetLogger()->Info(MESSAGE_ORIGIN, "Constructing ZDC geometry from ROOT file %s",
        // fileName.Data());
        LOG(info) << "Constructing ZDC geometry from ROOT file" << fileName.Data();
        ConstructRootGeometry();
    }
 
    FairGeoLoader* geoLoad = FairGeoLoader::Instance();
    FairGeoInterface* geoFace = geoLoad->getGeoInterface();
    BmnZdcGeo* zdcGeo = new BmnZdcGeo();
    zdcGeo->setGeomFile(GetGeometryFileName());
    geoFace->addGeoModule(zdcGeo);
 
    Bool_t rc = geoFace->readSet(zdcGeo);
    if (rc)
        zdcGeo->create(geoLoad->getGeoBuilder());
    TList* volList = zdcGeo->getListOfVolumes();
 
    // store geo parameter
    FairRun* fRun = FairRun::Instance();
    FairRuntimeDb* rtdb = FairRun::Instance()->GetRuntimeDb();
    BmnZdcGeoPar* par = (BmnZdcGeoPar*)(rtdb->getContainer("BmnZdcGeoPar"));
    TObjArray* fSensNodes = par->GetGeoSensitiveNodes();
    TObjArray* fPassNodes = par->GetGeoPassiveNodes();
 
    TListIter iter(volList);
    FairGeoNode* node = NULL;
    FairGeoVolume* aVol = NULL;
 
    while ((node = (FairGeoNode*)iter.Next())) {
        aVol = dynamic_cast<FairGeoVolume*>(node);
 
        if (node->isSensitive()) {
            fSensNodes->AddLast(aVol);
        } else {
            fPassNodes->AddLast(aVol);
        }
    }
    par->setChanged();
    par->setInputVersion(fRun->GetRunId(), 1);
 
    ProcessNodes(volList);
}
 
// Check if Sensitive-----------------------------------------------------------
Bool_t BmnZdc::CheckIfSensitive(std::string name)
{
    TString tsname = name;
    // if (tsname.Contains("zdc01s") || tsname.Contains("ScH")) {
    if (tsname.Contains("zdc01s") || tsname.Contains("zdc01s_after") || tsname.Contains("ScH")) {
        return kTRUE;
    }
    return kFALSE;
}
 
// -----   Private method AddHit   --------------------------------------------
// BmnZdcPoint* BmnZdc::AddHit(Int_t trackID, Int_t module_groupID, Int_t copyNo, Int_t copyNoMother,
BmnZdcPoint* BmnZdc::AddHit(Int_t trackID,
                            Int_t detID,
                            Int_t copyNo,
                            Int_t copyNoMother,
                            TVector3 pos,
                            TVector3 mom,
                            Double_t time,
                            Double_t length,
                            Double_t eLoss)
{
    TClonesArray& clref = *fZdcCollection;
    Int_t size = clref.GetEntriesFast();
    return new (clref[size]) BmnZdcPoint(trackID, detID, copyNo, copyNoMother, pos, mom, time, length, eLoss);
}