From 312f398d5eb6eed3c0ba35b9789b16e41a7ac297 Mon Sep 17 00:00:00 2001 From: Oliver Lantwin Date: Tue, 5 Nov 2024 11:15:48 +0100 Subject: [PATCH 1/4] Remove last vestiges of muflux --- passive/CMakeLists.txt | 1 - passive/MufluxTargetStation.cxx | 289 --- passive/MufluxTargetStation.h | 79 - passive/PassiveLinkDef.h | 1 - python/MufluxDigi.py | 204 --- python/MufluxDigiReco.py | 2016 --------------------- python/MufluxMuonTaggerDrifttubesCombo.py | 83 - python/MufluxMuonTaggerPatRec.py | 241 --- python/MufluxPatRec.py | 508 ------ 9 files changed, 3422 deletions(-) delete mode 100644 passive/MufluxTargetStation.cxx delete mode 100644 passive/MufluxTargetStation.h delete mode 100644 python/MufluxDigi.py delete mode 100644 python/MufluxDigiReco.py delete mode 100644 python/MufluxMuonTaggerDrifttubesCombo.py delete mode 100644 python/MufluxMuonTaggerPatRec.py delete mode 100644 python/MufluxPatRec.py diff --git a/passive/CMakeLists.txt b/passive/CMakeLists.txt index 8592889c47..025f58df26 100644 --- a/passive/CMakeLists.txt +++ b/passive/CMakeLists.txt @@ -23,7 +23,6 @@ ShipGeoCave.cxx ShipMagnet.cxx ShipChamber.cxx ShipTargetStation.cxx -MufluxTargetStation.cxx ShipMuonShield.cxx ShipPassiveContFact.cxx ShipTAUMagneticSpectrometer.cxx diff --git a/passive/MufluxTargetStation.cxx b/passive/MufluxTargetStation.cxx deleted file mode 100644 index 2b49197a16..0000000000 --- a/passive/MufluxTargetStation.cxx +++ /dev/null @@ -1,289 +0,0 @@ -#include "MufluxTargetStation.h" - -#include "TGeoManager.h" -#include "FairRun.h" // for FairRun -#include "FairRuntimeDb.h" // for FairRuntimeDb -#include // for ostream -#include "TList.h" // for TListIter, TList (ptr only) -#include "TObjArray.h" // for TObjArray -#include "TString.h" // for TString -#include "TGeoBBox.h" -#include "TGeoCompositeShape.h" -#include "TGeoShapeAssembly.h" -#include "TGeoTube.h" -#include "TGeoMaterial.h" -#include "FairGeoMedia.h" -#include "FairGeoBuilder.h" -#include "TGeoMedium.h" -#include // for NULL -#include // for operator<<, basic_ostream, etc - -using std::cout; -using std::endl; - - -MufluxTargetStation::~MufluxTargetStation() -{ -} -MufluxTargetStation::MufluxTargetStation() - : FairModule("MufluxTargetStation", "") -{ -} - -MufluxTargetStation::MufluxTargetStation(const char* name, const Double_t tl,const Double_t al,const Double_t tz, - const Double_t az, const int nS, const Double_t sl, const char* Title ) - : FairModule(name ,Title) -{ - fTargetLength = tl; - fAbsorberLength = al; - fAbsorberZ = az; - fTargetZ = tz; - fnS = nS; - fsl = sl; -} - -MufluxTargetStation::MufluxTargetStation(const char* name, const Double_t tl,const Double_t tz, - const int nS, const Double_t sl, const char* Title ) - : FairModule(name ,Title) -{ - fTargetLength = tl; - fAbsorberLength = 0; - fAbsorberZ = 0; - fTargetZ = tz; - fnS = nS; - fsl = sl; -} - -// ----- Private method InitMedium -Int_t MufluxTargetStation::InitMedium(const char* name) -{ - static FairGeoLoader *geoLoad=FairGeoLoader::Instance(); - static FairGeoInterface *geoFace=geoLoad->getGeoInterface(); - static FairGeoMedia *media=geoFace->getMedia(); - static FairGeoBuilder *geoBuild=geoLoad->getGeoBuilder(); - - FairGeoMedium *ShipMedium=media->getMedium(name); - - if (!ShipMedium) - { - Fatal("InitMedium","Material %s not defined in media file.", name); - return -1111; - } - TGeoMedium* medium=gGeoManager->GetMedium(name); - if (medium!=NULL) - return ShipMedium->getMediumIndex(); - return geoBuild->createMedium(ShipMedium); -} - -void MufluxTargetStation::SetIronAbsorber(Double_t absorber_x, Double_t absorber_y) -{ - fabsorber_x=absorber_x; - fabsorber_y=absorber_y; - -} - - -void MufluxTargetStation::SetAbsorberCutout(Double_t absorbercutout_x, Double_t absorbercutout_y) -{ - fabsorbercutout_x=absorbercutout_x; - fabsorbercutout_y=absorbercutout_y; -} - -void MufluxTargetStation::SetIronShield(Double_t ironshield_x, Double_t ironshield_y, Double_t ironshield_z) -{ - fironshield_x=ironshield_x; - fironshield_y=ironshield_y; - fironshield_z=ironshield_z; -} - -void MufluxTargetStation::SetConcreteShield(Double_t concreteshield_x, Double_t concreteshield_y, Double_t concreteshield_z) -{ - fconcreteshield_x=concreteshield_x; - fconcreteshield_y=concreteshield_y; - fconcreteshield_z=concreteshield_z; -} - -void MufluxTargetStation::SetAboveTargetShield(Double_t abovetargetshield_x, Double_t abovetargetshield_y,Double_t abovetargetshield_z) -{ - fabovetargetshield_x=abovetargetshield_x; - fabovetargetshield_y=abovetargetshield_y; - fabovetargetshield_z=abovetargetshield_z; -} - -void MufluxTargetStation::SetAboveAbsorberShield(Double_t aboveabsorbershield_x, Double_t aboveabsorbershield_y, Double_t aboveabsorbershield_z) -{ - faboveabsorbershield_x=aboveabsorbershield_x; - faboveabsorbershield_y=aboveabsorbershield_y; - faboveabsorbershield_z=aboveabsorbershield_z; -} - -void MufluxTargetStation::SetAboveAboveTargetShield(Double_t aboveabovetargetshield_y) -{ - - faboveabovetargetshield_y=aboveabovetargetshield_y; - -} - -void MufluxTargetStation::SetFloor(Double_t floor_x, Double_t floor_y, Double_t floor_z) -{ - ffloor_x=floor_x; - ffloor_y=floor_y; - ffloor_z=floor_z; -} - -void MufluxTargetStation::SetFloorT34(Double_t floorT34_x, Double_t floorT34_y, Double_t floorT34_z) -{ - ffloorT34_x=floorT34_x; - ffloorT34_y=floorT34_y; - ffloorT34_z=floorT34_z; -} - -void MufluxTargetStation::SetFloorRPC(Double_t floorRPC_x, Double_t floorRPC_y, Double_t floorRPC_z) -{ - ffloorRPC_x=floorRPC_x; - ffloorRPC_y=floorRPC_y; - ffloorRPC_z=floorRPC_z; -} - - -void MufluxTargetStation::ConstructGeometry() -{ - TGeoVolume *top=gGeoManager->GetTopVolume(); - InitMedium("tungstenmisis"); - TGeoMedium *tungstenmisis =gGeoManager->GetMedium("tungstenmisis"); - InitMedium("Concrete"); - TGeoMedium *concrete =gGeoManager->GetMedium("Concrete"); - InitMedium("iron"); - TGeoMedium *iron =gGeoManager->GetMedium("iron"); - InitMedium("H2O"); - TGeoMedium *water =gGeoManager->GetMedium("H2O"); - InitMedium("molybdenummisis"); - TGeoMedium *momisis =gGeoManager->GetMedium("molybdenummisis"); - InitMedium("air"); - TGeoMedium *air = gGeoManager->GetMedium("air"); - InitMedium("PlasticBase"); - TGeoMedium *plasticbase = gGeoManager->GetMedium("PlasticBase"); - InitMedium("steel"); - TGeoMedium *steel = gGeoManager->GetMedium("steel"); - - TGeoVolume *tTarget = new TGeoVolumeAssembly("TargetArea"); - //Double_t zPos = -2.09; //for front end endcap - Double_t zPos = 0.; //for front end endcap - Int_t slots = fnS; - slots = slots-1; - - TGeoVolume *frontcap; - TGeoVolume *endcap; - TGeoVolume *targettube; - frontcap = gGeoManager->MakeTube("FrontCap", steel, 0, 15./2., 0.5/2.); - frontcap->SetLineColor(20); - endcap = gGeoManager->MakeTube("EndCap", steel, 0, 15./2., 0.5/2.); - endcap->SetLineColor(20); - targettube = gGeoManager->MakeTube("TargetTube", steel, fDiameter/2.+0.5, fDiameter/2.+1.1, fTargetLength/2.-0.5); //subtract 1 (front&endcaps) - targettube->SetLineColor(20); - - TGeoVolume *target; - TGeoVolume *slit; - - for (Int_t i=0; iMakeTube(nmi, material, 0., fDiameter/2., fL.at(i)/2.); - } - else { - target = gGeoManager->MakeBox(nmi, material, fDiameter/2., fDiameter/2., fL.at(i)/2.); - } - if (fM.at(i)=="molybdenummisis") { - target->SetLineColor(28); - } else {target->SetLineColor(38);}; // silver/blue - tTarget->AddNode(target, 1, new TGeoTranslation(0, 0, zPos + fL.at(i)/2.) ); - if (i < slots){ - if(fnS == 18) { - slit = gGeoManager->MakeTube(sm, plasticbase, 0., fDiameter/2., fsl/2.); - } - else { - slit = gGeoManager->MakeBox(sm, water, fDiameter/2., fDiameter/2., fsl/2.); - } - slit->SetLineColor(7); // cyan - if (i==0) { - tTarget->AddNode(frontcap, 11, new TGeoTranslation(0, 0, zPos-0.2501) ); - } - tTarget->AddNode(slit, 1, new TGeoTranslation(0, 0, zPos+fL.at(i)+fsl/2.) ); - zPos+=fL.at(i)+fsl; - } else { - zPos+=fL.at(i); - if (fnS==18) { - tTarget->AddNode(targettube, 13, new TGeoTranslation(0, 0, fTargetLength/2-0.5)); - tTarget->AddNode(endcap, 12, new TGeoTranslation(0, 0, fTargetLength-0.7499) ); - } - zPos+=1.27; - } - } - - TGeoVolume *absorber; - TGeoVolume *absorbercutout; - TGeoVolume *aboveabsorbershield; - TGeoVolume *ironshield; - TGeoVolume *concreteshield; - TGeoVolume *abovetargetshield; - TGeoVolume *aboveabovetargetshield; - TGeoVolume *floor; - TGeoVolume *floorT34; - TGeoVolume *floorRPC; - - zPos = fTargetZ - fTargetLength/2.; - // Absorber made of iron - absorber = gGeoManager->MakeBox("Absorber", iron, fabsorber_x, fabsorber_y, fAbsorberLength/2.); - absorbercutout = gGeoManager->MakeBox("AbsorberCO", iron, fabsorbercutout_x, fabsorbercutout_y, fAbsorberLength/2.); - - ironshield = gGeoManager->MakeBox("IronShield", iron, fironshield_x, fironshield_y, fironshield_z); - concreteshield = gGeoManager->MakeBox("ConcreteShield", concrete, fconcreteshield_x, fconcreteshield_y, fconcreteshield_z); - abovetargetshield = gGeoManager->MakeBox("AboveTargetShield", concrete, fabovetargetshield_x, fabovetargetshield_y, fabovetargetshield_z); - aboveabsorbershield = gGeoManager->MakeBox("AboveBsorberShield", concrete, faboveabsorbershield_x, faboveabsorbershield_y, faboveabsorbershield_z); - aboveabovetargetshield = gGeoManager->MakeBox("AboveAboveTargetShield", concrete, fAbsorberLength, faboveabovetargetshield_y, fAbsorberLength); - floor = gGeoManager->MakeBox("Floor", concrete, ffloor_x, ffloor_y, ffloor_z); - floorT34 = gGeoManager->MakeBox("FloorT34", concrete, ffloorT34_x, ffloorT34_y, ffloorT34_z); - floorRPC = gGeoManager->MakeBox("FloorRPC", concrete, ffloorRPC_x, ffloorRPC_y, ffloorRPC_z); - - ironshield->SetLineColor(42); // brown / light red - concreteshield->SetLineColor(kGray); // gray - floor->SetLineColor(kGray); // gray - floorT34->SetLineColor(kGray); // gray - floorRPC->SetLineColor(kGray); // gray - aboveabsorbershield->SetLineColor(kGray); // gray - abovetargetshield->SetLineColor(kGray); // gray - aboveabovetargetshield->SetLineColor(kGray); // gray - frontcap->SetLineColor(29); // blue-ish - endcap->SetLineColor(29); // blue-ish - absorber->SetLineColor(42); // brown / light red - - tTarget->AddNode(absorbercutout, 11, new TGeoTranslation(18., 97.5, fAbsorberZ-zPos)); - tTarget->AddNode(absorber, 1, new TGeoTranslation(0., -27.5, fAbsorberZ-zPos)); - - //TGeoNode *node = GetNode("Absorber"); - //tTarget->ReplaceNode(node,new TGeoTube(1.,2.,10.),new TGeoTranslation(0,-100.2,12),iron); // - - tTarget->AddNode(floor, 7, new TGeoTranslation(0., -217., 500.)); - tTarget->AddNode(floorT34, 8, new TGeoTranslation(0., -120., fAbsorberZ-zPos+771.125)); - tTarget->AddNode(floorRPC, 9, new TGeoTranslation(0., -94.49, fAbsorberZ-zPos+1075.)); - - TGeoShapeAssembly* asmb = dynamic_cast(tTarget->GetShape()); - Double_t totLength = asmb->GetDZ(); - - top->AddNode(tTarget, 1, new TGeoTranslation(0, 0,fTargetZ - fTargetLength/2. + totLength)); - top->AddNode(ironshield, 2, new TGeoTranslation(-50., -47.5,fTargetZ - 83.49)); //3.6 - top->AddNode(concreteshield, 3, new TGeoTranslation(85.,-47.5,fTargetZ - 83.49)); - top->AddNode(abovetargetshield, 4, new TGeoTranslation(50.,77.5,fTargetZ - 83.49)); - top->AddNode(aboveabovetargetshield, 5, new TGeoTranslation(-50.,165.0,fAbsorberZ-zPos-624.5)); - top->AddNode(aboveabsorbershield, 6, new TGeoTranslation(0.,165.0,-26.48)); - - cout << "target and absorber positioned at " << fTargetZ <<" "<< fAbsorberZ << " m"<< endl; - -} diff --git a/passive/MufluxTargetStation.h b/passive/MufluxTargetStation.h deleted file mode 100644 index 85b607754e..0000000000 --- a/passive/MufluxTargetStation.h +++ /dev/null @@ -1,79 +0,0 @@ -#ifndef MufluxTargetStation_H -#define MufluxTargetStation_H - -#include "FairModule.h" // for FairModule -#include "Rtypes.h" // for MufluxTargetStation::Class, Bool_t, etc -#include // for string -#include - -class MufluxTargetStation : public FairModule -{ - public: - MufluxTargetStation(const char* name,const Double_t tl,const Double_t al,const Double_t tz,const Double_t az, - const int nS, const Double_t sl, const char* Title="MufluxTargetStation"); - MufluxTargetStation(const char* name,const Double_t tl,const Double_t tz, - const int nS, const Double_t sl, const char* Title="MufluxTargetStation"); - MufluxTargetStation(); - virtual ~MufluxTargetStation(); - - void SetIronAbsorber(Double_t absorber_x, Double_t absorber_y); - void SetAbsorberCutout(Double_t absorbercutout_x ,Double_t absorbercutout_y); - void SetIronShield(Double_t ironshield_x, Double_t ironshield_y, Double_t ironshield_z); - void SetConcreteShield(Double_t concreteshield_x, Double_t concreteshield_y, Double_t concreteshield_z); - void SetAboveTargetShield(Double_t abovetargetshield_x, Double_t abovetargetshield_y, Double_t abovetargetshield_z); - void SetAboveAbsorberShield(Double_t aboveabsorbershield_x, Double_t aboveabsorbershield_y, Double_t aboveabsorbershield_z); - void SetAboveAboveTargetShield(Double_t aboveabovetargetshield_y); - void SetFloor(Double_t floor_x, Double_t floor_y, Double_t floor_z); - void SetFloorT34(Double_t floorT34_x, Double_t floorT34_y, Double_t floorT34_z); - void SetFloorRPC(Double_t floorRPC_x, Double_t floorRPC_y, Double_t floorRPC_z); - - void ConstructGeometry(); - void SetLayerPosMat(Float_t d, std::vector L, std::vector M){ - fDiameter = d; - fL = L; - fM = M; - } - ClassDef(MufluxTargetStation,4) - - private: - Double_t fabsorber_x; - Double_t fabsorber_y; - Double_t fabsorbercutout_x; - Double_t fabsorbercutout_y; - Double_t fironshield_x; - Double_t fironshield_y; - Double_t fironshield_z; - Double_t fconcreteshield_x; - Double_t fconcreteshield_y; - Double_t fconcreteshield_z; - Double_t fabovetargetshield_x; - Double_t fabovetargetshield_y; - Double_t fabovetargetshield_z; - Double_t faboveabsorbershield_x; - Double_t faboveabsorbershield_y; - Double_t faboveabsorbershield_z; - Double_t faboveabovetargetshield_y; - Double_t ffloor_x; - Double_t ffloor_y; - Double_t ffloor_z; - Double_t ffloorT34_x; - Double_t ffloorT34_y; - Double_t ffloorT34_z; - Double_t ffloorRPC_x; - Double_t ffloorRPC_y; - Double_t ffloorRPC_z; - - protected: - - Double_t fTargetLength; // - Double_t fAbsorberLength; // - Double_t fAbsorberZ; // - Double_t fTargetZ; // - Double_t fDiameter; // - std::vector fL; // - std::vector fM; // - Int_t InitMedium(const char* name); - Int_t fnS; - Double_t fsl; -}; -#endif //MufluxTargetStation_H diff --git a/passive/PassiveLinkDef.h b/passive/PassiveLinkDef.h index 222d8186e5..12729761ae 100644 --- a/passive/PassiveLinkDef.h +++ b/passive/PassiveLinkDef.h @@ -10,7 +10,6 @@ #pragma link C++ class ShipCave+; #pragma link C++ class ShipChamber+; #pragma link C++ class ShipTargetStation+; -#pragma link C++ class MufluxTargetStation+; #pragma link C++ class ShipMuonShield+; #pragma link C++ class ShipGeoCave; #pragma link C++ class ShipPassiveContFact; diff --git a/python/MufluxDigi.py b/python/MufluxDigi.py deleted file mode 100644 index db3a40231d..0000000000 --- a/python/MufluxDigi.py +++ /dev/null @@ -1,204 +0,0 @@ -from builtins import range -import os -import ROOT -import global_variables -import shipunit as u -from array import array - -stop = ROOT.TVector3() -start = ROOT.TVector3() - -deadChannelsForMC = [10112001, 20112003, 30002041, 30012026, 30102025, 30112013, 30112018, 40012014] -ineffiency = {1:0.1,2:0.1,3:0.1,4:0.1,5:0.1} - -# defining constants for rpc properties -STRIP_XWIDTH = 0.8625 # internal STRIP V, WIDTH, in cm -EXT_STRIP_XWIDTH_L = 0.9625 # nominal (R&L) and Left measured external STRIP V, WIDTH, in cm (beam along z, out from the V plane) -EXT_STRIP_XWIDTH_R = 0.86 # measured Right external STRIP V, WIDTH,in cm (beam along z, out from the V plane) -V_STRIP_OFF = 0.200 -NR_VER_STRIPS = 184 -total_width = (NR_VER_STRIPS - 2) * STRIP_XWIDTH + EXT_STRIP_XWIDTH_L + EXT_STRIP_XWIDTH_R + (NR_VER_STRIPS - 1) * V_STRIP_OFF -reduced_width = total_width - (EXT_STRIP_XWIDTH_L + EXT_STRIP_XWIDTH_R) - -# function for calculating the strip number from a coordinate, for MuonTagger / RPC -def StripX(x): - if x < -total_width/2. or x > total_width/2.: - print("WARNING: x coordinate outside sensitive volume!",x) - if x < -total_width/2. + EXT_STRIP_XWIDTH_L + V_STRIP_OFF/2. : strip_x = 184 - elif x > total_width/2. - EXT_STRIP_XWIDTH_R - V_STRIP_OFF/2. : strip_x = 1 - else: - x_start = x - total_width/2. + EXT_STRIP_XWIDTH_R - strip_x = -int(x_start/reduced_width*182.)+1 - if not (0 < strip_x <= NR_VER_STRIPS-1): - print("WARNING: X strip outside range!",x,strip_x) - strip_x = 0 - return int(strip_x) - -def StripY(y): - STRIP_YWIDTH = 0.8625 # internal STRIP H, WIDTH, in cm - EXT_STRIP_YWIDTH = 0.3 # measured external STRIP H, WIDTH, in cm - H_STRIP_OFF = 0.1983 - NR_HORI_STRIPS = 116 - total_height = (NR_HORI_STRIPS - 2) * STRIP_YWIDTH + 2 * EXT_STRIP_YWIDTH + (NR_HORI_STRIPS - 1) * H_STRIP_OFF - y_start = total_height / 2. - strip_y = (y_start - EXT_STRIP_YWIDTH + 1.5 * STRIP_YWIDTH + H_STRIP_OFF - y)//(STRIP_YWIDTH + H_STRIP_OFF) - if not (0 < strip_y <= NR_HORI_STRIPS): - print("WARNING: Y strip outside range!") - strip_y = 0 - return int(strip_y) - -class MufluxDigi: - " convert FairSHiP MC hits / digitized hits to measurements" - def __init__(self,fout): - - self.iEvent = 0 - - outdir=os.getcwd() - outfile=outdir+"/"+fout - self.fn = ROOT.TFile(fout,'update') - self.sTree = self.fn.cbmsim - - # event header - self.header = ROOT.FairEventHeader() - self.eventHeader = self.sTree.Branch("ShipEventHeader",self.header,32000,1) - self.digiMufluxSpectrometer = ROOT.TClonesArray("MufluxSpectrometerHit") - self.digiMufluxSpectrometerBranch = self.sTree.Branch("Digi_MufluxSpectrometerHits",self.digiMufluxSpectrometer,32000,1) - self.digiLateMufluxSpectrometer = ROOT.TClonesArray("MufluxSpectrometerHit") - self.digiLateMufluxSpectrometerBranch = self.sTree.Branch("Digi_LateMufluxSpectrometerHits",self.digiMufluxSpectrometer,32000,1) - #muon taggger - if self.sTree.GetBranch("MuonTaggerPoint"): - self.digiMuonTagger = ROOT.TClonesArray("MuonTaggerHit") - self.digiMuonTaggerBranch = self.sTree.Branch("Digi_MuonTaggerHits", self.digiMuonTagger, 32000, 1) - # setup random number generator - ROOT.gRandom.SetSeed() - self.PDG = ROOT.TDatabasePDG.Instance() - # for the digitizing and reconstruction step - self.v_drift = global_variables.ShipGeo['MufluxSpectrometer']['v_drift'] - self.sigma_spatial = global_variables.ShipGeo['MufluxSpectrometer']['sigma_spatial'] - self.viewangle = global_variables.ShipGeo['MufluxSpectrometer']['ViewvAngle'] - - self.gMan = ROOT.gGeoManager - - def digitize(self): - - self.sTree.t0 = ROOT.gRandom.Rndm()*1*u.microsecond - self.header.SetEventTime( self.sTree.t0 ) - self.header.SetRunId( self.sTree.MCEventHeader.GetRunID() ) - self.header.SetMCEntryNumber( self.sTree.MCEventHeader.GetEventID() ) # counts from 1 - self.eventHeader.Fill() - self.digiMufluxSpectrometer.Delete() - self.digiLateMufluxSpectrometer.Delete() - self.digitizeMufluxSpectrometer() - self.digiMufluxSpectrometerBranch.Fill() - self.digiLateMufluxSpectrometerBranch.Fill() - self.digiMuonTagger.Delete() # produces a lot of warnings, rpc station 0 - self.digitizeMuonTagger() - self.digiMuonTaggerBranch.Fill() - - def digitizeMuonTagger(self, fake_clustering=False): - - station = 0 - strip = 0 - nav = ROOT.gGeoManager.GetCurrentNavigator() - DetectorID = set() # set of detector ids - already deduplicated - for MuonTaggerHit in self.sTree.MuonTaggerPoint: - # getting rpc nodes, name and matrix - detID = MuonTaggerHit.GetDetectorID() - s = str(detID//10000) - nav.cd('/VMuonBox_1/VSensitive'+s+'_'+s) - # translation from top to MuonBox_1 - point = array('d', [ - MuonTaggerHit.GetX(), - MuonTaggerHit.GetY(), - MuonTaggerHit.GetZ(), 1]) - # translation to local frame - point_local = array('d', [0, 0, 0, 1]) - nav.MasterToLocal(point, point_local) - - xcoord = point_local[0] - ycoord = point_local[1] - - # identify individual rpcs - station = detID//10000 - if station not in range(1, 6): # limiting the range of rpcs - print("WARNING: Invalid RPC number, something's wrong with the geometry ",station) - - # calculate strip - # x gives vertical direction - direction = 1 - strip = StripX(xcoord) - if not strip: - continue - # sampling number of strips around the exact strip for emulating clustering - detectorid = station*10000 + direction*1000 + strip - DetectorID.add(detectorid) - if fake_clustering: - s = ROOT.gRandom.Poisson(2) - if ROOT.gRandom.Rndm() < 0.5: strip = strip - int(s/2) - else: strip = strip + int(s/2) - for i in range(0, s): - detectorid = station*10000 + direction*1000 + strip + i - DetectorID.add(detectorid) - - # y gives horizontal direction - direction = 0 - strip = StripY(ycoord) - if not strip: - continue - # sampling number of strips around the exact strip for emulating clustering - detectorid = station*10000 + direction*1000 + strip - DetectorID.add(detectorid) - if fake_clustering: - s = ROOT.gRandom.Poisson(2) - if ROOT.gRandom.Rndm() < 0.5: strip = strip - int(s/2) - else: strip = strip + int(s/2) - for i in range(0, s): - detectorid = station*10000 + direction*1000 + strip + i - DetectorID.add(detectorid) - - self.digiMuonTagger.Expand(len(DetectorID)) - for index, detID in enumerate(DetectorID): - hit = ROOT.MuonTaggerHit(detID, 0) - self.digiMuonTagger[index] = hit - - if fake_clustering: - # cluster size loop - plotting the cluster size distribution - cluster_size = list() - DetectorID_list = list(DetectorID) # turn set into list to allow indexing - DetectorID_list.sort() # sorting the list - if len(DetectorID_list) > 1: - clusters = [[DetectorID_list[0]]] - for x in DetectorID_list[1:]: - if abs(x - clusters[-1][-1]) <= 1: - clusters[-1].append(x) - else: - clusters.append([x]) - cluster_size = [len(x) for x in clusters] - - def digitizeMufluxSpectrometer(self): - - # digitize FairSHiP MC hits - index = 0 - hitsPerDetId = {} - - for aMCPoint in self.sTree.MufluxSpectrometerPoint: - aHit = ROOT.MufluxSpectrometerHit(aMCPoint,self.sTree.t0) - aHit.setValid() - if index>0 and self.digiMufluxSpectrometer.GetSize() == index: self.digiMufluxSpectrometer.Expand(index+1000) - self.digiMufluxSpectrometer[index]=aHit - detID = aHit.GetDetectorID() - if detID in hitsPerDetId: - if self.digiMufluxSpectrometer[hitsPerDetId[detID]].tdc() > aHit.tdc(): - # second hit with smaller tdc - self.digiMufluxSpectrometer[hitsPerDetId[detID]].setInvalid() - hitsPerDetId[detID] = index - else: - hitsPerDetId[detID] = index - if aMCPoint.GetDetectorID() in deadChannelsForMC: aHit.setInvalid() - station = int(aMCPoint.GetDetectorID()//10000000) - if ROOT.gRandom.Rndm() < ineffiency[station]: aHit.setInvalid() - index+=1 - - def finish(self): - print('finished writing tree') - self.sTree.Write() diff --git a/python/MufluxDigiReco.py b/python/MufluxDigiReco.py deleted file mode 100644 index 77eed944e1..0000000000 --- a/python/MufluxDigiReco.py +++ /dev/null @@ -1,2016 +0,0 @@ -from builtins import range -import global_variables -from global_variables import h -import os -import ROOT -import MufluxPatRec -import shipunit as u -import rootUtils as ut - -from array import array - -import matplotlib -matplotlib.use('pdf') -import matplotlib.pyplot as plt -import numpy as np - -stop = ROOT.TVector3() -start = ROOT.TVector3() - - -# function for calculating the strip number from a coordinate, for MuonTagger / RPC -def StripX(x): - # defining constants for rpc properties - STRIP_XWIDTH = 0.8625 # internal STRIP V, WIDTH, in cm - EXT_STRIP_XWIDTH_L = 0.9625 # nominal (R&L) and Left measured external STRIP V, WIDTH, in cm (beam along z, out from the V plane) - EXT_STRIP_XWIDTH_R = 0.86 # measured Right external STRIP V, WIDTH,in cm (beam along z, out from the V plane) - V_STRIP_OFF = 0.200 - NR_VER_STRIPS = 184 - total_width = (NR_VER_STRIPS - 2) * STRIP_XWIDTH + EXT_STRIP_XWIDTH_L + EXT_STRIP_XWIDTH_R + (NR_VER_STRIPS - 1) * V_STRIP_OFF - x_start = (total_width - EXT_STRIP_XWIDTH_R + EXT_STRIP_XWIDTH_L) / 2 - # calculating strip as an integer - strip_x = (x_start - EXT_STRIP_XWIDTH_L + 1.5 * STRIP_XWIDTH + V_STRIP_OFF - x)//(STRIP_XWIDTH + V_STRIP_OFF) - if not (0 < strip_x <= NR_VER_STRIPS): - print("WARNING: X strip outside range!") - strip_x = 0 - return int(strip_x) - - -def StripY(y): - STRIP_YWIDTH = 0.8625 # internal STRIP H, WIDTH, in cm - EXT_STRIP_YWIDTH = 0.3 # measured external STRIP H, WIDTH, in cm - H_STRIP_OFF = 0.1983 - NR_HORI_STRIPS = 116 - total_height = (NR_HORI_STRIPS - 2) * STRIP_YWIDTH + 2 * EXT_STRIP_YWIDTH + (NR_HORI_STRIPS - 1) * H_STRIP_OFF - y_start = total_height / 2 - strip_y = (y_start - EXT_STRIP_YWIDTH + 1.5 * STRIP_YWIDTH + H_STRIP_OFF - y)//(STRIP_YWIDTH + H_STRIP_OFF) - if not (0 < strip_y <= NR_HORI_STRIPS): - print("WARNING: Y strip outside range!") - strip_y = 0 - return int(strip_y) - -class MufluxDigiReco: - " convert FairSHiP MC hits / digitized hits to measurements" - def __init__(self,fout): - - self.iEvent = 0 - - outdir=os.getcwd() - outfile=outdir+"/"+fout - self.fn = ROOT.TFile(fout,'update') - self.sTree = self.fn.cbmsim - - if self.sTree.GetBranch("FitTracks"): - print("remove RECO branches and rerun reconstruction") - self.fn.Close() - # make a new file without reco branches - f = ROOT.TFile(fout) - sTree = f.cbmsim - if sTree.GetBranch("FitTracks"): sTree.SetBranchStatus("FitTracks",0) - if sTree.GetBranch("Particles"): sTree.SetBranchStatus("Particles",0) - if sTree.GetBranch("fitTrack2MC"): sTree.SetBranchStatus("fitTrack2MC",0) - if sTree.GetBranch("FitTracks_PR"): sTree.SetBranchStatus("FitTracks_PR",0) - if sTree.GetBranch("Particles_PR"): sTree.SetBranchStatus("Particles_PR",0) - if sTree.GetBranch("fitTrack2MC_PR"): sTree.SetBranchStatus("fitTrack2MC_PR",0) - if sTree.GetBranch("MufluxSpectrometerPoint"): - if sTree.GetBranch("Digi_MufluxSpectrometerHits"): - sTree.SetBranchStatus("Digi_MufluxSpectrometerHits",0) - - rawFile = fout.replace("_rec.root","_raw.root") - recf = ROOT.TFile(rawFile,"recreate") - newTree = sTree.CloneTree(0) - for n in range(sTree.GetEntries()): - sTree.GetEntry(n) - rc = newTree.Fill() - sTree.Clear() - newTree.AutoSave() - f.Close() - recf.Close() - os.system('cp '+rawFile +' '+fout) - self.fn = ROOT.TFile(fout,'update') - self.sTree = self.fn.cbmsim - - # prepare for output - if simulation: - # event header - self.header = ROOT.FairEventHeader() - self.eventHeader = self.sTree.Branch("ShipEventHeader",self.header,32000,-1) - self.fitTrack2MC = ROOT.std.vector('int')() - self.mcLink = self.sTree.Branch("fitTrack2MC" + global_variables.realPR, self.fitTrack2MC, 32000, -1) - self.digiMufluxSpectrometer = ROOT.TClonesArray("MufluxSpectrometerHit") - self.digiMufluxSpectrometerBranch = self.sTree.Branch("Digi_MufluxSpectrometerHits",self.digiMufluxSpectrometer,32000,-1) - #muon taggger - if self.sTree.GetBranch("MuonTaggerPoint"): - self.digiMuonTagger = ROOT.TClonesArray("MuonTaggerHit") - self.digiMuonTaggerBranch = self.sTree.Branch("Digi_MuonTagger", self.digiMuonTagger, 32000, -1) - # setup random number generator - self.random = ROOT.TRandom() - ROOT.gRandom.SetSeed() - # fitted tracks - self.fGenFitArray = ROOT.TClonesArray("genfit::Track") - self.fGenFitArray.BypassStreamer(ROOT.kFALSE) - self.fitTracks = self.sTree.Branch("FitTracks" + global_variables.realPR, self.fGenFitArray, 32000, -1) - - self.PDG = ROOT.TDatabasePDG.Instance() - # for the digitizing and reconstruction step - self.v_drift = global_variables.modules["MufluxSpectrometer"].TubeVdrift() - self.sigma_spatial = global_variables.modules["MufluxSpectrometer"].TubeSigmaSpatial() - self.viewangle = global_variables.modules["MufluxSpectrometer"].ViewAngle() - - # access ShipTree - self.sTree.GetEvent(0) - self.geoMat = ROOT.genfit.TGeoMaterialInterface() - # init geometry and mag. field - self.gMan = ROOT.gGeoManager - self.bfield = ROOT.genfit.FairShipFields() - self.fM = ROOT.genfit.FieldManager.getInstance() - self.fM.init(self.bfield) - - ROOT.genfit.MaterialEffects.getInstance().init(self.geoMat) - - # init fitter, to be done before importing shipPatRec - self.fitter = ROOT.genfit.DAF() - self.fitter.setMaxIterations(20) - - if global_variables.debug: - - # self.fitter.setDebugLvl(1) # produces lot of printout - output_dir = 'pics/' - if os.path.exists(output_dir): - print('The directory already exists. It is OK.') - else: - os.mkdir(output_dir) - - def reconstruct(self): - ntracks = self.findTracks() - - def digitize(self): - - self.sTree.t0 = self.random.Rndm()*1*u.microsecond - self.header.SetEventTime( self.sTree.t0 ) - self.header.SetRunId( self.sTree.MCEventHeader.GetRunID() ) - self.header.SetMCEntryNumber( self.sTree.MCEventHeader.GetEventID() ) # counts from 1 - self.eventHeader.Fill() - self.digiMufluxSpectrometer.Delete() - self.digitizeMufluxSpectrometer() - self.digiMufluxSpectrometerBranch.Fill() - #self.digiMuonTagger.Delete() # produces a lot of warnings, rpc station 0 - #self.digitizeMuonTagger() - #self.digiMuonTaggerBranch.Fill() - - def digitizeMuonTagger(self, fake_clustering=True): - - station = 0 - strip = 0 - DetectorID = set() # set of detector ids - already deduplicated - for MuonTaggerHit in self.sTree.MuonTaggerPoint: - # getting rpc nodes, name and matrix - rpc_box = self.gMan.FindNode( - MuonTaggerHit.GetX(), - MuonTaggerHit.GetY(), - MuonTaggerHit.GetZ()) - rpc = rpc_box.GetName() - master_matrix = rpc_box.GetMatrix() - - # getting muon box volume (lower level) - muon_box = self.gMan.GetTopVolume().FindNode("VMuonBox_1") - muonbox_matrix = muon_box.GetMatrix() - - # translation from top to MuonBox_1 - point = array('d', [ - MuonTaggerHit.GetX(), - MuonTaggerHit.GetY(), - MuonTaggerHit.GetZ(), 1]) - point_muonbox = array('d', [0, 0, 0, 1]) - muonbox_matrix.MasterToLocal(point, point_muonbox) - - # translation to local frame - point_local = array('d', [0, 0, 0, 1]) - master_matrix.MasterToLocal(point_muonbox, point_local) - - xcoord = point_local[0] - ycoord = point_local[1] - - # identify individual rpcs - station = int(rpc[-1]) - if station not in range(1, 6): # limiting the range of rpcs - print("WARNING: Invalid RPC number, something's wrong with the geometry ",station) - - # calculate strip - # x gives vertical direction - direction = 1 - strip = StripX(xcoord) - if not strip: - continue - # sampling number of strips around the exact strip for emulating clustering - if fake_clustering: - s = np.random.poisson(3) - strip = strip - int(s/2) - for i in range(0, s): - detectorid = station*10000 + direction*1000 + strip + i - DetectorID.add(detectorid) - else: - detectorid = station*10000 + direction*1000 + strip - DetectorID.add(detectorid) - - # y gives horizontal direction - direction = 0 - strip = StripY(ycoord) - if not strip: - continue - # sampling number of strips around the exact strip for emulating clustering - if fake_clustering: - s = np.random.poisson(3) - strip = strip - int(s/2) - for i in range(0, s): - detectorid = station*10000 + direction*1000 + strip + i - DetectorID.add(detectorid) - else: - detectorid = station*10000 + direction*1000 + strip - DetectorID.add(detectorid) - - self.digiMuonTagger.Expand(len(DetectorID)) - for index, detID in enumerate(DetectorID): - hit = ROOT.MuonTaggerHit(detID, 0) - self.digiMuonTagger[index] = hit - - if fake_clustering: - # cluster size loop - plotting the cluster size distribution - cluster_size = list() - DetectorID_list = list(DetectorID) # turn set into list to allow indexing - DetectorID_list.sort() # sorting the list - if len(DetectorID_list) > 1: - clusters = [[DetectorID_list[0]]] - for x in DetectorID_list[1:]: - if abs(x - clusters[-1][-1]) <= 1: - clusters[-1].append(x) - else: - clusters.append([x]) - cluster_size = [len(x) for x in clusters] - for i in cluster_size: - h['muontagger_clusters'].Fill(i) - - - def digitizeMufluxSpectrometer(self): - - # digitize FairSHiP MC hits - index = 0 - hitsPerDetId = {} - - for aMCPoint in self.sTree.MufluxSpectrometerPoint: - aHit = ROOT.MufluxSpectrometerHit(aMCPoint,self.sTree.t0) - if self.digiMufluxSpectrometer.GetSize() == index: self.digiMufluxSpectrometer.Expand(index+1000) - self.digiMufluxSpectrometer[index]=aHit - detID = aHit.GetDetectorID() - if detID in hitsPerDetId: - if self.digiMufluxSpectrometer[hitsPerDetId[detID]].tdc() > aHit.tdc(): - # second hit with smaller tdc - self.digiMufluxSpectrometer[hitsPerDetId[detID]].setInvalid() - hitsPerDetId[detID] = index - else: - hitsPerDetId[detID] = index - index+=1 - - T1_entries_px = {} - T4_entries_px = {} - nMufluxHits = self.sTree.MufluxSpectrometerPoint.GetEntriesFast() - for i in range(nMufluxHits): - MufluxHit = self.sTree.MufluxSpectrometerPoint[i] - detector = self.gMan.FindNode(MufluxHit.GetX(),MufluxHit.GetY(),MufluxHit.GetZ()).GetName() - MufluxTrackId = MufluxHit.GetTrackID() - pid = MufluxHit.PdgCode() - xcoord = MufluxHit.GetX() - ycoord = MufluxHit.GetY() - if abs(pid)==13: - if (detector[0:8]=="gas_12_1"): - rc=h['hits-T1'].Fill(xcoord,ycoord) - if (detector[0:9]=="gas_12_10"): - rc=h['hits-T1x'].Fill(xcoord,ycoord) - if (detector[0:9]=="gas_12_11"): - rc=h['hits-T1u'].Fill(xcoord,ycoord) - if (detector[0:8]=="gas_12_2"): - rc=h['hits-T2'].Fill(xcoord,ycoord) - if (detector[0:9]=="gas_12_20"): - rc=h['hits-T2v'].Fill(xcoord,ycoord) - if (detector[0:9]=="gas_12_21"): - rc=h['hits-T2x'].Fill(xcoord,ycoord) - if (detector[0:5]=="gas_3"): - rc=h['hits-T3'].Fill(xcoord,ycoord) - if (detector[0:5]=="gas_4"): - rc=h['hits-T4'].Fill(xcoord,ycoord) - - if (detector[0:9]=="gas_12_10"): - if MufluxTrackId in T1_entries_px: - continue - else: - if abs(pid)==13 : - T1_entries_px[MufluxTrackId]=[MufluxHit.GetPx()] - - if (detector[0:5]=="gas_4"): - if MufluxTrackId in T4_entries_px: - continue - else: - pid = MufluxHit.PdgCode() - if abs(pid)==13 : - T4_entries_px[MufluxTrackId]=[MufluxHit.GetPx()] - - for i in range(nMufluxHits): - MufluxHit = self.sTree.MufluxSpectrometerPoint[i] - MufluxTrackId = MufluxHit.GetTrackID() - if (T1_entries_px.get(MufluxTrackId) is None or T4_entries_px.get(MufluxTrackId) is None) : - continue - else: - rc=h['pt-kick'].Fill(T1_entries_px.get(MufluxTrackId)[0]-T4_entries_px.get(MufluxTrackId)[0]) - - def withT0Estimate(self): - # loop over all straw tdcs and make average, correct for ToF - n = 0 - t0 = 0. - key = -1 - SmearedHits = [] - v_drift = global_variables.modules["MufluxSpectrometer"].TubeVdrift() - z1 = stop.z() - for aDigi in self.digiMufluxSpectrometer: - key+=1 - if not aDigi.isValid: continue - detID = aDigi.GetDetectorID() - # don't use hits from straw veto - station = int(detID/10000000) - if station > 4 : continue - aDigi.MufluxSpectrometerEndPoints(start,stop) - # MufluxSpectrometerHit::MufluxSpectrometerEndPoints(TVector3 &vbot, TVector3 &vtop) - delt1 = (start[2]-z1)/u.speedOfLight - t0+=aDigi.GetDigi()-delt1 - SmearedHits.append( {'digiHit':key,'xtop':stop.x(),'ytop':stop.y(),'z':stop.z(),'xbot':start.x(),'ybot':start.y(),'dist':aDigi.GetDigi(), 'detID':detID} ) - n+=1 - if n>0: - t0 = t0/n - 73.2*u.ns - print("t0 ",t0) - for s in SmearedHits: - delt1 = (s['z']-z1)/u.speedOfLight - s['dist'] = (s['dist'] -delt1 -t0)*v_drift - print("s['dist']",s['dist']) - return SmearedHits - - def smearHits(self,no_amb=None): - # smear strawtube points - SmearedHits = [] - key = -1 - for ahit in self.sTree.MufluxSpectrometerPoint: - key+=1 - detID = ahit.GetDetectorID() - top = ROOT.TVector3() - bot = ROOT.TVector3() - global_variables.modules["MufluxSpectrometer"].TubeEndPoints(detID, bot, top) - # MufluxSpectrometerHit::MufluxSpectrometerEndPoints(TVector3 &vbot, TVector3 &vtop) - #distance to wire, and smear it. - dw = ahit.dist2Wire() - smear = dw - if not no_amb: - smear = abs(self.random.Gaus(dw,self.sigma_spatial)) - - SmearedHits.append( {'digiHit':key,'xtop':top.x(),'ytop':top.y(),'z':top.z(),'xbot':bot.x(),'ybot':bot.y(),'dist':smear, 'detID':detID} ) - # Note: top.z()==bot.z() unless misaligned, so only add key 'z' to smearedHit - - if abs(top.y())==abs(bot.y()): h['disty'].Fill(dw) - if abs(top.y())>abs(bot.y()): h['distu'].Fill(dw) - if abs(top.y())1: - print("something wrong with detector id",detid) - pnb = 0 - return statnb,vnb,pnb,lnb,view - - def correctAlignment(self, hit): - # Taken from charmdet/drifttubeMonitoring.py - sqrt2 = ROOT.TMath.Sqrt(2.) - cos30 = ROOT.TMath.Cos(30./180.*ROOT.TMath.Pi()) - sin30 = ROOT.TMath.Sin(30./180.*ROOT.TMath.Pi()) - #delX=[[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]] - delX=[[-1.2,1.4,-0.9,0.7],[-1.69-1.31,0.3+2.82,0.6-0.4,-0.8-0.34],[-0.9,1.14,-0.80,1.1],[-0.8,1.29,0.15,2.0]] - #delUV = [[1.8,-2.,0.8,-1.],[1.2,-1.1,-0.5,0.3]] - delUV = [[0,0,0,0],[0,0,0,0]] - vtop = ROOT.TVector3() - vbot = ROOT.TVector3() - rc = hit.MufluxSpectrometerEndPoints(vbot,vtop) - s,v,p,l,view = self.stationInfo(hit) - if view=='_x': - cor = delX[s-1][2*l+p] - vbot[0]=vbot[0]+cor - vtop[0]=vtop[0]+cor - else: - if view=='_u': cor = delUV[0][2*l+p] - elif view=='_v': cor = delUV[1][2*l+p] - vbot[0]=vbot[0]+cor*cos30 - vtop[0]=vtop[0]+cor*cos30 - vbot[1]=vbot[1]+cor*sin30 - vtop[1]=vtop[1]+cor*sin30 - #tmp = vbot[0]*ROOT.TMath.Cos(rotUV)-vbot[1]*ROOT.TMath.Sin(rotUV) - #vbot[1]=vbot[1]*ROOT.TMath.Cos(rotUV)+vbot[0]*ROOT.TMath.Sin(rotUV) - #vbot[0]=tmp - #tmp = vtop[0]*ROOT.TMath.Cos(rotUV)-vtop[1]*ROOT.TMath.Sin(rotUV) - #vtop[1]=vtop[1]*ROOT.TMath.Cos(rotUV)+vtop[0]*ROOT.TMath.Sin(rotUV) - #vbot[0]=tmp - return vbot,vtop - - - def RT(self, s,t,function='parabola'): - # Taken from charmdet/drifttubeMonitoring.py - # rt relation, drift time to distance, drift time? - tMinAndTmax = {1:[587,1860],2:[587,1860],3:[610,2300],4:[610,2100]} - R = global_variables.ShipGeo.MufluxSpectrometer.InnerTubeDiameter/2. # = 3.63*u.cm - # parabola - if function == 'parabola' or 'rtTDC'+str(s)+'000_x' not in h: - p1p2 = {1:[688.,7.01],2:[688.,7.01],3:[923.,4.41],4:[819.,0.995]} - t0_corr = max(0,t-tMinAndTmax[s][0]) - tmp1 = ROOT.TMath.Sqrt(p1p2[s][0]**2+4.*p1p2[s][1]*t0_corr) - r = (2*tmp1-2*p1p2[s][0]+p1p2[s][0]*ROOT.TMath.Log(-p1p2[s][0]+2*tmp1)-p1p2[s][0]*ROOT.TMath.Log(p1p2[s][0]) )/(8*p1p2[s][1]) - else: - if t > tMinAndTmax[s][1]: r = R - elif t< tMinAndTmax[s][0]: r = 0 - else: r = h['rtTDC'+str(s)+'000_x'].Eval(t) - # h['TDC2R'].Fill(t,r) - return r - - - def smearHits_realData(self): - - # smear strawtube points - SmearedHits = [] - key = -1 - for ahit in self.sTree.Digi_MufluxSpectrometerHits: - key+=1 - detID = ahit.GetDetectorID() - top = ROOT.TVector3() - bot = ROOT.TVector3() - bot, top = self.correctAlignment(ahit) - # global_variables.modules["MufluxSpectrometer"].TubeEndPoints(detID, bot, top) - # ahit.MufluxSpectrometerEndPoints(bot,top) - # MufluxSpectrometerHit::MufluxSpectrometerEndPoints(TVector3 &vbot, TVector3 &vtop) - # distance to wire. - # dist = ahit.GetDigi() * 3.7 / (2000. * 2.) - tdc = ahit.GetDigi() - dist = self.RT(ahit.GetDetectorID()/10000000, tdc) - - SmearedHits.append( {'digiHit':key,'xtop':top.x(),'ytop':top.y(),'z':top.z(),'xbot':bot.x(),'ybot':bot.y(),'dist':dist, 'detID':detID} ) - - if abs(top.y())==abs(bot.y()): h['disty'].Fill(dist) - if abs(top.y())>abs(bot.y()): h['distu'].Fill(dist) - if abs(top.y()) 0: - frac = float(track[tmax]) / float(nh) - - return frac,tmax - - - def recognition_metrics(self, track_hit_ids, mode='all'): - - n_tracks_y12 = 0 - n_recognized_y12 = 0 - n_clones_y12 = 0 - n_ghosts_y12 = 0 - n_others_y12 = 0 - - n_tracks_stereo12 = 0 - n_recognized_stereo12 = 0 - n_clones_stereo12 = 0 - n_ghosts_stereo12 = 0 - n_others_stereo12 = 0 - - n_tracks_34 = 0 - n_recognized_34 = 0 - n_clones_34 = 0 - n_ghosts_34 = 0 - n_others_34 = 0 - - true_track_ids_y12 = [] - true_track_ids_stereo12 = [] - true_track_ids_34 = [] - - n_true_track_hits_y12 = {} - n_true_track_hits_stereo12 = {} - n_true_track_hits_34 = {} - - n_true_track_hits_1 = {} - n_true_track_hits_2 = {} - n_true_track_hits_3 = {} - n_true_track_hits_4 = {} - - for ahit in self.sTree.MufluxSpectrometerPoint: - - track_id = ahit.GetTrackID() - - detID = ahit.GetDetectorID() - statnb = detID // 10000000 - vnb = (detID - statnb * 10000000) // 1000000 - - is_y12 = (statnb == 1) * (vnb == 0) + (statnb == 2) * (vnb == 1) - is_stereo12 = (statnb == 1) * (vnb == 1) + (statnb == 2) * (vnb == 0) - is_34 = (statnb == 3) + (statnb == 4) - - if statnb == 1: - if track_id in n_true_track_hits_1: - n_true_track_hits_1[track_id] += 1 - else: - n_true_track_hits_1[track_id] = 1 - - if statnb == 2: - if track_id in n_true_track_hits_2: - n_true_track_hits_2[track_id] += 1 - else: - n_true_track_hits_2[track_id] = 1 - - if statnb == 3: - if track_id in n_true_track_hits_3: - n_true_track_hits_3[track_id] += 1 - else: - n_true_track_hits_3[track_id] = 1 - - if statnb == 4: - if track_id in n_true_track_hits_4: - n_true_track_hits_4[track_id] += 1 - else: - n_true_track_hits_4[track_id] = 1 - - if is_y12: - if track_id in n_true_track_hits_y12: - n_true_track_hits_y12[track_id] += 1 - else: - n_true_track_hits_y12[track_id] = 1 - - if is_stereo12: - if track_id in n_true_track_hits_stereo12: - n_true_track_hits_stereo12[track_id] += 1 - else: - n_true_track_hits_stereo12[track_id] = 1 - - if is_34: - if track_id in n_true_track_hits_34: - n_true_track_hits_34[track_id] += 1 - else: - n_true_track_hits_34[track_id] = 1 - - - if mode == 'all': - min_hits = 1 - elif mode == '3hits': - min_hits = 3 - - if mode == 'all' or mode == '3hits': - for key in n_true_track_hits_y12.keys(): - if n_true_track_hits_y12[key] >= min_hits: - true_track_ids_y12.append(key) - - for key in n_true_track_hits_stereo12.keys(): - if n_true_track_hits_stereo12[key] >= min_hits: - true_track_ids_stereo12.append(key) - - for key in n_true_track_hits_34.keys(): - if n_true_track_hits_34[key] >= min_hits: - true_track_ids_34.append(key) - - if mode == 'Tr4': - min_hits = 1 - for key in n_true_track_hits_y12.keys(): - if key in n_true_track_hits_1 and key in n_true_track_hits_2: - if key in n_true_track_hits_3 and key in n_true_track_hits_4: - if n_true_track_hits_y12[key] >= min_hits: - true_track_ids_y12.append(key) - - for key in n_true_track_hits_stereo12.keys(): - if key in n_true_track_hits_1 and key in n_true_track_hits_2: - if key in n_true_track_hits_3 and key in n_true_track_hits_4: - if n_true_track_hits_stereo12[key] >= min_hits: - true_track_ids_stereo12.append(key) - - for key in n_true_track_hits_34.keys(): - if key in n_true_track_hits_1 and key in n_true_track_hits_2: - if key in n_true_track_hits_3 and key in n_true_track_hits_4: - if n_true_track_hits_34[key] >= min_hits: - true_track_ids_34.append(key) - - - - n_tracks_y12 = len(true_track_ids_y12) - n_tracks_stereo12 = len(true_track_ids_stereo12) - n_tracks_34 = len(true_track_ids_34) - - found_track_ids_y12 = [] - found_track_ids_stereo12 = [] - found_track_ids_34 = [] - - for i_track in track_hit_ids.keys(): - - atrack = track_hit_ids[i_track] - atrack_y12 = atrack['y12'] - atrack_stereo12 = atrack['stereo12'] - atrack_34 = atrack['34'] - atrack_y_tagger = atrack['y_tagger'] - - if len(atrack_y12) > 0: - reco_hit_ids_y12 = [] - for i_hit in atrack_y12: - ahit = self.sTree.MufluxSpectrometerPoint[i_hit['digiHit']] - reco_hit_ids_y12.append(ahit.GetTrackID()) - frac_y12, tmax_y12 = self.fracMCsame(reco_hit_ids_y12) - if frac_y12 >= 0.7: - if tmax_y12 in true_track_ids_y12 and tmax_y12 not in found_track_ids_y12: - n_recognized_y12 += 1 - found_track_ids_y12.append(tmax_y12) - elif tmax_y12 in true_track_ids_y12 and tmax_y12 in found_track_ids_y12: - n_clones_y12 += 1 - elif tmax_y12 not in true_track_ids_y12: - n_others_y12 += 1 - else: - n_ghosts_y12 += 1 - - - if len(atrack_stereo12) > 0: - reco_hit_ids_stereo12 = [] - for i_hit in atrack_stereo12: - ahit = self.sTree.MufluxSpectrometerPoint[i_hit['digiHit']] - reco_hit_ids_stereo12.append(ahit.GetTrackID()) - frac_stereo12, tmax_stereo12 = self.fracMCsame(reco_hit_ids_stereo12) - if frac_stereo12 >= 0.7: - if tmax_stereo12 in true_track_ids_stereo12 and tmax_stereo12 not in found_track_ids_stereo12: - n_recognized_stereo12 += 1 - found_track_ids_stereo12.append(tmax_stereo12) - elif tmax_stereo12 in true_track_ids_stereo12 and tmax_stereo12 in found_track_ids_stereo12: - n_clones_stereo12 += 1 - elif tmax_stereo12 not in true_track_ids_stereo12: - n_others_stereo12 += 1 - else: - n_ghosts_stereo12 += 1 - - - if len(atrack_34) > 0: - reco_hit_ids_34 = [] - for i_hit in atrack_34: - ahit = self.sTree.MufluxSpectrometerPoint[i_hit['digiHit']] - reco_hit_ids_34.append(ahit.GetTrackID()) - frac_34, tmax_34 = self.fracMCsame(reco_hit_ids_34) - if frac_34 >= 0.7: - if tmax_34 in true_track_ids_34 and tmax_34 not in found_track_ids_34: - n_recognized_34 += 1 - found_track_ids_34.append(tmax_34) - elif tmax_34 in true_track_ids_34 and tmax_34 in found_track_ids_34: - n_clones_34 += 1 - elif tmax_34 not in true_track_ids_34: - n_others_34 += 1 - else: - n_ghosts_34 += 1 - - output = (n_tracks_y12, n_recognized_y12, n_clones_y12, n_ghosts_y12, n_others_y12, - n_tracks_stereo12, n_recognized_stereo12, n_clones_stereo12, n_ghosts_stereo12, n_others_stereo12, - n_tracks_34, n_recognized_34, n_clones_34, n_ghosts_34, n_others_34) - - return output - - - def target_metrics(self, track_hit_ids): - - n_tracks = 0 - n_recognized = 0 - n_clones = 0 - n_ghosts = 0 - n_others = 0 - - true_track_ids = [] - - n_true_track_hits_y12 = {} - n_true_track_hits_stereo12 = {} - n_true_track_hits_34 = {} - - for ahit in self.sTree.MufluxSpectrometerPoint: - - track_id = ahit.GetTrackID() - - detID = ahit.GetDetectorID() - statnb = detID // 10000000 - vnb = (detID - statnb * 10000000) // 1000000 - - is_y12 = (statnb == 1) * (vnb == 0) + (statnb == 2) * (vnb == 1) - is_stereo12 = (statnb == 1) * (vnb == 1) + (statnb == 2) * (vnb == 0) - is_34 = (statnb == 3) + (statnb == 4) - - if is_y12: - if track_id in n_true_track_hits_y12: - n_true_track_hits_y12[track_id] += 1 - else: - n_true_track_hits_y12[track_id] = 1 - - if is_stereo12: - if track_id in n_true_track_hits_stereo12: - n_true_track_hits_stereo12[track_id] += 1 - else: - n_true_track_hits_stereo12[track_id] = 1 - - if is_34: - if track_id in n_true_track_hits_34: - n_true_track_hits_34[track_id] += 1 - else: - n_true_track_hits_34[track_id] = 1 - - min_hits = 3 - for key in n_true_track_hits_y12.keys(): - if n_true_track_hits_y12[key] >= min_hits: - if key in n_true_track_hits_stereo12: - if n_true_track_hits_stereo12[key] >= min_hits: - if key in n_true_track_hits_34: - if n_true_track_hits_34[key] >= min_hits: - true_track_ids.append(key) - - n_tracks = len(true_track_ids) - - found_track_ids = [] - - for i_track in track_hit_ids.keys(): - - atrack = track_hit_ids[i_track] - atrack_y12 = atrack['y12'] - atrack_stereo12 = atrack['stereo12'] - atrack_34 = atrack['34'] - atrack_y_tagger = atrack['y_tagger'] - - if len(atrack_y12) >= min_hits and len(atrack_stereo12) >= min_hits and len(atrack_34) >= min_hits and len(atrack_y_tagger) >= withNTaggerHits: - - reco_hit_ids_y12 = [] - for i_hit in atrack_y12: - ahit = self.sTree.MufluxSpectrometerPoint[i_hit['digiHit']] - reco_hit_ids_y12.append(ahit.GetTrackID()) - frac_y12, tmax_y12 = self.fracMCsame(reco_hit_ids_y12) - - reco_hit_ids_stereo12 = [] - for i_hit in atrack_stereo12: - ahit = self.sTree.MufluxSpectrometerPoint[i_hit['digiHit']] - reco_hit_ids_stereo12.append(ahit.GetTrackID()) - frac_stereo12, tmax_stereo12 = self.fracMCsame(reco_hit_ids_stereo12) - - reco_hit_ids_34 = [] - for i_hit in atrack_34: - ahit = self.sTree.MufluxSpectrometerPoint[i_hit['digiHit']] - reco_hit_ids_34.append(ahit.GetTrackID()) - frac_34, tmax_34 = self.fracMCsame(reco_hit_ids_34) - - - min_eff = 0.7 - if tmax_y12 == tmax_stereo12 and tmax_y12 == tmax_34: - if frac_y12 >= min_eff and frac_stereo12 >= min_eff and frac_34 >= min_eff: - - if tmax_y12 in true_track_ids and tmax_y12 not in found_track_ids: - n_recognized += 1 - found_track_ids.append(tmax_y12) - elif tmax_y12 in true_track_ids and tmax_y12 in found_track_ids: - n_clones += 1 - elif tmax_y12 not in true_track_ids: - n_others += 1 - - else: - n_ghosts += 1 - else: - n_ghosts += 1 - - else: - n_ghosts += 1 - - output = (n_tracks, n_recognized, n_clones, n_ghosts, n_others) - - return output - - - def all_tracks_metrics(self, track_hit_ids): - - n_tracks = 0 - n_recognized = 0 - n_clones = 0 - n_ghosts = 0 - n_others = 0 - - true_track_ids = [] - for ahit in self.sTree.MufluxSpectrometerPoint: - track_id = ahit.GetTrackID() - if track_id not in true_track_ids: - true_track_ids.append(track_id) - - n_tracks = len(true_track_ids) - found_track_ids = [] - min_hits = 3 - - for i_track in track_hit_ids.keys(): - - atrack = track_hit_ids[i_track] - atrack_y12 = atrack['y12'] - atrack_stereo12 = atrack['stereo12'] - atrack_34 = atrack['34'] - atrack_y_tagger = atrack['y_tagger'] - - if len(atrack_y12) >= min_hits and len(atrack_stereo12) >= min_hits and len(atrack_34) >= min_hits and len(atrack_y_tagger) >= withNTaggerHits: - - reco_hit_ids_y12 = [] - for i_hit in atrack_y12: - ahit = self.sTree.MufluxSpectrometerPoint[i_hit['digiHit']] - reco_hit_ids_y12.append(ahit.GetTrackID()) - frac_y12, tmax_y12 = self.fracMCsame(reco_hit_ids_y12) - - reco_hit_ids_stereo12 = [] - for i_hit in atrack_stereo12: - ahit = self.sTree.MufluxSpectrometerPoint[i_hit['digiHit']] - reco_hit_ids_stereo12.append(ahit.GetTrackID()) - frac_stereo12, tmax_stereo12 = self.fracMCsame(reco_hit_ids_stereo12) - - reco_hit_ids_34 = [] - for i_hit in atrack_34: - ahit = self.sTree.MufluxSpectrometerPoint[i_hit['digiHit']] - reco_hit_ids_34.append(ahit.GetTrackID()) - frac_34, tmax_34 = self.fracMCsame(reco_hit_ids_34) - - min_eff = 0.7 - if tmax_y12 == tmax_stereo12 and tmax_y12 == tmax_34: - if frac_y12 >= min_eff and frac_stereo12 >= min_eff and frac_34 >= min_eff: - - if tmax_y12 in true_track_ids and tmax_y12 not in found_track_ids: - n_recognized += 1 - found_track_ids.append(tmax_y12) - elif tmax_y12 in true_track_ids and tmax_y12 in found_track_ids: - n_clones += 1 - elif tmax_y12 not in true_track_ids: - n_others += 1 - - else: - n_ghosts += 1 - else: - n_ghosts += 1 - - else: - n_ghosts += 1 - - output = (n_tracks, n_recognized, n_clones, n_ghosts, n_others) - - return output - - - def muon_metrics(self, track_hit_ids): - - n_tracks = 0 - n_recognized = 0 - n_clones = 0 - n_ghosts = 0 - n_others = 0 - - true_track_ids = [] - for ahit in self.sTree.MufluxSpectrometerPoint: - track_id = ahit.GetTrackID() - pdg = ahit.PdgCode() - if track_id not in true_track_ids: - if abs(pdg)==13: - true_track_ids.append(track_id) - - n_tracks = len(true_track_ids) - found_track_ids = [] - min_hits = 3 - - for i_track in track_hit_ids.keys(): - - atrack = track_hit_ids[i_track] - atrack_y12 = atrack['y12'] - atrack_stereo12 = atrack['stereo12'] - atrack_34 = atrack['34'] - atrack_y_tagger = atrack['y_tagger'] - - if len(atrack_y12) >= min_hits and len(atrack_stereo12) >= min_hits and len(atrack_34) >= min_hits and len(atrack_y_tagger) >= withNTaggerHits: - - reco_hit_ids_y12 = [] - for i_hit in atrack_y12: - ahit = self.sTree.MufluxSpectrometerPoint[i_hit['digiHit']] - reco_hit_ids_y12.append(ahit.GetTrackID()) - frac_y12, tmax_y12 = self.fracMCsame(reco_hit_ids_y12) - - reco_hit_ids_stereo12 = [] - for i_hit in atrack_stereo12: - ahit = self.sTree.MufluxSpectrometerPoint[i_hit['digiHit']] - reco_hit_ids_stereo12.append(ahit.GetTrackID()) - frac_stereo12, tmax_stereo12 = self.fracMCsame(reco_hit_ids_stereo12) - - reco_hit_ids_34 = [] - for i_hit in atrack_34: - ahit = self.sTree.MufluxSpectrometerPoint[i_hit['digiHit']] - reco_hit_ids_34.append(ahit.GetTrackID()) - frac_34, tmax_34 = self.fracMCsame(reco_hit_ids_34) - - min_eff = 0.7 - if tmax_y12 == tmax_stereo12 and tmax_y12 == tmax_34: - if frac_y12 >= min_eff and frac_stereo12 >= min_eff and frac_34 >= min_eff: - - if tmax_y12 in true_track_ids and tmax_y12 not in found_track_ids: - n_recognized += 1 - found_track_ids.append(tmax_y12) - elif tmax_y12 in true_track_ids and tmax_y12 in found_track_ids: - n_clones += 1 - elif tmax_y12 not in true_track_ids: - n_others += 1 - - else: - n_ghosts += 1 - else: - n_ghosts += 1 - - else: - n_ghosts += 1 - - output = (n_tracks, n_recognized, n_clones, n_ghosts, n_others) - - return output - - - def muon_with_tagger_metrics(self, track_hit_ids): - - n_tracks = 0 - n_recognized = 0 - n_clones = 0 - n_ghosts = 0 - n_others = 0 - - true_tagger_track_ids = [] - for ahit in self.sTree.MuonTaggerPoint: - track_id = ahit.GetTrackID() - pdg = ahit.PdgCode() - if track_id not in true_tagger_track_ids: - if abs(pdg)==13: - true_tagger_track_ids.append(track_id) - - true_track_ids = [] - for ahit in self.sTree.MufluxSpectrometerPoint: - track_id = ahit.GetTrackID() - pdg = ahit.PdgCode() - if track_id not in true_track_ids: - if abs(pdg)==13: - if track_id in true_tagger_track_ids: - true_track_ids.append(track_id) - - n_tracks = len(true_track_ids) - found_track_ids = [] - min_hits = 3 - - for i_track in track_hit_ids.keys(): - - atrack = track_hit_ids[i_track] - atrack_y12 = atrack['y12'] - atrack_stereo12 = atrack['stereo12'] - atrack_34 = atrack['34'] - atrack_y_tagger = atrack['y_tagger'] - - if len(atrack_y12) >= min_hits and len(atrack_stereo12) >= min_hits and len(atrack_34) >= min_hits and len(atrack_y_tagger) >= withNTaggerHits: - - reco_hit_ids_y12 = [] - for i_hit in atrack_y12: - ahit = self.sTree.MufluxSpectrometerPoint[i_hit['digiHit']] - reco_hit_ids_y12.append(ahit.GetTrackID()) - frac_y12, tmax_y12 = self.fracMCsame(reco_hit_ids_y12) - - reco_hit_ids_stereo12 = [] - for i_hit in atrack_stereo12: - ahit = self.sTree.MufluxSpectrometerPoint[i_hit['digiHit']] - reco_hit_ids_stereo12.append(ahit.GetTrackID()) - frac_stereo12, tmax_stereo12 = self.fracMCsame(reco_hit_ids_stereo12) - - reco_hit_ids_34 = [] - for i_hit in atrack_34: - ahit = self.sTree.MufluxSpectrometerPoint[i_hit['digiHit']] - reco_hit_ids_34.append(ahit.GetTrackID()) - frac_34, tmax_34 = self.fracMCsame(reco_hit_ids_34) - - min_eff = 0.7 - if tmax_y12 == tmax_stereo12 and tmax_y12 == tmax_34: - if frac_y12 >= min_eff and frac_stereo12 >= min_eff and frac_34 >= min_eff: - - if tmax_y12 in true_track_ids and tmax_y12 not in found_track_ids: - n_recognized += 1 - found_track_ids.append(tmax_y12) - elif tmax_y12 in true_track_ids and tmax_y12 in found_track_ids: - n_clones += 1 - elif tmax_y12 not in true_track_ids: - n_others += 1 - - else: - n_ghosts += 1 - else: - n_ghosts += 1 - - else: - n_ghosts += 1 - - output = (n_tracks, n_recognized, n_clones, n_ghosts, n_others) - - return output - - - def muon_metrics_vs_p(self, track_hit_ids): - - n_tracks = 0 - n_recognized = 0 - n_clones = 0 - n_ghosts = 0 - n_others = 0 - - true_track_ids = [] - true_track_p = {} - for ahit in self.sTree.MufluxSpectrometerPoint: - track_id = ahit.GetTrackID() - pdg = ahit.PdgCode() - if track_id not in true_track_ids: - if abs(pdg)==13: - true_track_ids.append(track_id) - if track_id not in true_track_p: - Ptruth,Ptruthx,Ptruthy,Ptruthz = self.getPtruthFirst(track_id) - true_track_p[track_id] = Ptruth - h['True_all_tracks_vs_p_true'].Fill(Ptruth) - if abs(pdg)==13: - h['True_muons_vs_p_true'].Fill(Ptruth) - - n_tracks = len(true_track_ids) - found_track_ids = [] - min_hits = 3 - - for i_track in track_hit_ids.keys(): - - atrack = track_hit_ids[i_track] - atrack_y12 = atrack['y12'] - atrack_stereo12 = atrack['stereo12'] - atrack_34 = atrack['34'] - atrack_y_tagger = atrack['y_tagger'] - - reco_hit_ids_y12 = [] - for i_hit in atrack_y12: - ahit = self.sTree.MufluxSpectrometerPoint[i_hit['digiHit']] - reco_hit_ids_y12.append(ahit.GetTrackID()) - frac_y12, tmax_y12 = self.fracMCsame(reco_hit_ids_y12) - - reco_hit_ids_stereo12 = [] - for i_hit in atrack_stereo12: - ahit = self.sTree.MufluxSpectrometerPoint[i_hit['digiHit']] - reco_hit_ids_stereo12.append(ahit.GetTrackID()) - frac_stereo12, tmax_stereo12 = self.fracMCsame(reco_hit_ids_stereo12) - - reco_hit_ids_34 = [] - for i_hit in atrack_34: - ahit = self.sTree.MufluxSpectrometerPoint[i_hit['digiHit']] - reco_hit_ids_34.append(ahit.GetTrackID()) - frac_34, tmax_34 = self.fracMCsame(reco_hit_ids_34) - - if len(atrack_y12) >= min_hits and len(atrack_stereo12) >= min_hits and len(atrack_34) >= min_hits and len(atrack_y_tagger) >= withNTaggerHits: - - min_eff = 0.7 - if tmax_y12 == tmax_stereo12 and tmax_y12 == tmax_34: - if frac_y12 >= min_eff and frac_stereo12 >= min_eff and frac_34 >= min_eff: - - if tmax_y12 in true_track_ids and tmax_y12 not in found_track_ids: - n_recognized += 1 - found_track_ids.append(tmax_y12) - h['Reco_muons_vs_p_true'].Fill(true_track_p[tmax_y12]) - elif tmax_y12 in true_track_ids and tmax_y12 in found_track_ids: - n_clones += 1 - elif tmax_y12 not in true_track_ids: - n_others += 1 - - else: - n_ghosts += 1 - h['Ghosts_muons_vs_p_true'].Fill(true_track_p[tmax_y12]) - else: - n_ghosts += 1 - h['Ghosts_muons_vs_p_true'].Fill(true_track_p[tmax_y12]) - - else: - n_ghosts += 1 - h['Ghosts_muons_vs_p_true'].Fill(true_track_p[tmax_y12]) - - output = (n_tracks, n_recognized, n_clones, n_ghosts, n_others) - - return output - - - - def findTracks(self): - - hitPosLists = {} - hitPosLists_noT4 = {} - stationCrossed = {} - stationCrossed_noT4 = {} - trackDigiHits = {} - trackDigiHits_noT4 = {} - trackMomentums = {} - trackMomentums_noT4 = {} - trackCandidates = [] - trackCandidates_noT4 = [] - trackCandidates_all = [] - nTrack = -1 - - self.fGenFitArray.Delete() - self.fitTrack2MC.clear() - - # hit smearing - if self.sTree.GetBranch("MufluxSpectrometerPoint"): - if global_variables.withT0: - self.SmearedHits = self.withT0Estimate() - self.TaggerHits = self.muonTaggerHits_mc() - else: - self.SmearedHits = self.smearHits(global_variables.withNoStrawSmearing) - self.TaggerHits = self.muonTaggerHits_mc() - elif self.sTree.GetBranch("Digi_MufluxSpectrometerHits"): - self.SmearedHits = self.smearHits_realData() - self.TaggerHits = self.muonTaggerHits_realData() - else: - print('No branches "MufluxSpectrometer" or "Digi_MufluxSpectrometerHits".') - return nTrack - - if global_variables.realPR: - - # Do real PatRec - track_hits = MufluxPatRec.execute(self.SmearedHits, self.TaggerHits, withNTaggerHits, withDist2Wire) - - # Create hitPosLists for track fit - for i_track in track_hits.keys(): - - atrack = track_hits[i_track] - atrack_y12 = atrack['y12'] - atrack_stereo12 = atrack['stereo12'] - atrack_34 = atrack['34'] - atrack_smeared_hits = list(atrack_y12) + list(atrack_stereo12) + list(atrack_34) - atrack_p = np.abs(atrack['p']) - atrack_y_in_magnet = atrack['y_in_magnet'] - atrack_x_in_magnet = atrack['x_in_magnet'] - - for sm in atrack_smeared_hits: - - # detID = self.digiMufluxSpectrometer[sm['digiHit']].GetDetectorID() - detID = sm['detID'] - station = int(detID/10000000) - trID = i_track - - # T1-4 for track fit - if trID not in hitPosLists: - hitPosLists[trID] = ROOT.std.vector('TVectorD')() - stationCrossed[trID] = {} - trackDigiHits[trID] = [] - trackMomentums[trID] = atrack_p - m = array('d',[sm['xtop'],sm['ytop'],sm['z'],sm['xbot'],sm['ybot'],sm['z'],sm['dist']]) - hitPosLists[trID].push_back(ROOT.TVectorD(7,m)) - if station not in stationCrossed[trID]: - stationCrossed[trID][station]=0 - stationCrossed[trID][station]+=1 - trackDigiHits[trID].append(sm['digiHit']) - - # T1-3 for track fit - if (int(detID/1000000)!=40): - if trID not in hitPosLists_noT4: - hitPosLists_noT4[trID] = ROOT.std.vector('TVectorD')() - stationCrossed_noT4[trID] = {} - trackDigiHits_noT4[trID] = [] - trackMomentums_noT4[trID] = atrack_p - m_noT4 = array('d',[sm['xtop'],sm['ytop'],sm['z'],sm['xbot'],sm['ybot'],sm['z'],sm['dist']]) - hitPosLists_noT4[trID].push_back(ROOT.TVectorD(7,m_noT4)) - if station not in stationCrossed_noT4[trID]: - stationCrossed_noT4[trID][station]=0 - stationCrossed_noT4[trID][station]+=1 - trackDigiHits_noT4[trID].append(sm['digiHit']) - - # Do fake PatRec - else: - - if not self.sTree.GetBranch("MufluxSpectrometerPoint"): - print('Fake PatRec is impossible. No "MufluxSpectrometerPoint" branch.') - return nTrack - - track_hits = {} - for sm in self.SmearedHits: - - # detID = self.digiMufluxSpectrometer[sm['digiHit']].GetDetectorID() - detID = sm['detID'] - station = int(detID/10000000) - vnb = (detID - station * 10000000) // 1000000 - is_y12 = (station == 1) * (vnb == 0) + (station == 2) * (vnb == 1) - is_stereo12 = (station == 1) * (vnb == 1) + (station == 2) * (vnb == 0) - is_34 = (station == 3) + (station == 4) - - trID = self.sTree.MufluxSpectrometerPoint[sm['digiHit']].GetTrackID() - - # PatRec - if trID not in track_hits: - atrack = {'y12': [], 'stereo12': [], '34': []} - track_hits[trID] = atrack - if is_y12: - track_hits[trID]['y12'].append(sm) - if is_stereo12: - track_hits[trID]['stereo12'].append(sm) - if is_34: - track_hits[trID]['34'].append(sm) - - # T1-4 for track fit - if trID not in hitPosLists: - hitPosLists[trID] = ROOT.std.vector('TVectorD')() - stationCrossed[trID] = {} - trackDigiHits[trID] = [] - trackMomentums[trID] = 3. - m = array('d',[sm['xtop'],sm['ytop'],sm['z'],sm['xbot'],sm['ybot'],sm['z'],sm['dist']]) - hitPosLists[trID].push_back(ROOT.TVectorD(7,m)) - if station not in stationCrossed[trID]: - stationCrossed[trID][station]=0 - stationCrossed[trID][station]+=1 - trackDigiHits[trID].append(sm['digiHit']) - - # T1-3 for track fit - if (int(detID/1000000)!=40): - if trID not in hitPosLists_noT4: - hitPosLists_noT4[trID] = ROOT.std.vector('TVectorD')() - stationCrossed_noT4[trID] = {} - trackDigiHits_noT4[trID] = [] - trackMomentums_noT4[trID] = 3. - m_noT4 = array('d',[sm['xtop'],sm['ytop'],sm['z'],sm['xbot'],sm['ybot'],sm['z'],sm['dist']]) - hitPosLists_noT4[trID].push_back(ROOT.TVectorD(7,m_noT4)) - if station not in stationCrossed_noT4[trID]: - stationCrossed_noT4[trID][station]=0 - stationCrossed_noT4[trID][station]+=1 - trackDigiHits_noT4[trID].append(sm['digiHit']) - - - # All tracks fit - for atrack in hitPosLists: - - if atrack < 0: - continue # these are hits not assigned to MC track because low E cut - pdg = 13 # only muons - if not abs(pdg)==13: - continue # only keep muons - - meas = hitPosLists[atrack] - mom_init = trackMomentums[atrack] - nM = meas.size() - - #if nM < 12 : continue # not enough hits to make a good trackfit - #comment for hits in t1-3 - #if len(stationCrossed[atrack]) < 4 : - # continue # not enough stations crossed to make a good trackfit - - charge = self.PDG.GetParticle(pdg).Charge()/(3.) - posM = ROOT.TVector3(0, 0, 0) - momM = ROOT.TVector3(0,0,mom_init * u.GeV) - # approximate covariance - covM = ROOT.TMatrixDSym(6) - resolution = self.sigma_spatial - if global_variables.withT0: - resolution = resolution*1.4 # worse resolution due to t0 estimate - for i in range(3): - covM[i][i] = resolution*resolution - covM[0][0]=resolution*resolution*100. - for i in range(3,6): - covM[i][i] = ROOT.TMath.Power(resolution / nM / ROOT.TMath.Sqrt(3), 2) - # trackrep - rep = ROOT.genfit.RKTrackRep(pdg) - # smeared start state - stateSmeared = ROOT.genfit.MeasuredStateOnPlane(rep) - rep.setPosMomCov(stateSmeared, posM, momM, covM) - # create track - seedState = ROOT.TVectorD(6) - seedCov = ROOT.TMatrixDSym(6) - rep.get6DStateCov(stateSmeared, seedState, seedCov) - theTrack = ROOT.genfit.Track(rep, seedState, seedCov) - hitCov = ROOT.TMatrixDSym(7) - hitCov[6][6] = resolution*resolution - for m in meas: - tp = ROOT.genfit.TrackPoint(theTrack) # note how the point is told which track it belongs to - measurement = ROOT.genfit.WireMeasurement(m,hitCov,1,6,tp) # the measurement is told which trackpoint it belongs to - measurement.setMaxDistance(1.85*u.cm) - # measurement.setLeftRightResolution(-1) - tp.addRawMeasurement(measurement) # package measurement in the TrackPoint - theTrack.insertPoint(tp) # add point to Track - trackCandidates_all.append([theTrack,atrack]) - - - # T1-4 track fit - for atrack in hitPosLists: - - if atrack < 0: - continue # these are hits not assigned to MC track because low E cut - pdg = 13 # only muons - if not abs(pdg)==13: - continue # only keep muons - - meas = hitPosLists[atrack] - mom_init = trackMomentums[atrack] - nM = meas.size() - - #if nM < 12 : continue # not enough hits to make a good trackfit - #comment for hits in t1-3 - if len(stationCrossed[atrack]) < 4 : - continue # not enough stations crossed to make a good trackfit - - charge = self.PDG.GetParticle(pdg).Charge()/(3.) - posM = ROOT.TVector3(0, 0, 0) - momM = ROOT.TVector3(0,0,mom_init * u.GeV) - # approximate covariance - covM = ROOT.TMatrixDSym(6) - resolution = self.sigma_spatial - if global_variables.withT0: - resolution = resolution*1.4 # worse resolution due to t0 estimate - for i in range(3): - covM[i][i] = resolution*resolution - covM[0][0]=resolution*resolution*100. - for i in range(3,6): - covM[i][i] = ROOT.TMath.Power(resolution / nM / ROOT.TMath.Sqrt(3), 2) - # trackrep - rep = ROOT.genfit.RKTrackRep(pdg) - # smeared start state - stateSmeared = ROOT.genfit.MeasuredStateOnPlane(rep) - rep.setPosMomCov(stateSmeared, posM, momM, covM) - # create track - seedState = ROOT.TVectorD(6) - seedCov = ROOT.TMatrixDSym(6) - rep.get6DStateCov(stateSmeared, seedState, seedCov) - theTrack = ROOT.genfit.Track(rep, seedState, seedCov) - hitCov = ROOT.TMatrixDSym(7) - hitCov[6][6] = resolution*resolution - for m in meas: - tp = ROOT.genfit.TrackPoint(theTrack) # note how the point is told which track it belongs to - measurement = ROOT.genfit.WireMeasurement(m,hitCov,1,6,tp) # the measurement is told which trackpoint it belongs to - measurement.setMaxDistance(1.85*u.cm) - # measurement.setLeftRightResolution(-1) - tp.addRawMeasurement(measurement) # package measurement in the TrackPoint - theTrack.insertPoint(tp) # add point to Track - trackCandidates.append([theTrack,atrack]) - - # T1-3 track fit - for atrack in hitPosLists_noT4: - - if atrack < 0: - continue # these are hits not assigned to MC track because low E cut - pdg = 13 # only muons - if not abs(pdg)==13: - continue # only keep muons - meas = hitPosLists_noT4[atrack] - mom_init = trackMomentums_noT4[atrack] - nM = meas.size() - - #if nM < 6 : continue # not enough hits to make a good trackfit - #comment for hits in t1-3 - if len(stationCrossed[atrack]) < 3 : - continue # not enough stations crossed to make a good trackfit - - charge = self.PDG.GetParticle(pdg).Charge()/(3.) - posM = ROOT.TVector3(0, 0, 0) - momM = ROOT.TVector3(0,0,mom_init * u.GeV) - # approximate covariance - covM = ROOT.TMatrixDSym(6) - resolution = self.sigma_spatial - if global_variables.withT0: - resolution = resolution*1.4 # worse resolution due to t0 estimate - for i in range(3): - covM[i][i] = resolution*resolution - covM[0][0]=resolution*resolution*100. - for i in range(3,6): - covM[i][i] = ROOT.TMath.Power(resolution / nM / ROOT.TMath.Sqrt(3), 2) - # trackrep - rep = ROOT.genfit.RKTrackRep(pdg) - # smeared start state - stateSmeared = ROOT.genfit.MeasuredStateOnPlane(rep) - rep.setPosMomCov(stateSmeared, posM, momM, covM) - # create track - seedState = ROOT.TVectorD(6) - seedCov = ROOT.TMatrixDSym(6) - rep.get6DStateCov(stateSmeared, seedState, seedCov) - theTrack = ROOT.genfit.Track(rep, seedState, seedCov) - hitCov = ROOT.TMatrixDSym(7) - hitCov[6][6] = resolution*resolution - for m in meas: - tp = ROOT.genfit.TrackPoint(theTrack) # note how the point is told which track it belongs to - measurement = ROOT.genfit.WireMeasurement(m,hitCov,1,6,tp) # the measurement is told which trackpoint it belongs to - measurement.setMaxDistance(1.85*u.cm) - # measurement.setLeftRightResolution(-1) - tp.addRawMeasurement(measurement) # package measurement in the TrackPoint - theTrack.insertPoint(tp) # add point to Track - trackCandidates_noT4.append([theTrack,atrack]) - - if global_variables.debug: - - z_hits_y = [] - x_hits_y = [] - for ahit in self.SmearedHits: - detID = ahit['detID'] - station = int(detID/10000000) - vnb = (detID - station * 10000000) // 1000000 - is_y12 = (station == 1) * (vnb == 0) + (station == 2) * (vnb == 1) - is_stereo12 = (station == 1) * (vnb == 1) + (station == 2) * (vnb == 0) - is_34 = (station == 3) + (station == 4) - if is_y12 or is_34: - z_hits_y.append(ahit['z']) - x_hits_y.append(ahit['xtop']) - plt.figure(figsize=(22, 6)) - plt.subplot(1, 2, 1) - plt.scatter(z_hits_y, x_hits_y, color='b') - - z_tagger_hits_y = [] - x_tagger_hits_y = [] - for ahit in self.TaggerHits: - detID = ahit['detID'] - station = detID // 10000 - direction = (detID - station * 10000) // 1000 - if direction == 1 or detID < 10: - z_tagger_hits_y.append(ahit['z']) - x_tagger_hits_y.append((ahit['xtop'] + ahit['xbot']) / 2.) - plt.scatter(z_tagger_hits_y, x_tagger_hits_y, color='b') - - for i_track in track_hits: - atrack = track_hits[i_track] - atrack_y12 = atrack['y12'] - atrack_stereo12 = atrack['stereo12'] - atrack_34 = atrack['34'] - atrack_tagger = atrack['y_tagger'] - atrack_y = list(atrack_y12) + list(atrack_34) + list(atrack_tagger) - z_atrack_y = [ahit['z'] for ahit in atrack_y] - x_atrack_y = [ahit['xtop'] for ahit in atrack_y] - plt.scatter(z_atrack_y, x_atrack_y, color=np.random.rand(3,)) - - plt.xlabel('z') - plt.ylabel('x') - plt.xlim(0, 1200) - plt.ylim(-150, 150) - - - plt.subplot(1, 2, 2) - for ahit in self.SmearedHits: - detID = ahit['detID'] - station = int(detID/10000000) - vnb = (detID - station * 10000000) // 1000000 - is_y12 = (station == 1) * (vnb == 0) + (station == 2) * (vnb == 1) - is_stereo12 = (station == 1) * (vnb == 1) + (station == 2) * (vnb == 0) - is_34 = (station == 3) + (station == 4) - if is_y12 or is_stereo12: - xbot, xtop = ahit['xbot'], ahit['xtop'] - ybot, ytop = ahit['ybot'], ahit['ytop'] - plt.plot([ybot, ytop], [xbot, xtop], color='b') - - for i_track in track_hits: - atrack = track_hits[i_track] - atrack_y12 = atrack['y12'] - atrack_stereo12 = atrack['stereo12'] - atrack_34 = atrack['34'] - atrack_tagger = atrack['y_tagger'] - atrack_12 = list(atrack_y12) + list(atrack_stereo12) - color = np.random.rand(3,) - for ahit in atrack_12: - xbot, xtop = ahit['xbot'], ahit['xtop'] - ybot, ytop = ahit['ybot'], ahit['ytop'] - plt.plot([ybot, ytop], [xbot, xtop], color=color) - - plt.xlabel('y') - plt.ylabel('x') - plt.xlim(-50, 50) - plt.ylim(-50, 50) - - plt.savefig('pics/event_'+str(self.iEvent)+'.pdf', format='pdf', bbox_inches='tight') - plt.clf() - plt.close() - - - - - - # Metrics - if self.sTree.GetBranch("MufluxSpectrometerPoint"): - - - (n_tracks, n_recognized, n_clones, n_ghosts, n_others) = self.target_metrics(track_hits) - h['Reco_target'].Fill("N total", n_tracks) - h['Reco_target'].Fill("N recognized tracks", n_recognized) - h['Reco_target'].Fill("N clones", n_clones) - h['Reco_target'].Fill("N ghosts", n_ghosts) - h['Reco_target'].Fill("N others", n_others) - - (n_tracks, n_recognized, n_clones, n_ghosts, n_others) = self.muon_metrics(track_hits) - h['Reco_muon'].Fill("N total", n_tracks) - h['Reco_muon'].Fill("N recognized tracks", n_recognized) - h['Reco_muon'].Fill("N clones", n_clones) - h['Reco_muon'].Fill("N ghosts", n_ghosts) - h['Reco_muon'].Fill("N others", n_others) - - (n_tracks, n_recognized, n_clones, n_ghosts, n_others) = self.muon_with_tagger_metrics(track_hits) - h['Reco_muon_with_tagger'].Fill("N total", n_tracks) - h['Reco_muon_with_tagger'].Fill("N recognized tracks", n_recognized) - h['Reco_muon_with_tagger'].Fill("N clones", n_clones) - h['Reco_muon_with_tagger'].Fill("N ghosts", n_ghosts) - h['Reco_muon_with_tagger'].Fill("N others", n_others) - - (n_tracks, n_recognized, n_clones, n_ghosts, n_others) = self.muon_metrics_vs_p(track_hits) - - (n_tracks, n_recognized, n_clones, n_ghosts, n_others) = self.all_tracks_metrics(track_hits) - h['Reco_all_tracks'].Fill("N total", n_tracks) - h['Reco_all_tracks'].Fill("N recognized tracks", n_recognized) - h['Reco_all_tracks'].Fill("N clones", n_clones) - h['Reco_all_tracks'].Fill("N ghosts", n_ghosts) - h['Reco_all_tracks'].Fill("N others", n_others) - - - (n_tracks_y12, n_recognized_y12, n_clones_y12, n_ghosts_y12, n_others_y12, - n_tracks_stereo12, n_recognized_stereo12, n_clones_stereo12, n_ghosts_stereo12, n_others_stereo12, - n_tracks_34, n_recognized_34, n_clones_34, n_ghosts_34, n_others_34) = self.recognition_metrics(track_hits) - - h['NTrueTracks'].Fill("Stations 1&2, Y views", n_tracks_y12) - h['NTrueTracks'].Fill("Stations 1&2, Stereo views", n_tracks_stereo12) - h['NTrueTracks'].Fill("Stations 3&4", n_tracks_34) - - h['Reco_y12'].Fill("N total", n_tracks_y12) - h['Reco_y12'].Fill("N recognized tracks", n_recognized_y12) - h['Reco_y12'].Fill("N clones", n_clones_y12) - h['Reco_y12'].Fill("N ghosts", n_ghosts_y12) - h['Reco_y12'].Fill("N others", n_others_y12) - - h['Reco_stereo12'].Fill("N total", n_tracks_stereo12) - h['Reco_stereo12'].Fill("N recognized tracks", n_recognized_stereo12) - h['Reco_stereo12'].Fill("N clones", n_clones_stereo12) - h['Reco_stereo12'].Fill("N ghosts", n_ghosts_stereo12) - h['Reco_stereo12'].Fill("N others", n_others_stereo12) - - h['Reco_34'].Fill("N total", n_tracks_34) - h['Reco_34'].Fill("N recognized tracks", n_recognized_34) - h['Reco_34'].Fill("N clones", n_clones_34) - h['Reco_34'].Fill("N ghosts", n_ghosts_34) - h['Reco_34'].Fill("N others", n_others_34) - - - (n_tracks_y12, n_recognized_y12, n_clones_y12, n_ghosts_y12, n_others_y12, - n_tracks_stereo12, n_recognized_stereo12, n_clones_stereo12, n_ghosts_stereo12, n_others_stereo12, - n_tracks_34, n_recognized_34, n_clones_34, n_ghosts_34, n_others_34) = self.recognition_metrics(track_hits, mode='3hits') - - h['NTrueTracks_3hits'].Fill("Stations 1&2, Y views", n_tracks_y12) - h['NTrueTracks_3hits'].Fill("Stations 1&2, Stereo views", n_tracks_stereo12) - h['NTrueTracks_3hits'].Fill("Stations 3&4", n_tracks_34) - - h['Reco_y12_3hits'].Fill("N total", n_tracks_y12) - h['Reco_y12_3hits'].Fill("N recognized tracks", n_recognized_y12) - h['Reco_y12_3hits'].Fill("N clones", n_clones_y12) - h['Reco_y12_3hits'].Fill("N ghosts", n_ghosts_y12) - h['Reco_y12_3hits'].Fill("N others", n_others_y12) - - h['Reco_stereo12_3hits'].Fill("N total", n_tracks_stereo12) - h['Reco_stereo12_3hits'].Fill("N recognized tracks", n_recognized_stereo12) - h['Reco_stereo12_3hits'].Fill("N clones", n_clones_stereo12) - h['Reco_stereo12_3hits'].Fill("N ghosts", n_ghosts_stereo12) - h['Reco_stereo12_3hits'].Fill("N others", n_others_stereo12) - - h['Reco_34_3hits'].Fill("N total", n_tracks_34) - h['Reco_34_3hits'].Fill("N recognized tracks", n_recognized_34) - h['Reco_34_3hits'].Fill("N clones", n_clones_34) - h['Reco_34_3hits'].Fill("N ghosts", n_ghosts_34) - h['Reco_34_3hits'].Fill("N others", n_others_34) - - - (n_tracks_y12, n_recognized_y12, n_clones_y12, n_ghosts_y12, n_others_y12, - n_tracks_stereo12, n_recognized_stereo12, n_clones_stereo12, n_ghosts_stereo12, n_others_stereo12, - n_tracks_34, n_recognized_34, n_clones_34, n_ghosts_34, n_others_34) = self.recognition_metrics(track_hits, mode='Tr4') - - h['NTrueTracks_Tr4'].Fill("Stations 1&2, Y views", n_tracks_y12) - h['NTrueTracks_Tr4'].Fill("Stations 1&2, Stereo views", n_tracks_stereo12) - h['NTrueTracks_Tr4'].Fill("Stations 3&4", n_tracks_34) - - h['Reco_y12_Tr4'].Fill("N total", n_tracks_y12) - h['Reco_y12_Tr4'].Fill("N recognized tracks", n_recognized_y12) - h['Reco_y12_Tr4'].Fill("N clones", n_clones_y12) - h['Reco_y12_Tr4'].Fill("N ghosts", n_ghosts_y12) - h['Reco_y12_Tr4'].Fill("N others", n_others_y12) - - h['Reco_stereo12_Tr4'].Fill("N total", n_tracks_stereo12) - h['Reco_stereo12_Tr4'].Fill("N recognized tracks", n_recognized_stereo12) - h['Reco_stereo12_Tr4'].Fill("N clones", n_clones_stereo12) - h['Reco_stereo12_Tr4'].Fill("N ghosts", n_ghosts_stereo12) - h['Reco_stereo12_Tr4'].Fill("N others", n_others_stereo12) - - h['Reco_34_Tr4'].Fill("N total", n_tracks_34) - h['Reco_34_Tr4'].Fill("N recognized tracks", n_recognized_34) - h['Reco_34_Tr4'].Fill("N clones", n_clones_34) - h['Reco_34_Tr4'].Fill("N ghosts", n_ghosts_34) - h['Reco_34_Tr4'].Fill("N others", n_others_34) - - - - n_true_track_hits_y12 = {} - n_true_track_hits_stereo12 = {} - n_true_track_hits_34 = {} - for ahit in self.sTree.MufluxSpectrometerPoint: - - track_id = ahit.GetTrackID() - - detID = ahit.GetDetectorID() - statnb = detID // 10000000 - vnb = (detID - statnb * 10000000) // 1000000 - - is_y12 = (statnb == 1) * (vnb == 0) + (statnb == 2) * (vnb == 1) - is_stereo12 = (statnb == 1) * (vnb == 1) + (statnb == 2) * (vnb == 0) - is_34 = (statnb == 3) + (statnb == 4) - - if is_y12: - if track_id in n_true_track_hits_y12: - n_true_track_hits_y12[track_id] += 1 - else: - n_true_track_hits_y12[track_id] = 1 - if is_stereo12: - if track_id in n_true_track_hits_stereo12: - n_true_track_hits_stereo12[track_id] += 1 - else: - n_true_track_hits_stereo12[track_id] = 1 - if is_34: - if track_id in n_true_track_hits_34: - n_true_track_hits_34[track_id] += 1 - else: - n_true_track_hits_34[track_id] = 1 - - for key in n_true_track_hits_y12.keys(): - n = n_true_track_hits_y12[key] - h['NHits_true_y12'].Fill(n) - for key in n_true_track_hits_stereo12.keys(): - n = n_true_track_hits_stereo12[key] - h['NHits_true_stereo12'].Fill(n) - for key in n_true_track_hits_34.keys(): - n = n_true_track_hits_34[key] - h['NHits_true_34'].Fill(n) - - for i_track in track_hits.keys(): - atrack = track_hits[i_track] - atrack_y12 = atrack['y12'] - atrack_stereo12 = atrack['stereo12'] - atrack_34 = atrack['34'] - - if len(atrack_y12) > 0: - h['NHits_reco_y12'].Fill(len(atrack_y12)) - if len(atrack_stereo12) > 0: - h['NHits_reco_stereo12'].Fill(len(atrack_stereo12)) - if len(atrack_34) > 0: - h['NHits_reco_34'].Fill(len(atrack_34)) - - - - - for entry in trackCandidates: - #check - #print "fitting with stereo" - atrack = entry[1] - theTrack = entry[0] - if not theTrack.checkConsistency(): - print('Problem with track before fit, not consistent',atrack,theTrack) - continue - # do the fit - try: self.fitter.processTrack(theTrack) # processTrackWithRep(theTrack,rep,True) - except: - print("genfit failed to fit track") - continue - #check - if not theTrack.checkConsistency(): - #print 'Problem with track after fit, not consistent',atrack,theTrack - continue - fitStatus = theTrack.getFitStatus() - nmeas = fitStatus.getNdf() - chi2 = fitStatus.getChi2()/nmeas - pvalue = fitStatus.getPVal() - #if pvalue < 0.05: - # print "P value too low. Rejecting track." - # continue - h['nmeas'].Fill(nmeas) - h['chi2'].Fill(chi2) - h['p-value'].Fill(pvalue) - try: - fittedState = theTrack.getFittedState() - fittedMom = fittedState.getMomMag() - h['p-fittedtracks'].Fill(fittedMom) - h['1/p-fittedtracks'].Fill(1./fittedMom) - Px,Py,Pz = fittedState.getMom().x(),fittedState.getMom().y(),fittedState.getMom().z() - P = fittedMom - Pt = ROOT.TMath.Sqrt(Px*Px+Py*Py) - h['p/pt'].Fill(P,Pt) - h['pt-fittedtracks'].Fill(Pt) - h['1/pt-fittedtracks'].Fill(1./Pt) - - if self.sTree.GetBranch("MufluxSpectrometerPoint"): - - atrack_ids = [] - for digi_hit in trackDigiHits[atrack]: - ahit = self.sTree.MufluxSpectrometerPoint[digi_hit] - atrack_ids.append(ahit.GetTrackID()) - frac, tmax = self.fracMCsame(atrack_ids) - Ptruth,Ptruthx,Ptruthy,Ptruthz = self.getPtruthFirst(tmax) # MC Truth - Pgun,Pgunx,Pguny,Pgunz = self.getPtruthAtOrigin(tmax) # MC Truth - Pttruth = ROOT.TMath.Sqrt(Ptruthx*Ptruthx+Ptruthy*Ptruthy) - h['p/pt_truth'].Fill(Ptruth,Pttruth) - perr = (P - Ptruth) / Ptruth - pterr = (Pt - Pttruth) / Pttruth - h['p_rel_error'].Fill(perr) - h['pt_rel_error'].Fill(pterr) - - mom_init = trackMomentums[atrack] - print("P_precalc, P_truth: ", mom_init, Ptruth) - - if Pz !=0: - pxpzfitted = Px/Pz - pypzfitted = Py/Pz - if Ptruthz !=0: - pxpztrue = Ptruthx/Ptruthz - pypztrue = Ptruthy/Ptruthz - h['Px/Pzfitted'].Fill(pxpzfitted) - h['Py/Pzfitted'].Fill(pypzfitted) - h['Px/Pztrue'].Fill(pxpztrue) - h['Py/Pztrue'].Fill(pypztrue) - h['Px/Pzfitted-Px/Pztruth'].Fill(Ptruth,pxpzfitted-pxpztrue) - h['Py/Pzfitted-Py/Pztruth'].Fill(Ptruth,pypzfitted-pypztrue) - h['ptruth'].Fill(Ptruth) - delPOverP = (P/Ptruth)-1 - invdelPOverP = (Ptruth/P)-1 - if 1==0: - if invdelPOverP < -0.8: - print("invdelPOverP = ",invdelPOverP) - print("Ptruth =",Ptruth," Pfitted =",P) - for n in range(hitPosLists[atrack].size()): - print("hit=",n," x(top) ",hitPosLists[atrack][n][0]," y(top) ",hitPosLists[atrack][n][1]," z ",hitPosLists[atrack][n][2]," x(bot) ",hitPosLists[atrack][n][3]," y(bot) ", hitPosLists[atrack][n][4], " dist ", hitPosLists[atrack][n][6]) - nMufluxHits = self.sTree.MufluxSpectrometerPoint.GetEntriesFast() - for i in range(nMufluxHits): - MufluxHit = self.sTree.MufluxSpectrometerPoint[i] - if ((hitPosLists[atrack][n][0]+1.8 > MufluxHit.GetX()) or(hitPosLists[atrack][n][3]+1.8 > MufluxHit.GetX())) and ((hitPosLists[atrack][n][0]-1.8MufluxHit.GetZ()) and (hitPosLists[atrack][n][2]-1. 200: - print("Too many hits.") - return track_hits - - min_hits = 3 - max_shared_hits = 0 - - # Split hits into ortogonal planes - TaggerHits_x, TaggerHits_y = split_hits(TaggerHits) - - # PatRec in zx plane - tracks_zx = pat_rec_plane(TaggerHits_x, coord='x', min_hits=min_hits, max_shared_hits=max_shared_hits) - - # PatRec in zy plane - tracks_zy = pat_rec_plane(TaggerHits_y, coord='y', min_hits=min_hits, max_shared_hits=max_shared_hits) - - # Combine tracks - track_combinations = combine_tracks(tracks_zx, tracks_zy) - - # Generate output - for i_track in range(len(track_combinations)): - - atrack = track_combinations[i_track] - - if len(atrack['hits_x']) >= min_hits and len(atrack['hits_y']) >= min_hits: - track_hits[i_track] = {} - track_hits[i_track]['hits_x'] = sort_hits(atrack['hits_x']) - track_hits[i_track]['hits_y'] = sort_hits(atrack['hits_y']) - - - if debug: - print("From MuonTaggerPatRec.py: ") - print("tracks_zx: ") - print(len(tracks_zx)) - print("tracks_zy: ") - print(len(tracks_zy)) - print("track_combinations: ") - print(len(track_combinations)) - print("Recognized tracks: ") - for i_track in track_hits.keys(): - print("Track ", i_track) - print("Hits_x: ", [hit['digiHit'] for hit in track_hits[i_track]['hits_x']]) - print("Hits_y: ", [hit['digiHit'] for hit in track_hits[i_track]['hits_y']]) - - return track_hits - - - -def finalize(): - pass - - -################################################################################################### - - -def split_hits(TaggerHits): - - TaggerHits_x = [] - TaggerHits_y = [] - for ahit in TaggerHits: - if abs(ahit['xtop'] - ahit['xbot']) < 10: - TaggerHits_x.append(ahit) - if abs(ahit['ytop'] - ahit['ybot']) < 10: - TaggerHits_y.append(ahit) - - return TaggerHits_x, TaggerHits_y - - -def pat_rec_plane(TaggerHits, coord='x', min_hits=3, max_shared_hits=2): - - long_tracks = [] - - # Take 2 hits as a track seed - for ahit1 in TaggerHits: - for ahit2 in TaggerHits: - - if ahit1['z'] >= ahit2['z']: - continue - if ahit1['detID'] == ahit2['detID']: - continue - - x1 = 0.5 * (ahit1[coord+'top'] + ahit1[coord+'bot']) - x2 = 0.5 * (ahit2[coord+'top'] + ahit2[coord+'bot']) - z1 = ahit1['z'] - z2 = ahit2['z'] - layer1 = ahit1['detID'] // 10000 - layer2 = ahit2['detID'] // 10000 - - k_bin = 1. * (x2 - x1) / (z2 - z1) - b_bin = x1 - k_bin * z1 - - if abs(k_bin) > 1: - continue - - atrack = {} - atrack['hits_'+coord] = [ahit1, ahit2] - atrack[coord+'_'+coord] = [x1, x2] - atrack['z_'+coord] = [z1, z2] - atrack['layer'] = [layer1, layer2] - - # Add new hits to the seed - for ahit3 in TaggerHits: - - if ahit3['detID'] == ahit1['detID'] or ahit3['detID'] == ahit2['detID']: - continue - - x3 = 0.5 * (ahit3[coord+'top'] + ahit3[coord+'bot']) - z3 = ahit3['z'] - layer3 = ahit3['detID'] // 10000 - - if layer3 in atrack['layer']: - continue - #pass - - in_bin = hit_in_window(z3, x3, k_bin, b_bin, window_width=5.0) - if in_bin: - atrack['hits_'+coord].append(ahit3) - atrack['z_'+coord].append(z3) - atrack[coord+'_'+coord].append(x3) - atrack['layer'].append(layer3) - - if len(atrack['hits_'+coord]) >= min_hits: - long_tracks.append(atrack) - - # Reduce number of clones - short_tracks = reduce_clones(long_tracks, coord, min_hits, max_shared_hits) - - # Fit tracks y12 - for atrack in short_tracks: - [atrack['k_'+coord], atrack['b_'+coord]] = np.polyfit(atrack['z_'+coord], atrack[coord+'_'+coord], deg=1) - - return short_tracks - - -def combine_tracks(tracks_zx, tracks_zy): - - track_combinations = [] - for atrack_zx in tracks_zx: - for atrack_zy in tracks_zy: - atrack_comb = atrack_zx.copy() - atrack_comb.update(atrack_zy) - track_combinations.append(atrack_comb) - - return track_combinations - - - -def reduce_clones(long_tracks, coord='x', min_hits=3, max_shared_hits=2): - - # Remove clones - used_hits = [] - short_tracks = [] - n_hits = [len(atrack['hits_'+coord]) for atrack in long_tracks] - - for i_track in np.argsort(n_hits)[::-1]: - - atrack = long_tracks[i_track] - n_shared_hits = 0 - - for i_hit in range(len(atrack['hits_'+coord])): - ahit = atrack['hits_'+coord][i_hit] - if ahit['digiHit'] in used_hits: #digiHit - n_shared_hits += 1 - - if len(atrack['hits_'+coord]) >= min_hits and n_shared_hits <= max_shared_hits: - short_tracks.append(atrack) - for ahit in atrack['hits_'+coord]: - used_hits.append(ahit['digiHit']) #digiHit - - return short_tracks - - -def hit_in_window(x, y, k_bin, b_bin, window_width=1.): - """ - Counts hits in a bin of track parameter space (b, k). - - Parameters - --------- - x : array-like - Array of x coordinates of hits. - y : array-like - Array of x coordinates of hits. - k_bin : float - Track parameter: y = k_bin * x + b_bin - b_bin : float - Track parameter: y = k_bin * x + b_bin - - Return - ------ - track_inds : array-like - Hit indexes of a track: [ind1, ind2, ...] - """ - - - y_approx = k_bin * x + b_bin - - flag = False - if np.abs(y_approx - y) <= window_width: - flag = True - - return flag - - -def sort_hits(hits): - - sorted_hits = [] - hits_z = [ahit['z'] for ahit in hits] - sort_index = np.argsort(hits_z) - for i_hit in sort_index: - sorted_hits.append(hits[i_hit]) - - return sorted_hits diff --git a/python/MufluxPatRec.py b/python/MufluxPatRec.py deleted file mode 100644 index c369cc7a46..0000000000 --- a/python/MufluxPatRec.py +++ /dev/null @@ -1,508 +0,0 @@ -__author__ = 'Mikhail Hushchyn' - -import numpy as np - -# Globals -ReconstructibleMCTracks = [] -theTracks = [] - - -def initialize(fgeo): - pass - - -def execute(SmearedHits, TaggerHits, withNTaggerHits, withDist2Wire, debug=0): - """ - Main function of track pattern recognition. - - Parameters: - ----------- - SmearedHits : list - Smeared hits. SmearedHits = [{'digiHit': key, - 'xtop': xtop, 'ytop': ytop, 'z': ztop, - 'xbot': xbot, 'ybot': ybot, - 'dist': dist2wire, 'detID': detID}, {...}, ...] - """ - - # withDist2Wire = False - # withNTaggerHits = 0 - # TaggerHits = [] - - - fittedtrackids = [] - track_hits = {} - if len(SmearedHits) > 100: - print("Too large hits in the event!") - return track_hits - - min_hits = 3 - max_shared_hits = 2 - - #### Separate hits - SmearedHits_y12, SmearedHits_stereo12, SmearedHits_34 = hits_split(SmearedHits) - - - #### PatRec in 12y - short_tracks_y12 = pat_rec_y_views(SmearedHits_y12, min_hits, max_shared_hits) - - - #### PatRec in stereo12 - short_tracks_12 = pet_rec_stereo_views(SmearedHits_stereo12, short_tracks_y12, min_hits) - - - #### PatRec in 34 - short_tracks_34 = pat_rec_y_views(SmearedHits_34, min_hits, max_shared_hits) - - - #### Combine tracks - z_center=350.75 - track_combinations = combine_tracks_before_and_after_the_magnet(short_tracks_12, short_tracks_34, z_center) - - - #### MuID is in other script - - - # Prepare output of PatRec - track_hits = {} - for i_track in range(len(track_combinations)): - - atrack = track_combinations[i_track] - - hits_y12 = atrack['hits_y12'] - hits_stereo12 = atrack['hits_stereo12'] - hits_34 = atrack['hits_y34'] - p = atrack['p'] - - [k, b] = np.polyfit(atrack['z_y12'], atrack['x_y12'], deg=1) - x_in_magnet = k * z_center + b - - [k, b] = np.polyfit(atrack['z_stereo12'], atrack['y_stereo12'], deg=1) - y_in_magnet = k * z_center + b - - if len(hits_y12) >= min_hits and len(hits_stereo12) >= min_hits and len(hits_34) >= min_hits: - - atrack = {'y12': sort_hits(hits_y12), - 'stereo12': sort_hits(hits_stereo12), - '34': sort_hits(hits_34), - 'y_tagger': [], - 'p': p, - 'x_in_magnet': x_in_magnet, - 'y_in_magnet': y_in_magnet} - track_hits[i_track] = atrack - - - if debug: - print("Recognized tracks:") - for i_track in track_hits.keys(): - atrack = track_hits[i_track] - print("Track ", i_track) - print("Z_y12", [str(np.around(hit['z'], 2)) for hit in atrack['y12']]) - print("X_y12", [str(np.around(hit['xtop'], 2)) for hit in atrack['y12']]) - print("Z_stereo12", [str(np.around(hit['z'], 2)) for hit in atrack['stereo12']]) - print("X_stereo12", [str(np.around(hit['xtop'], 2)) for hit in atrack['stereo12']]) - print("Z_34", [str(np.around(hit['z'], 2)) for hit in atrack['34']]) - print("X_34", [str(np.around(hit['xtop'], 2)) for hit in atrack['34']]) - - return track_hits - - -def finalize(): - pass - -################################################################################################### - -def hits_split(SmearedHits): - - # Separate hits - SmearedHits_12y = [] - SmearedHits_12stereo = [] - SmearedHits_34 = [] - for i_hit in range(len(SmearedHits)): - ahit = SmearedHits[i_hit] - detID = ahit['detID'] - statnb, vnb, pnb, lnb, snb = decodeDetectorID(detID) - is_y12 = (statnb == 1) * (vnb == 0) + (statnb == 2) * (vnb == 1) - is_stereo12 = (statnb == 1) * (vnb == 1) + (statnb == 2) * (vnb == 0) - is_34 = (statnb == 3) + (statnb == 4) - if is_y12: - SmearedHits_12y.append(ahit) - if is_stereo12: - SmearedHits_12stereo.append(ahit) - if is_34: - SmearedHits_34.append(ahit) - - return SmearedHits_12y, SmearedHits_12stereo, SmearedHits_34 - - -def pat_rec_y_views(SmearedHits, min_hits=3, max_shared_hits=2): - - long_tracks = [] - - # Take 2 hits as a track seed - for ahit1 in SmearedHits: - for ahit2 in SmearedHits: - - if ahit1['z'] >= ahit2['z']: - continue - if ahit1['detID'] == ahit2['detID']: - continue - - x1 = ahit1['xtop'] - x2 = ahit2['xtop'] - z1 = ahit1['z'] - z2 = ahit2['z'] - layer1 = ahit1['detID'] // 10000 - layer2 = ahit2['detID'] // 10000 - - k_bin = 1. * (x2 - x1) / (z2 - z1) - b_bin = x1 - k_bin * z1 - - if abs(k_bin) > 1: - continue - - atrack = {} - atrack['hits_y'] = [ahit1, ahit2] - atrack['x_y'] = [x1, x2] - atrack['z_y'] = [z1, z2] - atrack['layer'] = [layer1, layer2] - - # Add new hits to the seed - for ahit3 in SmearedHits: - - if ahit3['detID'] == ahit1['detID'] or ahit3['detID'] == ahit2['detID']: - continue - - x3 = ahit3['xtop'] - z3 = ahit3['z'] - layer3 = ahit3['detID'] // 10000 - - if layer3 in atrack['layer']: - continue - - in_bin = hit_in_window(z3, x3, k_bin, b_bin, window_width=3.0)#3 - if in_bin: - atrack['hits_y'].append(ahit3) - atrack['z_y'].append(z3) - atrack['x_y'].append(x3) - atrack['layer'].append(layer3) - - if len(atrack['hits_y']) >= min_hits: - long_tracks.append(atrack) - - # Reduce number of clones - short_tracks = reduce_clones(long_tracks, min_hits, max_shared_hits) - - # Fit tracks y12 - for atrack in short_tracks: - [atrack['k_y'], atrack['b_y']] = np.polyfit(atrack['z_y'], atrack['x_y'], deg=1) - - return short_tracks - - -def pet_rec_stereo_views(SmearedHits_stereo, short_tracks_y, min_hits=3): - - ### PatRec in stereo12 - long_tracks_stereo = [] - short_tracks_stereo = [] - used_hits = [] - # deg = np.deg2rad(ShipGeo['MufluxSpectrometer']['ViewAngle']) - - for i_track_y in range(len(short_tracks_y)): - - atrack_y = short_tracks_y[i_track_y] - k_y = atrack_y['k_y'] - b_y = atrack_y['b_y'] - - # Get hit zx projections - for ahit in SmearedHits_stereo: - y_center = get_zy_projection(ahit['z'], ahit['xtop'], ahit['ytop'], ahit['xbot'], ahit['ybot'], k_y, b_y) - ahit['zy_projection'] = y_center - - - temp_tracks_stereo = [] - - for ahit1 in SmearedHits_stereo: - for ahit2 in SmearedHits_stereo: - - if ahit1['z'] >= ahit2['z']: - continue - if ahit1['detID'] == ahit2['detID']: - continue - if ahit1['digiHit'] in used_hits: - continue - if ahit2['digiHit'] in used_hits: - continue - - y1_center = ahit1['zy_projection'] - y2_center = ahit2['zy_projection'] - - if abs(y1_center ) > 70 or abs(y2_center ) > 70: - continue - - y1 = y1_center - y2 = y2_center - z1 = ahit1['z'] - z2 = ahit2['z'] - layer1 = ahit1['detID'] // 10000 - layer2 = ahit2['detID'] // 10000 - - k_bin = 1. * (y2 - y1) / (z2 - z1) - b_bin = y1 - k_bin * z1 - - atrack = {} - atrack['hits_stereo'] = [ahit1, ahit2] - atrack['y_stereo'] = [y1, y2] - atrack['z_stereo'] = [z1, z2] - atrack['layer'] = [layer1, layer2] - - for ahit3 in SmearedHits_stereo: - - if ahit3['digiHit'] == ahit1['digiHit'] or ahit3['digiHit'] == ahit2['digiHit']: - continue - if ahit3['digiHit'] in used_hits: - continue - - y3_center = ahit3['zy_projection'] - z3 = ahit3['z'] - layer3 = ahit3['detID'] // 10000 - - if abs(y3_center) > 70: - continue - - if layer3 in atrack['layer']: - continue - - in_bin = hit_in_window(z3, y3_center, k_bin, b_bin, window_width=10.0)#10.0 - if in_bin: - atrack['hits_stereo'].append(ahit3) - atrack['z_stereo'].append(z3) - atrack['y_stereo'].append(y3_center) - atrack['layer'].append(layer3) - - if len(atrack['hits_stereo']) >= min_hits: - temp_tracks_stereo.append(atrack) - long_tracks_stereo.append(atrack) - - # Remove clones - max_track = None - max_n_hits = -999 - - for atrack in temp_tracks_stereo: - if len(atrack['hits_stereo']) > max_n_hits: - max_track = atrack - max_n_hits = len(atrack['hits_stereo']) - - if max_track is not None: - atrack = {} - atrack['hits_y'] = atrack_y['hits_y'] - atrack['z_y'] = atrack_y['z_y'] - atrack['x_y'] = atrack_y['x_y'] - atrack['k_y'] = atrack_y['k_y'] - atrack['b_y'] = atrack_y['b_y'] - atrack['hits_stereo'] = max_track['hits_stereo'] - atrack['z_stereo'] = max_track['z_stereo'] - atrack['y_stereo'] = max_track['y_stereo'] - - short_tracks_stereo.append(atrack) - for ahit in max_track['hits_stereo']: - #used_hits.append(ahit['digiHit']) - pass - - return short_tracks_stereo - - -def combine_tracks_before_and_after_the_magnet(short_tracks_12, short_tracks_34, z_center): - - # Combine track y12 and 34 and momentum calculation - i_track_12 = [] - i_track_34 = [] - deltas_y = [] - momentums = [] - for i_12 in range(len(short_tracks_12)): - atrack_12 = short_tracks_12[i_12] - y_center_12 = atrack_12['k_y'] * z_center + atrack_12['b_y'] - alpha_12 = np.arctan(atrack_12['k_y']) - for i_34 in range(len(short_tracks_34)): - atrack_34 = short_tracks_34[i_34] - y_center_34 = atrack_34['k_y'] * z_center + atrack_34['b_y'] - alpha_34 = np.arctan(atrack_34['k_y']) - i_track_12.append(i_12) - i_track_34.append(i_34) - deltas_y.append(abs(y_center_12 - y_center_34)) - momentums.append(1.03 / (alpha_12 - alpha_34)) - - max_dy = 50 - used_12 = [] - used_34 = [] - track_combinations = [] - for i in np.argsort(deltas_y): - dy = deltas_y[i] - mom = momentums[i] - i_12 = i_track_12[i] - i_34 = i_track_34[i] - if dy < max_dy: - if i_12 not in used_12: - if i_34 not in used_34: - atrack = {} - for key in short_tracks_12[i_12].keys(): - atrack[key+'12'] = short_tracks_12[i_12][key] - for key in short_tracks_34[i_34].keys(): - atrack[key+'34'] = short_tracks_34[i_34][key] - atrack['p'] = abs(mom) - track_combinations.append(atrack) - #used_12.append(i_12) - #used_34.append(i_34) - - return track_combinations - - -def reduce_clones(long_tracks, min_hits=3, max_shared_hits=2): - - # Remove clones in y12 - used_hits = [] - short_tracks = [] - n_hits = [len(atrack['hits_y']) for atrack in long_tracks] - - for i_track in np.argsort(n_hits)[::-1]: - - atrack = long_tracks[i_track] - n_shared_hits = 0 - - for i_hit in range(len(atrack['hits_y'])): - ahit = atrack['hits_y'][i_hit] - if ahit['digiHit'] in used_hits: #digiHit - n_shared_hits += 1 - - if len(atrack['hits_y']) >= min_hits and n_shared_hits <= max_shared_hits: - short_tracks.append(atrack) - for ahit in atrack['hits_y']: - used_hits.append(ahit['digiHit']) #digiHit - - return short_tracks - - -################################################################################################### - - -def decodeDetectorID(detID): - """ - Decodes detector ID. - - Parameters - ---------- - detID : int or array-like - Detector ID values. - - Returns - ------- - statnb : int or array-like - Station numbers. - vnb : int or array-like - View numbers. - pnb : int or array-like - Plane numbers. - lnb : int or array-like - Layer numbers. - snb : int or array-like - Straw tube numbers. - """ - - statnb = detID // 10000000 - vnb = (detID - statnb * 10000000) // 1000000 - pnb = (detID - statnb * 10000000 - vnb * 1000000) // 100000 - lnb = (detID - statnb * 10000000 - vnb * 1000000 - pnb * 100000) // 10000 - snb = detID - statnb * 10000000 - vnb * 1000000 - pnb * 100000 - lnb * 10000 - 2000 - - return statnb, vnb, pnb, lnb, snb - - -def hit_in_bin(x, y, k_bin, b_bin, k_size, b_size): - """ - Counts hits in a bin of track parameter space (b, k). - - Parameters - --------- - x : array-like - Array of x coordinates of hits. - y : array-like - Array of x coordinates of hits. - k_bin : float - Track parameter: y = k_bin * x + b_bin - b_bin : float - Track parameter: y = k_bin * x + b_bin - - Return - ------ - track_inds : array-like - Hit indexes of a track: [ind1, ind2, ...] - """ - - - b_left = y - (k_bin - 0.5 * k_size) * x - b_right = y - (k_bin + 0.5 * k_size) * x - - sel = (b_left >= b_bin - 0.5 * b_size) * (b_right <= b_bin + 0.5 * b_size) + \ - (b_left <= b_bin + 0.5 * b_size) * (b_right >= b_bin - 0.5 * b_size) - - return sel - -def hit_in_window(x, y, k_bin, b_bin, window_width=1.): - """ - Counts hits in a bin of track parameter space (b, k). - - Parameters - --------- - x : array-like - Array of x coordinates of hits. - y : array-like - Array of x coordinates of hits. - k_bin : float - Track parameter: y = k_bin * x + b_bin - b_bin : float - Track parameter: y = k_bin * x + b_bin - - Return - ------ - track_inds : array-like - Hit indexes of a track: [ind1, ind2, ...] - """ - - - y_approx = k_bin * x + b_bin - - flag = False - if np.abs(y_approx - y) <= window_width: - flag = True - - return flag - - -def get_zy_projection(z, xtop, ytop, xbot, ybot, k_y, b_y): - - x = k_y * z + b_y - - k = (ytop - ybot) / (xtop - xbot + 10**-6) - b = ytop - k * xtop - y = k * x + b - - return y - - -def sort_hits(hits): - - sorted_hits = [] - hits_z = [ahit['z'] for ahit in hits] - sort_index = np.argsort(hits_z) - for i_hit in sort_index: - sorted_hits.append(hits[i_hit]) - - return sorted_hits - - - - - - - - - From dd8a2b0011f1ba5878a46a4e34c025865db530ab Mon Sep 17 00:00:00 2001 From: ukmartylee <92804063+ukmartylee@users.noreply.github.com> Date: Fri, 23 Aug 2024 15:04:48 +0200 Subject: [PATCH 2/4] Remove all plane-related codes --- macro/flux_map.py | 2 +- muonShieldOptimization/ana_ShipMuon.py | 11 +- python/shipDet_conf.py | 6 - python/shipDigiReco.py | 12 +- python/shipPatRec.py | 11 +- python/shipStrawTracking.py | 15 +- strawtubes/strawtubes.cxx | 245 +++++++++---------------- strawtubes/strawtubes.h | 7 +- strawtubes/strawtubesHit.cxx | 10 +- strawtubes/strawtubes_single.cxx | 59 ++---- 10 files changed, 138 insertions(+), 240 deletions(-) diff --git a/macro/flux_map.py b/macro/flux_map.py index 1261d48681..8639248984 100755 --- a/macro/flux_map.py +++ b/macro/flux_map.py @@ -143,7 +143,7 @@ def main(): pid = hit.PdgCode() assert pid not in [12, -12, 14, -14, 16, -16] detector_ID = hit.GetDetectorID() - station = detector_ID // 10000000 + station = detector_ID // 1000000 h['T{}_all'.format(station)].Fill(x, y, weight) if abs(pid) == 13: muon = True diff --git a/muonShieldOptimization/ana_ShipMuon.py b/muonShieldOptimization/ana_ShipMuon.py index 67aea5bc4e..35304a9511 100644 --- a/muonShieldOptimization/ana_ShipMuon.py +++ b/muonShieldOptimization/ana_ShipMuon.py @@ -446,16 +446,15 @@ def makeProd(): os.chdir('../'+prefix+str(i+1)) def strawEncoding(detid): - # statnb*10000000+vnb*1000000+pnb*100000+lnb*10000+1000+snb - # vnb=view number; pnb=plane number; lnb=layer number; snb=straw number + # statnb*1000000+vnb*100000+lnb*10000+1000+snb + # vnb=view number; lnb=layer number; snb=straw number # statnb = station number. 1,2,3,4 tracking stations, 5 veto station vnb = ROOT.Long() - pnb = ROOT.Long() lnb = ROOT.Long() snb = ROOT.Long() statnb = ROOT.Long() - modules['Strawtubes'].StrawDecode(detid,statnb,vnb,pnb,lnb,snb) - return [statnb,vnb,pnb,lnb,snb] + modules['Strawtubes'].StrawDecode(detid,statnb,vnb,lnb,snb) + return [statnb,vnb,lnb,snb] def detMap(): sGeo = ROOT.gGeoManager @@ -658,7 +657,7 @@ def executeOneFile(fn,output=None,pid=None): detID = ahit.GetDetectorID() if ahit.GetName() == 'strawtubesPoint': tmp = strawEncoding(detID) - # detName = str(tmp[0]*10000000+tmp[1]*1000000+tmp[2]*100000+tmp[3]*10000) + # detName = str(tmp[0]*1000000+tmp[1]*100000+tmp[2]*10000) detName = "strawstation_"+str(tmp[0]) x = ahit.GetX() y = ahit.GetY() diff --git a/python/shipDet_conf.py b/python/shipDet_conf.py index cf8bdc62a8..94bdc33373 100644 --- a/python/shipDet_conf.py +++ b/python/shipDet_conf.py @@ -685,10 +685,6 @@ def configure(run, ship_geo): tauHpt.GetNumberofTargets(ship_geo.NuTauTarget.target) detectorList.append(tauHpt) - # for backward compatibility - if not hasattr(ship_geo.strawtubes, "YPlaneOffset"): - ship_geo.strawtubes.YLayerOffset = ship_geo.strawtubes.StrawPitch / 2.0 - ship_geo.strawtubes.YPlaneOffset = ship_geo.strawtubes.StrawPitch / 4.0 if ship_geo.strawDesign > 1: # for backward compatibility if ship_geo.strawDesign == 10 and not hasattr( @@ -719,10 +715,8 @@ def configure(run, ship_geo): Strawtubes.SetStrawPitch( ship_geo.strawtubes.StrawPitch, ship_geo.strawtubes.YLayerOffset, - ship_geo.strawtubes.YPlaneOffset, ) Strawtubes.SetDeltazLayer(ship_geo.strawtubes.DeltazLayer) - Strawtubes.SetDeltazPlane(ship_geo.strawtubes.DeltazPlane) Strawtubes.SetStrawsPerLayer(ship_geo.strawtubes.StrawsPerLayer) Strawtubes.SetStereoAngle(ship_geo.strawtubes.ViewAngle) Strawtubes.SetWireThickness(ship_geo.strawtubes.WireThickness) diff --git a/python/shipDigiReco.py b/python/shipDigiReco.py index e137cbe9c6..4880fbadd1 100644 --- a/python/shipDigiReco.py +++ b/python/shipDigiReco.py @@ -760,14 +760,14 @@ def withT0Estimate(self): key = -1 SmearedHits = [] v_drift = global_variables.modules["Strawtubes"].StrawVdrift() - global_variables.modules["Strawtubes"].StrawEndPoints(10002001, start, stop) + global_variables.modules["Strawtubes"].StrawEndPoints(1002001, start, stop) z1 = stop.z() for aDigi in self.digiStraw: key+=1 if not aDigi.isValid(): continue detID = aDigi.GetDetectorID() # don't use hits from straw veto - station = int(detID//10000000) + station = int(detID//1000000) if station > 4 : continue global_variables.modules["Strawtubes"].StrawEndPoints(detID, start, stop) delt1 = (start[2]-z1)/u.speedOfLight @@ -785,14 +785,14 @@ def smearHits(self,no_amb=None): SmearedHits = [] key = -1 v_drift = global_variables.modules["Strawtubes"].StrawVdrift() - global_variables.modules["Strawtubes"].StrawEndPoints(10002001, start, stop) + global_variables.modules["Strawtubes"].StrawEndPoints(1002001, start, stop) z1 = stop.z() for aDigi in self.digiStraw: key+=1 if not aDigi.isValid(): continue detID = aDigi.GetDetectorID() # don't use hits from straw veto - station = int(detID//10000000) + station = int(detID//1000000) if station > 4 : continue global_variables.modules["Strawtubes"].StrawEndPoints(detID, start, stop) #distance to wire @@ -845,7 +845,7 @@ def findTracks(self): atrack_smeared_hits = list(atrack_y12) + list(atrack_stereo12) + list(atrack_y34) + list(atrack_stereo34) for sm in atrack_smeared_hits: detID = sm['detID'] - station = int(detID//10000000) + station = int(detID//1000000) trID = i_track # Collect hits for track fit if trID not in hitPosLists: @@ -861,7 +861,7 @@ def findTracks(self): else: # do fake pattern recognition for sm in self.SmearedHits: detID = self.digiStraw[sm['digiHit']].GetDetectorID() - station = int(detID//10000000) + station = int(detID//1000000) trID = self.sTree.strawtubesPoint[sm['digiHit']].GetTrackID() if trID not in hitPosLists: hitPosLists[trID] = ROOT.std.vector('TVectorD')() diff --git a/python/shipPatRec.py b/python/shipPatRec.py index c28aff6bb8..018b92740a 100644 --- a/python/shipPatRec.py +++ b/python/shipPatRec.py @@ -799,7 +799,7 @@ def hits_split(smeared_hits): ahit = smeared_hits[i_hit] detID = ahit['detID'] - statnb, vnb, pnb, lnb, snb = decodeDetectorID(detID) + statnb, vnb, lnb, snb = decodeDetectorID(detID) is_y12 = ((statnb == 1) + (statnb == 2)) * ((vnb == 0) + (vnb == 3)) is_stereo12 = ((statnb == 1) + (statnb == 2)) * ((vnb == 1) + (vnb == 2)) is_y34 = ((statnb == 3) + (statnb == 4)) * ((vnb == 0) + (vnb == 3)) @@ -932,8 +932,6 @@ def decodeDetectorID(detID): Station numbers. vnb : int or array-like View numbers. - pnb : int or array-like - Plane numbers. lnb : int or array-like Layer numbers. snb : int or array-like @@ -942,11 +940,10 @@ def decodeDetectorID(detID): statnb = detID // 10000000 vnb = (detID - statnb * 10000000) // 1000000 - pnb = (detID - statnb * 10000000 - vnb * 1000000) // 100000 - lnb = (detID - statnb * 10000000 - vnb * 1000000 - pnb * 100000) // 10000 - snb = detID - statnb * 10000000 - vnb * 1000000 - pnb * 100000 - lnb * 10000 - 2000 + lnb = (detID - statnb * 10000000 - vnb * 1000000) // 10000 + snb = detID - statnb * 10000000 - vnb * 1000000 - lnb * 10000 - 2000 - return statnb, vnb, pnb, lnb, snb + return statnb, vnb, lnb, snb diff --git a/python/shipStrawTracking.py b/python/shipStrawTracking.py index a2b53869ca..3e929446cf 100644 --- a/python/shipStrawTracking.py +++ b/python/shipStrawTracking.py @@ -213,7 +213,7 @@ def run_track_pattern_recognition(input_file, geo_file, output_file, method): hits[key] = numpy.array(hits[key]) # Decoding - statnb, vnb, pnb, lnb, snb = decodeDetectorID(hits['DetID']) + statnb, vnb, lnb, snb = decodeDetectorID(hits['DetID']) is_stereo = ((vnb == 1) + (vnb == 2)) is_y = ((vnb == 0) + (vnb == 3)) is_before = ((statnb == 1) + (statnb == 2)) @@ -553,8 +553,6 @@ def decodeDetectorID(detID): Station numbers. vnb : int or array-like View numbers. - pnb : int or array-like - Plane numbers. lnb : int or array-like Layer numbers. snb : int or array-like @@ -563,11 +561,10 @@ def decodeDetectorID(detID): statnb = detID // 10000000 vnb = (detID - statnb * 10000000) // 1000000 - pnb = (detID - statnb * 10000000 - vnb * 1000000) // 100000 - lnb = (detID - statnb * 10000000 - vnb * 1000000 - pnb * 100000) // 10000 - snb = detID - statnb * 10000000 - vnb * 1000000 - pnb * 100000 - lnb * 10000 - 2000 + lnb = (detID - statnb * 10000000 - vnb * 1000000) // 10000 + snb = detID - statnb * 10000000 - vnb * 1000000 - lnb * 10000 - 2000 - return statnb, vnb, pnb, lnb, snb + return statnb, vnb, lnb, snb ######################################################################################################################## @@ -632,10 +629,10 @@ def getReconstructibleTracks(iEvent, sTree, sGeo, ShipGeo): VetoStationEndZ = VetoStationZ + (ShipGeo.strawtubes.DeltazView + ShipGeo.strawtubes.OuterStrawDiameter) / 2 TStationz = ShipGeo.TrackStation1.z - Zpos = TStationz - 3. /2. * ShipGeo.strawtubes.DeltazView - 1. / 2. * ShipGeo.strawtubes.DeltazPlane - 1. / 2. * ShipGeo.strawtubes.DeltazLayer + Zpos = TStationz - 3. /2. * ShipGeo.strawtubes.DeltazView - 1. / 2. * ShipGeo.strawtubes.DeltazLayer TStation1StartZ = Zpos - ShipGeo.strawtubes.OuterStrawDiameter / 2 - Zpos = TStationz + 3. /2. * ShipGeo.strawtubes.DeltazView + 1. / 2. * ShipGeo.strawtubes.DeltazPlane + 1. / 2. * ShipGeo.strawtubes.DeltazLayer + Zpos = TStationz + 3. /2. * ShipGeo.strawtubes.DeltazView + 1. / 2. * ShipGeo.strawtubes.DeltazLayer TStation4EndZ = Zpos + ShipGeo.strawtubes.OuterStrawDiameter / 2 diff --git a/strawtubes/strawtubes.cxx b/strawtubes/strawtubes.cxx index d4530799c4..52618b52d7 100644 --- a/strawtubes/strawtubes.cxx +++ b/strawtubes/strawtubes.cxx @@ -220,11 +220,10 @@ void strawtubes::SetOuterStrawDiameter(Double_t outerstrawdiameter) } -void strawtubes::SetStrawPitch(Double_t strawpitch,Double_t layer_offset, Double_t plane_offset) +void strawtubes::SetStrawPitch(Double_t strawpitch, Double_t layer_offset) { fStraw_pitch = strawpitch; //! Distance (x) between straws in one layer fOffset_layer12 = layer_offset; - fOffset_plane12 = plane_offset; } void strawtubes::SetDeltazLayer(Double_t deltazlayer) @@ -232,11 +231,6 @@ void strawtubes::SetDeltazLayer(Double_t deltazlayer) fDeltaz_layer12 = deltazlayer; //! Distance (z) between layer 1&2 } -void strawtubes::SetDeltazPlane(Double_t deltazplane) -{ - fDeltaz_plane12 = deltazplane; //! Distance (z) between plane 1&2 -} - void strawtubes::SetStrawsPerLayer(Int_t strawsperlayer) { fStraws_per_layer = strawsperlayer; //! number of straws in one layer @@ -244,7 +238,7 @@ void strawtubes::SetStrawsPerLayer(Int_t strawsperlayer) void strawtubes::SetStereoAngle(Double_t stereoangle) { - fView_angle = stereoangle; //! Stereo angle of planes in a view + fView_angle = stereoangle; //! Stereo angle of layers in a view fcosphi=cos(TMath::Pi()*fView_angle/180.); fsinphi=sin(TMath::Pi()*fView_angle/180.); } @@ -351,8 +345,6 @@ void strawtubes::ConstructGeometry() Double_t framewidth = 40.; //width of view Double_t viewwidth = fDeltaz_view-eps; - //width of plane - Double_t planewidth = fOuter_Straw_diameter+fDeltaz_layer12-eps; //width of layer Double_t layerwidth = fOuter_Straw_diameter; @@ -443,7 +435,7 @@ void strawtubes::ConstructGeometry() //vetovac->SetTransparency(80); //Veto station - //vnb=view number; pnb=plane number; lnb=layer number; snb=straw number + //vnb=view number; lnb=layer number; snb=straw number TString nmveto = "Veto"; TStationz=fT0z; for (Int_t vnb=0; vnb<2; vnb++) { @@ -478,52 +470,36 @@ void strawtubes::ConstructGeometry() TGeoTranslation t5p; - for (Int_t pnb=0; pnb<1; pnb++) { - //plane loop - TString nmplane_veto = nmveto+"_plane_"; nmplane_veto += pnb; - //width of the planes: z distance between layers + outer straw diameter - TGeoBBox *plane_veto = new TGeoBBox("plane box", fStraw_length_veto+eps/2., fvetoydim+eps/2., planewidth/2.+3.*eps/2.); - TGeoVolume *planebox_veto = new TGeoVolume(nmplane_veto, plane_veto, med); - //the planebox sits in the viewframe - //hence z translate the plane wrt to the view - - t5.SetTranslation(0, 0,(vnb-1./2.)*fDeltaz_view+(pnb-1./2.)*fDeltaz_plane12); - TGeoCombiTrans d5(t5, r5); - TGeoHMatrix *j5 = new TGeoHMatrix(d5); - vetovac->AddNode(planebox_veto, statnb*10000000+vnb*1000000+pnb*100000,j5); - - for (Int_t lnb=0; lnb<2; lnb++) { - TString nmlayer_veto = nmplane_veto+"_layer_"; nmlayer_veto += lnb; - //width of the layer: (plane width-2eps)/2 - layer_veto = new TGeoBBox("layer box_veto", fStraw_length_veto+eps/4., fvetoydim+eps/4., layerwidth/2.+eps/4.); - TGeoVolume *layerbox_veto = new TGeoVolume(nmlayer_veto, layer_veto, med); - //z translate the layerbox wrt the plane box (which is already rotated) - planebox_veto->AddNode(layerbox_veto, statnb*10000000+vnb*1000000+pnb*100000+lnb*10000,new TGeoTranslation(0,0,(lnb-1./2.)*fDeltaz_layer12)); - //layer loop - TGeoRotation r6v; - TGeoTranslation t6v; - Int_t nr = statnb*10000000+vnb*1000000+pnb*100000+lnb*10000; - for (Int_t snb=1; snbAddNode(straw_veto,nr+1000+snb,h6v); - layerbox_veto->AddNode(gas_veto, nr+2000+snb,h6v); - layerbox_veto->AddNode(wire_veto, nr+3000+snb,h6v); - //end of straw loop - } - //end of layer loop + for (Int_t lnb=0; lnb<2; lnb++) { + TString nmlayer_veto = nmveto+"_layer_"; nmlayer_veto += lnb; + layer_veto = new TGeoBBox("layer box_veto", fStraw_length_veto+eps/4., fvetoydim+eps/4., layerwidth/2.+eps/4.); + TGeoVolume *layerbox_veto = new TGeoVolume(nmlayer_veto, layer_veto, med); + //the layerbox sits in the viewframe + //hence z translate the layer wrt the view + vetovac->AddNode(layerbox_veto, statnb*10000000+vnb*1000000+lnb*10000,new TGeoTranslation(0,0,(lnb-1./2.)*fDeltaz_layer12)); + //layer loop + TGeoRotation r6v; + TGeoTranslation t6v; + Int_t nr = statnb*10000000+vnb*1000000+lnb*10000; + for (Int_t snb=1; snbAddNode(straw_veto,nr+1000+snb,h6v); + layerbox_veto->AddNode(gas_veto, nr+2000+snb,h6v); + layerbox_veto->AddNode(wire_veto, nr+3000+snb,h6v); + //end of straw loop } - //end of plane loop + //end of layer loop } //end of view loop } } // end of veto station loop //Tracking stations - //statnb=station number; vnb=view number; pnb=plane number; lnb=layer number; snb=straw number + //statnb=station number; vnb=view number; lnb=layer number; snb=straw number //New scalable endpoints of vacuum boxes which cover rotated view frames @@ -620,49 +596,31 @@ void strawtubes::ConstructGeometry() vac_12->AddNode(viewframe_12, statnb*10000000+vnb*1000000,h5); viewframe_12->SetLineColor(kRed); - for (Int_t pnb=0; pnb<1; pnb++) { - //plane loop - TString nmplane_12 = nmview_12+"_plane_"; - nmplane_12 += pnb; - TGeoBBox *plane_12 = new TGeoBBox("plane box_12", fStraw_length_12+eps/2, ftr12ydim+eps/2, planewidth/2.+3.*eps/2); - TGeoVolume *planebox_12 = new TGeoVolume(nmplane_12, plane_12, med); - - //the planebox sits in the viewframe - //hence z translate the plane wrt to the view - TGeoTranslation t3; - t3.SetTranslation(0, 0,(vnb-3./2.)*(fDeltaz_view)+(pnb-1./2.)*fDeltaz_plane12); - TGeoCombiTrans d3(t3, r5); - TGeoHMatrix *j3 = new TGeoHMatrix(d3); - vac_12->AddNode(planebox_12, statnb*10000000+vnb*1000000+pnb*100000,j3); - - for (Int_t lnb=0; lnb<2; lnb++) { - - //width of the layer: (plane width-2*eps)/2 - - //z translate the layerbox wrt the plane box (which is already rotated) - TString nmlayer_12 = nmplane_12+"_layer_"; nmlayer_12 += lnb; - TGeoBBox *layer_12 = new TGeoBBox("layer box_12", fStraw_length_12+eps/4, ftr12ydim+eps/4, layerwidth/2.+eps/4); - TGeoVolume *layerbox_12 = new TGeoVolume(nmlayer_12, layer_12, med); - planebox_12->AddNode(layerbox_12, statnb*10000000+vnb*1000000+pnb*100000+lnb*10000,new TGeoTranslation(0,0,(lnb-1./2.)*fDeltaz_layer12)); - - //layer loop - TGeoRotation r6s; - TGeoTranslation t6s; - for (Int_t snb=1; snbAddNode(straw_12,statnb*10000000+vnb*1000000+pnb*100000+lnb*10000+1000+snb,h6s); - layerbox_12->AddNode(gas_12,statnb*10000000+vnb*1000000+pnb*100000+lnb*10000+2000+snb,h6s); - layerbox_12->AddNode(wire_12,statnb*10000000+vnb*1000000+pnb*100000+lnb*10000+3000+snb,h6s); - - //end of straw loop - } - //end of layer loop + for (Int_t lnb=0; lnb<2; lnb++) { + + //the layerbox sits in the viewframe + //hence z translate the layer wrt the view + TString nmlayer_12 = nmview_12+"_layer_"; nmlayer_12 += lnb; + TGeoBBox *layer_12 = new TGeoBBox("layer box_12", fStraw_length_12+eps/4, ftr12ydim+eps/4, layerwidth/2.+eps/4); + TGeoVolume *layerbox_12 = new TGeoVolume(nmlayer_12, layer_12, med); + vac_12->AddNode(layerbox_12, statnb*10000000+vnb*1000000+lnb*10000,new TGeoTranslation(0,0,(lnb-1./2.)*fDeltaz_layer12)); + + //layer loop + TGeoRotation r6s; + TGeoTranslation t6s; + for (Int_t snb=1; snbAddNode(straw_12,statnb*10000000+vnb*1000000+lnb*10000+1000+snb,h6s); + layerbox_12->AddNode(gas_12,statnb*10000000+vnb*1000000+lnb*10000+2000+snb,h6s); + layerbox_12->AddNode(wire_12,statnb*10000000+vnb*1000000+lnb*10000+3000+snb,h6s); + + //end of straw loop } - //end of plane loop + //end of layer loop } //end of view loop } @@ -717,49 +675,31 @@ void strawtubes::ConstructGeometry() vac->AddNode(viewframe, statnb*10000000+vnb*1000000,h5); viewframe->SetLineColor(kRed); - for (Int_t pnb=0; pnb<1; pnb++) { - //plane loop - TString nmplane = nmview+"_plane_"; - nmplane += pnb; - TGeoBBox *plane = new TGeoBBox("plane box", fStraw_length+eps/2, ftr34ydim+eps/2, planewidth/2.+3.*eps/2); - TGeoVolume *planebox = new TGeoVolume(nmplane, plane, med); - - //the planebox sits in the viewframe - //hence z translate the plane wrt to the view - TGeoTranslation t3; - t3.SetTranslation(0, 0,(vnb-3./2.)*(fDeltaz_view)+(pnb-1./2.)*fDeltaz_plane12); - TGeoCombiTrans d3(t3, r5); - TGeoHMatrix *j3 = new TGeoHMatrix(d3); - vac->AddNode(planebox, statnb*10000000+vnb*1000000+pnb*100000,j3); - - for (Int_t lnb=0; lnb<2; lnb++) { - - //width of the layer: (plane width-2*eps)/2 - - //z translate the layerbox wrt the plane box (which is already rotated) - TString nmlayer = nmplane+"_layer_"; nmlayer += lnb; - TGeoBBox *layer = new TGeoBBox("layer box", fStraw_length+eps/4, ftr34ydim+eps/4, layerwidth/2.+eps/4); - TGeoVolume *layerbox = new TGeoVolume(nmlayer, layer, med); - planebox->AddNode(layerbox, statnb*10000000+vnb*1000000+pnb*100000+lnb*10000,new TGeoTranslation(0,0,(lnb-1./2.)*fDeltaz_layer12)); - - //layer loop - TGeoRotation r6s; - TGeoTranslation t6s; - for (Int_t snb=1; snbAddNode(straw,statnb*10000000+vnb*1000000+pnb*100000+lnb*10000+1000+snb,h6s); - layerbox->AddNode(gas,statnb*10000000+vnb*1000000+pnb*100000+lnb*10000+2000+snb,h6s); - layerbox->AddNode(wire,statnb*10000000+vnb*1000000+pnb*100000+lnb*10000+3000+snb,h6s); - - //end of straw loop - } - //end of layer loop + for (Int_t lnb=0; lnb<2; lnb++) { + + //the layerbox sits in the viewframe + //hence z translate the layer wrt the view + TString nmlayer = nmview+"_layer_"; nmlayer += lnb; + TGeoBBox *layer = new TGeoBBox("layer box", fStraw_length+eps/4, ftr34ydim+eps/4, layerwidth/2.+eps/4); + TGeoVolume *layerbox = new TGeoVolume(nmlayer, layer, med); + vac->AddNode(layerbox, statnb*10000000+vnb*1000000+lnb*10000,new TGeoTranslation(0,0,(lnb-1./2.)*fDeltaz_layer12)); + + //layer loop + TGeoRotation r6s; + TGeoTranslation t6s; + for (Int_t snb=1; snbAddNode(straw,statnb*10000000+vnb*1000000+lnb*10000+1000+snb,h6s); + layerbox->AddNode(gas,statnb*10000000+vnb*1000000+lnb*10000+2000+snb,h6s); + layerbox->AddNode(wire,statnb*10000000+vnb*1000000+lnb*10000+3000+snb,h6s); + + //end of straw loop } - //end of plane loop + //end of layer loop } //end of view loop } @@ -773,13 +713,12 @@ void strawtubes::ConstructGeometry() } // ----- Public method StrawDecode ------------------------------------------- // ----- returns station layer ... numbers ----------------------------------- -void strawtubes::StrawDecode(Int_t detID,int &statnb,int &vnb,int &pnb,int &lnb, int &snb) +void strawtubes::StrawDecode(Int_t detID,int &statnb,int &vnb,int &lnb, int &snb) { statnb = detID/10000000; vnb = (detID - statnb*10000000)/1000000; - pnb = (detID - statnb*10000000 - vnb*1000000)/100000; - lnb = (detID - statnb*10000000 - vnb*1000000 - pnb*100000)/10000; - snb = detID - statnb*10000000 - vnb*1000000 - pnb*100000 - lnb*10000 - 2000; + lnb = (detID - statnb*10000000 - vnb*1000000)/10000; + snb = detID - statnb*10000000 - vnb*1000000 - lnb*10000 - 2000; } // ----- Public method StrawEndPoints ------------------------------------------- // ----- returns top(left) and bottom(right) coordinate of straw ----------------------------------- @@ -788,8 +727,7 @@ void strawtubes::StrawEndPoints(Int_t fDetectorID, TVector3 &vbot, TVector3 &vto { Int_t statnb = fDetectorID/10000000; Int_t vnb = (fDetectorID - statnb*10000000)/1000000; - Int_t pnb = (fDetectorID- statnb*10000000 - vnb*1000000)/100000; - Int_t lnb = (fDetectorID - statnb*10000000 - vnb*1000000 - pnb*100000)/10000; + Int_t lnb = (fDetectorID - statnb*10000000 - vnb*1000000)/10000; TString stat = "Tr";stat+=+statnb;stat+="_";stat+=statnb; if (statnb==5){stat="Veto_5";} TString view; @@ -813,14 +751,13 @@ void strawtubes::StrawEndPoints(Int_t fDetectorID, TVector3 &vbot, TVector3 &vto TString prefix = "Tr"; if (statnb==5){prefix="Veto";} else{prefix+=statnb;} - prefix+=view;prefix+="_plane_";prefix+=pnb;prefix+="_"; - TString plane = prefix;plane+=statnb;plane+=vnb;plane+=+pnb;plane+="00000"; - TString layer = prefix+"layer_";layer+=lnb;layer+="_";layer+=statnb;layer+=vnb;layer+=pnb;layer+=lnb;layer+="0000"; + prefix+=view;prefix+="_"; + TString layer = prefix+"layer_";layer+=lnb;layer+="_";layer+=statnb;layer+=vnb;layer+=lnb;layer+="0000"; TString wire = "wire_"; if (statnb==5){wire+="veto_";} wire+=(fDetectorID+1000); if (statnb<3){wire = "wire_12_";wire+=(fDetectorID+1000);} - TString path = "/";path+=stat;path+="/";path+=plane;path+="/";path+=layer;path+="/";path+=wire; + TString path = "/";path+=stat;path+="/";path+=layer;path+="/";path+=wire; Bool_t rc = nav->cd(path); if (not rc){ cout << "strawtubes::StrawDecode, TgeoNavigator failed "<cd(path); if (not rc){ cout << "strawtubes::StrawDecode, TgeoNavigator failed "<SetVisibility(kFALSE); //planebox->VisibleDaughters(kTRUE); - /* TGeoRotation r5l; + TGeoRotation r5l; TGeoTranslation t5l; for (Int_t lnb=0; lnb<2; lnb++) { TString nmlayer = nmplane+"_layer_"; nmlayer += lnb; @@ -430,15 +424,15 @@ void strawtubes::ConstructGeometry() //end of straw loop } //end of layer loop - }*/ + } //end of plane loop - } + } */ //end of view loop } // end of veto station loop //Tracking stations - //statnb=station number; vnb=view number; pnb=plane number; lnb=layer number; snb=straw number + //statnb=station number; vnb=view number; lnb=layer number; snb=straw number TGeoBBox *vacbox = new TGeoBBox("vacbox", fVacBox_x, fVacBox_y, 2.*fDeltaz_view+eps ); for (Int_t statnb=0;statnb<4;statnb++) { @@ -534,54 +528,39 @@ void strawtubes::ConstructGeometry() wire->SetLineColor(6); //wire->SetVisibility(kTRUE); - for (Int_t pnb=0; pnb<2; pnb++) { - //plane loop - TString nmplane = nmview+"_plane_"; nmplane += pnb; - //width of the planes: z distance between layers + outer straw diameter - TGeoBBox *plane = new TGeoBBox("plane box", fStraw_length-eps, fStraw_length-eps, planewidth/2.); - TGeoVolume *planebox = new TGeoVolume(nmplane, plane, med); - //the planebox sits in the viewframe - //hence z translate the plane wrt to the view - viewframe->AddNode(planebox, statnb*10000000+vnb*1000000+pnb*100000,new TGeoTranslation(0, 0,(pnb-1./2.)*fDeltaz_plane12)); - //planebox->SetVisibility(kFALSE); - //planebox->VisibleDaughters(kTRUE); - TGeoRotation r5l; - TGeoTranslation t5l; for (Int_t lnb=0; lnb<2; lnb++) { TString nmlayer = nmplane+"_layer_"; nmlayer += lnb; - //width of the layer: (plane width-2eps)/2 TGeoBBox *layer = new TGeoBBox("layer box", fStraw_length-2*eps, fStraw_length-2*eps, layerwidth/2.); TGeoVolume *layerbox = new TGeoVolume(nmlayer, layer, med); - //z translate the layerbox wrt the plane box (which is already rotated) - planebox->AddNode(layerbox, statnb*10000000+vnb*1000000+pnb*100000+lnb*10000,new TGeoTranslation(0,0,(lnb-1./2.)*fDeltaz_layer12)); + //the layerbox sits in the viewframe + //hence z translate the layer wrt the view + vac->AddNode(layerbox, statnb*10000000+vnb*1000000+lnb*10000,new TGeoTranslation(0,0,(lnb-1./2.)*fDeltaz_layer12)); //layerbox->SetVisibility(kTRUE); //layerbox->VisibleDaughters(kTRUE); - //std::cout <AddNode(straw,statnb*10000000+vnb*1000000+pnb*100000+lnb*10000+1000+snb,h6s); - layerbox->AddNode(gas,statnb*10000000+vnb*1000000+pnb*100000+lnb*10000+2000+snb,h6s); - layerbox->AddNode(wire,statnb*10000000+vnb*1000000+pnb*100000+lnb*10000+3000+snb,h6s); + layerbox->AddNode(straw,statnb*10000000+vnb*1000000+lnb*10000+1000+snb,h6s); + layerbox->AddNode(gas,statnb*10000000+vnb*1000000+lnb*10000+2000+snb,h6s); + layerbox->AddNode(wire,statnb*10000000+vnb*1000000+lnb*10000+3000+snb,h6s); //only the straws are sensitive //end of straw loop } //end of layer loop } - //end of plane loop - } //end of view loop } //end of station From dc64264eac0da14fd3e133a3862fa4297f1c72d5 Mon Sep 17 00:00:00 2001 From: ukmartylee <92804063+ukmartylee@users.noreply.github.com> Date: Mon, 26 Aug 2024 14:22:38 +0200 Subject: [PATCH 3/4] Add update information --- strawtubes/strawtubes.cxx | 2 ++ 1 file changed, 2 insertions(+) diff --git a/strawtubes/strawtubes.cxx b/strawtubes/strawtubes.cxx index 52618b52d7..4e091ffc47 100644 --- a/strawtubes/strawtubes.cxx +++ b/strawtubes/strawtubes.cxx @@ -2,6 +2,8 @@ // 7/10/2015 // E. van Herwijnen eric.van.herwijnen@cern.ch // Also contains (for the moment) the veto station +// +// 23/08/2024 updated by W.-C. Marty Lee #include "strawtubes.h" #include "strawtubesPoint.h" From 766bebead5696b01b37804f99d103e5a8a1f7ea8 Mon Sep 17 00:00:00 2001 From: wcmartylee Date: Mon, 14 Oct 2024 15:39:46 +0200 Subject: [PATCH 4/4] Move station & view number by one digit in detID (remove plane number digit completely) --- python/shipPatRec.py | 8 +++---- python/shipStrawTracking.py | 8 +++---- strawtubes/strawtubes.cxx | 43 ++++++++++++++++++------------------ strawtubes/strawtubesHit.cxx | 6 ++--- 4 files changed, 33 insertions(+), 32 deletions(-) diff --git a/python/shipPatRec.py b/python/shipPatRec.py index 018b92740a..8d688ce500 100644 --- a/python/shipPatRec.py +++ b/python/shipPatRec.py @@ -938,10 +938,10 @@ def decodeDetectorID(detID): Straw tube numbers. """ - statnb = detID // 10000000 - vnb = (detID - statnb * 10000000) // 1000000 - lnb = (detID - statnb * 10000000 - vnb * 1000000) // 10000 - snb = detID - statnb * 10000000 - vnb * 1000000 - lnb * 10000 - 2000 + statnb = detID // 1000000 + vnb = (detID - statnb * 1000000) // 100000 + lnb = (detID - statnb * 1000000 - vnb * 100000) // 10000 + snb = detID - statnb * 1000000 - vnb * 100000 - lnb * 10000 - 2000 return statnb, vnb, lnb, snb diff --git a/python/shipStrawTracking.py b/python/shipStrawTracking.py index 3e929446cf..60848e81db 100644 --- a/python/shipStrawTracking.py +++ b/python/shipStrawTracking.py @@ -559,10 +559,10 @@ def decodeDetectorID(detID): Straw tube numbers. """ - statnb = detID // 10000000 - vnb = (detID - statnb * 10000000) // 1000000 - lnb = (detID - statnb * 10000000 - vnb * 1000000) // 10000 - snb = detID - statnb * 10000000 - vnb * 1000000 - lnb * 10000 - 2000 + statnb = detID // 1000000 + vnb = (detID - statnb * 1000000) // 100000 + lnb = (detID - statnb * 1000000 - vnb * 100000) // 10000 + snb = detID - statnb * 1000000 - vnb * 100000 - lnb * 10000 - 2000 return statnb, vnb, lnb, snb diff --git a/strawtubes/strawtubes.cxx b/strawtubes/strawtubes.cxx index 4e091ffc47..be5bb04b4a 100644 --- a/strawtubes/strawtubes.cxx +++ b/strawtubes/strawtubes.cxx @@ -3,7 +3,8 @@ // E. van Herwijnen eric.van.herwijnen@cern.ch // Also contains (for the moment) the veto station // -// 23/08/2024 updated by W.-C. Marty Lee +// The compact straw geometry +// updated by W.-C. Marty Lee 23.08.2024 #include "strawtubes.h" #include "strawtubesPoint.h" @@ -467,7 +468,7 @@ void strawtubes::ConstructGeometry() r5.SetAngles(angle,0,0); TGeoCombiTrans c5(t5, r5); TGeoHMatrix *h5 = new TGeoHMatrix(c5); - vetovac->AddNode(viewframe_veto, statnb*10000000+vnb*1000000,h5); + vetovac->AddNode(viewframe_veto, statnb*1000000+vnb*100000,h5); viewframe_veto->SetLineColor(kRed); TGeoTranslation t5p; @@ -478,11 +479,11 @@ void strawtubes::ConstructGeometry() TGeoVolume *layerbox_veto = new TGeoVolume(nmlayer_veto, layer_veto, med); //the layerbox sits in the viewframe //hence z translate the layer wrt the view - vetovac->AddNode(layerbox_veto, statnb*10000000+vnb*1000000+lnb*10000,new TGeoTranslation(0,0,(lnb-1./2.)*fDeltaz_layer12)); + vetovac->AddNode(layerbox_veto, statnb*1000000+vnb*100000+lnb*10000,new TGeoTranslation(0,0,(lnb-1./2.)*fDeltaz_layer12)); //layer loop TGeoRotation r6v; TGeoTranslation t6v; - Int_t nr = statnb*10000000+vnb*1000000+lnb*10000; + Int_t nr = statnb*1000000+vnb*100000+lnb*10000; for (Int_t snb=1; snbAddNode(viewframe_12, statnb*10000000+vnb*1000000,h5); + vac_12->AddNode(viewframe_12, statnb*1000000+vnb*100000,h5); viewframe_12->SetLineColor(kRed); for (Int_t lnb=0; lnb<2; lnb++) { @@ -605,7 +606,7 @@ void strawtubes::ConstructGeometry() TString nmlayer_12 = nmview_12+"_layer_"; nmlayer_12 += lnb; TGeoBBox *layer_12 = new TGeoBBox("layer box_12", fStraw_length_12+eps/4, ftr12ydim+eps/4, layerwidth/2.+eps/4); TGeoVolume *layerbox_12 = new TGeoVolume(nmlayer_12, layer_12, med); - vac_12->AddNode(layerbox_12, statnb*10000000+vnb*1000000+lnb*10000,new TGeoTranslation(0,0,(lnb-1./2.)*fDeltaz_layer12)); + vac_12->AddNode(layerbox_12, statnb*1000000+vnb*100000+lnb*10000,new TGeoTranslation(0,0,(lnb-1./2.)*fDeltaz_layer12)); //layer loop TGeoRotation r6s; @@ -616,9 +617,9 @@ void strawtubes::ConstructGeometry() r6s.SetAngles(90,90,0); TGeoCombiTrans c6s(t6s, r6s); TGeoHMatrix *h6s = new TGeoHMatrix(c6s); - layerbox_12->AddNode(straw_12,statnb*10000000+vnb*1000000+lnb*10000+1000+snb,h6s); - layerbox_12->AddNode(gas_12,statnb*10000000+vnb*1000000+lnb*10000+2000+snb,h6s); - layerbox_12->AddNode(wire_12,statnb*10000000+vnb*1000000+lnb*10000+3000+snb,h6s); + layerbox_12->AddNode(straw_12,statnb*1000000+vnb*100000+lnb*10000+1000+snb,h6s); + layerbox_12->AddNode(gas_12,statnb*1000000+vnb*100000+lnb*10000+2000+snb,h6s); + layerbox_12->AddNode(wire_12,statnb*1000000+vnb*100000+lnb*10000+3000+snb,h6s); //end of straw loop } @@ -674,7 +675,7 @@ void strawtubes::ConstructGeometry() TGeoCombiTrans c5(t5, r5); TGeoHMatrix *h5 = new TGeoHMatrix(c5); - vac->AddNode(viewframe, statnb*10000000+vnb*1000000,h5); + vac->AddNode(viewframe, statnb*1000000+vnb*100000,h5); viewframe->SetLineColor(kRed); for (Int_t lnb=0; lnb<2; lnb++) { @@ -684,7 +685,7 @@ void strawtubes::ConstructGeometry() TString nmlayer = nmview+"_layer_"; nmlayer += lnb; TGeoBBox *layer = new TGeoBBox("layer box", fStraw_length+eps/4, ftr34ydim+eps/4, layerwidth/2.+eps/4); TGeoVolume *layerbox = new TGeoVolume(nmlayer, layer, med); - vac->AddNode(layerbox, statnb*10000000+vnb*1000000+lnb*10000,new TGeoTranslation(0,0,(lnb-1./2.)*fDeltaz_layer12)); + vac->AddNode(layerbox, statnb*1000000+vnb*100000+lnb*10000,new TGeoTranslation(0,0,(lnb-1./2.)*fDeltaz_layer12)); //layer loop TGeoRotation r6s; @@ -695,9 +696,9 @@ void strawtubes::ConstructGeometry() r6s.SetAngles(90,90,0); TGeoCombiTrans c6s(t6s, r6s); TGeoHMatrix *h6s = new TGeoHMatrix(c6s); - layerbox->AddNode(straw,statnb*10000000+vnb*1000000+lnb*10000+1000+snb,h6s); - layerbox->AddNode(gas,statnb*10000000+vnb*1000000+lnb*10000+2000+snb,h6s); - layerbox->AddNode(wire,statnb*10000000+vnb*1000000+lnb*10000+3000+snb,h6s); + layerbox->AddNode(straw,statnb*1000000+vnb*100000+lnb*10000+1000+snb,h6s); + layerbox->AddNode(gas,statnb*1000000+vnb*100000+lnb*10000+2000+snb,h6s); + layerbox->AddNode(wire,statnb*1000000+vnb*100000+lnb*10000+3000+snb,h6s); //end of straw loop } @@ -717,19 +718,19 @@ void strawtubes::ConstructGeometry() // ----- returns station layer ... numbers ----------------------------------- void strawtubes::StrawDecode(Int_t detID,int &statnb,int &vnb,int &lnb, int &snb) { - statnb = detID/10000000; - vnb = (detID - statnb*10000000)/1000000; - lnb = (detID - statnb*10000000 - vnb*1000000)/10000; - snb = detID - statnb*10000000 - vnb*1000000 - lnb*10000 - 2000; + statnb = detID/1000000; + vnb = (detID - statnb*1000000)/100000; + lnb = (detID - statnb*1000000 - vnb*100000)/10000; + snb = detID - statnb*1000000 - vnb*100000 - lnb*10000 - 2000; } // ----- Public method StrawEndPoints ------------------------------------------- // ----- returns top(left) and bottom(right) coordinate of straw ----------------------------------- void strawtubes::StrawEndPoints(Int_t fDetectorID, TVector3 &vbot, TVector3 &vtop) // method to get end points from TGeoNavigator { - Int_t statnb = fDetectorID/10000000; - Int_t vnb = (fDetectorID - statnb*10000000)/1000000; - Int_t lnb = (fDetectorID - statnb*10000000 - vnb*1000000)/10000; + Int_t statnb = fDetectorID/1000000; + Int_t vnb = (fDetectorID - statnb*1000000)/100000; + Int_t lnb = (fDetectorID - statnb*1000000 - vnb*100000)/10000; TString stat = "Tr";stat+=+statnb;stat+="_";stat+=statnb; if (statnb==5){stat="Veto_5";} TString view; diff --git a/strawtubes/strawtubesHit.cxx b/strawtubes/strawtubesHit.cxx index 41a88a3726..5c42c9ec9f 100644 --- a/strawtubes/strawtubesHit.cxx +++ b/strawtubes/strawtubesHit.cxx @@ -45,9 +45,9 @@ strawtubesHit::strawtubesHit(strawtubesPoint* p, Double_t t0) } void strawtubesHit::StrawEndPoints(TVector3 &vbot, TVector3 &vtop) { - Int_t statnb = fDetectorID/10000000; - Int_t vnb = (fDetectorID - statnb*10000000)/1000000; - Int_t lnb = (fDetectorID - statnb*10000000 - vnb*1000000)/10000; + Int_t statnb = fDetectorID/1000000; + Int_t vnb = (fDetectorID - statnb*1000000)/100000; + Int_t lnb = (fDetectorID - statnb*1000000 - vnb*100000)/10000; TString stat = "Tr";stat+=+statnb;stat+="_";stat+=statnb; if (statnb==5){stat="Veto_5";} TString view;