create_stsdigipar_GEMS.C

Go to the documentation of this file.

/******************************************************************************
 ** Creation of STS digitisation parameter file in ASCII format.
 **
 ** @file create_stsdigipar_v12.C
 ** @author Volker Friese <v.friese@gsi.de>
 ** @since 21 June 2012
 ** @date 21.06.2012
 ** @Modified by Alexander Zinchenko for BM@N GEM STS
 ** @date 21.12.2013
 **
 ** This macro is compatible with geometries created with create_geometry_v12.C
 **
 ** The digitisation properties are the same for all sensors.
 ** The numbering of modules (sectors) is consecutively in each station,
 ** in the order given by the geometry. This has to be the same in
 ** CbmStsDigiScheme::Init().
 *****************************************************************************/
 
#include "TString.h"
 
#include <cstdio>
#include <iostream>
 
using namespace std;
 
// -------------   Steering variables       -----------------------------------
 
// ---> Sensor digitisation type.
// ---> 1 = pixel sensor, 2 = double-sided strip sensor,
// ---> 3 = double-sided strip sensor with double metal layer
Int_t gkSensorType = 2;   // 3;
 
// ---> Size of guard ring. It defines an inactive border on each side of the
// ---> sensor. MCPoints in this area will not be processed.
Double_t gkGuardRingWidth = 0.13;
 
// ---> Strip pitch on front and back side
Double_t gkStripPitchFront = 0.0400;
Double_t gkStripPitchBack = 0.0400;
// Double_t gkStripPitchFront = 0.0058;
// Double_t gkStripPitchBack  = 0.0058;
 
// ---> Strip stereo angle on front and back side
Double_t gkStereoFront = 0.;
Double_t gkStereoBack = -15.;
 
// ---> Geometry and digitisation tags
TString gkGeoTag = "v16a";
TString gkDigiTag = "gem";
// ----------------------------------------------------------------------------
 
// -------------   Other global variables   -----------------------------------
// ---> TGeoManager (too lazy to write out 'Manager' all the time
TGeoManager* geoMan = NULL;   // will be set later
// ----------------------------------------------------------------------------
 
// ============================================================================
// ======                         Main function                           =====
// ============================================================================
void create_stsdigipar_GEMS(const char* geoTag = gkGeoTag, const char* digiTag = gkDigiTag)
{
 
    // ---> Open parameter file
    TString parFileName = "sts_";
    parFileName = parFileName + geoTag + "_" + digiTag + ".digi.par";
    FILE* parFile;
    parFile = fopen(parFileName.Data(), "w");
    if (!parFile) {
        cout << "-E- Could not open parameter " << parFileName << endl;
        exit(0);
    }
    fprintf(parFile, "%s\n", "[CbmStsDigiPar]");
 
    // ---> Read geometry from file
    TString geoPath = getenv("VMCWORKDIR");
    // TString geoFileName = geoPath + "/geometry/sts/sts_";
    // geoFileName = geoFileName + geoTag + ".geo.root";
    TString geoFileName = "data/geofile_full.root";
    Bool_t init = CreateGeoFromFile(geoFileName.Data());
    if (!init) {
        cout << "-E- Could not read STS geometry!" << endl;
        exit(0);
    }
    TGeoManager* geoMan = gGeoManager;
 
    // ---> Top node
    geoMan->CdTop();
    TGeoNode* cave = geoMan->GetCurrentNode();
    cout << "Top node: " << cave->GetName() << ", daughters: " << cave->GetNdaughters() << endl;
 
    // ---> GEM STS
    if (!geoMan->cd("GEMS_0")) {
        cout << "-E- Mother volume GEMS not found " << endl;
        exit(0);
    }
 
    TGeoNode* sts = geoMan->GetCurrentNode();
    if (!TString(sts->GetName()).Contains("GEMS")) {
        cout << "-E- Illegal GEMS node name " << sts->GetName() << endl;
        exit(0);
    }
    cout << "STS node: " << sts->GetName() << ", daughters: " << sts->GetNdaughters()
         << ", stations: " << GetNofDaughters(sts, "Station") << endl;
 
    // ---> Stations
    Int_t statNr = 0;   // AZ
    for (Int_t iCaveD = 0; iCaveD < sts->GetNdaughters(); iCaveD++) {
        geoMan->CdDown(iCaveD);
        TGeoNode* station = geoMan->GetCurrentNode();
        if (!TString(station->GetName()).Contains("Station")) {
            geoMan->CdUp();
            continue;
        }
        // AZ Int_t statNr = station->GetNumber();
        ++statNr;
        if (statNr > 4)
            gkStripPitchFront = gkStripPitchBack = 0.08;   // 800 um
        Int_t nModules = GetNofModules(station);
        cout << "Station Nr. " << statNr << ", node " << station->GetName() << ", modules: " << nModules << endl;
        fprintf(parFile, "%d   %d   %d\n", statNr, 0, nModules * 2);
        Int_t moduleNr = 0;   // Running module number
        TGeoMatrix* matrixSt = geoMan->GetCurrentMatrix();
        Double_t* rot = matrixSt->GetRotationMatrix();
        Double_t phiStat = -TMath::ASin(rot[1]) * TMath::RadToDeg();
        gkStereoBack *= -1;
 
        // ---> Modules
        for (Int_t iStatD = 0; iStatD < nModules; iStatD++) {
            geoMan->CdDown(iStatD);
 
            TGeoNode* module = geoMan->GetCurrentNode();
            if (!TString(module->GetName()).Contains("Module")) {
                geoMan->CdUp();
                continue;
            }
            // moduleNr++;
            // fprintf(parFile, "%4d %4d\n", moduleNr, GetNofDaughters(module, "Senso"));
 
            // ---> Sensors
            for (Int_t iModD = 0; iModD < module->GetNdaughters(); iModD++) {
                geoMan->CdDown(iModD);
                TGeoNode* sensor = geoMan->GetCurrentNode();
                // if ( ! TString(sensor->GetName()).Contains("Sensor") ) {
                if (!TString(sensor->GetName()).Contains("Senso")) {
                    geoMan->CdUp();
                    continue;
                }
 
                moduleNr++;
                fprintf(parFile, "%4d %4d\n", moduleNr, 1);
                // ---> Sensor number
                Int_t sensNr = 1;
                // Int_t sensNr = sensor->GetNumber();
 
                // ---> Get sensor position
                TGeoMatrix* matrix = geoMan->GetCurrentMatrix();
                Double_t* xTrans = matrix->GetTranslation();
 
                // ---> Get sensor angle.
                // Sensors in lower half of station will be read out to bottom.
                Int_t phi = 0.;
                if (xTrans[1] < 0.)
                    phi = 180;
 
                // ---> Get sensor size
                TGeoBBox* shape = (TGeoBBox*)sensor->GetVolume()->GetShape();
                Double_t dx = 2. * shape->GetDX();
                Double_t dy = 2. * shape->GetDY();
                Double_t dz = 2. * shape->GetDZ();
                // AZ
                Double_t* rot = matrix->GetRotationMatrix();
                // if (TMath::Abs(rot[1]) > 0.9) {
                if (iStatD % 2) {
                    // Rotation by 90 deg
                    // Double_t dy = 2. * shape->GetDX();
                    // Double_t dx = 2. * shape->GetDY();
                }
                // AZ
 
                // ---> Subtract guard rings
                // dx -= 2 * gkGuardRingWidth;
                // dy -= 2 * gkGuardRingWidth;
 
                // ---> Write sensor information
                fprintf(parFile,
                        //"%6d %6d %9.3f %9.3f %9.3f %8d %8.2f %8.2f %8.2f %8.4f %8.4f %8.2f %8.2f\n",
                        "%6d %6d %9.3f %9.3f %9.3f %8.2f %8.2f %8.2f %8.2f %8.4f %8.4f %8.2f %8.2f\n",
                        // AZ sensNr, gkSensorType, xTrans[0], xTrans[1], xTrans[2], phi,
                        sensNr, gkSensorType, xTrans[0], xTrans[1], xTrans[2], phiStat, dx, dy, dz, gkStripPitchFront,
                        gkStripPitchBack, gkStereoFront, gkStereoBack);
 
                geoMan->CdUp();   // back to module
            }   // sensor loop
            geoMan->CdUp();   // back to module
        }   // module loop
        geoMan->CdUp();   // back to sts
    }   // station loop
 
    delete geoMan;
    fclose(parFile);
    cout << "-I- Parameters successfully written to " << parFileName << endl;
}
// ============================================================================
// ======                   End of main function                          =====
// ============================================================================
 
/** ===========================================================================
 ** CreateGeoFromFile
 **
 ** Create geometry, read STS node from file and add it to geometry.
 **/
Bool_t CreateGeoFromFile(const char* fileName)
{
 
    // ---> Create TGeoManager and top node (cave)
    geoMan = new TGeoManager("CBM STS", "STS Geometry");
    TGeoVolume* cave = geoMan->MakeBox("Cave", NULL, 2000., 2000., 2000.);
    geoMan->SetTopVolume(cave);
 
    // ---> Open geometry file
    cout << "Reading geometry from file " << fileName << endl;
    TFile* geoFile = new TFile(fileName, "READ");
    if (!geoFile->IsOpen()) {
        cout << "-E- Could not open file!" << endl;
        return kFALSE;
    }
 
    // ---> Get STS node from file
    TList* keyList = geoFile->GetListOfKeys();
    TIter next(keyList);
    TKey* key = NULL;
    TGeoVolume* top = NULL;
    while ((key = (TKey*)next())) {
        TString className(key->GetClassName());
        if (className.BeginsWith("TGeoVolume")) {
            top = dynamic_cast<TGeoVolume*>(key->ReadObj());
            std::cout << "Found volume " << top->GetName() << endl;
            break;
        }
        // AZ
        if (className.BeginsWith("TGeoManager")) {
            // File with full geometry
            // geoMan = dynamic_cast<TGeoManager*> (key->ReadObj() );
            key->ReadObj();
            std::cout << "Found GeoManager " << gGeoManager->GetName() << endl;
            return kTRUE;
        }
        // AZ
    }
    if (!top) {
        cout << "-E- Could not find volume object in file" << endl;
        return kFALSE;
    }
    TGeoVolume* sts = top->GetNode(0)->GetVolume();
    TString stsName = sts->GetName();
    if (!stsName.BeginsWith("STS")) {
        cout << "-E- No STS volume found on second level" << endl;
        return kFALSE;
    }
    TString stsVersion = stsName;
    stsVersion.Remove(0, 4);
    cout << "Geometry version " << stsVersion << endl;
    top->GetNode(0)->InspectNode();
 
    // ---> Add STS node to cave
    TGeoMatrix* stsmatrix = top->GetNode(0)->GetMatrix();
    cave->AddNode(sts, 0, stsmatrix);
 
    // ---> Finish geometry and close geometry file
    geoMan->CloseGeometry();
    geoFile->Close();
 
    return kTRUE;
}
/** ======================================================================= **/
 
Int_t GetNofModules(TGeoNode* station)
{
 
    Int_t nModules = station->GetNdaughters();
 
    // --- Modules
    /*
    for (Int_t iModule = 0; iModule < station->GetNdaughters(); iModule++) {
      TGeoNode* module = station->GetDaughter(iModule);
      if ( TString(module->GetName()).Contains("Module") ) nModules++;
    }
    */
 
    return nModules;
}
 
Int_t GetNofDaughters(TGeoNode* node, const char* name)
{
    Int_t nDaughters = 0;
    for (Int_t iNode = 0; iNode < node->GetNdaughters(); iNode++)
        if (TString(node->GetDaughter(iNode)->GetName()).Contains(name))
            nDaughters++;
    return nDaughters;
}