Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • azhemch/spdroot
  • Tkachenko/spdroot
  • AndreiMaltsev/spdroot
  • mzhabitsky/spdroot
  • gag/spdroot
  • Vlandr57cms/spdroot
  • polianskivv/spdroot
  • Kurbatov/spdroot
  • jemtchou/spdroot
  • ElinaK/spdroot
  • vvpopov/spdroot
  • azelenov/spdroot
  • shtejer/spdroot
  • ybedfer/spdroot
14 results
Show changes
Commits on Source (45)
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 1207 additions and 887 deletions
.gitlab-ci.yml
View file @ e7a31252
......@@ -121,7 +121,9 @@ deploy_spdroot:
- rm -rf build
- mkdir build
- cd build
- cat SetEnv.sh
- source ../SetEnv.sh
- gcc --version
- echo "Installing in /cvmfs/nica.jinr.ru/spd/$CI_ENVIRONMENT_NAME/spdroot/$CI_COMMIT_REF_NAME"
- cmake -DCMAKE_INSTALL_PREFIX=/cvmfs/nica.jinr.ru/spd/$CI_ENVIRONMENT_NAME/spdroot/$CI_COMMIT_REF_NAME ..
- make -j 5
......
CMakeLists.txt
View file @ e7a31252
......@@ -370,7 +370,7 @@ Install(DIRECTORY macro/run macro/primgen macro/analysis
DESTINATION share/examples
PATTERN ".svn" EXCLUDE)
Install(DIRECTORY macro/examples/K0decay macro/examples/chic macro/examples/ecal macro/examples/jpsi-ee macro/examples/jpsi-mumu
Install(DIRECTORY macro/examples/K0decay macro/examples/chic macro/examples/ecal macro/examples/jpsi-mumu
DESTINATION share/examples/fullchain/
PATTERN ".svn" EXCLUDE)
......
bbc/SpdBbc.cxx
View file @ e7a31252
......@@ -44,7 +44,9 @@ SpdBbc::SpdBbc():
SpdDetector("B-B counter (qsl)", kSpdBbc, kTRUE),
fTrackID(-1),
fVolumeID(-1),
fModuleID(0),
fLocalModuleID(0),
fRingID(0),
fSectorID(0),
fPos(),fMom(),
fPosOut(),fMomOut(),
fTime(-1.),fTimeOut(-1.),
......@@ -69,7 +71,9 @@ SpdBbc::SpdBbc(const char* name, Bool_t active):
SpdDetector(name, kSpdBbc, active),
fTrackID(-1),
fVolumeID(-1),
fModuleID(0),
fLocalModuleID(0),
fRingID(),
fSectorID(),
fPos(),fMom(),
fPosOut(),fMomOut(),
fTime(-1.),fTimeOut(-1.),
......@@ -165,13 +169,18 @@ Bool_t SpdBbc::ProcessHits(FairVolume* vol)
fSegmentLength = gMC->TrackLength() - fLength;
fTimeOut = gMC->TrackTime() * 1.e9; // s -> ns
fVolumeID = kSpdBbc;
fVolumePath = gMC->CurrentVolPath();
if (fVolumePath.Contains("zplus")) fModuleID = 1;
else if (fVolumePath.Contains("zminus")) fModuleID = -1;
else fModuleID = -10; //???
fNodePath = gMC->CurrentVolPath();
SpdGeopathParser parser;
parser.ParsePath(fNodePath);
//cave: 1, Bbcmodule: 2, BBcRing:3, Bbcsector:4
fVolumeID = kSpdBbc;
fLocalModuleID = parser.Num(2, true);
fRingID = parser.Num(3, true);
if (fRingID==1) fSectorID =0;
else fSectorID = parser.Num(4, true);
gMC->TrackPosition(fPosOut);
gMC->TrackMomentum(fMomOut);
......@@ -198,7 +207,7 @@ void SpdBbc::AddHit()
TClonesArray& clref = *fPointCollection;
Int_t size = clref.GetEntriesFast();
new(clref[size]) SpdBbcPoint(fTrackID,fVolumeID,fModuleID,
new(clref[size]) SpdBbcPoint(fTrackID,fVolumeID,fLocalModuleID,fRingID,fSectorID,
TVector3(fPos.X(),fPos.Y(),fPos.Z()),
TVector3(fMom.Px(),fMom.Py(),fMom.Pz()),
TVector3(fPosOut.X(),fPosOut.Y(),fPosOut.Z()),
......@@ -266,40 +275,154 @@ void SpdBbc::ConstructDetector()
cout << "-E- <SpdBbc::BuildModule> Unknown material: " << mat << endl;
return;
}
// Double_t mthick = mapper->GetModuleThickness();
// Double_t msize = mapper->GetModuleSize();
// Double_t mwidth = mapper->GetModuleWidth();
Double_t mthick = mapper->GetModuleThickness();
Double_t msize = mapper->GetModuleSize();
Double_t mwidth = mapper->GetModuleWidth();
fModule1 = gGeoManager->MakeTube(mapper->AddPrefix("zplus"),material,msize-mwidth,msize,0.5*mthick);
if (!fModule1) {
cout << "-E- <SpdBbc::ConstructGeometry> Module (1) cannot be created " << endl;
return;
}
fModule2 = gGeoManager->MakeTube(mapper->AddPrefix("zminus"),material,msize-mwidth,msize,0.5*mthick);
if (!fModule2) {
cout << "-E- <SpdBbc::ConstructGeometry> Module (2) cannot be created " << endl;
return;
}
Double_t bbcthickness = mapper->GetModuleThickness(); //cm
//**************
TGeoVolume *fRing1 = new TGeoVolumeAssembly(mapper->AddPrefix("Ring"));
TGeoVolume *fSector1 = gGeoManager->MakeTube(mapper->AddPrefix("Sector"),material,4.25,9.5,0.5*bbcthickness);
fSector1->SetFillColor(kMagenta);
fSector1->SetLineColor(kMagenta);
fSector1->SetTransparency(15);
TGeoVolume *fRing2 = new TGeoVolumeAssembly(mapper->AddPrefix("Ring"));
TGeoVolume *fSector2 = gGeoManager->MakeTrd1(mapper->AddPrefix("Sector"), material, 3.5/2, 8.9/2 ,0.5*bbcthickness, 13.7/2 );
fSector2->SetFillColor(kMagenta);
fSector2->SetLineColor(kMagenta);
fSector2->SetTransparency(15);
TGeoVolume *fRing3 = new TGeoVolumeAssembly(mapper->AddPrefix("Ring"));
TGeoVolume *fSector3 = gGeoManager->MakeTrd1(mapper->AddPrefix("Sector"), material, 9.3/2, 14.8/2 ,0.5*bbcthickness, 13.7/2 );
fSector3->SetFillColor(kMagenta);
fSector3->SetLineColor(kMagenta);
fSector3->SetTransparency(15);
TGeoVolume *fRing4 = new TGeoVolumeAssembly(mapper->AddPrefix("Ring"));
TGeoVolume *fSector4 = gGeoManager->MakeTrd1(mapper->AddPrefix("Sector"), material, 15.2/2, 20.6/2 ,0.5*bbcthickness, 13.7/2 );
fSector4->SetFillColor(kMagenta);
fSector4->SetLineColor(kMagenta);
fSector4->SetTransparency(15);
TGeoVolume *fRing5 = new TGeoVolumeAssembly(mapper->AddPrefix("Ring"));
TGeoVolume *fSector5 = gGeoManager->MakeTrd1(mapper->AddPrefix("Sector"), material, 21./2, 26.5/2 ,0.5*bbcthickness, 13.7/2);
fSector5->SetFillColor(kMagenta);
fSector5->SetLineColor(kMagenta);
fSector5->SetTransparency(15);
TGeoVolume *fRing6 = new TGeoVolumeAssembly(mapper->AddPrefix("Ring"));
TGeoVolume *fSector6 = gGeoManager->MakeTrd1(mapper->AddPrefix("Sector"), material, 26.9/2, 32.3/2 ,0.5*bbcthickness, 13.7/2 );
fSector6->SetFillColor(kMagenta);
fSector6->SetLineColor(kMagenta);
fSector6->SetTransparency(15);
AddSensitiveVolume(fSector1); //ATTENTION FIXME ATTENTION
AddSensitiveVolume(fSector2); //ATTENTION FIXME ATTENTION
AddSensitiveVolume(fSector3); //ATTENTION FIXME ATTENTION
AddSensitiveVolume(fSector4); //ATTENTION FIXME ATTENTION
AddSensitiveVolume(fSector5); //ATTENTION FIXME ATTENTION
AddSensitiveVolume(fSector6); //ATTENTION FIXME ATTENTION
for (int j=0; j<6; j++){
Int_t inode = 0;
if (j==0){
TGeoRotation rot0;
TGeoTranslation trans0;
trans0.SetTranslation(0,0,0);
fRing1->AddNode(fSector1, 1, new TGeoCombiTrans(trans0, rot0));
}
else{
for (int i = 0; i<16; i++){
TGeoRotation rot1;
TGeoTranslation trans1;
rot1.RotateX(-90);
Double_t angle = 22.5/2+i*22.5;
rot1.RotateZ(-angle);
Double_t h = 13.7;
Double_t rho = 10.0+13.7/2 + (j-1)*(h+1.);
// Double_t rho = 10.0+13.7/2 + 0*(h+1.);
Double_t xd = rho*sin(angle*TMath::Pi()/180.);
Double_t yd = rho*cos(angle*TMath::Pi()/180.);
trans1.SetTranslation(xd, yd, 0);
if(j==1) fRing2->AddNode(fSector2, ++inode, new TGeoCombiTrans(trans1, rot1));
if(j==2) fRing3->AddNode(fSector3, ++inode, new TGeoCombiTrans(trans1, rot1));
if(j==3) fRing4->AddNode(fSector4, ++inode, new TGeoCombiTrans(trans1, rot1));
if(j==4) fRing5->AddNode(fSector5, ++inode, new TGeoCombiTrans(trans1, rot1));
if(j==5) fRing6->AddNode(fSector6, ++inode, new TGeoCombiTrans(trans1, rot1));
}
}
fModule1->SetFillColor(kMagenta);
fModule1->SetLineColor(kMagenta);
fModule1->SetTransparency(30);
}
fModule2->SetFillColor(kMagenta);
fModule2->SetLineColor(kMagenta);
fModule2->SetTransparency(30);
TGeoRotation rot0;
TGeoTranslation trans0;
trans0.SetTranslation(0,0,0);
// Double_t shift = 171.6;
Double_t shift = mapper->GetMinDistance();
//TGeoVolume *
fModule1 = new TGeoVolumeAssembly(mapper->AddPrefix("zplus"));
AddSensitiveVolume(fModule1); //ATTENTION FIXME ATTENTION
AddSensitiveVolume(fModule2); //ATTENTION FIXME ATTENTION
fModule1->AddNode(fRing1, 1, new TGeoCombiTrans(trans0, rot0));
fModule1->AddNode(fRing2, 2, new TGeoCombiTrans(trans0, rot0));
fModule1->AddNode(fRing3, 3, new TGeoCombiTrans(trans0, rot0));
fModule1->AddNode(fRing4, 4, new TGeoCombiTrans(trans0, rot0));
fModule1->AddNode(fRing5, 5, new TGeoCombiTrans(trans0, rot0));
fModule1->AddNode(fRing6, 6, new TGeoCombiTrans(trans0, rot0));
TGeoRotation rot1;
// rot1.RotateX(0);
TGeoTranslation trans1;
trans1.SetTranslation(0,0,shift);
fMasterVolume->AddNode(fModule1, 1, new TGeoCombiTrans(trans1, rot1));
// TGeoVolume *
fModule2 = new TGeoVolumeAssembly(mapper->AddPrefix("zminus"));
fModule2->AddNode(fRing1, 1, new TGeoCombiTrans(trans0, rot0));
fModule2->AddNode(fRing2, 2, new TGeoCombiTrans(trans0, rot0));
fModule2->AddNode(fRing3, 3, new TGeoCombiTrans(trans0, rot0));
fModule2->AddNode(fRing4, 4, new TGeoCombiTrans(trans0, rot0));
fModule2->AddNode(fRing5, 5, new TGeoCombiTrans(trans0, rot0));
fModule2->AddNode(fRing6, 6, new TGeoCombiTrans(trans0, rot0));
Double_t shift = mapper->GetMinDistance() + 0.5*mthick;
TGeoRotation rot2;
// rot1.RotateX(0);
TGeoTranslation trans2;
trans2.SetTranslation(0,0,-shift);
fMasterVolume->AddNode(fModule2, 2, new TGeoCombiTrans(trans2, rot2));
fMasterVolume->AddNode(fModule1,1,new TGeoTranslation(0,0, shift));
fMasterVolume->AddNode(fModule2,1,new TGeoTranslation(0,0,-shift));
}
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
......
bbc/SpdBbc.h
View file @ e7a31252
......@@ -63,7 +63,10 @@ private:
Int_t fTrackID; //! track index
Int_t fVolumeID; //! volume id (Detector Id)
Int_t fModuleID; //! module id (Endcap z+, z-)
Int_t fLocalModuleID; //! module id (z+, z-)
Int_t fRingID;
Int_t fSectorID;
TString fNodePath;
TString fVolumePath; //! full path to current volume
TLorentzVector fPos; //! position at entrance
TLorentzVector fMom; //! momentum at entrance
......
bbc/SpdBbcPoint.cxx
View file @ e7a31252
......@@ -11,18 +11,20 @@ using std::endl;
ClassImp(SpdBbcPoint)
//_____________________________________________________________________________________
SpdBbcPoint::SpdBbcPoint():FairMCPoint(),fModuleId(0),fTimeOut(0),fSegmentLength(0)
SpdBbcPoint::SpdBbcPoint():FairMCPoint(),fLocalModId(0),fTimeOut(0),fSegmentLength(0)
{
fELoss = -1.;
}
//_____________________________________________________________________________________
SpdBbcPoint::SpdBbcPoint(Int_t trackID, Int_t detID, Int_t modID,
SpdBbcPoint::SpdBbcPoint(Int_t trackID, Int_t detID, Int_t modID, Int_t ringID, Int_t sectorID,
TVector3 pos, TVector3 mom, TVector3 posOut, TVector3 momOut,
Double_t tof, Double_t tofOut,
Double_t length, Double_t segmentlength, Double_t ELoss)
:FairMCPoint(trackID, detID, pos, mom, tof, length, ELoss),
fModuleId(modID),
fLocalModId(modID),
fRingID(ringID),
fSectorId(sectorID),
fPosOut(posOut),fMomOut(momOut),fTimeOut(tofOut),
fSegmentLength(segmentlength)
{
......@@ -41,7 +43,7 @@ void SpdBbcPoint::Print(const Option_t* opt) const
cout << "<SpdBbcPoint::Print> " << endl;
cout << "\n";
cout << " Track: " << fTrackID << endl;
cout << " DetectorID/ModuleID: " << fDetectorID << "/" << fModuleId << endl;
cout << " DetectorID/ModuleID: " << fDetectorID << "/" << fLocalModId << endl;
cout << " Position (in,out): " << "(" << fX << ", " << fY << ", " << fZ << "); "
<< "(" << fPosOut.X() << ", " << fPosOut.Y() << ", " << fPosOut.Z() << ") "
<< "[cm]" << endl;
......
bbc/SpdBbcPoint.h
View file @ e7a31252
......@@ -34,7 +34,7 @@ public:
SpdBbcPoint();
SpdBbcPoint(Int_t trackID, Int_t detID, Int_t modID,
SpdBbcPoint(Int_t trackID, Int_t detID, Int_t localmodID, Int_t ringID, Int_t sectorID,
TVector3 pos, TVector3 mom, TVector3 posOut, TVector3 momOut,
Double_t tof, Double_t tofOut,
Double_t length, Double_t segmentlength, Double_t ELoss);
......@@ -43,7 +43,9 @@ public:
virtual ~SpdBbcPoint();
inline Int_t GetModuleId() const { return fModuleId; }
inline Int_t GetLocalModId() const { return fLocalModId; }
inline Int_t GetRingId() const { return fRingID; }
inline Int_t GetSectorId() const { return fSectorId; }
inline TVector3 GetPosIn() const { return TVector3(fX,fY,fZ); } // cm
inline TVector3 GetMomIn() const { return TVector3(fPx,fPy,fPz); } // GeV/c
......@@ -57,7 +59,9 @@ public:
private:
Int_t fModuleId; // -1 = endcap z-; +1 = endcap z+
Int_t fLocalModId; // 1 = endcap z-; 2 = endcap z+
Int_t fRingID; // Ring: 1-6
Int_t fSectorId; // Sector in ring: 1-16
TVector3 fPosOut; // Track coordinates at exit to active volume, [cm]
TVector3 fMomOut; // Track coordinates at exit to active volume, [cm]
......
common/SpdCommonGeoMapper.cxx
View file @ e7a31252
......@@ -36,7 +36,7 @@ TGeoVolume* SpdCommonGeoMapper::theMasterVolume = 0;
/*==================================================================*/
Int_t SpdCommonGeoMapper::theNGeoSectors = 8; // > 2 and < 13
Double_t SpdCommonGeoMapper::theSectorClearance = 0.2; // cm
Double_t SpdCommonGeoMapper::theSectorClearance = 0.4; // cm
/*============================= PIPE ===============================*/
......@@ -67,18 +67,18 @@ Double_t SpdCommonGeoMapper::theItsMaxLength = 200.; // cm
/*===================== BEAM-BEAM COUNTER (BBC) ======================*/
Int_t SpdCommonGeoMapper::theBbcDefGeoType = 1;
TString SpdCommonGeoMapper::theBbcBaseMaterial = "air";
Double_t SpdCommonGeoMapper::theBbcThickness = 5.; //2.; // cm
TString SpdCommonGeoMapper::theBbcBaseMaterial = "silicon"; // "air"
Double_t SpdCommonGeoMapper::theBbcThickness = 0.5; //2.; // cm
Double_t SpdCommonGeoMapper::theBbcSize = 84.7; // cm
Double_t SpdCommonGeoMapper::theBbcWidth = 78.7; // cm
Double_t SpdCommonGeoMapper::theBbcMinDist = 171.6; //167.8; // cm
Double_t SpdCommonGeoMapper::theBbcMinDist = 186.6; // cm
/*===================== TIME-OF-FLIGHT SYSTEM (BARREL) ======================*/
Int_t SpdCommonGeoMapper::theTofBDefGeoType = 1;
TString SpdCommonGeoMapper::theTofBBaseMaterial = "air";
Double_t SpdCommonGeoMapper::theTofBLength = 371.2; // cm
Double_t SpdCommonGeoMapper::theTofBSize = 106.4; // cm
Double_t SpdCommonGeoMapper::theTofBSize = 112.0; // cm
Double_t SpdCommonGeoMapper::theTofBWidth = 18.7; // cm
/*===================== TIME-OF-FLIGHT SYSTEM (ENDCAPS) =====================*/
......@@ -86,18 +86,18 @@ Double_t SpdCommonGeoMapper::theTofBWidth = 18.7; // cm
Int_t SpdCommonGeoMapper::theTofECDefGeoType = 1;
TString SpdCommonGeoMapper::theTofECBaseMaterial = "air";
Double_t SpdCommonGeoMapper::theTofECThickness = 6.; //28.; // cm
Double_t SpdCommonGeoMapper::theTofECSize = 84.7; // cm
Double_t SpdCommonGeoMapper::theTofECWidth = 69.7; // cm
Double_t SpdCommonGeoMapper::theTofECMinDist = 179.6; //171.6; // cm
Double_t SpdCommonGeoMapper::theTofECSize = 89.2; // cm
Double_t SpdCommonGeoMapper::theTofECWidth = 74.2; // cm
Double_t SpdCommonGeoMapper::theTofECMinDist = 194.0; // cm
/*===================== AEROGEL (AEG) ======================*/
Int_t SpdCommonGeoMapper::theAegDefGeoType = 1;
TString SpdCommonGeoMapper::theAegBaseMaterial = "air";
Double_t SpdCommonGeoMapper::theAegThickness = 16.0; //2.; // cm
Double_t SpdCommonGeoMapper::theAegThickness = 30.0; // cm
Double_t SpdCommonGeoMapper::theAegSize = 84.7; // cm
Double_t SpdCommonGeoMapper::theAegWidth = 78.7; // cm
Double_t SpdCommonGeoMapper::theAegMinDist = 152.6; //199.6; // cm
Double_t SpdCommonGeoMapper::theAegMinDist = 153.6; //199.6; // cm
/*==================================================================*/
/*= ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ =*/
......@@ -286,13 +286,13 @@ TString SpdCommonGeoMapper::theSolMagnetTwoCryoMaterial = "steel";
Int_t SpdCommonGeoMapper::theSolMagnetTwoNSections = 10;
Double_t SpdCommonGeoMapper::theSolMagnetTwoTotalLen = 397.6; // cm // changed from 380
Double_t SpdCommonGeoMapper::theSolMagnetTwoTotalSize = 188.4; // cm
Double_t SpdCommonGeoMapper::theSolMagnetTwoTotalWidth = 23.; // cm
Double_t SpdCommonGeoMapper::theSolMagnetTwoTotalLen = 408.0; // cm //changed from 397.6
Double_t SpdCommonGeoMapper::theSolMagnetTwoTotalSize = 200.4; // cm //changed from 188.4
Double_t SpdCommonGeoMapper::theSolMagnetTwoTotalWidth = 27.; // cm //changed from 23
Double_t SpdCommonGeoMapper::theSolMagnetTwoLen = 39.76; // cm
Double_t SpdCommonGeoMapper::theSolMagnetTwoSize = 188.4; // cm
Double_t SpdCommonGeoMapper::theSolMagnetTwoWidth = 23.; // cm
Double_t SpdCommonGeoMapper::theSolMagnetTwoLen = 40.8; // cm //changed from 39.76
Double_t SpdCommonGeoMapper::theSolMagnetTwoSize = 200.4; // cm //changed from 188.4
Double_t SpdCommonGeoMapper::theSolMagnetTwoWidth = 27.; // cm //changed from 23
Double_t SpdCommonGeoMapper::theSolMagnetTwoInnerCryoW = 0.5; // cm
Double_t SpdCommonGeoMapper::theSolMagnetTwoOuterCryoW = 0.5; // cm
......@@ -1020,8 +1020,8 @@ void SpdCommonGeoMapper::DefineQslGeometrySet()
theTsTECDefGeoType = 3; // TsTEC
theTsTECMinDist3 = 121.; // cm
theTsTECLength3 = 28.6; // cm
theTsTECMinDist3 = 126.6; // cm
theTsTECLength3 = 24.0; // cm
theTsTECSize = 80.; // cm
theTsTECWidth = 80.-15.; // cm
......@@ -1040,9 +1040,9 @@ void SpdCommonGeoMapper::DefineQslGeometrySet()
theEcalTBBaseMaterial = "air";
theEcalTBLen2 = 371.2; // (510,TOR) cm
theEcalTBSize2 = 163.4; // (265,TOR) cm
theEcalTBWidth2 = 55.0; // (60,TOR) cm
theEcalTBLen2 = 372.0; // (510,TOR) cm
theEcalTBSize2 = 168.4; // (265,TOR) cm
theEcalTBWidth2 = 57.0; // (60,TOR) cm
/*===== ECALTEC (TOROIDAL ELECTROMAGNETIC CALORIMETER ENDCAPS) =====*/
......@@ -1050,7 +1050,7 @@ void SpdCommonGeoMapper::DefineQslGeometrySet()
theEcalTECBaseMaterial = "air";
theEcalTECMinDist1 = 188.6; //201.6; // (320,HYB) cm
theEcalTECMinDist1 = 204.0; // (320,HYB) cm
theEcalTECSize = 155.1; // 163.4; // (265,HYB) cm
theEcalTECWidth = 155.4 - 32.14/2; // (255,HYB) cm
......@@ -1060,18 +1060,18 @@ void SpdCommonGeoMapper::DefineQslGeometrySet()
theRsTBDefGeoType = 2;
theRsTBLen2 = 493.2; // (510,TOR) cm
theRsTBSize2 = 331.4; // (340,TOR) cm
theRsTBLen2 = 523.2; // 493.2 cm
theRsTBSize2 = 341.4; // 331.4 cm
theRsTBWidth2 = 139.0; // (60,TOR) cm
/*============ RSTEC (TOROIDAL RANGE SYSTEM ENDCAPS) ===============*/
theRsTECDefGeoType = 1;
theRsTECMinDist1 = 246.6; // (400,HYB) cm
theRsTECMinDist1 = 261.6; // 246.6 cm
theRsTECSize = 331.4; // (419,HYB) cm
theRsTECWidth = 304.4; // (409,HYB) cm
theRsTECSize = 341.4; // 331.4 cm
theRsTECWidth = 331.4; // 304.4 cm
theRsTECThickness = 139.0; // (110.5,HYB) cm
}
......
ecalt/barrel/SpdEcalB.cxx
View file @ e7a31252
......@@ -303,7 +303,6 @@ void SpdEcalB::ConstructGeometryTDR() {
Double_t barrelRadius = mapper->GetBarrelRadius();
Double_t basketPhiClearance = mapper->GetBasketPhiClearance();
Double_t basketZClearance = mapper->GetBasketZClearance();
Double_t sectorClearance = mapper->GetSectorClearance();
Int_t nSectors = mapper->GetNSectors();
Int_t nBasketsZ = mapper->GetNZBaskets();
......@@ -311,32 +310,15 @@ void SpdEcalB::ConstructGeometryTDR() {
// assuming clearance is not only between baskets, but also at the edges of barrel in Z direction
Double_t basketZSize = (barrelLength - basketZClearance*nBasketsZ)/nBasketsZ;
Double_t basketAngle = TMath::DegToRad()*360./nSectors;
//the small (< 180°) angle at the trapezoid base
Double_t basketAngle = TMath::DegToRad()*360./nSectors - basketPhiClearance/barrelRadius;
//the small (< 180 degrees) angle at the trapezoid base
Double_t angleBasketBase = (TMath::Pi() - basketAngle)/2;
//distance between adjacent baskets in direction parallel to the bases
//in other words, this is subtracted from basket φ dimensions
Double_t sectorClearanceZProj = (sectorClearance + basketPhiClearance)/TMath::Sin(angleBasketBase);
Double_t basketPhiInnerSize = barrelRadius*TMath::Tan(basketAngle/2)*2 - sectorClearanceZProj;
Double_t basketPhiOuterSize = (barrelRadius + cellThickness)*TMath::Tan(basketAngle/2)*2 - sectorClearanceZProj;
Double_t basketVertices[16];
basketVertices[0] = basketPhiInnerSize/2; basketVertices[1] = barrelRadius;
basketVertices[2] = -basketPhiInnerSize/2; basketVertices[3] = barrelRadius;
basketVertices[4] = -basketPhiOuterSize/2; basketVertices[5] = barrelRadius + cellThickness;
basketVertices[6] = basketPhiOuterSize/2; basketVertices[7] = barrelRadius + cellThickness;
basketVertices[8] = basketPhiInnerSize/2; basketVertices[9] = barrelRadius;
basketVertices[10] = -basketPhiInnerSize/2; basketVertices[11] = barrelRadius;
basketVertices[12] = -basketPhiOuterSize/2; basketVertices[13] = barrelRadius + cellThickness;
basketVertices[14] = basketPhiOuterSize/2; basketVertices[15] = barrelRadius + cellThickness;
fBasket = gGeoManager->MakeArb8("EcalBasket", Air, basketZSize/2, basketVertices);
fBasket = gGeoManager->MakeTubs("EcalBasket", Air, barrelRadius, barrelRadius + cellThickness, basketZSize/2, 90.0 - TMath::RadToDeg()*basketAngle/2, 90.0 + TMath::RadToDeg()*basketAngle/2); //by convention, starting with x=0, y=barrelRadius
fBasket->SetFillColor(kBlue);
fBasket->SetLineColor(kBlue);
fBasket->SetTransparency(50);
Double_t cellSizeTheta = mapper->GetCellZSize();
Double_t cellInnerSizePhi = mapper->GetCellInnerSizePhi();
......@@ -361,9 +343,9 @@ void SpdEcalB::ConstructGeometryTDR() {
mapper->SetBasketNModulesPhi(nModulesPhi);
cout << "-I- <SpdEcalB::ConstructGeometryTDR> : cell Z size set to " << cellSizeTheta << " cm" << endl;
cout << "-I- <SpdEcalB::ConstructGeometryTDR> : cell φ inner size set to " << cellInnerSizePhi << " cm" << endl;
cout << "-I- <SpdEcalB::ConstructGeometryTDR> : cell φ outer size set to " << cellOuterSizePhi << " cm" << endl;
cout << "-I- <SpdEcalB::ConstructGeometryTDR> : cell φ inner size corresponds to angle of " << TMath::RadToDeg()*2*TMath::ATan((cellInnerSizePhi/2)/barrelRadius) << " degrees" << endl;
cout << "-I- <SpdEcalB::ConstructGeometryTDR> : cell phi inner size set to " << cellInnerSizePhi << " cm" << endl;
cout << "-I- <SpdEcalB::ConstructGeometryTDR> : cell phi outer size set to " << cellOuterSizePhi << " cm" << endl;
cout << "-I- <SpdEcalB::ConstructGeometryTDR> : cell phi inner size corresponds to angle of " << TMath::RadToDeg()*2*TMath::ATan((cellInnerSizePhi/2)/barrelRadius) << " degrees" << endl;
vector<Double_t> rotAnglesModules;
Double_t currAngle = -basketAngle/2 + moduleAngle/2;
......@@ -462,7 +444,6 @@ void SpdEcalB::ConstructGeometryTDR() {
fModule->SetFillColor(kBlue);
fModule->SetLineColor(kBlue);
for (int iz = 0; iz < nModulesZ; ++iz) {
for (int iphi = 0; iphi < nModulesPhi; ++iphi) {
TGeoRotation * rot = new TGeoRotation();
......@@ -502,7 +483,6 @@ void SpdEcalB::ConstructGeometryTDR() {
}
}
}
......
external/GenFit2/fitters/src/KalmanFitterRefTrack.cc
View file @ e7a31252
......@@ -706,6 +706,14 @@ bool KalmanFitterRefTrack::prepareTrack(Track* tr, const AbsTrackRep* rep, bool
referenceState->setForwardNoiseMatrix(FNoiseMatrix_);
referenceState->setForwardDeltaState(forwardDeltaState_);
TVectorD st = referenceState->getState();
for (int j=0; j<st.GetNrows(); j++) {
if (!TMath::Finite(st[j])) {
Exception exc("KalmanFitterRefTrack::prepareTrack ==> nan or inf propageted t referenceState ",__LINE__,__FILE__);
exc.setFatal();
throw exc;
}
}
referenceState->resetBackward();
fitterInfo->setReferenceState(referenceState);
......
external/KFParticle/KFParticle/KFParticleBase.cxx
View file @ e7a31252
......@@ -3010,7 +3010,7 @@ float KFParticleBase::GetDeviationFromParticle( const KFParticleBase &p ) const
float mP1[8], mC1[36];
float mP2[8], mC2[36];
Transport(ds[0], dsdr[0], mP1, mC1, dsdr[1], F1, F2);
Transport(ds[0], dsdr[0], mP1, mC1, dsdr[1], F1, F2);
p.Transport(ds[1], dsdr[3], mP2, mC2, dsdr[2], F4, F3);
MultQSQt(F2, p.fC, V0Tmp, 6);
......@@ -3019,6 +3019,8 @@ float KFParticleBase::GetDeviationFromParticle( const KFParticleBase &p ) const
for(int iC=0; iC<6; iC++)
mC1[iC] += V0Tmp[iC] + mC2[iC] + V1Tmp[iC];
InvertCholetsky3(mC1);
float d[3]={ mP2[0]-mP1[0], mP2[1]-mP1[1], mP2[2]-mP1[2]};
return ( ( mC1[0]*d[0] + mC1[1]*d[1] + mC1[3]*d[2])*d[0]
......@@ -3452,4 +3454,3 @@ void KFParticleBase::MultQSQt( const float Q[], const float S[], float SOut[], c
// 72-charachters line to define the printer border
//3456789012345678901234567890123456789012345678901234567890123456789012
macro/README.md 0 → 100644
View file @ e7a31252
This file describes the most important files in this directory.
# `rootlogon.C` and `rootlogoff.C`.
These files are needed to load necessary dynamic libraries and headers before running your script or properly exiting the ROOT session. The `spdroot.py` script is a wrapper that runs `rootlogon.C` first, then your scripts, and finally `rootlogoff.C`. You can copy them to your working directory and use `root' instead.
# Simulation
Historically, the reference simulation file is called `SimQslPy8.C`. This is the only simulation file that is kept up to date. You can find a Pythia6 usage example in `SimuQslPy6.C`, but the detector description there is outdated.
# Reconstruction
Historically, the reference reconstruction file is `analysis/RecoEventFull.C`. Similar to `SimuQslPy8.C`, this is the only file that is kept up to date.
# Performance tests
The `performance-tests` folder contains maintained simulation files and reference numbers for detector performance.
# Other Files
The other files can be useful for various purposes, such as visualizing detectors, checking radiation length, etc., but they are mostly not maintained. You can use them as a starting point for your studies.
\ No newline at end of file
macro/SimuQslPy8.C
View file @ e7a31252
This diff is collapsed.
macro/analysis/RecoEventFull.C
View file @ e7a31252
......@@ -304,14 +304,57 @@ void RecoEventFull()
// Output: rc-vertices + vertices-fit-parameters
SpdRCVerticesFinder* vtxs_finder = new SpdRCVerticesFinder();
//vtxs_finder->SaveVertices(false);
vtxs_finder->SetMinTrackPVchi2(0.);
vtxs_finder->SetFitSecondaries(false); // FIT PRIMARY VERTECES ONLY
vtxs_finder->SetVerboseLevel(1);
Run->AddTask(vtxs_finder);
// You may use this code to reconstruct vertices based on MC-truth as "PID".
// See K0 decay for and example standalone use of KFParticle (to be used for
// reaslistic PID)
//=======================================================================//
// [REAL SECONDARY VERTEX FINDER (FITTED MC-TRACKS)]
// Input: mc-event, mc-particles,
// mc-tracks + track-fit-parameters,
// primary rc-vertex + primary rc-vertex-fit-parameters
// Output: secondary rc-vertices + secondary vertices-fit-parameters
//---------------------
//
// SetMinItsHits(Int_t) - minimum number of hits in ITS for each track (default = 3)
// SetTrackSelOption(Int_t) - if option: > 0 : hard cut for track selection will be used (default = true)
// SetUsedTypeOfPV(Int_t) - type of used primary vertex (PV): sim vertex(0), reco vertex (1 = default)
// SetConstrainToPV(Int_t) - constrain with PV (1 = default), no constrain (0)
// SetMinPVTrackChi2(Double_t) - minimum chi2 track to PV (default = 0.5)
// SetMaxTwoTracksChi2(Double_t) - maximum chi2 between tracks (default = 10.)
// SetMinimumMass(Double_t) - minimum mass window edge (default = 0.0 GeV)
// SetMaximumMass(Double_t) - maximum mass window edge (default = 2.5 GeV)
//
/*
SpdRCKFpartV0Finder* v0_finder = new SpdRCKFpartV0Finder();
v0_finder->SetMinItsHits(3);
v0_finder->SetTrackSelOption(0);
v0_finder->SetUsedTypeOfPV(1);
v0_finder->SetConstrainToPV(1);
v0_finder->SetMinPVTrackChi2(2.0);
v0_finder->SetMaxTwoTracksChi2(2.0);
v0_finder->SetMinimumMass(0);
v0_finder->SetMaximumMass(10);
// search for decay mode of K0
v0_finder->AddVertexCandidate(-211,211);
v0_finder->SetVerboseLevel(1);
Run->AddTask(v0_finder);
*/
// **************************************************************************************************************
// **************************************************************************************************************
//--------------------------------------------------------//
// [REAL ECAL CLUSTERIZATION AND
......@@ -392,6 +435,16 @@ void RecoEventFull()
//========================================================================//
SpdMCAegParticleProducer* mcaeg_part = new SpdMCAegParticleProducer();
//
//mctof_part->SaveParticles(false);
mcaeg_part->SetVerboseLevel(1);
//
Run->AddTask(mcaeg_part);
cout << "Start ... " << endl;
Run->Initialize();
......
macro/examples/chic/analyze_chic.C
View file @ e7a31252
......@@ -19,19 +19,8 @@
void analyze_chic(int nevents_max = 1000, int chic_pdg = 20443)
{
TFile* fdata = new TFile( "reco_full.root" );
if (!fdata) {
cout << "Input file not found" << endl;
return;
}
TTree* datatree = (TTree*)fdata->Get("spdsim");
if (!datatree) {
cout << "No tree in the file" << endl;
return;
}
// The output histograms
// Book output histograms
TH1D* hangle = new TH1D("hangle", "Angle between cluster and closest MC photon;#alpha [rad];counts", 40, 0, 0.1);
TH1D* hE = new TH1D("hE", "#chi_{c} photon energy spectrum;E_{#gamma} [GeV];counts", 40, 0, 2);
TH1D* hdM = new TH1D("hdM", "M(#gamma#mu^{+}#mu^{-}) - M(#mu^{+}#mu^{-}) (#alpha<0.05);M(#gamma#mu^{+}#mu^{-}) - M(#mu^{+}#mu^{-}) [GeV];counts", 100, 0, 1);
......@@ -40,39 +29,37 @@ void analyze_chic(int nevents_max = 1000, int chic_pdg = 20443)
TH1D* hVtxZ = new TH1D( "hVtxZ", "Primary vertex rec. Z;Z [cm];counts", 500, -100, 100);
TH2D* hVtxXY = new TH2D( "hVtxXY", "Primary vertex rec. XY;X [cm];Y [cm]", 50, -1, 1, 50, -1, 1);
TClonesArray* mcparticles = nullptr;
TClonesArray* rcecalpartilces = nullptr;
TClonesArray* mctracks = nullptr; // - mc-tracks "MCTracks" TClonesArray[SpdTrackMC(+SpdTrackFitPar)]
TClonesArray* rcvertices = nullptr;
// Attaching branches
datatree->SetBranchAddress("MCParticles", &mcparticles) ;
datatree->SetBranchAddress("RCEcalParticles", &rcecalpartilces) ;
datatree->SetBranchAddress("MCTracks", &mctracks) ;
datatree->SetBranchAddress("RCVertices", &rcvertices) ;
// Prepare event reader
SpdMCDataIterator IT; // Do I need to read pars?
IT.AddSourceFile("reco.root");
IT.ActivateBranch("all");
IT.Init();
// aliases for frequently used brances
const TClonesArray* mcparticles = IT.GetParticles();
const TClonesArray* mctracks = IT.GetTracks();
const TClonesArray* rcecalpartilces = IT.GetEcalParticlesRC();
const TClonesArray* rcvertices = IT.GetVerticesRC();
SpdTrackPropagatorGF trkProp;
// The main event loop
int nevents = datatree->GetEntries();
cout << "The mumber of events in the tree is " << nevents << endl;
for (int iev = 0; iev < std::min(nevents, nevents_max); iev++ ) {
datatree->GetEntry( iev );
// Loop over collection of SpdMCParticle to find
// the photon from the chic_c decay
int ne_total(0);
while ( (IT.NextEvent()) && (ne_total < nevents_max) )
{
// Find photon candidate for chic -> gamma J/psi as the closest
// to true photon from SpdEmcParticles.
// First, fined the refrence photon from MC-truth.
SpdMCParticle *etac_photon = nullptr;
for (int ipart=0; ipart<std::min(mcparticles->GetEntriesFast(), 100); ipart++) {
SpdMCParticle* part = (SpdMCParticle*) mcparticles->At(ipart);
//cout << part->GetPdgCode() << "\t" << part->GetId() << "\t" << part->GetMotherId()
// << "\t" << part->GetGeneration() << "\t" << part->GetSpecialMotherId() << endl;
if ( part->GetPdgCode() != 22 ) continue;
int mid = part->GetMotherId();
if ( mid >= 0 ) {
SpdMCParticle* mother_part = (SpdMCParticle*) mcparticles->At( mid );
if (mother_part->GetPdgCode() == chic_pdg) {
etac_photon = part;
//cout << "Found eta_c photon: " << ipart << endl;
}
}
}
......@@ -95,10 +82,11 @@ void analyze_chic(int nevents_max = 1000, int chic_pdg = 20443)
break;
}
}
if (!found_vrec) continue;
hVtxZ->Fill( vrec.Z() );
hVtxXY->Fill( vrec.X(), vrec.Y() );
if (!found_vrec) continue;
// Now find the closet EMC shower (particle)
double min_angle = 1e10;
......@@ -108,7 +96,7 @@ void analyze_chic(int nevents_max = 1000, int chic_pdg = 20443)
SpdEcalRCParticle* ecal_part = (SpdEcalRCParticle*) rcecalpartilces->At(ishower);
if (!ecal_part) continue;
TVector3 shower_pos = ecal_part->GetPosition();
if (ecal_part->GetEnergy() < 0.1) {
if (ecal_part->GetEnergy() < 0.4) {
continue;
}
double angle = (shower_pos - vrec).Angle(ptrue);
......@@ -140,13 +128,9 @@ void analyze_chic(int nevents_max = 1000, int chic_pdg = 20443)
if (!mctrack) continue; // why??
if ( TMath::Abs( mctrack->GetParticlePdg() ) == 13 ) {
// Get the corresponding paritlce and check id it
// originates from J/psi
int partid = mctrack->GetParticleId();
SpdMCParticle* part = (SpdMCParticle*) mcparticles->At( partid );
//if ( part->GetSpecialMotherId() != 5 ) continue;
// Extrapolate to Ecal
SpdTrackFitPar* fitpars = mctrack->GetFitPars();
if ( !fitpars ) continue;
......@@ -170,6 +154,7 @@ void analyze_chic(int nevents_max = 1000, int chic_pdg = 20443)
TLorentzVector pjpsi = pmup + pmum;
hdM->Fill( (pjpsi + p4phot).Mag() - pjpsi.Mag() );
}
ne_total++;
} // events
gStyle->SetOptStat(11);
......@@ -190,7 +175,7 @@ void analyze_chic(int nevents_max = 1000, int chic_pdg = 20443)
c->cd(5);
hdM->Draw();
hdM->Fit("gaus");
hdM->Fit("gaus", "L", "", 0.3, 0.5);
c->cd(3);
hVtxZ->Draw();
......
macro/examples/jpsi-ee/reco_event_jpsi.C → macro/examples/chic/reco.C
View file @ e7a31252
//______________________________________________________________________________
void reco_event_jpsi()
void reco()
{
Int_t nEvents = 0;
TString inFile = "run.root"; // Input file (MC data)
TString inFile = "simu.root"; // Input file (MC data)
TString parFile = "params.root"; // Input file with parameters
TString outFile = "reco_full.root"; // Output file
TString outFile = "reco.root"; // Output file
SpdRunAna* Run = new SpdRunAna(); // Main task
......@@ -49,10 +49,15 @@ void reco_event_jpsi()
// - its-points "SpdItsPoint" TClonesArray[SpdItsPoint]
// - ts-barrel-points "SpdTsTB2Point" TClonesArray[SpdTsTB2Point]
// - ts-endcaps-points "SpdTsTEC2Point" TClonesArray[SpdTsTEC2Point]
// - tof-barrel-points "SpdTofBPoint" TClonesArray[SpdTofBPoint]
// - tof-endcaps-points "SpdTofECPoint" TClonesArray[SpdTofECPoint]
// - ecal-barrel-points "SpdEcalTB2Point" TClonesArray[SpdEcalTB2Point]
// - ecal-endcaps-points "SpdEcalTEC2Point" TClonesArray[SpdEcalTEC2Point]
// - rs-barrel-points "SpdRsTB2Point" TClonesArray[SpdRsTB2Point]
// - rs-endcaps-points "SpdRsTEC2Point" TClonesArray[SpdRsTEC2Point]
// - bbc-points "SpdBbcPoint" TClonesArray[SpdBbcPoint]
// - aeg-points "SpdAegPoint" TClonesArray[SpdAegPoint]
// - zdc-points "SpdZdcPoint" TClonesArray[SpdZdcPoint]
//
//--------------------------------------------------------------------------------------------
//
......@@ -66,13 +71,23 @@ void reco_event_jpsi()
// - mc-vertices SpdMCEventMaker
// - mc-its-hits SpdItsMCHitProducer
// - mc-ts-hits SpdTsMCHitProducer
// - mc-tof-hits SpdTofMCHitProducer
// - mc-ecal-hits SpdEcalMCHitProducer
// - mc-rs-hits SpdRsMCHitProducer
// - mc-bbc-hits SpdBbcMCHitProducer
// - mc-aeg-hits SpdAegMCHitProducer
// - mc-zdc-hits SpdZdcMCHitProducer
// - mc-tracks SpdMCTracksFinder (+ fit pars., optionally)
// - vertices-fit-pars.(mc) SpdMCVerticesFitter (mc-vertices update)
// - rc-vertices SpdRCVerticesFinder (+ fit pars. )
// - rc-ecal-clusters SpdEcalRCMaker
// - rc-ecal-particles SpdEcalRCMaker
// - mc-ecal-clusters-info SpdEcalClusterMCInfoMaker
// - mc-ecal-particles SpdEcalClusterMCInfoMaker
// - mc-rs-clusters SpdRsMCClusterMaker
// - mc-rs-particles SpdRsMCClusterMaker
// - ts-particles SpdMCTsParticleProducer
// - tof-particles SpdMCTofParticleProducer
//
// (full data list): DATA BRANCH NAME: BRANCH OBJECT TYPE:
//
......@@ -83,11 +98,21 @@ void reco_event_jpsi()
// - mc-its-hits "ItsMCHits" TClonesArray[SpdMCSiliconHit]
// - mc-ts-hits "TsMCHits" TClonesArray[SpdMCStrawHit(1D/2D)]
// - mc-ecal-hits "EcalMCHits" TClonesArray[SpdEcalMCHit]
// - mc-tof-hits "TofMCHits" TClonesArray[SpdTofMCHit]
// - mc-rs-hits "RsMCHits" TClonesArray[SpdRsMCHit]
// - mc-bbc-hits "BbcMCHits" TClonesArray[SpdBbcMCHit]
// - mc-aeg-hits "AegMCHits" TClonesArray[SpdAegMCHit]
// - mc-zdc-hits "ZdcMCHits" TClonesArray[SpdZdcMCHit]
// - mc-tracks "MCTracks" TClonesArray[SpdTrackMC(+SpdTrackFitPar)]
// - rc-vertices "RCVertices" TClonesArray[SpdVertexRC(+SpdVertexFitPar)]
// - rc-ecal-clusters "RCEcalClusters" TClonesArray[SpdEcalRCCluster]
// - rc-ecal-particles "RCEcalParticles" TClonesArray[SpdEcalRCParticle]
// - mc-ecal-clusters-info "MCEcalClustersInfo" TClonesArray[SpdEcalClusterMCInfo]
// - mc-ecal-particles "MCEcalParticles" TClonesArray[SpdEcalMCParticle]
// - mc-rs-clusters "MCRsClusters" TClonesArray[SpdRsMCCluster]
// - mc-rs-particles "MCRsParticles" TClonesArray[SpdRsMCParticle]
// - ts-particles "TsParticles" TClonesArray[SpdTsParticle]
// - tof-particles "TofParticles" TClonesArray[SpdTofParticle]
//
//--------------------------------------------------------------------------------------------
//
......@@ -151,6 +176,20 @@ void reco_event_jpsi()
Run->AddTask(ts_hits_producer);
//--------------------------------------------------------//
// [TIME-OF-FLIGHT (TOF) SYSTEM (BARREL+CAPS) HIT PRODUCER]
// Input: mc-event, mc-particles,
// tof-barrel-points and/or tof-endcaps-points
// Output: mc-tof-hits
SpdTofMCHitProducer* tof_hits_producer = new SpdTofMCHitProducer();
//tof_hits_producer->SaveHits(false); // default: true
tof_hits_producer->SetVerboseLevel(1);
Run->AddTask(tof_hits_producer);
//--------------------------------------------------------//
// [ELECTROMAGNETIC CALORIMETER (BARREL+CAPS) HIT PRODUCER]
// Input: mc-event, mc-particles,
......@@ -173,12 +212,53 @@ void reco_event_jpsi()
SpdRsMCHitProducer* rs_hits_producer = new SpdRsMCHitProducer();
//rs_hits_producer->MakeStripHits(true);
//rs_hits_producer->SaveHits(false); // default: true
rs_hits_producer->SetVerboseLevel(1);
Run->AddTask(rs_hits_producer);
//--------------------------------------------------------//
// [BEAM_BEAM COUNTER HIT PRODUCER]
// Input: mc-event, mc-particles, bbc-points
// Output: mc-bbc-hits
SpdBbcMCHitProducer* bbc_hits_producer = new SpdBbcMCHitProducer();
//bbc_hits_producer->SaveHits(false); // default: true
bbc_hits_producer->SetVerboseLevel(1);
Run->AddTask(bbc_hits_producer);
//--------------------------------------------------------//
// [AEROGEL HIT PRODUCER]
// Input: mc-event, mc-particles, aeg-points
// Output: mc-aeg-hits
SpdAegMCHitProducer* aeg_hits_producer = new SpdAegMCHitProducer();
//aeg_hits_producer->SaveHits(false); // default: true
aeg_hits_producer->SetVerboseLevel(1);
Run->AddTask(aeg_hits_producer);
//--------------------------------------------------------//
// [ZERO-DEGREE CALORIMETER HIT PRODUCER]
// Input: mc-event, mc-particles, zdc-points
// Output: mc-zdc-hits
SpdZdcMCHitProducer* zdc_hits_producer = new SpdZdcMCHitProducer();
//zdc_hits_producer->SaveHits(false); // default: true
zdc_hits_producer->SetVerboseLevel(1);
Run->AddTask(zdc_hits_producer);
//--------------------------------------------------------//
// [MC-TRACK FINDER]
// Input: mc-event, mc-particles,
......@@ -194,6 +274,11 @@ void reco_event_jpsi()
//track_fitter->SetFitterMaxTrials(2); // default: 2
//track_fitter->SetStartSeedMethod(1); // (0,1), default: 1
//track_fitter->StoreImPoints(true); // default: false
//track_fitter->StoreImMomentum(true); // default: true
//track_fitter->StoreImPosition(true); // default: false
//track_fitter->StoreImCovariance(true); // default: false
track_fitter->SetVerboseLevel(1);
track_finder->SetVerboseLevel(1);
......@@ -243,14 +328,84 @@ void reco_event_jpsi()
Run->AddTask(ecal_rc);
//--------------------------------------------------------//
// [MC-TRUTH FOR ELECTROMAGNETIC CALORIMETER CLUSTERS]
// Input: mc-events, mc-particles,
// mc-ecal-hits, rc-ecal-clusters
// Output: mc-ecal-clusters-info, mc-ecal-particles
SpdEcalClusterMCInfoMaker* ecal_mc_info = new SpdEcalClusterMCInfoMaker();
//ecal_mc_info->SaveInfo(false); // default = true
//ecal_mc_info->SaveEcalParticles(false); // default = true
//ecal_mc_info->CreateEcalParticles(false); // default = true
//ecal_mc_info->RemoveHitMCTruth(true); // default = false
ecal_mc_info->SetVerboseLevel(1);
Run->AddTask(ecal_mc_info);
//--------------------------------------------------------//
// [MC-CLUSTERING AND MC-RECONSTRUCTION FOR RANGE SYSTEM]
// Input: mc-event, mc-particles,
// mc-rs-hits
// Output: mc-rs-clusters, mc-rs-particles
SpdRsMCClusterMaker* rs_cl = new SpdRsMCClusterMaker();
//rs_cl->SaveClusters(false); // default: true
//rs_cl->SaveRsParticles(false); // default: true
//rs_cl->CreateRsParticles(false); // default = true
//rs_cl->RemoveHitMCTruth(true); // default = false
rs_cl->SetVerboseLevel(1);
Run->AddTask(rs_cl);
// ///--------------------------------------------------------//
// // [MC-PRODUCER FOR TS-PARTICLES]
// // Input: mc-event, mc-particles, mc-tracks(+ fit pars.)
// // Output: ts-particles
//
// SpdMCTsParticleProducer* mcts_part = new SpdMCTsParticleProducer();
//
// //mcts_part->SaveParticles(false);
//
// mcts_part->SetVerboseLevel(1);
//
// Run->AddTask(mcts_part);
//
// ///--------------------------------------------------------//
// // [MC-PRODUCER FOR TOF-PARTICLES]
// // Input: mc-event, mc-particles, mc-tracks(+ fit pars.), mc-tof-hits
// // Output: tof-particles
//
// SpdMCTofParticleProducer* mctof_part = new SpdMCTofParticleProducer();
//
// //mctof_part->SaveParticles(false);
//
// mctof_part->SetVerboseLevel(1);
//
// Run->AddTask(mctof_part);
//========================================================================//
SpdMCAegParticleProducer* mcaeg_part = new SpdMCAegParticleProducer();
//
//mctof_part->SaveParticles(false);
mcaeg_part->SetVerboseLevel(1);
//
Run->AddTask(mcaeg_part);
cout << "Start ... " << endl;
Run->Initialize();
//-------------------------------------//
cout << "Run ... " << endl;
......@@ -260,6 +415,8 @@ void reco_event_jpsi()
//-------------------------------------//
//OutputFile->Close();
timer.Stop();
Double_t rtime = timer.RealTime();
......@@ -271,6 +428,9 @@ void reco_event_jpsi()
cout << endl;
//SpdMCEventHelper::PrintMCProcesses();
//SpdMCEventHelper::Instance()->PrintParticles();
//SpdMCEventHelper::Instance()->PrintVertices();
}
macro/examples/chic/rootlogoff.C deleted 100644 → 0
View file @ 3e5006bc
void rootlogoff()
{
//cout << "::rootlogoff::" << endl;
if (gROOT) {
if (gROOT->GetVersionInt() > 61000) gROOT->~TROOT();
}
}
\ No newline at end of file
macro/examples/chic/rootlogon.C deleted 100644 → 0
View file @ 3e5006bc
#define LOAD_MACROPATH
#define LOAD_SPDLIBRARIES
//#define LOAD_SPDINCLUDES
void SetMacroPath();
void LoadLibraries();
void LoadIncludes();
void rootlogon()
{
//cout << "::rootlogon::" << endl;
if (!gROOT) return;
#ifdef LOAD_MACROPATH
SetMacroPath();
#endif
#ifdef LOAD_SPDLIBRARIES
LoadLibraries();
#endif
#ifdef LOAD_SPDINCLUDES
LoadIncludes();
#endif
}
//______________________________________________________________
void SetMacroPath()
{
//cout << ":SetMacroPath:" << endl;
TString PROJECT = gSystem->Getenv("VMCWORKDIR");
if (PROJECT == "") PROJECT = ".";
//cout << "Project install directory: " << PROJECT <<endl;
TString macropath;
//macropath = gROOT->GetMacroPath();
//cout << "macropath: " << macropath << endl;
TString mpath = gSystem->ConcatFileName(PROJECT,"macro");
macropath = ".:" + mpath + ":";
gROOT->SetMacroPath(macropath);
// add macro subdirectories
//macropath += mpath + "/geom" + ":";
//macropath += mpath + "/analysis" + ":";
//macropath += mpath + "/my" + ":";
//cout << "macropath: " << macropath << endl;
gROOT->SetMacroPath(macropath);
//return;
// print macro path
int id(0),idx(0);
cout << "path to ROOT macros: " << endl;
while ((idx = macropath.Index(":",id)) != -1) {
cout << macropath(id,idx-id) << endl;
id = idx+1;
}
}
//______________________________________________________________
void LoadLibraries()
{
//cout << ":LoadLibraries:" << endl;
cout << "load libraries ... ";
gSystem->Load("libSpdData");
gSystem->Load("libSpdCommon");
gSystem->Load("libSpdField");
gSystem->Load("libSpdGen");
gSystem->Load("libSpdGeometry");
gSystem->Load("libSpdDisplay");
gSystem->Load("libClustering");
gSystem->Load("libSpdPassive");
gSystem->Load("libSpdIts");
gSystem->Load("libSpdTst");
gSystem->Load("libSpdEcalt");
gSystem->Load("libSpdRst");
gSystem->Load("libSpdTest");
// gSystem->Load("libSpdTss");
// gSystem->Load("libSpdEcals");
// gSystem->Load("libSpdRss");
// gSystem->Load("libMpdEmc");
// gSystem->Load("libSpdEcal");
// tracking
gSystem->Load("libGenFit2");
gSystem->Load("libTrackingGF");
gSystem->Load("libSpdReco");
gSystem->Load("libKFParticle");
gSystem->Load("libVertexKF");
//gSystem->ListLibraries();
cout << "Ok " << endl;
}
//______________________________________________________________
void LoadIncludes()
{
//cout << ":LoadIncludes:" << endl;
cout << "load includes ... ";
gROOT->ProcessLine(".include $SIMPATH/include/");
gROOT->ProcessLine(".include $FAIRROOTPATH/include/");
gROOT->ProcessLine(".include ../spddata/");
gROOT->ProcessLine(".include ../spddata/params");
gROOT->ProcessLine(".include ../spddata/mcdata");
gROOT->ProcessLine(".include ../spddata/IdealTrackData");
gROOT->ProcessLine(".include ../common/");
gROOT->ProcessLine(".include ../common/geometry");
gROOT->ProcessLine(".include ../common/region");
gROOT->ProcessLine(".include ../common/checks");
gROOT->ProcessLine(".include ../field/");
gROOT->ProcessLine(".include ../spdgenerators/");
gROOT->ProcessLine(".include ../spdgeometry/");
gROOT->ProcessLine(".include ../spdgeometry/sol/");
gROOT->ProcessLine(".include ../spdgeometry/its/");
gROOT->ProcessLine(".include ../spdgeometry/tst/");
gROOT->ProcessLine(".include ../spdgeometry/ecalt/");
gROOT->ProcessLine(".include ../spdgeometry/rst/");
gROOT->ProcessLine(".include ../spddisplay/");
gROOT->ProcessLine(".include ../proc/");
gROOT->ProcessLine(".include ../procmc/");
gROOT->ProcessLine(".include ../reco/");
gROOT->ProcessLine(".include ../passive");
gROOT->ProcessLine(".include ../its");
gROOT->ProcessLine(".include ../tst/");
gROOT->ProcessLine(".include ../tst/barrel");
gROOT->ProcessLine(".include ../tst/ecps");
gROOT->ProcessLine(".include ../ecalt/");
gROOT->ProcessLine(".include ../ecalt/barrel");
gROOT->ProcessLine(".include ../ecalt/ecps");
gROOT->ProcessLine(".include ../rst/");
gROOT->ProcessLine(".include ../rst/barrel");
gROOT->ProcessLine(".include ../rst/ecps");
gROOT->ProcessLine(".include ../tss/");
gROOT->ProcessLine(".include ../tss/barrel");
gROOT->ProcessLine(".include ../tss/ecps");
gROOT->ProcessLine(".include ../ecals/");
gROOT->ProcessLine(".include ../ecals/barrel");
gROOT->ProcessLine(".include ../ecals/ecps");
gROOT->ProcessLine(".include ../rss/");
gROOT->ProcessLine(".include ../rss/barrel");
gROOT->ProcessLine(".include ../rss/ecps");
gROOT->ProcessLine(".include ../test/tracker");
gROOT->ProcessLine(".include ../test/ecalcell");
gROOT->ProcessLine(".include ../test/sts/");
gROOT->ProcessLine(".include ../mpdgenerators/");
gROOT->ProcessLine(".include ../emc/");
gROOT->ProcessLine(".include ../ecal/");
gROOT->ProcessLine(".include ../proc/clustering");
gROOT->ProcessLine(".include ../proc/tracking");
gROOT->ProcessLine(".include ../proc/vertex");
gROOT->ProcessLine(".include ../external/KFParticle/KFParticle");
gROOT->ProcessLine(".include ../external/KFParticle/KFParticleTest");
cout << " Ok " << endl;
}
macro/examples/chic/run_all.sh
View file @ e7a31252
#/bin/bash
root -l -b -q simu_chic.C
root -l -b -q reco_event_chic.C
root -l -b -q analyze_chic.C
spdroot.py -l -b -q "simu.C(1000)"
spdroot.py -l -b -q reco.C
spdroot.py -l -b -q analyze_chic.C
macro/examples/chic/simu.C 0 → 100644
View file @ e7a31252
/*
This is a default simulation file for the second stage
of the SPD experiment operation. To change the central tracker
to MCT please see the exmple at
macro/performance-tests/track-fitting/its-straw
In case of issues please email to iden AT jinr.ru.
*/
// #define _COMPILE_MACRO_
#if defined(_COMPILE_MACRO_)
#include <TRint.h>
#include <TStopwatch.h>
#include "FairParRootFileIo.h"
#include "SpdRunSim.h"
#include "SpdRootFileSink.h"
#include "SpdMCEventHeader.h"
#include "SpdFields.hh"
#include "SpdGenerators.hh"
#include "SpdCommonGeoMapper.h"
#include "SpdCave.h"
#include "SpdPipe.h"
#include "SpdHybMagnet.h"
#include "SpdIts.h"
#include "SpdTsTB.h"
#include "SpdTsTEC.h"
#include "SpdEcalTB2.h"
#include "SpdEcalTEC2.h"
#include "SpdRsTB2.h"
#include "SpdRsTEC2.h"
#include "SpdTofB.h"
#include "SpdTofEC.h"
#include "SpdZdc.h"
#include "SpdItsGeoMapperX.h"
#include "SpdTsTBGeoMapper.h"
#include "SpdTsTBGeoBuilder.h"
#include "SpdTsTECGeoMapper.h"
#endif
void CustomIts(Int_t option);
void CustomTsB(Double_t deg);
//_________________________________________________________________________
void simu(Int_t nEvents=1, int seed=0)
{
TString outFile = "simu.root"; // Output data file name
TString parFile = "params.root"; // Output parameters file name
SpdRunSim *run = new SpdRunSim();
run->SetName("TGeant4");
// gSystem->Setenv("GEOMPATH","/home/artur/Projects/1/");
run->SetMaterials("media.geo");
SpdRootFileSink *OutputFile = new SpdRootFileSink(outFile);
run->SetSink(OutputFile);
TString decayer_config = "DecayConfigPythia8.C";
run->SetPythiaDecayer(decayer_config);
run->SetMCEventHeader(new SpdMCEventHeader);
/*--------------------------------------------------*/
/* +++++++++ GEOMETRY (QSOLENOID) ++++++++ */
/* +++++++++ [hystorical name] ++++++++ */
/*--------------------------------------------------*/
SpdCommonGeoMapper::Instance()->DefineQslGeometrySet();
/* ++++++++++++++++++ CAVE ++++++++++++++++++ */
FairModule *cave = new SpdCave("CAVE");
// cave->SetGeometryFileName("cave.geo");
// cave->SetGeometryFileName("cave_precise.geo");
run->AddModule(cave);
/* ++++++++++++++++++ PIPE ++++++++++++++++++ */
SpdPipe *Pipe = new SpdPipe();
run->AddModule(Pipe);
/* +++++++++ MAGNET +++++++++ */
SpdSolMagnet2 *Magnet = new SpdSolMagnet2();
run->AddModule(Magnet);
/* +++++++++++++++++ DETECTORS ++++++++++++++ */
SpdTsTB *ts_barrel = new SpdTsTB(); /* +++++++++ TS (BARREL) ++++++++++++ */
SpdTsTEC *ts_ecps = new SpdTsTEC(); /* +++++++++ TS (ENDCAPS) +++++++++++ */
SpdTofB *tof_barrel = new SpdTofB(); /* +++++++++ TOF (BARREL) +++++++++++ */
SpdTofEC *tof_ecps = new SpdTofEC(); /* +++++++++ TOF (ENDCAPS) ++++++++++ */
SpdEcalB *ecal_barrel = new SpdEcalB(); /* +++++++++ ECAL (BARREL) ++++++++++ */
SpdEcalEC *ecal_ecps = new SpdEcalEC(); /* +++++++++ ECAL (ENDCAPS) +++++++++ */
SpdRsTB2 *rs_barrel = new SpdRsTB2(); /* +++++++++ RS (BARREL) ++++++++++++ */
SpdRsTEC2 *rs_ecps = new SpdRsTEC2(); /* +++++++++ RS (ENDCAPS) +++++++++++ */
SpdBbc *bbc = new SpdBbc(); /* +++++++++ BEAM-BEAM COUNTER +++++++ */
SpdAeg *aeg = new SpdAeg(); /* +++++++++ AEROGEL +++++++++++++++++ */
SpdZdc *zdc = new SpdZdc(); /* +++++++++ Z0 CALORIMETER ++++++++++ */
run->AddModule(ts_barrel);
run->AddModule(ts_ecps);
run->AddModule(tof_barrel);
run->AddModule(tof_ecps);
run->AddModule(ecal_barrel);
run->AddModule(ecal_ecps);
run->AddModule(rs_barrel);
run->AddModule(rs_ecps);
run->AddModule(bbc);
run->AddModule(aeg);
run->AddModule(zdc);
/* ===== Vertex detector ===== */
SpdIts *its = new SpdIts();
run->AddModule(its);
// its->SaveEmptyHits();
/*--------------------------------------------------*/
/* ++++++++++++ CUSTOMIZE GEOMETRY +++++++++++++ */
/*--------------------------------------------------*/
CustomIts(1);
CustomTsB(-1.);
/*--------------------------------------------------*/
/* ++++++++++++ CUSTOMIZE MATERIALS +++++++++++++ */
/*--------------------------------------------------*/
// SpdParameter* par;
/* --- alter TsTB straw active media (gas) --- */
// par = ts_barrel->GetMapper()->AddParameter("TsTBBaseStrawMaterial");
// if (par) *par = "arco28020x1p5"; // ArCO2, 80/20 %, density x 1.5
// if (par) *par = "arco28020x2p0"; // ArCO2, 80/20 %, density x 2.0
// if (par) *par = "arco27030x1p5"; // ArCO2, 70/30 %, density x 1.5
// if (par) *par = "arco27030x2p0"; // ArCO2, 70/30 %, density x 2.0
/* --- alter TsTEC straw active media (gas) --- */
// par = ts_ecps->GetMapper()->AddParameter("TsTECBaseStrawMaterial");
// if (par) *par = "arco28020x1p5"; // ArCO2, 80/20 %, density x 1.5
// if (par) *par = "arco28020x2p0"; // ArCO2, 80/20 %, density x 2.0
// if (par) *par = "arco27030x1p5"; // ArCO2, 70/30 %, density x 1.5
// if (par) *par = "arco27030x2p0"; // ArCO2, 70/30 %, density x 2.0
/* --- reset all module's materials ---*/
// ts_barrel->GetMapper()->UnsetMaterials("air"); // "air" or "vacuum"
// ts_ecps->GetMapper()->UnsetMaterials("air"); // "air" or "vacuum"
// its->GetMapper()->UnsetMaterials("air"); // "air" or "vacuum"
/*--------------------------------------------------*/
/* ++++++++++++++ MAGNETIC FIELD ++++++++++++++ */
/*--------------------------------------------------*/
SpdFieldMap1_8 *MagField = new SpdFieldMap1_8("full_map");
MagField->InitData("field_full1_8.bin");
SpdRegion *reg = MagField->CreateFieldRegion("box");
reg->SetBoxRegion(-330, 330, -330, 330, -386, 386); // (X,Y,Z)_(min,max), cm
run->SetField(MagField);
MagField->Print();
/*--------------------------------------------------*/
/* ++++++++++ DEFINE PRIMARY GENERATORS +++++++++++ */
/*--------------------------------------------------*/
SpdPrimaryGenerator *primGen = new SpdPrimaryGenerator();
//============================
// PYTHIA 8 GENERATOR
//============================
SpdPythia8Generator *P8gen = new SpdPythia8Generator();
P8gen -> SetSeed(seed);
P8gen->SetBeam(2212, 2212, 27.); // pdg(A), pdg(B), E_cms (CM energy of collision)
P8gen -> SetParameters("Charmonium:gg2ccbar(3PJ)[3PJ(1)]g = off,on,off");
P8gen -> SetParameters("Charmonium:qg2ccbar(3PJ)[3PJ(1)]q = off,on,off");
P8gen -> SetParameters("Charmonium:qqbar2ccbar(3PJ)[3PJ(1)]g = off,on,off");
P8gen -> SetParameters("Charmonium:gg2ccbar(3PJ)[3S1(8)]g = off,on,off");
P8gen -> SetParameters("Charmonium:qg2ccbar(3PJ)[3S1(8)]q = off,on,off");
P8gen -> SetParameters("Charmonium:qqbar2ccbar(3PJ)[3S1(8)]g = off,on,off");
// chi_c decays to gamma + J/psi
P8gen -> SetParameters("20443:oneChannel = 1 1 0 22 443");
// Jpsi decays to mu+ mu-
P8gen->SetParameters("443:onMode = off");
P8gen->SetParameters("443:onIfAny = 13 -13");
// P8gen->SetVerboseLevel(-1);
primGen->AddGenerator(P8gen);
run->SetGenerator(primGen);
primGen->SetBeam(0., 0., 0.1, 0.1); // (X,Y)- position, (X,Y)- (2*delta or sigma) [cm]
// primGen->SmearVertexXY(kTRUE);
primGen->SmearGausVertexXY(kTRUE);
primGen->SetTarget(0., 30.);// beam Z-position and smearing (2*delta or sigma), [cm]
// primGen->SmearVertexZ(kTRUE);
primGen->SmearGausVertexZ(kTRUE);
primGen->SetVerboseLevel(-10);
primGen->SetVerbose(0);
/* ------------------------------------------------ */
/* +++++++++++++++ GLOBAL OPTIONS +++++++++++++++++ */
/* ------------------------------------------------ */
// run->SetStoreTraj(kTRUE);
// SpdCommonGeoMapper::Instance()->SaveEmptyHits();
// run->ForceParticleLifetime(-211, 26.033/5.); // pdg, life time [ns]
// run->ForceParticleLifetime( 211, 26.033/5.); // pdg, life time [ns]
// SpdCommonGeoMapper::Instance()->UnsetMaterials(false);
/* ----------------------------------------------------------------------- */
/* >>>>>>>>>>>>>>>>>>>>>>>>>>> INITALIZE RUN <<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
/* ----------------------------------------------------------------------- */
cout << "\n\t" << "++++++++++++++++++++++++++++++++++++++++++++" << endl;
cout << "\t" << "+ +" << endl;
cout << "\t" << "+ Init Run (start) +" << endl;
cout << "\t" << "+ +" << endl;
cout << "\t" << "++++++++++++++++++++++++++++++++++++++++++++\n" << endl;
run->Init();
cout << "\n\t" << "++++++++++++++++++++++++++++++++++++++++++++" << endl;
cout << "\t" << "+ +" << endl;
cout << "\t" << "+ Init Run (finish) +" << endl;
cout << "\t" << "+ +" << endl;
cout << "\t" << "++++++++++++++++++++++++++++++++++++++++++++\n" << endl;
/*--------------------------------------------------*/
/* +++++++++++++ CREATE RUN PARAMETERS +++++++++++ */
/*--------------------------------------------------*/
Bool_t MergePars = kFALSE;
FairParRootFileIo *parOut = new FairParRootFileIo(MergePars);
if (MergePars)
parOut->open(parFile.Data());
else
parOut->open(parFile.Data(), "RECREATE");
FairRuntimeDb *rtdb = run->GetRuntimeDb();
rtdb->setOutput(parOut);
/* ----------------------------------------------------------------------- */
/* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> RUN <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
/* ----------------------------------------------------------------------- */
TStopwatch timer;
timer.Start();
run->Run(nEvents);
timer.Stop();
/*--------------------------------------------------*/
/* ++++++++++++ SAVE RUN PARAMETERS +++++++++++++ */
/*--------------------------------------------------*/
rtdb->saveOutput();
/*--------------------------------------------------*/
Double_t rtime = timer.RealTime();
Double_t ctime = timer.CpuTime();
//-------------------------------------------------------------
cout << endl
<< endl;
cout << "Macro finished succesfully." << endl;
cout << "Output file is " << outFile << endl;
cout << "Parameter file is " << parFile << endl;
cout << "Real time " << rtime << " s, CPU time " << ctime << "s" << endl;
cout << endl;
/*--------------------------------------------------*/
// SpdCommonGeoMapper::Instance()->PrintGeometryList();
SpdCommonGeoMapper::Instance()->PrintGeometry();
// SpdItsGeoMapperX::Instance()->PrintVolPars();
// SpdItsGeoMapperX::Instance()->PrintGeoTable();
// SpdItsGeoMapperX::Instance()->Print("");
/*--------------------------------------------------*/
gApplication->Terminate();
}
//___________________________________________________________________________
void CustomIts(Int_t option)
{
if (option < 1)
return;
Double_t mm = 0.1;
SpdItsGeoMapperX *mapper = SpdItsGeoMapperX::Instance();
if (option == 1) // MAPS
{
//----------------------------------------------
// MAPS 4 layers, no EndCap
//----------------------------------------------
mapper->SetGeometryPars(1, 1);
mapper->EnableEndcaps(0);
return;
}
if (option == 2) // DSSD
{
//----------------------------------------------
// DSSD 3 layers w EndCap
//----------------------------------------------
mapper->SetGeometryPars(3, 2); // see https://git.jinr.ru/nica/spdroot/-/wikis/Simulation
mapper->EnableEndcaps(1);
return;
}
}
//________________________________________________________________________________________________
void CustomTsB(Double_t deg)
{
if (deg < 0 || deg > 90)
return;
SpdTsTBGeoMapper *mapper = SpdTsTBGeoMapper::Instance();
Double_t d = 1.0; // straw diameter;
Double_t gap1 = 0.; // gap between straws in the sublayer
Double_t gap2 = 0.; // gap between the sublayers in the double layer
Double_t gap3 = 0.; // gap between the double layers in the module
Double_t w = SpdTsTBGeoMapper::GetLayerWidthFromStrawDiameter(d, gap1, gap2);
mapper->SetStrawModulePars(1, w, gap3);
mapper->SetStrawLayerPars(1, 0., d, gap1, gap2); // z layer
mapper->SetStrawLayerPars(1, +deg, d, gap1, gap2); // u layer
mapper->SetStrawLayerPars(1, -deg, d, gap1, gap2); // v layer
}
macro/examples/jpsi-ee/analyze_jpsi.C deleted 100644 → 0
View file @ 3e5006bc
/*
This is a physics analysis example using SpdMCDataIterator class.
It shows access to
-- mc truth,
-- fitted tracks,
-- EMC particles ("separated" showers in EMC),
-- track extrapolation (to match the proper MC particle).
This macro analyzes J/psi -> e+ e- events and provides a set of histograms as an output.
It macro takes the "reco_full.root" as an input. This file is produced by the default
reco_event_jpsi.C (="macro/analysis/RecoEventFull.C") macro.
*/
#include "fit_dimu.C"
void analyze_jpsi(int events_max = 1000, int lepton_pdg = 11)
{
double lepton_mass = 0.10566;
if (lepton_pdg == 11) {
lepton_mass = 0.511e-3;
}
SpdMCDataIterator* IT = 0;
const SpdSetParSet* Pars_ = 0;
//----------------------------------
// pointers to current data objects
const TClonesArray* mcparticles = 0; //! List of mc-particles
const TClonesArray* mctracks = 0; //! List of mc-tracks
const TClonesArray* rcecalpartilces = 0; //! List of ECAL rs-particles
IT = new SpdMCDataIterator();
IT->AddSourceFile("reco_full.root");
IT->ActivateBranch("MCParticles");
IT->ActivateBranch("RCEcalParticles");
IT->ActivateBranch("MCTracks");
IT->Init();
//----------------------------------------------------------
// Get data pointers
Pars_ = IT->GetParameters();
mcparticles = IT->GetParticles();
mctracks = IT->GetTracks();
rcecalpartilces = IT->GetEcalParticlesRC();
TH1D* hpmup = new TH1D("hpmup", "Momentum of e^{+}; p(e^{+}) [GeV]", 50, 0, 5);
TH1D* hMmumu = new TH1D("hMmumu", "Dimuon mass; M_{#mu#mu} [GeV]", 50, 2.6, 3.6);
TH1D* hjpsi_pT = new TH1D("hjpsi_pT", "J/#psi p_{T} spectrum; p_{T} [GeV]", 100, 0, 5);
TH1D* hMuonEcal = new TH1D("hMuonEcal", "Lepton shower in ecal (d<10 cm);E [GeV]", 100, 0, 5);
TH1D* hd = new TH1D("hd", "Distance from extr. track to the closest shower;d [cm]", 100, 0, 100);
TH2D* hCosMupCosMum = new TH2D("hCosMupCosMum", "cos#theta_{e^{+}} Vs. cos#theta_{e^{-}};cos#theta_{e^{+}};cos#theta_{e^{-}}", 50, -1, 1, 50, -1, 1);
SpdTrackPropagatorGF trkProp;
trkProp.Init();
trkProp.SetMaterialEffects(false);
Int_t ne_total(0);
while ( (IT->NextEvent()) && (ne_total < events_max) ) {
TLorentzVector pmum;
TLorentzVector pmup;
bool fitted_pmum = false;
bool fitted_pmup = false;
for (int imctrack=0; imctrack<mctracks->GetEntriesFast(); imctrack++) {
SpdTrackMC* mctrack = (SpdTrackMC*) mctracks->At(imctrack);
if (!mctrack) continue; // why??
if ( TMath::Abs( mctrack->GetParticlePdg() ) == lepton_pdg ) {
// Get the corresponding lepton and check id it
// originates from J/psi
int partid = mctrack->GetParticleId();
SpdMCParticle* part = (SpdMCParticle*) mcparticles->At( partid );
//if ( part->GetSpecialMotherId() != 5 ) continue; // may not work for all events?
// Instead
int mid = part->GetMotherId();
if ( mid < 0 ) continue;
SpdMCParticle* mother_part = (SpdMCParticle*) mcparticles->At( mid );
if (mother_part->GetPdgCode() != 443) continue;
// Extrapolate to Ecal
SpdTrackFitPar* fitpars = mctrack->GetFitPars();
if ( !fitpars ) continue;
SpdTrackState* trklaststate = fitpars->GetLastState();
if (!trklaststate) continue;
//cout << "The muon last state pos: ";
TVector3 lpos = trklaststate->GetPosition();
//lpos.Print();
// Get momentum form the final state fit parameters
SpdTrackState* trkfinalstate = fitpars->GetFinalState();
if (!trkfinalstate) continue;
TVector3 p3 = trkfinalstate->GetMomentum();
//cout << "Muon momentum: ";
//p3.Print();
int pdg = mctrack->GetParticlePdg();
if ( pdg == lepton_pdg ) {
pmum.SetXYZM(p3.X(), p3.Y(), p3.Z(), lepton_mass);
fitted_pmum = true;
}
if ( pdg == -lepton_pdg ) {
pmup.SetXYZM(p3.X(), p3.Y(), p3.Z(), lepton_mass);
fitted_pmup = true;
hpmup->Fill( p3.Mag() );
}
trkProp.InitTrack( pdg );
// Loop over SpdEcalRCParticle, find the closest
double dmin = 1e10;
SpdEcalRCParticle* closest_emc_particle = nullptr;
for (int ishower=0; ishower<rcecalpartilces->GetEntriesFast(); ishower++) {
SpdEcalRCParticle* ecal_part = (SpdEcalRCParticle*) rcecalpartilces->At(ishower);
if (!ecal_part) continue;
//cout << " ---------------- EMC particle! -----------------" << endl;
TVector3 pos = ecal_part->GetPosition();
//cout << "EMC energy " << ecal_part->GetEnergy() << endl;
if (ecal_part->GetEnergy() < 0.1) {
//cout << "Skipping shower" << endl;
continue;
}
//cout << "The shower pos: ";
//pos.Print();
SpdTrackState trkextstate;
//cout << "Extrapolation distance: ";
trkProp.ExtrapolateToPoint( pos, *trklaststate, trkextstate);
TVector3 extpos = trkextstate.GetPosition();
//cout << "The muon extrapolated pos: ";
//extpos.Print();
double d = (pos - extpos).Mag();
if ( d < dmin ) {
dmin = d;
closest_emc_particle = ecal_part;
cout << " *********** d = " << d << " ************* " << endl;
cout << "Charged particle momentum: ";
p3.Print();
cout << " Shower postion: ";
pos.Print();
cout << " Track extr. postion (poca): ";
extpos.Print();
}
} // ishower
if ( closest_emc_particle ) {
//cout << "dmin = " << dmin << endl;
hd->Fill(dmin);
if ( dmin < 10 ) {
hMuonEcal->Fill( closest_emc_particle->GetEnergy() );
}
}
}
} // MCtracks
if (fitted_pmum && fitted_pmup) {
TLorentzVector pjpsi = pmup + pmum;
hMmumu->Fill( pjpsi.Mag() );
hjpsi_pT->Fill( pjpsi.Pt() );
hCosMupCosMum->Fill( pmup.CosTheta(), pmum.CosTheta() );
}
ne_total++;
}
cout << "Events processed: " << ne_total << endl;
gStyle->SetOptStat(11);
TCanvas* c = new TCanvas("c", "", 1200, 800);
c->Divide(3,2);
c->cd(1);
hpmup->Draw();
c->cd(2);
hMmumu->Draw();
FitDimuSpectrum( hMmumu );
c->cd(3);
hjpsi_pT->Draw();
c->cd(4);
hCosMupCosMum->Draw("colz");
c->cd(5);
hd->Draw();
c->cd(6);
hMuonEcal->Draw();
c->SaveAs("jpsi.pdf");
}