CbmRoot
CbmLitFieldQa.cxx
Go to the documentation of this file.
1 
6 #include "CbmLitFieldQa.h"
7 #include "CbmHistManager.h"
8 #include "CbmLitFieldQaReport.h"
9 #include "CbmRichDetectorData.h" // for CbmRichPmtData, CbmRichPixelData
10 #include "CbmRichDigiMapManager.h"
11 #include "CbmRichGeoManager.h"
12 #include "CbmUtils.h"
13 
14 #include "FairField.h"
15 #include "FairRunAna.h"
16 #include "FairRuntimeDb.h"
17 
18 #include "TGraph.h"
19 #include "TGraph2D.h"
20 #include "TVector3.h"
21 #include <TFile.h>
22 
23 #include <boost/assign/list_of.hpp>
24 #include <cmath>
25 #include <sstream>
26 
27 using boost::assign::list_of;
28 using Cbm::ToString;
29 using namespace std;
30 
32  : fField(NULL)
33  , fNofSlices(0)
34  , fZSlicePosition()
35  , fXSlicePosition()
36  , fYSlicePosition()
37  , fAlongZAngles()
38  , fAlongZXY()
39  , fZMin(-10.)
40  , fZMax(300.)
41  , fZStep(1.)
42  ,
43  //fAcceptanceAngleX(35.),
44  //fAcceptanceAngleY(35.),
45  fNofBinsX(100)
46  , fNofBinsY(100)
47  , fMinZFieldIntegral(171.)
48  , fMaxZFieldIntegral(330.)
49  , fHM(NULL)
50  , fOutputDir("./") {}
51 
53 
54 InitStatus CbmLitFieldQa::Init() {
55  fNofSlices = fZSlicePosition.size();
56 
57  // Calculate (X, Y) window for each slice
60  for (Int_t i = 0; i < fNofSlices; i++) {
61  // Double_t tanXangle = tan(fAcceptanceAngleX * 3.14159265 / 180); //
62  // Double_t tanYangle = tan(fAcceptanceAngleY * 3.14159265 / 180); //
63  // fXSlicePosition[i] = fZSlicePosition[i] * tanXangle;
64  // fYSlicePosition[i] = fZSlicePosition[i] * tanYangle;
65  fXSlicePosition[i] = 250.;
66  fYSlicePosition[i] = 250.;
67  }
68 
69  vector<Double_t> tmp = list_of(0.)(10.)(20.)(30.);
70  fAlongZAngles = tmp;
71 
72  fAlongZXY.push_back(make_pair(0., 0.));
73  fAlongZXY.push_back(make_pair(100., 0.));
74  fAlongZXY.push_back(make_pair(0., 100.));
75 
76 
77  fField = FairRunAna::Instance()->GetField();
78 
79  fHM = new CbmHistManager();
80 
81  CreateHistos();
82  FillBHistos();
83 
85 
87  report->Create(fHM, fOutputDir);
88  delete report;
89  TDirectory* oldir = gDirectory;
90  TFile* outFile = FairRootManager::Instance()->GetOutFile();
91  if (outFile != NULL) {
92  outFile->cd();
93  fHM->WriteToFile();
94  }
95  gDirectory->cd(oldir->GetPath());
96 
97  // print some field values at RICH entrance
98  Double_t B[3];
99  vector<vector<Double_t>> pos = {
100  {0., 0., 170.}, {0., 80., 170.}, {50., 0., 170.}, {0., 0., 250.}};
101  for (Int_t i = 0; i < pos.size(); i++) {
102  fField->GetFieldValue(&pos[i][0], B);
103  Double_t magB = sqrt(B[0] * B[0] + B[1] * B[1] + B[2] * B[2]);
104  cout << "B at (" << pos[i][0] << ", " << pos[i][1] << ", " << pos[i][2]
105  << ") = " << magB << " kGauss" << endl;
106  }
107 
108  cout << "CbmLitFieldQa::Init() end" << endl;
109 
110  return kSUCCESS;
111 }
112 
113 void CbmLitFieldQa::Exec(Option_t* opt) {}
114 
116 
118  string names[] = {"Bx", "By", "Bz", "Mod"};
119  string zTitle[] = {
120  "B_{x} [kGauss]", "B_{y} [kGauss]", "B_{z} [kGauss]", "|B| [kGauss]"};
121  for (Int_t v = 0; v < 4; v++) {
122  for (Int_t i = 0; i < fNofSlices; i++) {
123  TGraph2D* graph = new TGraph2D();
124  graph->SetNpx(200);
125  graph->SetNpy(200);
126  string name =
127  "hmf_" + names[v] + "_Graph2D_" + ToString<Int_t>(fZSlicePosition[i]);
128  string title = name + ";X [cm];Y [cm];" + zTitle[v];
129  graph->SetNameTitle(name.c_str(), title.c_str());
130  fHM->Add(name, graph);
131  }
132  }
133 
134  for (Int_t v = 0; v < 4; v++) {
135  for (Int_t i = 0; i < fAlongZAngles.size(); i++) {
136  TGraph* graph = new TGraph();
137  string name = "hmf_" + names[v] + "AlongZAngle_Graph_"
138  + ToString<Int_t>(fAlongZAngles[i]);
139  string title = name + ";Z [cm];B [kGauss]";
140  graph->SetNameTitle(name.c_str(), title.c_str());
141  fHM->Add(name, graph);
142  }
143  for (Int_t i = 0; i < fAlongZXY.size(); i++) {
144  TGraph* graph = new TGraph();
145  string name = "hmf_" + names[v] + "AlongZXY_Graph_"
146  + ToString<Int_t>(fAlongZXY[i].first) + "_"
147  + ToString<Int_t>(fAlongZXY[i].second);
148  string title = name + ";Z [cm];B [kGauss]";
149  graph->SetNameTitle(name.c_str(), title.c_str());
150  fHM->Add(name, graph);
151  }
152  for (Int_t i = 0; i < fAlongZXY.size(); i++) {
153  TGraph* graph = new TGraph();
154  string name = "hmf_" + names[v] + "AlongZXYIntegral_Graph_"
155  + ToString<Int_t>(fAlongZXY[i].first) + "_"
156  + ToString<Int_t>(fAlongZXY[i].second);
157  string title = name + ";Z [cm];B_{Int_t} [kGauss*m]";
158  graph->SetNameTitle(name.c_str(), title.c_str());
159  fHM->Add(name, graph);
160  }
161  }
162 
163  for (Int_t p = 0; p < 2; p++) {
164  for (Int_t v = 0; v < 4; v++) {
165  TGraph2D* graph = new TGraph2D();
166  graph->SetNpx(200);
167  graph->SetNpy(200);
168  string name =
169  "hmf_" + names[v] + "_Rich_Pmt_" + ((p == 0) ? "up" : "down");
170  string title = name + ";X [cm];Y [cm];" + zTitle[v];
171  graph->SetNameTitle(name.c_str(), title.c_str());
172  fHM->Add(name, graph);
173  }
174  }
175 }
176 
178  string names[] = {"Bx", "By", "Bz", "Mod"};
179  // Fill graphs for magnetic field for each (X, Y) slice
180  for (UInt_t iSlice = 0; iSlice < fNofSlices; iSlice++) { // loop over slices
181  Double_t Z = fZSlicePosition[iSlice];
182 
183  Int_t cnt = 0;
184  Double_t HX =
185  2 * fXSlicePosition[iSlice] / fNofBinsX; // step size for X position
186  Double_t HY =
187  2 * fYSlicePosition[iSlice] / fNofBinsY; // step size for Y position
188  for (Int_t iX = 0; iX < fNofBinsX; iX++) { // loop over x position
189  Double_t X = -fXSlicePosition[iSlice] + (iX + 0.5) * HX;
190  for (Int_t iY = 0; iY < fNofBinsY; iY++) { // loop over y position
191  Double_t Y = -fYSlicePosition[iSlice] + (iY + 0.5) * HY;
192 
193  // Get field value
194  Double_t pos[3] = {X, Y, Z};
195  Double_t B[4];
196  fField->GetFieldValue(pos, B);
197 
198  B[3] = sqrt(B[0] * B[0] + B[1] * B[1] + B[2] * B[2]);
199  for (Int_t v = 0; v < 4; v++) {
200  string name = "hmf_" + names[v] + "_Graph2D_"
201  + ToString<Int_t>(fZSlicePosition[iSlice]);
202  fHM->G2(name)->SetPoint(cnt, X, Y, B[v]);
203  }
204  cnt++;
205  }
206  }
207  }
208 
209  // Fill histograms for magnetic field along Z for different angles
210  for (Int_t i = 0; i < fAlongZAngles.size(); i++) {
211  Int_t nofSteps = Int_t((fZMax - fZMin) / fZStep);
212  for (Int_t istep = 0; istep < nofSteps; istep++) {
213  Double_t Z = fZMin + istep * fZStep;
214  Double_t tanXangle = tan(fAlongZAngles[i] * 3.14159265 / 180); //
215  Double_t tanYangle = tan(fAlongZAngles[i] * 3.14159265 / 180); //
216  Double_t X = Z * tanXangle;
217  Double_t Y = Z * tanYangle;
218 
219  // Get field value
220  Double_t pos[3] = {X, Y, Z};
221  Double_t B[4];
222  fField->GetFieldValue(pos, B);
223 
224  B[3] = sqrt(B[0] * B[0] + B[1] * B[1] + B[2] * B[2]);
225  for (Int_t v = 0; v < 4; v++) {
226  string name = "hmf_" + names[v] + "AlongZAngle_Graph_"
227  + ToString<Int_t>(fAlongZAngles[i]);
228  fHM->G1(name)->SetPoint(istep, Z, B[v]);
229  }
230  }
231  }
232  // Fill histograms for magnetic field along Z for different X position
233  for (Int_t i = 0; i < fAlongZXY.size(); i++) {
234  Int_t nofSteps = Int_t((fZMax - fZMin) / fZStep);
235  Double_t integralB[4] = {0., 0., 0., 0.};
236  for (Int_t istep = 0; istep < nofSteps; istep++) {
237  Double_t Z = fZMin + istep * fZStep;
238  Double_t X = fAlongZXY[i].first;
239  Double_t Y = fAlongZXY[i].second;
240 
241  // Get field value
242  Double_t pos[3] = {X, Y, Z};
243  Double_t B[4];
244  fField->GetFieldValue(pos, B);
245 
246  B[3] = sqrt(B[0] * B[0] + B[1] * B[1] + B[2] * B[2]);
247 
248  for (Int_t v = 0; v < 4; v++) {
249  string name = "hmf_" + names[v] + "AlongZXY_Graph_"
250  + ToString<Int_t>(fAlongZXY[i].first) + "_"
251  + ToString<Int_t>(fAlongZXY[i].second);
252  fHM->G1(name)->SetPoint(istep, Z, B[v]);
253  // Calculate field integral in the RICH detector
254  if (Z >= fMinZFieldIntegral && Z <= fMaxZFieldIntegral) {
255  integralB[v] += 0.01 * fZStep * fabs(B[v]); // in kGauss * meter
256  string name = "hmf_" + names[v] + "AlongZXYIntegral_Graph_"
257  + ToString<Int_t>(fAlongZXY[i].first) + "_"
258  + ToString<Int_t>(fAlongZXY[i].second);
259  fHM->G1(name)->SetPoint(istep, Z, integralB[v]);
260  fHM->G1(name)->SetMaximum(1.1 * integralB[v]);
261  }
262  }
263  }
264  }
265 }
266 
268  string names[] = {"Bx", "By", "Bz", "Mod"};
269  vector<Int_t> pixels =
271  if (pixels.size() == 0) return;
272 
273  Double_t maxModB = 0.;
274  int iS = -1;
275  for (Int_t i = 0; i < pixels.size(); i++) {
276  TVector3 inPos;
277  CbmRichPixelData* pixelData =
279  inPos.SetXYZ(pixelData->fX, pixelData->fY, pixelData->fZ);
280  string ud = "up";
281  if (inPos.Y() < 0) ud = "down";
282  TVector3 outPos;
283  CbmRichGeoManager::GetInstance().RotatePoint(&inPos, &outPos);
284 
285  // Get field value
286  // Take B-field values from real PMT position
287  Double_t pos[3] = {inPos.X(), inPos.Y(), inPos.Z()};
288 
289  Double_t B[4];
290  fField->GetFieldValue(pos, B);
291 
292  B[3] = sqrt(B[0] * B[0] + B[1] * B[1] + B[2] * B[2]);
293  if (maxModB < B[3]) maxModB = B[3];
294  for (Int_t v = 0; v < 4; v++) {
295  string name = "hmf_" + names[v] + "_Rich_Pmt_" + ud;
296  // Take B-field values from real PMT position, but display them for rotated X and Y
297  fHM->G2(name)->SetPoint(
298  fHM->G2(name)->GetN(), outPos.X(), outPos.Y(), B[v]);
299  }
300  }
301 
302  cout << "Maximum Bmod onto RICH PMT:" << maxModB << " kGauss" << endl;
303 }
304 
CbmLitFieldQa::fField
FairField * fField
Definition: CbmLitFieldQa.h:88
CbmLitFieldQa::~CbmLitFieldQa
virtual ~CbmLitFieldQa()
Destructor.
Definition: CbmLitFieldQa.cxx:52
CbmRichPixelData
Definition: CbmRichDetectorData.h:17
CbmLitFieldQa::fZMin
Double_t fZMin
Definition: CbmLitFieldQa.h:110
CbmLitFieldQa::fNofBinsX
Int_t fNofBinsX
Definition: CbmLitFieldQa.h:102
CbmRichDigiMapManager::GetPixelAddresses
std::vector< Int_t > GetPixelAddresses()
Definition: CbmRichDigiMapManager.cxx:280
CbmLitFieldQa::fAlongZAngles
vector< Double_t > fAlongZAngles
Definition: CbmLitFieldQa.h:107
CbmLitFieldQa::fNofBinsY
Int_t fNofBinsY
Definition: CbmLitFieldQa.h:103
CbmHistManager::WriteToFile
void WriteToFile()
Write all histograms to current opened file.
Definition: core/base/CbmHistManager.cxx:103
CbmRichGeoManager::GetInstance
static CbmRichGeoManager & GetInstance()
Definition: CbmRichGeoManager.h:29
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
CbmLitFieldQa::fXSlicePosition
vector< Double_t > fXSlicePosition
Definition: CbmLitFieldQa.h:96
CbmLitFieldQa::Init
virtual InitStatus Init()
Inherited from FairTask.
Definition: CbmLitFieldQa.cxx:54
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
CbmRichDigiMapManager.h
CbmLitFieldQa::FillBHistos
void FillBHistos()
Fill graphs and histos for field map for each field component (Bx, By, Bz, |B|).
Definition: CbmLitFieldQa.cxx:177
CbmLitFieldQa
Field map QA.
Definition: CbmLitFieldQa.h:29
CbmLitFieldQa::fNofSlices
Int_t fNofSlices
Definition: CbmLitFieldQa.h:99
CbmHistManager::G1
TGraph * G1(const std::string &name) const
Return pointer to TGraph.
Definition: CbmHistManager.h:223
CbmLitFieldQa::FillRichPmtPlaneBHistos
void FillRichPmtPlaneBHistos()
Fill B-field histograms for RICH PMT plane.
Definition: CbmLitFieldQa.cxx:267
CbmLitFieldQa.h
Field map QA.
CbmHistManager.h
Histogram manager.
CbmLitFieldQa::fMaxZFieldIntegral
Double_t fMaxZFieldIntegral
Definition: CbmLitFieldQa.h:105
CbmLitFieldQa::fYSlicePosition
vector< Double_t > fYSlicePosition
Definition: CbmLitFieldQa.h:97
CbmLitFieldQa::fZStep
Double_t fZStep
Definition: CbmLitFieldQa.h:112
CbmLitFieldQa::fHM
CbmHistManager * fHM
Definition: CbmLitFieldQa.h:114
CbmRichGeoManager.h
CbmLitFieldQa::CbmLitFieldQa
CbmLitFieldQa()
Constructor.
Definition: CbmLitFieldQa.cxx:31
CbmRichDetectorData.h
CbmHistManager
Histogram manager.
Definition: CbmHistManager.h:41
CbmLitFieldQa::fZSlicePosition
vector< Double_t > fZSlicePosition
Definition: CbmLitFieldQa.h:95
CbmRichPixelData::fY
Double_t fY
Definition: CbmRichDetectorData.h:21
CbmLitFieldQa::fZMax
Double_t fZMax
Definition: CbmLitFieldQa.h:111
CbmHistManager::G2
TGraph2D * G2(const std::string &name) const
Return pointer to TGraph2D.
Definition: CbmHistManager.h:243
CbmLitFieldQa::Finish
virtual void Finish()
Inherited from FairTask.
Definition: CbmLitFieldQa.cxx:115
CbmSimulationReport::Create
void Create(CbmHistManager *histManager, const std::string &outputDir)
Main function which creates report data.
CbmRichPixelData::fX
Double_t fX
Definition: CbmRichDetectorData.h:20
CbmRichGeoManager::RotatePoint
void RotatePoint(TVector3 *inPos, TVector3 *outPos, Bool_t noTilting=false)
Definition: CbmRichGeoManager.cxx:407
ClassImp
ClassImp(CbmLitFieldQa)
first
bool first
Definition: LKFMinuit.cxx:143
CbmUtils.h
CbmLitFieldQa::CreateHistos
void CreateHistos()
Create histograms.
Definition: CbmLitFieldQa.cxx:117
CbmLitFieldQaReport.h
Creates field QA report.
CbmLitFieldQa::fAlongZXY
vector< std::pair< Double_t, Double_t > > fAlongZXY
Definition: CbmLitFieldQa.h:109
v
__m128 v
Definition: L1/vectors/P4_F32vec4.h:1
CbmLitFieldQa::fMinZFieldIntegral
Double_t fMinZFieldIntegral
Definition: CbmLitFieldQa.h:104
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:60
pos
TVector3 pos
Definition: CbmMvdSensorDigiToHitTask.cxx:60
CbmSimulationReport
Base class for simulation reports.
Definition: CbmSimulationReport.h:28
CbmRichDigiMapManager::GetInstance
static CbmRichDigiMapManager & GetInstance()
Definition: CbmRichDigiMapManager.h:29
CbmLitFieldQa::fOutputDir
string fOutputDir
Definition: CbmLitFieldQa.h:117
CbmLitFieldQaReport
Creates field QA report.
Definition: CbmLitFieldQaReport.h:20
Cbm::ToString
std::string ToString(const T &value)
Definition: CbmUtils.h:16
CbmLitFieldQa::Exec
virtual void Exec(Option_t *opt)
Inherited from FairTask.
Definition: CbmLitFieldQa.cxx:113
CbmRichDigiMapManager::GetPixelDataByAddress
CbmRichPixelData * GetPixelDataByAddress(Int_t address)
Definition: CbmRichDigiMapManager.cxx:267
CbmRichPixelData::fZ
Double_t fZ
Definition: CbmRichDetectorData.h:22
CbmHistManager::Add
void Add(const std::string &name, TNamed *object)
Add new named object to manager.
Definition: CbmHistManager.h:58