FairRoot/PandaRoot
anaplaneclust.C
Go to the documentation of this file.
1 // root macro to analyze the clusterization output
2 {
3  int nEvents = 1000;
4  bool verbose = false;
5 
6  // ----- Load libraries ------------------------------------------------
7  gROOT->LoadMacro("$VMCWORKDIR/gconfig/basiclibs.C");
8  basiclibs();
9  gSystem->Load("libGeoBase");
10  gSystem->Load("libParBase");
11  gSystem->Load("libBase");
12  gSystem->Load("libPndData");
13  gSystem->Load("libField");
14  gSystem->Load("libGen");
15  gSystem->Load("libPassive");
16 
17  gSystem->Load("libgenfit");
18  gSystem->Load("libtrackrep");
19  gSystem->Load("libtpc");
20  gSystem->Load("libtpcreco");
21  gSystem->Load("librecotasks");
22 
23  gSystem->Load("libMvd");
24  gSystem->Load("libMvdReco");
25 
26  // ----- Timer --------------------------------------------------------
27  TStopwatch timer;
28  timer.Start();
29  // ------------------------------------------------------------------------
30 
31 
32  TFile* f = new TFile("../data/testMCCluster.root"); // the sim file you want to analyse
33  TTree *t=(TTree *) f->Get("pndsim") ;
34  TClonesArray* hit_array=new TClonesArray("MvdCluster");
35  t->SetBranchAddress("MVDCluster",&hit_array);//Branch names
36 
37  TFile* F = new TFile("testMC.root"); // the sim file you want to analyse
38  TTree *T=(TTree *) F->Get("pndsim") ;
39  TClonesArray* mc_array=new TClonesArray("MvdMCPoint");
40  T->SetBranchAddress("MVDPoint",&mc_array);//Branch names
41 
42  TGeoManager *geoMan = (TGeoManager*) gDirectory->Get("FAIRGeom");
43 
44  TH2D* hisxy = new TH2D("hisxy","",400,-15.,15.,400,-15.,15.);
45  hisxy->SetTitle("MVD MC Cluster, xy view;x / cm;y / cm");
46  TH2D* hisrz = new TH2D("hisrz","",400,-20.,20.,400,-15.,25.);
47  hisrz->SetTitle("MVD MC Cluster, rz view;z / cm;r/ cm");
48  TH2D* hisSensxy = new TH2D("hisSensxy","",400,-15.,15.,400,-15.,15.);
49  hisSensxy->SetTitle("MVD MC Sensor, xy view;x / cm;y / cm");
50  TH2D* hisSensrz = new TH2D("hisSensrz","",400,-20.,20.,400,-15.,25.);
51  hisSensrz->SetTitle("MVD MC Sensor, rz view;z / cm;r/ cm");
52 
53  TH2D* hisDiffXY = new TH2D("hisdiffxy","",100,-0.05,0.05,100,-0.05,0.05);
54  hisDiffXY->SetTitle("(MC - RECO) Hit coordinates xy view;#Deltax / cm;#Deltay / cm");
55  TH2D* hisDiffRZ = new TH2D("hisdiffrz","",100,-0.05,0.05,100,-0.00,0.10);
56  hisDiffRZ->SetTitle("(MC - RECO) Hit coordinates rz view;#Deltaz / cm;#Deltar / cm");
57 
58  TH1D* hisDiffX = new TH1D("hisdiffx","",100,-0.05,0.05);
59  hisDiffX->SetTitle("(MC - RECO) Hit coordinate x;#Deltax / cm;");
60  TH1D* hisDiffY = new TH1D("hisdiffy","",100,-0.05,0.05);
61  hisDiffY->SetTitle("(MC - RECO) Hit coordinate y;#Deltay / cm;");
62  TH1D* hisDiffZ = new TH1D("hisdiffz","",100,-0.05,0.05);
63  hisDiffZ->SetTitle("(MC - RECO) Hit coordinate z;#Deltaz / cm;");
64 
65  TH1D* hisDiffLocalX = new TH1D("hisdiffLocalx","",100,-0.05,0.05);
66  hisDiffLocalX->SetTitle("(MC - RECO) Local Hit coordinate x;#Deltax_{L} / cm;");
67  TH1D* hisDiffLocalY = new TH1D("hisdiffLocaly","",100,-0.05,0.05);
68  hisDiffLocalY->SetTitle("(MC - RECO) Local Hit coordinate y;#Deltay_{L} / cm;");
69  TH1D* hisDiffLocalZ = new TH1D("hisdiffLocalz","",100,-0.05,0.05);
70  hisDiffLocalZ->SetTitle("(MC - RECO) Local Hit coordinate z;#Deltaz_{L} / cm;");
71  TH2D* hisDiffLocalXY = new TH2D("hisdiffLocalxy","",100,-0.05,0.05,100,-0.05,0.05);
72  hisDiffLocalXY->SetTitle("(MC - RECO) Local Hit XY;#Deltax_{L} / cm;#Deltay_{L} / cm");
73  TH2D* hisDiffLocalRZ = new TH2D("hisdiffLocalrz","",100,-0.05,0.05,100,-0.00,0.10);
74  hisDiffLocalRZ->SetTitle("(MC - RECO) Local Hit RZ;#Deltaz_{L} / cm;#Deltar_{L} / cm");
75 
76  TH2D* hisLocalXY = new TH2D("hisLocalxy","",100,-5.,5.,100,-5.,5.);
77  hisLocalXY->SetTitle("(MC - RECO) Local Hit XY;x_{L} / cm;y_{L} / cm");
78  TH2D* hisLocalRZ = new TH2D("hisLocalrz","",100,-5.,5.,100,-0.,10.);
79  hisLocalRZ->SetTitle("(MC - RECO) Local Hit RZ;z_{L} / cm;r_{L} / cm");
80  TH2D* hisLocalYZ = new TH2D("hisLocalyz","",100,-5.,5.,100,-5.,5.);
81  hisLocalYZ->SetTitle("(MC - RECO) Local Hit YZ;z_{L} / cm;y_{L} / cm");
82  TH2D* hisLocalXZ = new TH2D("hisLocalxz","",100,-5.,5.,100,-5.,5.);
83  hisLocalXZ->SetTitle("(MC - RECO) Local Hit XZ;z_{L} / cm;x_{L} / cm");
84 
85  TH1D* hisPLUV = new TH1D("hispluv","sensor unit vectors",100,0.9999,1.0001);
86 
87  TVector3 vecs, vecmc, vecdiff;
88  Double_t tmpx,tmpy,tmpz;
89  TVector2 locals, localmc, localdiff;
90 // TGeoHMatrix* currentTransMat;
91  for (Int_t j=0; j<nEvents && j<t->GetEntriesFast(); j++)
92  {
93  t->GetEntry(j);
94  T->GetEntry(j);
95  if(verbose) cout<<"Event No "<<j<<endl;
96  for (Int_t i=0; i<hit_array->GetEntriesFast(); i++)
97  {
98  if(verbose) cout<<"Point No "<<i<<endl;
99  PndSdsCluster *hit=(PndSdsCluster*)hit_array->At(i);
100  PndSdsMCPoint *point=(PndSdsMCPoint*)mc_array->At(i);
101  int mcpdg = -1;
102 
103  geoMan->cd(point->GetDetName());
104 // currentTransMat = geoMan->GetCurrentMatrix();
105 
106  vecs.SetXYZ(hit->GetX(), hit->GetY(), hit->GetZ());
107  vecmc.SetXYZ(point->GetX(),point->GetY(),point->GetZ());
108  vecdiff = vecmc - vecs;
109 
110  Int_t layer = Int_t(10.*vecs->Mag());
111  if(verbose) cout<<vecs.x()<<"\t"<<vecs.y()<<"\t"<<vecs.z()<<endl;
112  hisxy->Fill(vecs.x(),vecs.y());
113  if(vecs.y() > 0.) hisrz->Fill(vecs.z(),vecs.Perp());
114  else hisrz->Fill(vecs.z(),-1.*vecs.Perp());
115 
116  hisDiffXY->Fill(vecdiff.x(),vecdiff.y());
117  hisDiffRZ->Fill(vecdiff.z(),vecdiff.Perp());
118  hisDiffX->Fill(vecdiff.x());
119  hisDiffY->Fill(vecdiff.y());
120  hisDiffZ->Fill(vecdiff.z());
121 
122  // --- Now move vectors to the local sensor system
123  DetPlane plane = hit->getDetPlane();
124  TVector3 vnorm = plane->getNormal(); //vnorm->SetMag(1.0);
125  TVector3 plO = plane->getO();
126  TVector3 plU = plane->getU(); //plU->SetMag(1.0);
127  TVector3 plV = plane->getV(); //plV->SetMag(1.0);
128  hisSensxy->Fill(plO.x(),plO.y());
129  if(plO.y()>0.)hisSensrz->Fill(plO.z(),plO.Perp());
130  else hisSensrz->Fill(plO.z(),-1.*plO.Perp());
131 
132  vecs -= plO;
133  tmpz = vecs.Dot(plU);
134  tmpx = vecs.Dot(plV);
135  tmpy = vecs.Dot(vnorm);
136  vecs.SetXYZ(tmpx,tmpy,tmpz);
137  hisLocalXY->Fill(vecs.x(),vecs.y());
138  hisLocalXZ->Fill(vecs.z(),vecs.x());
139  hisLocalYZ->Fill(vecs.z(),vecs.y());
140  hisLocalRZ->Fill(vecs.z(),vecs.Perp());
141 
142  vecmc -= plO;
143  tmpz = vecmc.Dot(plU);
144  tmpx = vecmc.Dot(plV);
145  tmpy = vecmc.Dot(vnorm);
146  vecmc.SetXYZ(tmpx,tmpy,tmpz);
147  vecdiff = vecmc - vecs;
148 
149  tmpx = plU.Dot(plU);tmpy = plU.Dot(plV);tmpz = plU.Dot(vnorm);
150  hisPLUV->Fill(tmpx);hisPLUV->Fill(tmpy);hisPLUV->Fill(tmpz);
151 
152  tmpx = plV.Dot(plU);tmpy = plV.Dot(plV);tmpz = plV.Dot(vnorm);
153  hisPLUV->Fill(tmpx);hisPLUV->Fill(tmpy);hisPLUV->Fill(tmpz);
154 
155  tmpx = vnorm.Dot(plU);tmpy = vnorm.Dot(plV);tmpz = vnorm.Dot(vnorm);
156  hisPLUV->Fill(tmpx);hisPLUV->Fill(tmpy);hisPLUV->Fill(tmpz);
157 
158  hisDiffLocalX->Fill(vecdiff.X());
159  hisDiffLocalY->Fill(vecdiff.Y());
160  hisDiffLocalZ->Fill(vecdiff.Z());
161  hisDiffLocalXY->Fill(vecdiff.X(),vecdiff.Y());
162  hisDiffLocalRZ->Fill(vecdiff.Z(),vecdiff.Perp());
163 
164  }//end for i (points in event)
165  }// end for j (events)
166 
167  gStyle->SetOptFit();
168 
169  Double_t par[6];
170 
171  TF1 *g1 = new TF1("g1","gaus",-0.007,0.007);
172  TF1 *g2 = new TF1("g2","gaus",-0.03,0.03);
173  TF1 *total = new TF1("total","gaus(0)+gaus(3)",-0.03,0.03);
174 
175  total->SetLineColor(kRed);
176  total->SetLineWidth(1);
177  g1->SetLineWidth(1);
178  g2->SetLineWidth(1);
179 
180  hisDiffX->Fit(g1,"R");
181  hisDiffX->Fit(g2,"R+");
182  g1->GetParameters(&par[0]);
183  g2->GetParameters(&par[3]);
184  total->SetParameters(par);
185  hisDiffX->Fit(total,"R");//"R+"
186 
187  hisDiffY->Fit(g1,"R");
188  hisDiffY->Fit(g2,"R+");
189  g1->GetParameters(&par[0]);
190  g2->GetParameters(&par[3]);
191  total->SetParameters(par);
192  hisDiffY->Fit(total,"R");//"R+"
193 
194  hisDiffZ->Fit(g1,"R");
195  hisDiffZ->Fit(g2,"R+");
196  g1->GetParameters(&par[0]);
197  g2->GetParameters(&par[3]);
198  total->SetParameters(par);
199  hisDiffZ->Fit(total,"R");//"R+"
200 
201  hisDiffLocalX->Fit("gaus");
202  hisDiffLocalX->GetFunction("gaus")->SetLineColor(kRed);
203  hisDiffLocalX->GetFunction("gaus")->SetLineWidth(1);
204 
205  hisDiffLocalY->Fit("gaus");
206  hisDiffLocalY->GetFunction("gaus")->SetLineColor(kRed);
207  hisDiffLocalY->GetFunction("gaus")->SetLineWidth(1);
208 
209  hisDiffLocalZ->Fit("gaus");
210  hisDiffLocalZ->GetFunction("gaus")->SetLineColor(kRed);
211  hisDiffLocalZ->GetFunction("gaus")->SetLineWidth(1);
212 
213 Int_t a = 5, b = 4;
214 
215 TCanvas* can1 = new TCanvas("MvdTestPlot","MCHit view in MVD",0,0,a*250,b*250);
216 can1->Divide(a,b);
217 
218 can1->cd(1);hisxy->DrawCopy("col");
219 can1->cd(2);hisrz->DrawCopy("col");
220 can1->cd(3);hisSensxy->DrawCopy("BOX");
221 can1->cd(4);hisSensrz->DrawCopy("BOX");
222 can1->cd(5);hisPLUV->DrawCopy();
223 
224 can1->cd(6);hisDiffX->DrawCopy("");
225 can1->cd(7);hisDiffY->DrawCopy("");
226 can1->cd(8);hisDiffZ->DrawCopy("");
227 can1->cd(9);hisDiffXY->DrawCopy("col");
228 can1->cd(10);hisDiffRZ->DrawCopy("col");
229 
230 can1->cd(11);hisDiffLocalX->DrawCopy("");
231 can1->cd(12);hisDiffLocalY->DrawCopy("");
232 can1->cd(13);hisDiffLocalZ->DrawCopy("");
233 can1->cd(14);hisDiffLocalXY->DrawCopy("col");
234 can1->cd(15);hisDiffLocalRZ->DrawCopy("col");
235 
236 can1->cd(16);hisLocalXY->DrawCopy("col");
237 can1->cd(17);hisLocalXZ->DrawCopy("col");
238 can1->cd(18);hisLocalYZ->DrawCopy("col");
239 can1->cd(19);hisLocalRZ->DrawCopy("col");
240 can1->cd(20);
241 
242 
243 can1->Print("testoutput.ps");
244 
245 
246  // ----- Finish -------------------------------------------------------
247  timer.Stop();
248  Double_t rtime = timer.RealTime();
249  Double_t ctime = timer.CpuTime();
250  cout << endl << endl;
251  cout << "Macro finished succesfully." << endl;
252  cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << endl;
253  cout << endl;
254  // ------------------------------------------------------------------------
255 
256 }
hisDiffXY
TH2D * hisDiffXY
Definition: anaLmdReco.C:57
basiclibs
basiclibs()
i
Int_t i
Definition: run_full.C:25
b
TTree * b
Definition: GammaSpectraAnalysis.C:24
PndSdsCluster
Class to store the Digis which belong to one cluster This class holds the information which Digi belo...
Definition: PndSdsCluster.h:19
hisDiffX
TH1F * hisDiffX
Definition: anaLmdReco.C:62
verbose
#define verbose
Definition: PndRadMapBoxMesh.cxx:15
geoMan
TGeoManager * geoMan
Definition: createRootGeoFileBarrel_2018v1.C:21
hisSensxy
TH2D * hisSensxy
Definition: anaplaneclust.C:48
hisDiffLocalZ
TH1D * hisDiffLocalZ
Definition: anaplaneclust.C:69
par
Double_t par[3]
Definition: anaLmdCluster.C:146
localdiff
TVector2 localdiff
Definition: anaLmdDigi.C:66
hisDiffLocalXY
TH2D * hisDiffLocalXY
Definition: anaplaneclust.C:71
mc_array
TClonesArray * mc_array
Definition: anaLmdCluster.C:25
hisDiffZ
TH1F * hisDiffZ
Definition: anaLmdReco.C:66
g2
TF1 * g2
Definition: anaplaneclust.C:172
hisDiffLocalX
TH1D * hisDiffLocalX
Definition: anaplaneclust.C:65
g1
TF1 * g1
Definition: anaplaneclust.C:171
hit_array
TClonesArray * hit_array
Definition: anal_hit_digi_cluster_fwendcap.C:21
hisLocalRZ
TH2D * hisLocalRZ
Definition: anaplaneclust.C:78
vecs
TVector3 vecs
Definition: GammaSpectraAnalysis.C:43
hisDiffRZ
TH2D * hisDiffRZ
Definition: anaLmdReco.C:59
T
TTree * T
Definition: anaLmdReco.C:32
a
Int_t a
Definition: anaLmdDigi.C:126
localmc
TVector2 localmc
Definition: anaLmdDigi.C:66
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
locals
TVector2 locals
Definition: anaLmdDigi.C:66
F
TFile * F
Definition: anaLmdReco.C:31
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
hisLocalXZ
TH2D * hisLocalXZ
Definition: anaplaneclust.C:82
total
TF1 * total
Definition: anaLmdCluster.C:147
vecmc
TVector3 vecmc
Definition: anaLmdCluster.C:52
hisxy
TH2D * hisxy
Definition: anaLmdDigi.C:38
hisLocalYZ
TH2D * hisLocalYZ
Definition: anaplaneclust.C:80
hisDiffLocalRZ
TH2D * hisDiffLocalRZ
Definition: anaplaneclust.C:73
mcpdg
int mcpdg
Definition: GammaSpectraAnalysis.C:44
ctime
Double_t ctime
Definition: hit_dirc.C:114
hisrz
TH2D * hisrz
Definition: anaLmdDigi.C:41
hisPLUV
TH1D * hisPLUV
Definition: anaplaneclust.C:85
tmpx
Double_t tmpx
Definition: anaLmdDigi.C:65
layer
Int_t layer
Definition: reco_muo.C:36
t
TTree * t
Definition: bump_analys.C:13
hisSensrz
TH2D * hisSensrz
Definition: anaplaneclust.C:50
hisDiffY
TH1F * hisDiffY
Definition: anaLmdReco.C:64
hit
PndSdsMCPoint * hit
Definition: anasim.C:70
can1
TCanvas * can1
Definition: anaLmdCluster.C:151
tmpy
Double_t tmpy
Definition: anaLmdDigi.C:65
rtime
Double_t rtime
Definition: hit_dirc.C:113
vecdiff
TVector3 vecdiff
Definition: anaLmdReco.C:77
hisLocalXY
TH2D * hisLocalXY
Definition: anaplaneclust.C:76
hisDiffLocalY
TH1D * hisDiffLocalY
Definition: anaplaneclust.C:67
tmpz
Double_t tmpz
Definition: anaLmdDigi.C:65
point
PndSdsMCPoint * point
Definition: anaLmdCluster.C:72
Generated on Wed Apr 3 2019 08:43:49 for FairRoot/PandaRoot by   doxygen 1.8.5