FairRoot/PandaRoot
simpleMCandRECmatch.C
Go to the documentation of this file.
1 //################################################################
2 //# Macro for study difference between MC and REC tracks (with simple case: 1 trk per event)
3 //# author: Anastasia Karavdina
4 //# date: 14/04/2013
5 //#
6 //# to compile it add in "# install #" part of CMakeLists.txt lines:
7 //# add_executable(simple_mc_rec_match simpleMCandRECmatch.C)
8 //# target_link_libraries(simple_mc_rec_match ${ROOT_LIBRARIES} Lmd GeoBase ParBase Geom PndData TrkBase VMC EG GeomPainter generalTools FairTools LmdReco LmdTrk Geane trackrep RecoHits genfitAdapters genfit SdsReco Sds Stt Fts Proof MathMore Minuit FairDB Base)
9 //#
10 //# to run it (and see options):
11 //# ${PANDAROOT}/build/bin/./simple_mc_rec_match --help
12 //################################################################
13 
14 
15 #include <iostream>
16 #include <sstream>
17 
18 #include<TApplication.h>
19 #include<TCanvas.h>
20 #include<TROOT.h>
21 #include<TString.h>
22 #include<TChain.h>
23 #include<TFile.h>
24 #include<TClonesArray.h>
25 #include<TSystem.h>
26 #include<TH1.h>
27 #include<TH2.h>
28 #include<TF1.h>
29 #include<TRotation.h>
30 #include<TVector3.h>
31 #include<TMath.h>
32 #include<TGaxis.h>
33 #include<TNtuple.h>
34 #include<TMatrixD.h>
35 #include<TStopwatch.h>
36 #include<TNtuple.h>
37 #include<PndMCTrack.h>
38 #include<PndTrack.h>
39 #include<PndSdsMCPoint.h>
40 #include "PndLmdDim.h"
41 
42 //lmd track
43 #include<PndLinTrack.h>
44 #include<PndTrackCand.h>
45 #include<PndTrackCandHit.h>
46 #include<PndSdsHit.h>
47 #include<PndSdsClusterStrip.h>
48 #include<PndSdsDigiStrip.h>
49 #include<PndSdsDigiPixel.h>
50 #include<PndSdsClusterPixel.h>
51 #include<PndSdsMergedHit.h>
52 
53 // needed for geane backtracking
54 #include<FairRunAna.h>
55 #include<FairRootManager.h>
56 #include<FairGeane.h>
57 #include<FairRtdbRun.h>
58 #include<FairRuntimeDb.h>
59 #include<FairParRootFileIo.h>
60 #include<FairLogger.h>
61 #include<FairTrackParH.h>
62 
63 #include<string>
64 
65 using namespace std;
66 int main(int __argc,char *__argv[]) {
67  //gROOT->Macro("/PANDA/pandaroot/macro/lmd/Style_Imported_Style.C");
68  // gROOT->SetStyle("Imported_Style");
69 
70  //TODO: read this like params!
71  // const int nEvents=500000;
72  int nEvents=1000;
73  int nMCtracks=1;
74  int startEvent=0;
75  TString storePath="/data/FAIRsorf/pandaroot/trunk/macro/lmd/tmpOutput";
76  double Plab=15.;
77  int verboseLevel=5;
78  int mh = 0;//if>0 : use merged hit collection
79  int dnu=0;//if>0 : use namespace for pixel
80  std::string startStr="", momStr="", nStr="", pathStr="", verbStr="" , mcTrkStr="", useNewDStr="",usedMHStr="";
81  // decode arguments
82  if( __argc>1 && ( strcmp( __argv[1], "-help" ) == 0
83  || strcmp( __argv[1], "--help" ) == 0 ) ){
84 
85  std::cout << "This is script for comparision reconstructed and simulated tracks with parameters\n"
86  <<"-s start event \n"
87  <<"-n Number of events \n"
88  <<"-mom Beam Momentum \n"
89  <<"-path path to the file(s) \n"
90  <<"-v verbose Level (if>0, print out some information) \n"
91  <<"Have fun! \n"
92  << std::endl;
93  return 0;
94  }
95  while ((optind < (__argc-1) ) && (__argv[optind][0]=='-')) {
96  bool found=false;
97  std::string sw = __argv[optind];
98  if (sw=="-s") {
99  optind++;
100  startStr = __argv[optind];
101  found=true;
102  }
103  if (sw=="-n"){
104  optind++;
105  nStr = __argv[optind];
106  found=true;
107  }
108  if (sw=="-path"){
109  optind++;
110  pathStr = __argv[optind];
111  found=true;
112  }
113  if (sw=="-mom"){
114  optind++;
115  momStr = __argv[optind];
116  found=true;
117  }
118  if (sw=="-v"){
119  optind++;
120  verbStr = __argv[optind];
121  found=true;
122  }
123  if (!found){
124  std::cout<< "Unknown switch: "
125  << __argv[optind] <<std::endl;
126  optind++;
127  }
128 
129  while ( (optind < __argc ) && __argv[optind][0]!='-' ) optind++;
130  }
131 
132  std::stringstream startSStr(startStr), momSStr(momStr), nSStr(nStr), pathSStr(pathStr), verbSStr(verbStr);
133  startSStr >> startEvent;
134  momSStr >> Plab;
135  nSStr >> nEvents;
136  pathSStr >> storePath;
137  verbSStr >> verboseLevel;
138  nMCtracks=1;
139  cout<<"====================================="<<endl;
140  cout<<"some INFO about this macro params: "<<endl;
141  cout<<"expected number of events: "<<nEvents<<endl;
142  cout<<"Pbeam: "<<Plab<<endl;
143  cout<<"Will be used Path: "<<storePath<<endl;
144  cout<<"====================================="<<endl;
145  // ---- Load libraries -------------------------------------------------
146  gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.C");
147  gSystem->Load("libLmdTrk");
148  // gROOT->LoadMacro("line3Dfit.C");
149  // ------------------------------------------------------------------------
150 
151  // ----- Timer --------------------------------------------------------
152  TStopwatch timer;
153  timer.Start();
154  // ------------------------------------------------------------------------
155 
156  // ---- Input files --------------------------------------------------------
157  TString simMC=storePath+"/Lumi_MC_";
158  simMC += startEvent;
159  simMC += ".root";
160  TChain tMC("pndsim");
161  tMC.Add(simMC);
162 
163  TString DigiFile = storePath+"/Lumi_digi_";
164  DigiFile += startEvent;
165  DigiFile += ".root";
166  TChain tdigiHits("pndsim");
167  tdigiHits.Add(DigiFile);
168 
169  TString recHit=storePath+"/Lumi_reco";
170  // if(dnu>0) recHit+="Merged";
171  recHit +="_";
172  recHit += startEvent;
173  recHit += ".root";
174  TChain tHits("pndsim");
175  tHits.Add(recHit);
176 
177  TString recHitmerged=storePath+"/Lumi_recoMerged_";
178  recHitmerged += startEvent;
179  recHitmerged += ".root";
180  TChain tHitsMerged("pndsim");
181  tHitsMerged.Add(recHitmerged);
182 
183 
184  TString trkCand = storePath+"/Lumi_TCand_";
185  trkCand += startEvent;
186  trkCand += ".root";
187  TChain tTrkCand("pndsim");
188  tTrkCand.Add(trkCand);
189 
190  TString recTrack;
191  TChain tTrkRec("pndsim");
192  recTrack = storePath+"/Lumi_Track_";
193  recTrack += startEvent;
194  recTrack += ".root";
195  tTrkRec.Add(recTrack);
196 
197  TString geaneFile = storePath+"/Lumi_Geane_";
198  geaneFile += startEvent;
199  geaneFile += ".root";
200  TChain tgeane("pndsim");
201  tgeane.Add(geaneFile);
202 
203  // ---------------------------------------------------------------------------------
204 
205  // ---- Output file ----------------------------------------------------------------
206  TString out=storePath+"/Lumi_compare_MC_and_REC_trks_";
207  out += startEvent;
208  out += ".root";
209  TFile *f = new TFile(out,"RECREATE");
210  // ---------------------------------------------------------------------------------
211 
212  //--- MC info -----------------------------------------------------------------
213  TClonesArray* true_tracks=new TClonesArray("PndMCTrack");
214  tMC.SetBranchAddress("MCTrack",&true_tracks); //True Track to compare
215 
216  TClonesArray* true_points=new TClonesArray("PndSdsMCPoint");
217  tMC.SetBranchAddress("LMDPoint",&true_points); //True Points to compare
218  //----------------------------------------------------------------------------------
219 
220 
221  //--- Digitization info ------------------------------------------------------------
222  TClonesArray* fStripClusterArray;
223  fStripClusterArray = new TClonesArray("PndSdsClusterPixel");
224  tHits.SetBranchAddress("LMDPixelClusterCand",&fStripClusterArray);
225 
226 
227  TClonesArray* fStripDigiArray;
228  fStripDigiArray = new TClonesArray("PndSdsDigiPixel");
229  tdigiHits.SetBranchAddress("LMDPixelDigis",&fStripDigiArray);
230  //----------------------------------------------------------------------------------
231 
232  //--- Real Hits --------------------------------------------------------------------
233  TClonesArray* rechit_array;
234  rechit_array = new TClonesArray("PndSdsMergedHit");
235  tHitsMerged.SetBranchAddress("LMDHitsMerged",&rechit_array); //Points for Tracks
236  //----------------------------------------------------------------------------------
237 
238 
239  //--- Track Candidate ---------------------------------------------------------------
240  TClonesArray* trkcand_array=new TClonesArray("PndTrackCand");
241  tTrkCand.SetBranchAddress("LMDTrackCand",&trkcand_array); //Points for Track Canidates
242  //-----------------------------------------------------------------------------------
243 
244  //--- Real tracks -------------------------------------------------------------------
245  // TClonesArray* rec_trk=new TClonesArray("PndLinTrack");
246  // tTrkRec.SetBranchAddress("LMDTrack",&rec_trk); //Tracks
247  TClonesArray* rec_trk=new TClonesArray("PndTrack");
248  tTrkRec.SetBranchAddress("LMDPndTrack",&rec_trk); //Tracks
249  //----------------------------------------------------------------------------------
250 
251  //--- Geane info ------------------------------------------------------------------
252  TClonesArray* geaneArray =new TClonesArray("FairTrackParH");
253  tgeane.SetBranchAddress("GeaneTrackFinal",&geaneArray); //Tracks with helix parametrisation
254 
255  cout<<"And we'll make some hists"<<endl;
256  //--- Output histogram -----------------------------------------------------
257  double thetarange[2]={0.001,0.01};
258  double thetam = thetarange[0];
259  if(Plab<5) thetam = thetarange[1];
260 
261  //Near IP
262  TH1 *hResMom = new TH1F("hResMom","P_{MC}-P_{REC};#deltaP,GeV/c",1e3,-1e-4,1e-4);
263  TH1 *hErrMom = new TH1F("hErrMom","#sigma_{P};#sigmaP,GeV/c",1e3,0,1e-3);
264  TH1 *hPullMom = new TH1F("hPullMom","(P_{MC}-P_{REC})/#sigma_{P};",1e3,-1e1,1e1);
265  // TH1 *hResTheta = new TH1F("hResTheta","#theta_{MC}-#theta_{REC};#delta#theta,rad",1e3,-1e-2,1e-2);
266  TH1 *hResTheta = new TH1F("hResTheta","#theta_{MC}-#theta_{REC};#delta#theta,rad",1e2,-thetam,thetam);//TEST
267  TH2 *hThetaResTheta = new TH2F("hThetaResTheta","#theta_{MC}-#theta_{REC};#theta,rad;#delta#theta,rad",2e1,0,1e-2,1e2,-thetam,thetam);//TEST
268  TH1 *hErrTheta = new TH1F("hErrTheta","#sigma(#theta_{REC});#sigma,rad",1e3,0,10*thetam);
269  TH1 *hPullTheta = new TH1F("hPullTheta","(#theta_{MC}-#theta_{REC})/#sigma_{#theta};",1e2,-10,10);
270  TH1 *hResPhi = new TH1F("hResPhi","#phi_{MC}-#phi_{REC};#delta#phi,rad",2e3,-1.,1.);
271  TH1 *hErrPhi = new TH1F("hErrPhi","#sigma(#phi_{REC});#sigma,rad",1e3,0,0.1);
272  TH1 *hPullPhi = new TH1F("hPullPhi","(#phi_{MC}-#phi_{REC})/#sigma_{#phi};",1e2,-10,10);
273 
274  TH1 *hResPointPx = new TH1F("hResPointPx","Px_{MC}-Px_{REC};#deltaPx, GeV/c",1e2,-0.01,0.01);
275  TH1 *hErrPointPx = new TH1F("hErrPointPx","#sigma_{Px};#sigmaPx, GeV/c",1e3,0,0.01);
276  TH1 *hPullPointPx = new TH1F("hPullPointPx","(Px_{MC}-Px_{REC})/#sigma_{Px};(Px_{MC}-Px_{REC})/#sigma_{Px}",1e2,-10,10);
277 
278  TH1 *hResPointPy = new TH1F("hResPointPy","Py_{MC}-Py_{REC};#deltaPy, GeV/c",1e2,-0.01,0.01);
279  TH1 *hErrPointPy = new TH1F("hErrPointPy","#sigma_{Py};#sigmaPy, GeV/c",1e3,0,0.01);
280  TH1 *hPullPointPy = new TH1F("hPullPointPy","(Py_{MC}-Py_{REC})/#sigma_{Py};(Py_{MC}-Py_{REC})/#sigma_{Py}",1e2,-10,10);
281 
282  TH1 *hResPointPz = new TH1F("hResPointPz","Pz_{MC}-Pz_{REC};#deltaPz, GeV/c",1e2,-1e-3,1e-3);
283  TH1 *hErrPointPz = new TH1F("hErrPointPz","#sigma_{Pz};#sigmaPz, GeV/c",1e3,0,1e-2);
284  TH1 *hPullPointPz = new TH1F("hPullPointPz","(Pz_{MC}-Pz_{REC})/#sigma_{Pz};(Pz_{MC}-Pz_{REC})/#sigma_{Pz}",1e2,-10,10);
285 
286  TH1 *hResPointX = new TH1F("hResPointX","X_{MC}-X_{REC};#deltaX,cm",1e2,-2.,2.);
287  TH1 *hPullPointX = new TH1F("hPullPointX","(X_{MC}-X_{REC})/#sigma_{X};(X_{MC}-X_{REC})/#sigma_{X}",1e2,-10,10);
288  TH1 *hResPointY = new TH1F("hResPointY","Y_{MC}-Y_{REC};#deltaY,cm",1e2,-2.,2.);
289  TH1 *hPullPointY = new TH1F("hPullPointY","(Y_{MC}-Y_{REC})/#sigma_{Y};(Y_{MC}-Y_{REC})/#sigma_{Y}",1e2,-10,10);
290  TH1 *hResPointZ = new TH1F("hResPointZ","Z_{MC}-Z_{REC};#deltaZ,cm",1e3,-0.15,0.15);
291  TH1 *hPullPointZ = new TH1F("hPullPointZ","(Z_{MC}-Z_{REC})/#sigma_{Z};(Z_{MC}-Z_{REC})/#sigma_{Z}",1e2,-10,10);
292 
293  //Near 1st LMD plane
294  TH1 *hhits = new TH1I("hhits","number of hits in trk",7,0,7);
295  TH1 *hchi2 = new TH1F("hchi2","#chi^2 for reconstructed tracks;#chi^2;",1.5e2,0,15.);
296  TH1 *hResLumiTrkMom = new TH1F("hResLumiTrkMom","P_{MC}-P_{REC}(near Lumi);#deltaP,GeV/c",1e3,-6e-7,6e-7);
297  TH1 *hResLumiTrkTheta = new TH1F("hResLumiTrkTheta","#theta_{MC}-#theta_{REC}(near Lumi);#delta#theta,rad",1e3,-6e-3,6e-3);
298  TH1 *hResLumiTrkPhi = new TH1F("hResLumiTrkPhi","#phi_{MC}-#phi_{REC}(near Lumi);#delta#phi,rad",2e3,-1e-1,1e-1);
299  TH1 *hResLumiTrkPointX = new TH1F("hResLumiTrkPointX","X_{MC}-X_{REC}(near Lumi);#deltaX,cm",1e2,-0.02,0.02);
300  TH1 *hResLumiTrkPointY = new TH1F("hResLumiTrkPointY","Y_{MC}-Y_{REC}(near Lumi);#deltaY,cm",1e2,-0.02,0.02);
301  TH1 *hResLumiTrkPointZ = new TH1F("hResLumiTrkPointZ","Z_{MC}-Z_{REC}(near Lumi);#deltaZ,cm",1e2,-0.02,0.02);
302 
303  TH1 *hResLumiTrkPointPx = new TH1F("hResLumiTrkPointPx","Px_{MC}-Px_{REC}(near Lumi);#deltaPx, GeV/c",1e2,-0.01,0.01);
304  TH1 *hResLumiTrkPointPy = new TH1F("hResLumiTrkPointPy","Py_{MC}-Py_{REC}(near Lumi);#deltaPy, GeV/c",1e2,-0.01,0.01);
305  TH1 *hResLumiTrkPointPz = new TH1F("hResLumiTrkPointPz","Pz_{MC}-Pz_{REC}(near Lumi);#deltaPz, GeV/c",1e2,-0.001,0.001);
306  // TH1 *hResLumiTrkPointP = new TH1F("hResLumiTrkPointP","P_{MC}-P_{REC}(near Lumi);#deltaP, GeV/c",1e2,-0.001,0.001);
307  // TH2 *hResLumiTrkPointPmcPrec = new TH2F("hResLumiTrkPointPmcPrec","P_{MC} vs. P_{REC}(near Lumi);P_{rec}, GeV/c;P_{mc}, GeV/c",1e2,Plab-0.1,Plab+0.1,1e2,Plab-0.1,Plab+0.1);
308 
309  TH1 *hResLumiTrkPointXErr = new TH1F("hResLumiTrkPointXErr","#sigma(X_{REC})(near Lumi);#sigma_{X},cm",1e2,0,0.02);
310  TH1 *hResLumiTrkPointYErr = new TH1F("hResLumiTrkPointYErr","#sigma(Y_{REC})(near Lumi);#sigma_{Y},cm",1e2,0,0.02);
311  TH1 *hResLumiTrkPointZErr = new TH1F("hResLumiTrkPointZErr","#sigma(Z_{REC})(near Lumi);#sigma_{Z},cm",1e2,0,0.02);
312  TH1 *hResLumiTrkPointPxErr = new TH1F("hResLumiTrkPointPxErr","#sigma(Px_{REC})(near Lumi);#sigma_{Px}, GeV/c",1e2,0,0.01);
313  TH1 *hResLumiTrkPointPyErr = new TH1F("hResLumiTrkPointPyErr","#sigma(Py_{REC})(near Lumi);#sigma_{Py}, GeV/c",1e2,0,0.01);
314  TH1 *hResLumiTrkPointPzErr = new TH1F("hResLumiTrkPointPzErr","#sigma(Pz_{REC})(near Lumi);#sigma_{Pz}, GeV/c",1e2,0,0.001);
315 
316  TH1 *hResLumiTrkPointXPull = new TH1F("hResLumiTrkPointXPull","(X_{MC}-X_{REC})/#sigma (near Lumi) ;(X_{MC}-X_{REC})/#sigma",1e2,-10.,10.);
317  TH1 *hResLumiTrkPointYPull = new TH1F("hResLumiTrkPointYPull","(Y_{MC}-Y_{REC})/#sigma (near Lumi);(Y_{MC}-Y_{REC})/#sigma",1e2,-10.,10.);
318  TH1 *hResLumiTrkPointZPull = new TH1F("hResLumiTrkPointZPull","(Z_{MC}-Z_{REC})/#sigma (near Lumi);(Z_{MC}-Z_{REC})/#sigma",1e3,-100.,100.);
319 
320  TH1 *hResLumiTrkPointPxPull = new TH1F("hResLumiTrkPointPxPull","(Px_{MC}-Px_{REC})/#sigma (near Lumi);(Px_{MC}-Px_{REC})/#sigma",1e2,-10,10.);
321  TH1 *hResLumiTrkPointPyPull = new TH1F("hResLumiTrkPointPyPull","(Py_{MC}-Py_{REC})/#sigma (near Lumi);(Py_{MC}-Py_{REC})/#sigma",1e2,-10,10);
322  TH1 *hResLumiTrkPointPzPull = new TH1F("hResLumiTrkPointPzPull","(Pz_{MC}-Pz_{REC})/#sigma (near Lumi);(Pz_{MC}-Pz_{REC})/#sigma",1e2,-10,10);
323  TH1 *hResLumiTrkThetaPull = new TH1F("hResLumiTrkThetaPull","(#theta_{MC}-#theta_{REC})/#sigma (near Lumi);#delta#theta, rad",1e2,-10,10);
324  TH1 *hResLumiTrkPhiPull = new TH1F("hResLumiTrkPhiPull","(#phi_{MC}-#phi_{REC})/#sigma (near Lumi);#delta#phi, rad",1e2,-10,10);
325 
326  // TH2 *hMCidRefID = new TH2I("hMCidRefID","; MCid; RefID",100,0,100,100,0,100);
327 
328  TNtuple *nBadTrks = new TNtuple("nBadTrks","Info about _bad_ rec.tracks ","xrec:yrec:zrec:pxrec:pyrec:pzrec:nrechits:xmc:ymc:zmc:pxmc:pymc:pzmc:mvID");
329 
330  int glBADGEANE=0;
331  // int glBadEv = 0;
332  // int glNoisehit = 0;// total number of noise hits
333  for (Int_t j=0; j<nEvents; j++){
334  // Read GEANE & MC info -----------------------------------------------------------------
335  // if(kf<1)
336  tTrkRec.GetEntry(j);
337  tgeane.GetEntry(j);
338  tMC.GetEntry(j);
339  tTrkCand.GetEntry(j);
340  tHits.GetEntry(j);
341  tdigiHits.GetEntry(j);
342  tHitsMerged.GetEntry(j);
343  const int nGeaneTrks = geaneArray->GetEntriesFast();
344  const int nParticles = true_tracks->GetEntriesFast();
345  const int numTrk = nGeaneTrks;
346  const int nRecHits = rechit_array->GetEntriesFast();
347  const int nMCHits = true_points->GetEntriesFast();
348 
349  const int nTrkCandidates = trkcand_array->GetEntriesFast();
350  const int nRecTrks = rec_trk->GetEntriesFast();
351  if(verboseLevel>0)
352  cout<<"Event #"<<j<<" has "<<nParticles<<" true particles, "<<" out of it "<<nRecHits<<" hits, "<<nTrkCandidates
353  <<" trk-cands, "<<numTrk<<" tracks and "<<nGeaneTrks<<" geane Trks!"<<endl;
354  // if(nParticles!=nMCtracks) continue;
355  for (Int_t iN=0; iN<nGeaneTrks; iN++){// loop over all reconstructed trks
356  FairTrackParH *fRes = (FairTrackParH*)geaneArray->At(iN);
357  Double_t lyambda = fRes->GetLambda();
358  if(lyambda==0){
359  cout<<"GEANE didn't propagate this trk!"<<endl;
360  glBADGEANE++;
361  nBadTrks->Fill(1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6);
362  }
363  if(lyambda==0) continue;
364 
366  TVector3 MomRecPCA = fRes->GetMomentum();
367  TVector3 MomRecPCAnotnorm = MomRecPCA;
368  MomRecPCA *= Plab/MomRecPCA.Mag();
369  TVector3 PosRecPCA = fRes->GetPosition();
370  Double_t errPx = fRes->GetDPx();
371  Double_t errPy = fRes->GetDPy();
372  Double_t errPz = fRes->GetDPz();
373  TVector3 errMomRecPCA(errPx,errPy,errPz);
374  Double_t errX = fRes->GetDX();
375  Double_t errY = fRes->GetDY();
376  Double_t errZ = fRes->GetDZ();
377  TVector3 errPosRecPCA(errX,errY,errZ);
378 
379  Double_t thetaBP = TMath::Pi()/2. - lyambda;
380  // Double_t err_lyambda = fRes->GetDLambda();
381  Double_t phiBP = fRes->GetPhi();
382  // Double_t err_phi = fRes->GetDPhi();
383 
384  //calculate theta & phi errors
385  double fLmPCA = TMath::ASin(MomRecPCA.Z()/MomRecPCA.Mag());
386  double cLmPCA= TMath::Cos(fLmPCA);
387  double sLmPCA= TMath::Sin(fLmPCA);
388  Double_t fPPCA =sqrt(MomRecPCA.X()*MomRecPCA.X()+MomRecPCA.Y()*MomRecPCA.Y()+MomRecPCA.Z()*MomRecPCA.Z());
389  Double_t fDPPCA= (2*MomRecPCA.X()*errMomRecPCA.X()+2*MomRecPCA.Y()*errMomRecPCA.Y()+2*MomRecPCA.Z()*errMomRecPCA.Z())/(2*fPPCA); //dp
390  Double_t err_lyambda = (-((MomRecPCA.Z()*fDPPCA)/pow(fPPCA,2)) + errMomRecPCA.Z()/fPPCA)/ TMath::Sqrt(1 - pow(MomRecPCA.Z(),2)/pow(fPPCA,2));
391  Double_t err_phi = (-((MomRecPCA.Y()*fDPPCA/cLmPCA)/pow(fPPCA,2)) + (errMomRecPCA.Y()/cLmPCA)/fPPCA +(MomRecPCA.Y()*err_lyambda*TMath::Tan(fLmPCA)/cLmPCA)/fPPCA) /TMath::Sqrt(1 - (pow(MomRecPCA.Y(),2)*pow(1/cLmPCA,2))/pow(fPPCA,2));
392 
393 
394  // Double_t CovGEANELAB[6][6];
395  // fRes->GetMARSCov(CovGEANELAB);
396  // Double_t errMomRecBP = fRes->GetDQp();
397  // ///get rid from most probably ghost track ----------
398  // //TODO: find reason for such trks in Kalman
399  // double pca_lim = 1.;//=10*sigma_Xpca~10*{0.093,0.11,0.12,0.22,0.55};
400  // if(Plab<5) pca_lim = 2.;
401  // if(Plab<2) pca_lim = 5.;
402  // if(fabs(PosRecPCA.X())>pca_lim && fabs(PosRecPCA.Y())>pca_lim) continue; // PCA_x and PCA_y should be < 10sigmaX
403  // ///get rid from most probably ghost track (END) ---
404  // ///------------------------------------------------------------------------------------
405 
407  PndTrack *trkpnd = (PndTrack*)rec_trk->At(iN);
408  double chi2 = trkpnd->GetChi2();
409  hchi2->Fill(chi2);
410  FairTrackParP fFittedTrkP = trkpnd->GetParamFirst();
411  TVector3 PosRecLMD(fFittedTrkP.GetX(),fFittedTrkP.GetY(),fFittedTrkP.GetZ());
412  TVector3 MomRecLMD(fFittedTrkP.GetPx(),fFittedTrkP.GetPy(),fFittedTrkP.GetPz());
413  MomRecLMD *=Plab/MomRecLMD.Mag();
414  double covMARS[6][6];
415  fFittedTrkP.GetMARSCov(covMARS);
416  TVector3 errMomRecLMD(sqrt(covMARS[0][0]),sqrt(covMARS[1][1]),sqrt(covMARS[2][2]));
417  TVector3 errPosRecLMD(sqrt(covMARS[3][3]),sqrt(covMARS[4][4]),sqrt(covMARS[5][5]));
418 
419  //calculate theta & phi errors
420  double fLm = TMath::ASin(MomRecLMD.Z()/MomRecLMD.Mag());
421  double cLm= TMath::Cos(fLm);
422  double sLm= TMath::Sin(fLm);
423  Double_t fP =sqrt(MomRecLMD.X()*MomRecLMD.X()+MomRecLMD.Y()*MomRecLMD.Y()+MomRecLMD.Z()*MomRecLMD.Z());
424  Double_t fDP= (2*MomRecLMD.X()*errMomRecLMD.X()+2*MomRecLMD.Y()*errMomRecLMD.Y()+2*MomRecLMD.Z()*errMomRecLMD.Z())/(2*fP); //dp
425  Double_t err_lyambdaLMD = (-((MomRecLMD.Z()*fDP)/pow(fP,2)) + errMomRecLMD.Z()/fP)/ TMath::Sqrt(1 - pow(MomRecLMD.Z(),2)/pow(fP,2));
426  Double_t err_phiLMD = (-((MomRecLMD.Y()*fDP/cLm)/pow(fP,2)) + (errMomRecLMD.Y()/cLm)/fP +(MomRecLMD.Y()*err_lyambdaLMD*TMath::Tan(fLm)/cLm)/fP) /TMath::Sqrt(1 - (pow(MomRecLMD.Y(),2)*pow(1/cLm,2))/pow(fP,2));
428 
429 
430 
431  //Matching between MC & Rec on 1st hit level-----------------------------------
432  int candID = trkpnd->GetRefIndex();
433  PndTrackCand *trkcand = (PndTrackCand*)trkcand_array->At(candID);
434  const int Ntrkcandhits= trkcand->GetNHits();
435  PndSdsMCPoint* MCPointHit;
436  int MCid;
437  bool hitmix = false;
438  if(Ntrkcandhits<4) continue; //require trks with hits on all planes
439  hhits->Fill(Ntrkcandhits);
440  for (Int_t iHit = 0; iHit < Ntrkcandhits; iHit++){ // loop over rec.hits
441  PndTrackCandHit candhit = (PndTrackCandHit)(trkcand->GetSortedHit(iHit));
442  Int_t hitID = candhit.GetHitId();
443  PndSdsHit* myHit = (PndSdsHit*)(rechit_array->At(hitID));
444 
445  //for pixel design
446  PndSdsClusterPixel* myCluster = (PndSdsClusterPixel*)(fStripClusterArray->At(myHit->GetClusterIndex()));
447  PndSdsDigiPixel* astripdigi = (PndSdsDigiPixel*)(fStripDigiArray->At(myCluster->GetDigiIndex(0)));
448  if (astripdigi->GetIndex(0) == -1)
449  continue;
450  PndSdsMCPoint* MCPoint = (PndSdsMCPoint*)(true_points->At(astripdigi->GetIndex(0)));
451  int MCidTOP = MCPoint->GetTrackID();
452  // hMCidRefID->Fill(MCidTOP,myHit->GetRefIndex());
453  if(iHit<1){
454  MCPointHit = MCPoint;
455  MCid = MCidTOP;
456  }
457  else
458  if(MCid!=MCidTOP){
459  cout<<"REC trk contains hits from different MC trks! Skip this event."<<endl;
460  hitmix = true;
461  }
462  }
463  if(hitmix) continue;
465 
466 
468 
470  PndMCTrack *mctrk =(PndMCTrack*) true_tracks->At(MCid);
471  Int_t mcID = mctrk->GetPdgCode();
472  int mvID = mctrk->GetMotherID();
473  TVector3 MomMCpca = mctrk->GetMomentum();
474  TVector3 PosMCpca = mctrk->GetStartVertex();
475  Double_t thetaMC = MomMCpca.Theta();
476  Double_t phiMC = MomMCpca.Phi();
478  hResPointPx->Fill(MomMCpca.X()-MomRecPCA.X());
479  hResPointPy->Fill(MomMCpca.Y()-MomRecPCA.Y());
480  hResPointPz->Fill(MomMCpca.Z()-MomRecPCA.Z());
481  hErrPointPx->Fill(errMomRecPCA.X());
482  hErrPointPy->Fill(errMomRecPCA.Y());
483  hErrPointPz->Fill(errMomRecPCA.Z());
484  hPullPointPx->Fill((MomMCpca.X()-MomRecPCA.X())/errMomRecPCA.X());
485  hPullPointPy->Fill((MomMCpca.Y()-MomRecPCA.Y())/errMomRecPCA.Y());
486  hPullPointPz->Fill((MomMCpca.Z()-MomRecPCA.Z())/errMomRecPCA.Z());
487  hResPointX->Fill(PosMCpca.X()-PosRecPCA.X());
488  hResPointY->Fill(PosMCpca.Y()-PosRecPCA.Y());
489  hResPointZ->Fill(PosMCpca.Z()-PosRecPCA.Z());
490  hPullPointX->Fill((PosMCpca.X()-PosRecPCA.X())/errPosRecPCA.X());
491  hPullPointY->Fill((PosMCpca.Y()-PosRecPCA.Y())/errPosRecPCA.Y());
492  hPullPointZ->Fill((PosMCpca.Z()-PosRecPCA.Z())/errPosRecPCA.Z());
493 
494  hResTheta->Fill(MomMCpca.Theta()-MomRecPCA.Theta());
495  hThetaResTheta->Fill(MomMCpca.Theta(),(MomMCpca.Theta()-MomRecPCA.Theta()));
496  hResPhi->Fill(MomMCpca.Phi()-MomRecPCA.Phi());
497  hPullTheta->Fill((MomMCpca.Theta()-MomRecPCA.Theta())/err_lyambda);
498  hPullPhi->Fill((MomMCpca.Phi()-MomRecPCA.Phi())/err_phi);
499  //TNtuple *nBadTrks = new TNtuple("nBadTrks","Info about _bad_ rec.tracks ","xrec:yrec:zrec:pxrec:pyrec:pzrec:nrechits:xmc:ymc:zmc:pxmc:pymc:pzmc");
500  if(fabs(MomRecPCAnotnorm.Mag()-Plab)>0.1*Plab){
501  nBadTrks->Fill(PosRecPCA.X(),PosRecPCA.Y(),PosRecPCA.Z(),MomRecPCAnotnorm.X(),MomRecPCAnotnorm.Y(),MomRecPCAnotnorm.Z(),Ntrkcandhits,PosMCpca.X(),PosMCpca.Y(),PosMCpca.Z(),MomMCpca.X(),MomMCpca.Y(),MomMCpca.Z(),mvID);
502  cout<<"Event #"<<j<<" contains BAD trk"<<endl;
503  }
504 
505  //Near 1st LMD plane
507  TVector3 PosMClmd = MCPointHit->GetPosition();
508  double pxTrue = MCPointHit->GetPx();
509  double pyTrue = MCPointHit->GetPy();
510  double pzTrue = MCPointHit->GetPz();
511  TVector3 MomMClmd(pxTrue,pyTrue,pzTrue);
512  TVector3 dirMClmd = MomMClmd;
513  dirMClmd *=1./MomMClmd.Mag();
514  double deltaZ = -PosMClmd.Z()+PosRecLMD.Z();
515  // double deltaZ = 0;
516  double xneu=PosMClmd.X()+dirMClmd.X()*deltaZ;
517  double yneu=PosMClmd.Y()+dirMClmd.Y()*deltaZ;
518  double zneu = PosMClmd.Z()+deltaZ;
519  PosMClmd.SetXYZ(xneu,yneu,zneu);
520  MomMClmd = dirMClmd*Plab;
521  // hResLumiTrkPointP->Fill((MomMClmd.Mag()-MomRecLMD.Mag()));
522  // hResLumiTrkPointPmcPrec->Fill(MomRecLMD.Mag(),MomMClmd.Mag());
523 
524 
525  // MomMClmd *=1./MomMClmd.Mag();//TEST
526  // MomRecLMD *=1./MomRecLMD.Mag();//TEST
527  // errMomRecLMD *=1./MomRecLMD.Mag();//TEST
528 
529  // cout<<"MomMClmd.Mag() = "<<MomMClmd.Mag()<<" MomRecLMD.Mag() = "<<MomRecLMD.Mag()<<endl;
530  // cout<<" MC - REC = "<<1e3*(MomMClmd.Mag()-MomRecLMD.Mag())<<" MeV"<<endl;
532  hResLumiTrkPointX->Fill(PosMClmd.X()-PosRecLMD.X());
533  hResLumiTrkPointY->Fill(PosMClmd.Y()-PosRecLMD.Y());
534  hResLumiTrkPointZ->Fill(PosMClmd.Z()-PosRecLMD.Z());
535  hResLumiTrkPointXPull->Fill((PosMClmd.X()-PosRecLMD.X())/errPosRecLMD.X());
536  hResLumiTrkPointYPull->Fill((PosMClmd.Y()-PosRecLMD.Y())/errPosRecLMD.Y());
537  hResLumiTrkPointZPull->Fill((PosMClmd.Z()-PosRecLMD.Z())/errPosRecLMD.Z());
538 
539  hResLumiTrkPointPx->Fill(MomMClmd.X()-MomRecLMD.X());
540  hResLumiTrkPointPy->Fill(MomMClmd.Y()-MomRecLMD.Y());
541  hResLumiTrkPointPz->Fill(MomMClmd.Z()-MomRecLMD.Z());
542  hResLumiTrkPointPxPull->Fill((MomMClmd.X()-MomRecLMD.X())/errMomRecLMD.X());
543  hResLumiTrkPointPyPull->Fill((MomMClmd.Y()-MomRecLMD.Y())/errMomRecLMD.Y());
544  hResLumiTrkPointPzPull->Fill((MomMClmd.Z()-MomRecLMD.Z())/errMomRecLMD.Z());
545  hResLumiTrkTheta->Fill(MomMClmd.Theta()-MomRecLMD.Theta());
546  hResLumiTrkPhi->Fill(MomMClmd.Phi()-MomRecLMD.Phi());
547  hResLumiTrkThetaPull->Fill((MomMClmd.Theta()-MomRecLMD.Theta())/err_lyambdaLMD);
548  hResLumiTrkPhiPull->Fill((MomMClmd.Phi()-MomRecLMD.Phi())/err_phiLMD);
549  hResLumiTrkPointXErr->Fill(errPosRecLMD.X());
550  hResLumiTrkPointYErr->Fill(errPosRecLMD.Y());
551  hResLumiTrkPointZErr->Fill(errPosRecLMD.Z());
552  hResLumiTrkPointPxErr->Fill(errMomRecLMD.X());
553  hResLumiTrkPointPyErr->Fill(errMomRecLMD.Y());
554  hResLumiTrkPointPzErr->Fill(errMomRecLMD.Z());
555 
557  }
558  }
559 
560  //Draw results on few canvases
561 
562  TCanvas *c1 = new TCanvas("pulls_before_bp");
563  c1->Divide(3,2);
564  c1->cd(1);
565  hResLumiTrkPointXPull->Draw();
566  c1->cd(2);
567  hResLumiTrkPointYPull->Draw();
568  c1->cd(3);
569  hResLumiTrkPointZPull->Draw();
570  c1->cd(4);
571  hResLumiTrkPointPxPull->Draw();
572  c1->cd(5);
573  hResLumiTrkPointPyPull->Draw();
574  c1->cd(6);
575  hResLumiTrkPointPzPull->Draw();
576  c1->Write();
577  c1->Close();
578 
579  TCanvas *c2 = new TCanvas("pulls_after_bp");
580  c2->Divide(3,2);
581  c2->cd(1);
582  hPullPointX->Draw();
583  c2->cd(2);
584  hPullPointY->Draw();
585  c2->cd(3);
586  hPullPointZ->Draw();
587  c2->cd(4);
588  hPullPointPx->Draw();
589  c2->cd(5);
590  hPullPointPy->Draw();
591  c2->cd(6);
592  hPullPointPz->Draw();
593  c2->Write();
594  c2->Close();
595 
596  hchi2->Write();
597  hResPointPx->Write();
598  hResPointPy->Write();
599  hResPointPz->Write();
600  hErrPointPx->Write();
601  hErrPointPy->Write();
602  hErrPointPz->Write();
603  hPullPointPx->Write();
604  hPullPointPy->Write();
605  hPullPointPz->Write();
606  hResPointX->Write();
607  hResPointY->Write();
608  hResPointZ->Write();
609  hPullPointX->Write();
610  hPullPointY->Write();
611  hPullPointZ->Write();
612  hResLumiTrkPointX->Write();
613  hResLumiTrkPointY->Write();
614  hResLumiTrkPointZ->Write();
615  hResLumiTrkPointXPull->Write();
616  hResLumiTrkPointYPull->Write();
617  hResLumiTrkPointZPull->Write();
618  hResLumiTrkPointPx->Write();
619  hResLumiTrkPointPy->Write();
620  hResLumiTrkPointPz->Write();
621  hResLumiTrkPointPxPull->Write();
622  hResLumiTrkPointPyPull->Write();
623  hResLumiTrkPointPzPull->Write();
624  hResTheta->Write();
625  hResPhi->Write();
626  hPullTheta->Write();
627  hPullPhi->Write();
628  hResLumiTrkTheta->Write();
629  hResLumiTrkPhi->Write();
630  hResLumiTrkThetaPull->Write();
631  hResLumiTrkPhiPull->Write();
632  hhits->Write();
633  hResLumiTrkPointXErr->Write();
634  hResLumiTrkPointYErr->Write();
635  hResLumiTrkPointZErr->Write();
636  hResLumiTrkPointPxErr->Write();
637  hResLumiTrkPointPyErr->Write();
638  hResLumiTrkPointPzErr->Write();
639  // hResLumiTrkPointP->Write();
640  // hResLumiTrkPointPmcPrec->Write();
641  // hMCidRefID->Write();
642  hThetaResTheta->Write();
643  nBadTrks->Write();
644  f->Close();
645  cout<<"Number of trks where GEANE failed: "<<glBADGEANE<<endl;
646 }
CAMath::ASin
static T ASin(const T &x)
verboseLevel
int verboseLevel
Definition: Lars/runMvdSim.C:7
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
CAMath::Sqrt
static T Sqrt(const T &x)
Definition: PndCAMath.h:37
PndMCTrack::GetPdgCode
Int_t GetPdgCode() const
Definition: PndMCTrack.h:73
PndSdsDigi::GetIndex
Int_t GetIndex(int i=0) const
Definition: PndSdsDigi.h:63
PndTrackCand::GetSortedHit
PndTrackCandHit GetSortedHit(UInt_t i)
Definition: PndTrackCand.h:54
CAMath::Sin
static T Sin(const T &x)
Definition: PndCAMath.h:42
MCPoint
TFile * MCPoint
Definition: anaLmdDigi.C:25
PndMCTrack::GetMomentum
TVector3 GetMomentum() const
Definition: PndMCTrack.h:78
startEvent
Int_t startEvent
Definition: emc/run_linkCreation.C:7
c2
c2
Definition: plot_dirc.C:39
CAMath::Tan
float Tan(float x)
Definition: PndCAMath.h:165
storePath
TString storePath
Definition: runGeaneTestSimBox.cxx:50
CAMath::Cos
static T Cos(const T &x)
Definition: PndCAMath.h:43
DigiFile
TString DigiFile
Definition: emc/run_linkCreation.C:14
PndSdsCluster::GetDigiIndex
Int_t GetDigiIndex(Int_t i) const
Definition: PndSdsCluster.h:40
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
nEvents
Int_t nEvents
Definition: hit_dirc.C:11
timer
TStopwatch timer
Definition: hit_dirc.C:51
f
TFile * f
Definition: bump_analys.C:12
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
out
TFile * out
Definition: reco_muo.C:20
c1
c1
Definition: plot_dirc.C:35
PndSdsMCPoint::GetPosition
TVector3 GetPosition() const
Definition: PndSdsMCPoint.h:90
PndTrack::GetChi2
Double_t GetChi2() const
Definition: PndTrack.h:34
PndTrackCand::GetNHits
UInt_t GetNHits() const
Definition: PndTrackCand.h:59
PndSdsDigiPixel
Data class to store the digi output of a pixel module.
Definition: PndSdsDigiPixel.h:25
PndSdsHit::GetClusterIndex
Int_t GetClusterIndex() const
Definition: PndSdsHit.h:94
PndMCTrack::GetStartVertex
TVector3 GetStartVertex() const
Definition: PndMCTrack.h:76
PndMCTrack::GetMotherID
Int_t GetMotherID() const
Definition: PndMCTrack.h:74
PndTrackCandHit::GetHitId
Int_t GetHitId() const
Definition: PndTrackCandHit.h:55
Pi
Double_t Pi
Definition: PndGeoHypGeCrystal.cxx:10
main
int main(int __argc, char *__argv[])
Definition: simpleMCandRECmatch.C:66
Generated on Sun Aug 18 2019 03:02:03 for FairRoot/PandaRoot by   doxygen 1.8.8