FairRoot/PandaRoot
PndGemMagneticFieldVsTrackParameters.cxx
Go to the documentation of this file.
1 //* $Id: */
2 
3 // -------------------------------------------------------------------------
4 // ----- PndGemMagneticFieldVsTrackParameters source file -----
5 // ----- Created 12/02/2009 by R. Karabowicz -----
6 // -------------------------------------------------------------------------
7 
8 // Includes from GEM
9 #include "PndGemMagneticFieldVsTrackParameters.h"
10 
11 // Includes from base
12 #include "FairRootManager.h"
13 #include "FairRunAna.h"
14 #include "FairRuntimeDb.h"
15 
16 // Includes from ROOT
17 #include "TClonesArray.h"
18 #include "TObjArray.h"
19 #include "TMath.h"
20 #include "TGeoManager.h"
21 #include "TGeoNode.h"
22 
23 #include "PndGemMCPoint.h"
24 #include "PndMCTrack.h"
25 #include "PndGemDigiPar.h"
26 
27 #include <iostream>
28 #include <iomanip>
29 #include <map>
30 
31 using std::cout;
32 using std::cerr;
33 using std::endl;
34 using std::pair;
35 using std::setw;
36 using std::left;
37 using std::right;
38 using std::fixed;
39 using std::setprecision;
40 using std::map;
41 
42 
43 
44 // ----- Default constructor ------------------------------------------
45 PndGemMagneticFieldVsTrackParameters::PndGemMagneticFieldVsTrackParameters() : FairTask("GEM MagneticFieldVsTrackParametersr", 1) {
46  fDigiPar = NULL;
47  fMCPointArray = NULL;
48  fMCTrackArray = NULL;
49  Reset();
50 }
51 // -------------------------------------------------------------------------
52 
53 
54 
55 // ----- Standard constructor ------------------------------------------
56 PndGemMagneticFieldVsTrackParameters::PndGemMagneticFieldVsTrackParameters(Int_t iVerbose)
57  : FairTask("GEM MagneticFieldVsTrackParametersr", iVerbose) {
58  fDigiPar = NULL;
59  fMCPointArray = NULL;
60  fMCTrackArray = NULL;
61  Reset();
62 }
63 // -------------------------------------------------------------------------
64 
65 
66 
67 // ----- Constructor with name -----------------------------------------
68 PndGemMagneticFieldVsTrackParameters::PndGemMagneticFieldVsTrackParameters(const char* name, Int_t iVerbose)
69  : FairTask(name, iVerbose) {
70  fDigiPar = NULL;
71  fMCPointArray = NULL;
72  fMCTrackArray = NULL;
73  Reset();
74 }
75 // -------------------------------------------------------------------------
76 
77 
78 
79 // ----- Destructor ----------------------------------------------------
80 PndGemMagneticFieldVsTrackParameters::~PndGemMagneticFieldVsTrackParameters() {
81  if ( fDigiPar) delete fDigiPar;
82  Reset();
83 }
84 // -------------------------------------------------------------------------
85 
86 // ----- Public method Exec --------------------------------------------
87 void PndGemMagneticFieldVsTrackParameters::Exec(Option_t*) {
88 
89  // Size of fMCTrackArray
90  //cout <<"# MC Tracks: "<< fMCTrackArray->GetEntriesFast() << endl;
91  // Size of fMCTrackArray
92  //cout <<"# MC Points: "<< fMCPointArray->GetEntriesFast() << endl;
93 
94  Fill1StationHistograms();
95  Fill2StationsHistograms();
96 
97  Reset();
98 
99 }
100 // -------------------------------------------------------------------------
101 
102 
103 
104 // ----- Private method SetParContainers -------------------------------
105 void PndGemMagneticFieldVsTrackParameters::SetParContainers() {
106 
107  // Get run and runtime database
108  FairRunAna* run = FairRunAna::Instance();
109  if ( ! run ) Fatal("SetParContainers", "No analysis run");
110 
111  FairRuntimeDb* db = run->GetRuntimeDb();
112  if ( ! db ) Fatal("SetParContainers", "No runtime database");
113 
114  // Get GEM digitisation parameter container
115  fDigiPar = (PndGemDigiPar*)(db->getContainer("PndGemDetectors"));
116 
117 }
118 // -------------------------------------------------------------------------
119 
120 
121 
122 // ----- Private method Init -------------------------------------------
123 InitStatus PndGemMagneticFieldVsTrackParameters::Init() {
124 
125  // Get input array
126  FairRootManager* ioman = FairRootManager::Instance();
127  if ( ! ioman ) Fatal("Init", "No FairRootManager");
128 
129  // Get MCPoint array
130  fMCPointArray = (TClonesArray*) ioman->GetObject("GEMPoint");
131  if( !fMCPointArray ) {
132  cout << "-E- "<< GetName() <<"::Init: No MCPoint array!"
133  << endl;
134  return kERROR;
135  }
136 
137  // Get MCTrack array
138  fMCTrackArray = (TClonesArray*) ioman->GetObject("MCTrack");
139  if( !fMCTrackArray ) {
140  cout << "-E- "<< GetName() <<"::Init: No MCTrack array!"
141  << endl;
142  return kERROR;
143  }
144 
145  fTrackFinderOnHits_ParThetaA = fDigiPar->GetTrackFinderOnHits_ParThetaA();
146  fTrackFinderOnHits_ParThetaB = fDigiPar->GetTrackFinderOnHits_ParThetaB();
147 
148  fTrackFinderOnHits_ParTheta0 = fDigiPar->GetTrackFinderOnHits_ParTheta0();
149  fTrackFinderOnHits_ParTheta1 = fDigiPar->GetTrackFinderOnHits_ParTheta1();
150  fTrackFinderOnHits_ParTheta2 = fDigiPar->GetTrackFinderOnHits_ParTheta2();
151  fTrackFinderOnHits_ParTheta3 = fDigiPar->GetTrackFinderOnHits_ParTheta3();
152 
153  fTrackFinderOnHits_ParRadPhi0 = fDigiPar->GetTrackFinderOnHits_ParRadPhi0();
154  fTrackFinderOnHits_ParRadPhi2 = fDigiPar->GetTrackFinderOnHits_ParRadPhi2();
155  for ( Int_t in = 0 ; in < 3 ; in++ ) {
156  fTrackFinderOnHits_ParMat0[in] = fDigiPar->GetTrackFinderOnHits_ParMat0(in);
157  fTrackFinderOnHits_ParMat1[in] = fDigiPar->GetTrackFinderOnHits_ParMat1(in);
158  }
159 
160  CreateHistos();
161 
162  cout << "-I- " << fName.Data() << "::Init(). There are " << fDigiPar->GetNStations() << " GEM stations." << endl;
163  cout << "-I- " << fName.Data() << "::Init(). Initialization succesfull." << endl;
164 
165  return kSUCCESS;
166 
167 }
168 // -------------------------------------------------------------------------
169 
170 
171 
172 // ----- Private method ReInit -----------------------------------------
173 InitStatus PndGemMagneticFieldVsTrackParameters::ReInit() {
174  return kSUCCESS;
175 }
176 // -------------------------------------------------------------------------
177 
178 
179 
180 // ----- Private method Reset ------------------------------------------
181 void PndGemMagneticFieldVsTrackParameters::Reset() {
182  fNofEvents = 0;
183 }
184 // -------------------------------------------------------------------------
185 
186 
187 // ---- Private function to fill histograms dependent on 1 station --------------------
188 Int_t PndGemMagneticFieldVsTrackParameters::Fill1StationHistograms() {
189  Int_t nofGemPoints = fMCPointArray->GetEntriesFast();
190  PndGemMCPoint* mcPoint;
191  PndMCTrack* mcTrack;
192 
193  for ( Int_t imcp = 0 ; imcp < nofGemPoints ; imcp++ ) {
194  mcPoint = (PndGemMCPoint*)fMCPointArray->At(imcp);
195 
196  Int_t stationNumber = fDigiPar->GetStationNr(mcPoint->GetSensorId());
197 
198  mcTrack = (PndMCTrack*)fMCTrackArray->At(mcPoint->GetTrackID());
199  TVector3 mcMomVec = mcTrack->GetMomentum();
200 
201  Double_t pRadius = TMath::Sqrt(mcPoint->GetX()*mcPoint->GetX()+mcPoint->GetY()*mcPoint->GetY());
202  Double_t pZ = mcPoint->GetZ();
203  Double_t tMom = mcMomVec.Mag();
204  Double_t tTheta = mcMomVec.Theta()*TMath::RadToDeg();
205 
206  fhThetaVsRadiusVsMomentum[stationNumber-1]->Fill(tMom,pRadius, tTheta);
207  fhThetaVsRadiusVsMomentumAll ->Fill(tMom,pRadius/pZ,tTheta);
208 
209  }
210  return 0;
211 }
212 // -------------------------------------------------------------------------
213 
214 // ---- Private function to fill histograms dependent on 2 stations --------------------
215 Int_t PndGemMagneticFieldVsTrackParameters::Fill2StationsHistograms() {
216  Int_t nofGemPoints = fMCPointArray->GetEntriesFast();
217  PndGemMCPoint* mcPoint1;
218  PndGemMCPoint* mcPoint2;
219  PndMCTrack* mcTrack;
220 
221  for ( Int_t imcp1 = 0 ; imcp1 < nofGemPoints ; imcp1++ ) {
222  mcPoint1 = (PndGemMCPoint*)fMCPointArray->At(imcp1);
223 
224  Double_t p1Radius = TMath::Sqrt(mcPoint1->GetX()*mcPoint1->GetX()+mcPoint1->GetY()*mcPoint1->GetY());
225  Double_t p1Z = mcPoint1->GetZ();
226  Double_t p1PhiAng = TMath::ACos(mcPoint1->GetX()/p1Radius);
227  if ( mcPoint1->GetY() < 0 )
228  p1PhiAng = 2.*TMath::Pi() - p1PhiAng;
229  p1PhiAng *= TMath::RadToDeg();
230 
231  Int_t mcTrackId = mcPoint1->GetTrackID();
232 
233  mcTrack = (PndMCTrack*)fMCTrackArray->At(mcTrackId);
234 
235  if ( mcTrack->GetMotherID() != -1 ) continue;
236 
237  TVector3 mcMomVec = mcTrack->GetMomentum();
238  Double_t tMom = mcMomVec.Mag();
239  Double_t tPt = mcMomVec.Pt();
240  Double_t pMom = TMath::RadToDeg()*mcMomVec.Phi() + (mcMomVec.Phi()>=0?0.:360.);
241 
242  for ( Int_t imcp2 = imcp1+1 ; imcp2 < nofGemPoints ; imcp2++ ) {
243  mcPoint2 = (PndGemMCPoint*)fMCPointArray->At(imcp2);
244  Double_t p2Z = mcPoint2->GetZ();
245  if ( TMath::Abs(p2Z-p1Z) < 5. ) continue;
246 
247  Double_t p2Radius = TMath::Sqrt(mcPoint2->GetX()*mcPoint2->GetX()+mcPoint2->GetY()*mcPoint2->GetY());
248 
249  Double_t p2PhiAng = TMath::ACos(mcPoint2->GetX()/p2Radius);
250  if ( mcPoint2->GetY() < 0 )
251  p2PhiAng = 2.*TMath::Pi() - p2PhiAng;
252  p2PhiAng *= TMath::RadToDeg();
253 
254  if ( p1PhiAng < 90. && p2PhiAng > 270. ) p2PhiAng -= 360.;
255  if ( p1PhiAng > 270. && p2PhiAng < 90. ) p2PhiAng += 360.;
256 
257  fhRadiusVsAngle->Fill((p2PhiAng-p1PhiAng)*p1Z/p2Z,p2Radius*p1Z/(p1Radius*p2Z));
258 
259  Int_t histNo = -1;
260  Int_t ihist = 0;
261  for ( Int_t istat1 = 0 ; istat1 < fDigiPar->GetNStations() ; istat1++ ) {
262  PndGemStation* stat1 = (PndGemStation*)fDigiPar->GetStation(istat1);
263  Double_t zStation1 = stat1->GetZ();
264  for ( Int_t istat2 = istat1+1 ; istat2 < fDigiPar->GetNStations() ; istat2++ ) {
265  PndGemStation* stat2 = (PndGemStation*)fDigiPar->GetStation(istat2);
266  Double_t zStation2 = stat2->GetZ();
267  if ( TMath::Abs(p1Z-zStation1) < 2. && TMath::Abs(p2Z-zStation2) < 2. ) histNo = ihist;
268  ihist++;
269  }
270  }
271  if ( histNo == -1 ) continue;
272  fhRadiusVsAnglePair[histNo]->Fill((p2PhiAng-p1PhiAng)*p1Z/p2Z,p2Radius*p1Z/(p1Radius*p2Z));
273 
274  // Double_t circRad = TMath::Sqrt( (gemHit2X-gemHit1X)*(gemHit2X-gemHit1X) + (gemHit2Y-gemHit1Y)*(gemHit2Y-gemHit1Y) ) / TMath::Sin(pangle-pangle2) / 2.;
275  fhMomTransVsHRadius[histNo]->Fill(TMath::Abs(TMath::Sqrt((mcPoint2->GetX()-mcPoint1->GetX())*(mcPoint2->GetX()-mcPoint1->GetX())+
276  (mcPoint2->GetY()-mcPoint1->GetY())*(mcPoint2->GetY()-mcPoint1->GetY())) /
277  TMath::Sin(TMath::DegToRad()*(p2PhiAng-p1PhiAng)) / 2.),
278  tPt);
279 
280  if ( mcTrackId != mcPoint2->GetTrackID() )
281  continue;
282 
283  fhMomentumVsPhiDiffVsRadius[histNo]->Fill(p1Radius,TMath::Abs(p2PhiAng-p1PhiAng),tMom);
284 
285  fhMomentumVsPhiDiffAll[histNo] ->Fill(TMath::Abs(p2PhiAng-p1PhiAng),tMom);
286  fhMomentumVsPhiDiff [histNo][Int_t(p1Radius/5.)]->Fill(TMath::Abs(p2PhiAng-p1PhiAng),tMom);
287 
288  fhTrackPhiVsHitPhis [histNo] ->Fill(p1PhiAng,(p1PhiAng-p2PhiAng)*p1Z/(p2Z-p1Z),pMom);
289  fhTrackPhiVsCalcPhi [histNo] ->Fill(p1PhiAng+(p1PhiAng-p2PhiAng)*p1Z/(p2Z-p1Z),pMom);
290 
291  }
292  }
293  return 0;
294 }
295 // -------------------------------------------------------------------------
296 
297 // ----- Private method CreateHistos --------------------------------------------
298 void PndGemMagneticFieldVsTrackParameters::CreateHistos() {
299  fHistoList = new TList();
300  // number of mc tracks, reco tracks, efficiency as function of MOMENTUM
301 
302  const Int_t nofMomBins = 19;
303  Double_t momBins[nofMomBins+1] = {0.05,0.15,0.25,0.35,0.45,
304  0.55,0.65,0.75,0.85,0.95,
305  1.05,1.15,1.25,1.75,2.25,
306  2.75,3.25,4.75,7.25,10.75};
307  const Int_t nofRadiusBins = 60;
308  Double_t minRadius = 0.;
309  Double_t maxRadius = 60.;
310  Double_t minRRadius = 0.; // radius/zdist
311  Double_t maxRRadius = 0.4; // radius/zdist
312  Double_t radBins[nofRadiusBins+1];
313  for ( Int_t itemp = 0 ; itemp < nofRadiusBins+1 ; itemp++ ) radBins[itemp] = minRadius + (Double_t)itemp*(maxRadius-minRadius)/(Double_t(nofRadiusBins));
314  Double_t rRadBins[nofRadiusBins+1];
315  for ( Int_t itemp = 0 ; itemp < nofRadiusBins+1 ; itemp++ ) rRadBins[itemp] = minRRadius + (Double_t)itemp*(maxRRadius-minRRadius)/(Double_t(nofRadiusBins));
316  const Int_t nofThetaBins = 60;
317  Double_t minTheta = 0.;
318  Double_t maxTheta = 30.;
319  Double_t theBins[nofThetaBins+1];
320  for ( Int_t itemp = 0 ; itemp < nofThetaBins+1 ; itemp++ ) theBins[itemp] = minTheta + (Double_t)itemp*(maxTheta-minTheta)/(Double_t(nofThetaBins));
321  const Int_t nofPhiBins = 160;
322  Double_t minPhi = 0.;
323  Double_t maxPhi = 40.;
324  Double_t phiBins[nofPhiBins+1];
325  for ( Int_t itemp = 0 ; itemp < nofPhiBins+1 ; itemp++ ) phiBins[itemp] = minPhi + (Double_t)itemp*(maxPhi-minPhi)/(Double_t(nofPhiBins));
326  const Int_t nofDetPhiBins = 400;
327  const Int_t nofDetMomBins = 100;
328  Double_t minMom = 0.;
329  Double_t maxMom = 10.;
330  const Int_t nofPhi1Bins = 360;
331  Double_t phi1Bins[nofPhi1Bins+1];
332  for ( Int_t itemp = 0 ; itemp < nofPhi1Bins+1 ; itemp++ ) phi1Bins[itemp] = itemp;
333  const Int_t nofPhi2Bins = 120;
334  Double_t minPhi2 = -60.;
335  Double_t maxPhi2 = 60.;
336  Double_t phi2Bins[nofPhi2Bins+1];
337  for ( Int_t itemp = 0 ; itemp < nofPhi2Bins+1 ; itemp++ ) phi2Bins[itemp] = minPhi2 + (Double_t)itemp*(maxPhi2-minPhi2)/(Double_t(nofPhi2Bins));
338 
339  fhThetaVsRadiusVsMomentumAll = new TH3F("fhThetaVsRadiusVsMomentumAll",
340  "Track Theta vs Point Radius/Z vs Track Momentum, all;p [GeV/c];r [cm];#theta [#circ]",
341  nofMomBins,momBins,
342  nofRadiusBins,rRadBins,
343  nofThetaBins,theBins);
344  for ( Int_t istat = 0 ; istat < fDigiPar->GetNStations() ; istat++ ) {
345  fhThetaVsRadiusVsMomentum[istat] = new TH3F(Form("fhThetaVsRadiusVsMomentum_s%d",istat+1),
346  Form("Track Theta vs Point Radius vs Track Momentum on station %d;p [GeV/c];r [cm];#theta [#circ]",istat+1),
347  nofMomBins,momBins,
348  nofRadiusBins,radBins,
349  nofThetaBins,theBins);
350  ffThetaVsRadiusVsMomentum[istat] = new TF2 (Form("ffThetaVsRadiusVsMomentum_s%d",istat+1),
351  "([0]+1/(x+[1]*[4]+[2]))*y/[4]+[3]",0.,20.,minRRadius,maxRRadius);
352  fHistoList->Add(fhThetaVsRadiusVsMomentum[istat]);
353  fHistoList->Add(ffThetaVsRadiusVsMomentum[istat]);
354  }
355  fHistoList->Add(fhThetaVsRadiusVsMomentumAll);
356 
357 
358 
359  fhRadiusVsAngle = new TH2F("fhRadiusVsAngle","Radius ratio/z ratio vs phi angle difference;#phi_{2}-#phi_{1} [#circ];#frac{r_{2}}{r_{1}}#frac{z_{1}}{z_{2}}",800,-40.,40.,200,0.,2.);
360  fHistoList->Add(fhRadiusVsAngle);
361  Int_t crHist = 0;
362  for ( Int_t istat1 = 0 ; istat1 < fDigiPar->GetNStations() ; istat1++ ) {
363  for ( Int_t istat2 = istat1+1 ; istat2 < fDigiPar->GetNStations() ; istat2++ ) {
364  fhRadiusVsAnglePair[crHist] = new TH2F(Form("fhRadiusVsAnglePair_s%d_s%d",istat1+1,istat2+1),
365  Form("Radius ratio/z ratio vs phi angle difference, stat%d vs stat%d;#phi_{2} - #phi_{1} [#circ];#frac{r_{2}}{r_{1}}#frac{z_{1}}{z_{2}}",istat1+1,istat2+1),
366  800,-40.,40.,200,0.,2.);
367  fHistoList->Add(fhRadiusVsAnglePair[crHist]);
368 
369  fhMomentumVsPhiDiffVsRadius[crHist] = new TH3F(Form("fhMomentumVsPhiDiffVsRadius_s%d_s%d",istat1+1,istat2+1),
370  Form("Momentum vs phi angle difference vs radius, stat%d vs stat%d;r [cm];#phi_{2} - #phi_{1} [#circ];p [GeV/c]",istat1+1,istat2+1),
371  nofRadiusBins,radBins,
372  nofPhiBins,phiBins,
373  nofMomBins,momBins);
374  fHistoList->Add(fhMomentumVsPhiDiffVsRadius[crHist]);
375 
376 
377  ffMomentumVsPhiDiffVsRadius[crHist] = new TF2(Form("ffMomentumVsPhiDiffVsRadius_s%d_s%d",istat1+1,istat2+1),
378  "([0]+[1]/x)/y",0.,100.,0.,40.);
379  fHistoList->Add(ffMomentumVsPhiDiffVsRadius[crHist]);
380 
381  fhMomentumVsPhiDiffAll[crHist] = new TH2F(Form("fhMomentumVsPhiDiffAll_s%d_s%d",istat1+1,istat2+1),
382  Form("Momentum vs phi angle difference, stat%d vs stat%d;#phi_{2} - #phi_{1} [#circ];p [GeV/c]",istat1+1,istat2+1),
383  nofDetPhiBins,minPhi,maxPhi,
384  nofDetMomBins,minMom,maxMom);
385  fHistoList->Add(fhMomentumVsPhiDiffAll[crHist]);
386 
387  ffMomentumVsPhiDiffAll[crHist] = new TF1(Form("ffMomentumVsPhiDiffAll_s%d_s%d",istat1+1,istat2+1),
388  "[0]/x",0.,40.);
389  fHistoList->Add(ffMomentumVsPhiDiffAll[crHist]);
390 
391  for ( Int_t irad = 0 ; irad < 20 ; irad++ ) {
392  fhMomentumVsPhiDiff[crHist][irad] = new TH2F(Form("fhMomentumVsPhiDiff_s%d_s%d_r%d",istat1+1,istat2+1,irad),
393  Form("Momentum vs phi angle difference, stat%d vs stat%d, rad %d to %d;#phi_{2} - #phi_{1} [#circ];p [GeV/c]",istat1+1,istat2+1,irad*5,(irad+1)*5),
394  nofDetPhiBins,minPhi,maxPhi,
395  nofDetMomBins,minMom,maxMom);
396  fHistoList->Add(fhMomentumVsPhiDiff[crHist][irad]);
397 
398  ffMomentumVsPhiDiff[crHist][irad] = new TF1(Form("ffMomentumVsPhiDiff_s%d_s%d_r%d",istat1+1,istat2+1,irad),
399  "[0]/x",0.,40.);
400  fHistoList->Add(ffMomentumVsPhiDiff[crHist][irad]);
401  }
402 
403  fhTrackPhiVsHitPhis [crHist] = new TH3F(Form("fhTrackPhiVsHitPhis_s%d_s%d",istat1+1,istat2+1),
404  Form("Track phi vs hit1 phi vs phi difference/z ratio, stat%d vs stat%d;#phi_{1} [#circ];(#phi_{2} - #phi_{1})*z__{2}/z_{1} [#circ];#phi_{T} [#circ]",istat1+1,istat2+1),
405  nofPhi1Bins,phi1Bins,
406  nofPhi2Bins,phi2Bins,
407  nofPhi1Bins,phi1Bins);
408  fhTrackPhiVsCalcPhi [crHist] = new TH2F(Form("fhTrackPhiVsCalcPhi_s%d_s%d",istat1+1,istat2+1),
409  Form("Track MC phi vs track calculated phi, stat%d vs stat%d;#phi_{calc} [#circ];#phi_{MC} [#circ]",istat1+1,istat2+1),
410  nofPhi1Bins,phi1Bins,
411  nofPhi1Bins,phi1Bins);
412  fHistoList->Add(fhTrackPhiVsHitPhis[crHist]);
413  fHistoList->Add(fhTrackPhiVsCalcPhi[crHist]);
414 
415 
416  fhMomTransVsHRadius[crHist] = new TH2F(Form("fhMomTransVsHRadius_s%d_s%d",istat1+1,istat2+1),
417  Form("Track MC momentum trans vs track calculated helix radius;radius_{radius} [cm];pt_{MC} [GeV/c]"),
418  nofDetMomBins,0.,800.,
419  nofDetMomBins,0.,maxMom/3.);
420  fHistoList->Add(fhMomTransVsHRadius[crHist]);
421 
422  crHist++;
423  }
424  }
425  for ( Int_t itemp = 0 ; itemp < 3 ; itemp++ ) {
426  ffRadiusVsAngle[itemp] = new TF1(Form("ffRadiusVsAngle%d",itemp),"pol2",-40.,40.);
427  fHistoList->Add(ffRadiusVsAngle[itemp]);
428  }
429 
430 
431 }
432 // -------------------------------------------------------------------------
433 
434 // ----- Private method Finish --------------------------------------------
435 void PndGemMagneticFieldVsTrackParameters::Finish() {
436 
437  for ( Int_t istat = 0 ; istat < fDigiPar->GetNStations() ; istat++ ) {
438  PndGemStation* station = (PndGemStation*)fDigiPar->GetStation(istat);
439  Double_t zStation = station->GetZ();
440 
441  ffThetaVsRadiusVsMomentum[istat]->SetParameters(fTrackFinderOnHits_ParTheta0,
442  fTrackFinderOnHits_ParTheta1,
443  fTrackFinderOnHits_ParTheta2,
444  fTrackFinderOnHits_ParTheta3,
445  zStation);
446  }
447  for ( Int_t itemp = 0 ; itemp < 3 ; itemp++ ) {
448  ffRadiusVsAngle[itemp]->SetParameters(fTrackFinderOnHits_ParRadPhi0+(-1+itemp)*0.05,0.,fTrackFinderOnHits_ParRadPhi2);
449  }
450  Int_t crHist = 0;
451  for ( Int_t istat1 = 0 ; istat1 < fDigiPar->GetNStations() ; istat1++ ) {
452  PndGemStation* stat1 = (PndGemStation*)fDigiPar->GetStation(istat1);
453  Double_t zStation1 = stat1->GetZ();
454  for ( Int_t istat2 = istat1+1 ; istat2 < fDigiPar->GetNStations() ; istat2++ ) {
455  PndGemStation* stat2 = (PndGemStation*)fDigiPar->GetStation(istat2);
456  Double_t zStation2 = stat2->GetZ();
457 
458  Double_t zCuDiff = zStation2*zStation2*zStation2-zStation1*zStation1*zStation1;
459  Double_t zSqDiff = zStation2*zStation2-zStation1*zStation1;
460  Double_t zDiff = zStation2-zStation1;
461 
462  Double_t par0_mom = fTrackFinderOnHits_ParMat0[0]*zCuDiff+fTrackFinderOnHits_ParMat0[1]*zSqDiff+fTrackFinderOnHits_ParMat0[2]*zDiff;
463  Double_t par1_mom = fTrackFinderOnHits_ParMat1[0]*zCuDiff+fTrackFinderOnHits_ParMat1[1]*zSqDiff+fTrackFinderOnHits_ParMat1[2]*zDiff;
464 
465  ffMomentumVsPhiDiffVsRadius[crHist]->SetParameters(par0_mom,par1_mom);
466  ffMomentumVsPhiDiffAll [crHist]->SetParameter(0,par0_mom);
467  for ( Int_t irad = 0 ; irad < 20 ; irad++ ) {
468  Double_t par = par0_mom + par1_mom/(Double_t(irad)*5.+2.5);
469  ffMomentumVsPhiDiff [crHist][irad]->SetParameter(0,par);
470  }
471 
472  crHist++;
473  }
474  }
475 
476  gDirectory->mkdir("GemMFVTP");
477  gDirectory->cd("GemMFVTP");
478  TIter next(fHistoList);
479  while ( TH1* histo = ((TH1*)next()) ) histo->Write();
480  gDirectory->cd("..");
481 }
482 // ------------------------------------------------------------
483 
484 ClassImp(PndGemMagneticFieldVsTrackParameters)
PndGemDigiPar::GetTrackFinderOnHits_ParTheta1
Double_t GetTrackFinderOnHits_ParTheta1()
Definition: PndGemDigiPar.h:68
run
Int_t run
Definition: autocutx.C:47
PndGemDigiPar::GetTrackFinderOnHits_ParThetaB
Double_t GetTrackFinderOnHits_ParThetaB()
Definition: PndGemDigiPar.h:65
PndGemDigiPar::GetTrackFinderOnHits_ParTheta3
Double_t GetTrackFinderOnHits_ParTheta3()
Definition: PndGemDigiPar.h:70
map
PndTransMap * map
Definition: sim_emc_apd.C:99
CAMath::Sqrt
static T Sqrt(const T &x)
Definition: PndCAMath.h:37
PndGemStation::GetZ
Double_t GetZ(Int_t it=0)
Definition: PndGemStation.cxx:102
CAMath::Sin
static T Sin(const T &x)
Definition: PndCAMath.h:42
par
Double_t par[3]
Definition: anaLmdCluster.C:146
PndMCTrack::GetMomentum
TVector3 GetMomentum() const
Definition: PndMCTrack.h:78
PndGemDigiPar::GetTrackFinderOnHits_ParMat1
Double_t GetTrackFinderOnHits_ParMat1(Int_t in)
Definition: PndGemDigiPar.h:76
PndGemDigiPar
Digitization Parameter Class for GEM part.
Definition: PndGemDigiPar.h:31
PndGemDigiPar::GetTrackFinderOnHits_ParTheta0
Double_t GetTrackFinderOnHits_ParTheta0()
Definition: PndGemDigiPar.h:67
PndGemDigiPar::GetTrackFinderOnHits_ParMat0
Double_t GetTrackFinderOnHits_ParMat0(Int_t in)
Definition: PndGemDigiPar.h:75
CAMath::Abs
static T Abs(const T &x)
Definition: PndCAMath.h:39
PndGemMCPoint::GetSensorId
Int_t GetSensorId() const
Definition: PndGemMCPoint.h:90
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
PndGemDigiPar::GetNStations
Int_t GetNStations()
Definition: PndGemDigiPar.h:45
name
TString name
Definition: createRootGeoFile.C:82
PndGemDigiPar::GetStationNr
Int_t GetStationNr(Int_t sensorId)
Definition: PndGemDigiPar.h:54
PndGemDigiPar::GetTrackFinderOnHits_ParRadPhi2
Double_t GetTrackFinderOnHits_ParRadPhi2()
Definition: PndGemDigiPar.h:73
ClassImp
ClassImp(PndAnaContFact)
iVerbose
Int_t iVerbose
Definition: Lars/runMvdDigi.C:7
PndMCTrack::GetMotherID
Int_t GetMotherID() const
Definition: PndMCTrack.h:74
next
static int next[96]
Definition: ranlxd.cxx:374
Pi
Double_t Pi
Definition: PndGeoHypGeCrystal.cxx:10
PndGemDigiPar::GetTrackFinderOnHits_ParThetaA
Double_t GetTrackFinderOnHits_ParThetaA()
Definition: PndGemDigiPar.h:64
PndGemDigiPar::GetTrackFinderOnHits_ParTheta2
Double_t GetTrackFinderOnHits_ParTheta2()
Definition: PndGemDigiPar.h:69
PndGemDigiPar::GetTrackFinderOnHits_ParRadPhi0
Double_t GetTrackFinderOnHits_ParRadPhi0()
Definition: PndGemDigiPar.h:72
PndGemDigiPar::GetStation
PndGemStation * GetStation(Int_t iStation)
Definition: PndGemDigiPar.cxx:176
Generated on Wed Apr 3 2019 08:43:56 for FairRoot/PandaRoot by   doxygen 1.8.5