FairRoot/PandaRoot
CreateHitsForAlignmentPixel.C
Go to the documentation of this file.
1 //################################################################
2 //# Macros to write hits point for alignment with Knossos
3 //# author: Anastasia Karavdina
4 //# date: Feb, 2013
5 //#
6 //# to compile it add in "# install #" part of CMakeLists.txt lines:
7 //# add_executable(hits_align_pixel CreateHitsForAlignmentPixel.C)
8 //# target_link_libraries(hits_align_pixel ${ROOT_LIBRARIES} genfit genfitAdapters trackrep RecoHits Stt Fts Lmd LmdReco LmdTrk Sds SdsReco Geane Geom Proof TrkBase FairTools PndData VMC EG GeomPainter generalTools Base GeoBase FairDB ParBase MathMore Minuit)
9 //#
10 //# to run it (and see options):
11 //# ${PANDAROOT}/build/bin/./hits_align_pixel --help
12 //################################################################
13 #include <iostream>
14 #include <sstream>
15 
16 #include<TApplication.h>
17 #include<TCanvas.h>
18 #include<TROOT.h>
19 #include<TString.h>
20 #include<TChain.h>
21 #include<TFile.h>
22 #include<TClonesArray.h>
23 #include<TSystem.h>
24 #include<TH1.h>
25 #include<TH2.h>
26 #include<TF1.h>
27 #include<TRotation.h>
28 #include<TVector3.h>
29 #include<TMath.h>
30 #include<TGaxis.h>
31 #include<TNtuple.h>
32 #include<TMatrixD.h>
33 #include<TStopwatch.h>
34 #include<TNtuple.h>
35 #include<PndMCTrack.h>
36 #include<PndSdsMCPoint.h>
37 
38 //lmd track
39 #include<PndLinTrack.h>
40 #include<PndTrackCand.h>
41 #include<PndTrackCandHit.h>
42 #include<PndSdsHit.h>
43 #include<PndSdsClusterPixel.h>
44 #include<PndSdsDigiPixel.h>
45 #include<PndTrack.h>
46 
47 // needed for geane backtracking
48 #include<FairRunAna.h>
49 #include<FairRootManager.h>
50 #include<FairGeane.h>
51 #include<FairRtdbRun.h>
52 #include<FairRuntimeDb.h>
53 #include<FairParRootFileIo.h>
54 #include<FairLogger.h>
55 #include<FairTrackParH.h>
56 
57 // needed for corditate from\out lumi frame calculation
58 #include<PndLmdDim.h>
59 
60 using namespace std;
61 
62 int main(int __argc,char *__argv[]) {
63  //TODO: read this like params!
64  // const int nEvents=500000;
65  ofstream output;
66  int nEvents=1000;
67  int startEvent=0;
68  TString storePath="/data/FAIRsorf/pandaroot/trunk/macro/lmd/tmpOutput";
69  double Plab=15.;
70  int verboseLevel=0;
71  int sectorPos=0;
72  std::string startStr="", momStr="", nStr="", pathStr="", verbStr="", outnStr="", outrootStr="",sectorStr="";
73  // decode arguments
74  if( __argc>1 && ( strcmp( __argv[1], "-help" ) == 0
75  || strcmp( __argv[1], "--help" ) == 0 ) ){
76 
77  std::cout << "This is script for alignment data creation with pixel sensors design \n"
78  << "with parameters\n"
79  <<"-s start event \n"
80  <<"-n Number of events \n"
81  <<"-path path to the file(s) \n"
82  <<"-out name for output file \n"
83  <<"-outhist name for output file with residuals histogram \n"
84  <<"-v verbose Level (if>0, print out some information) \n"
85  <<"-m sector#"
86  <<"Have fun! \n"
87  << std::endl;
88  return 0;
89  }
90  while ((optind < (__argc-1) ) && (__argv[optind][0]=='-')) {
91  bool found=false;
92  std::string sw = __argv[optind];
93  if (sw=="-s") {
94  optind++;
95  startStr = __argv[optind];
96  found=true;
97  }
98  if (sw=="-m") {
99  optind++;
100  sectorStr = __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(sw=="-out"){
124  optind++;
125  outnStr= __argv[optind];
126  found=true;
127  }
128  if(sw=="-outhist"){
129  optind++;
130  outrootStr= __argv[optind];
131  found=true;
132  }
133  if (!found){
134  std::cout<< "Unknown switch: "
135  << __argv[optind] <<std::endl;
136  optind++;
137  }
138 
139  while ( (optind < __argc ) && __argv[optind][0]!='-' ) optind++;
140  }
141 
142  std::stringstream startSStr(startStr), nSStr(nStr), pathSStr(pathStr), verbSStr(verbStr),sectorSStr(sectorStr);
143 
144  startSStr >> startEvent;
145  nSStr >> nEvents;
146  pathSStr >> storePath;
147  verbSStr >> verboseLevel;
148  sectorSStr>> sectorPos;
149  cout<<"For data files will be used path: "<<storePath<<endl;
150 
151  // ---- Open output file ------------------------------------------------
152  output.open(outnStr.c_str());
153  TFile *f = new TFile(outrootStr.c_str(),"RECREATE");
154  // ---- Load libraries -------------------------------------------------
155  gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.C");
156  gSystem->Load("libLmdTrk");
157  // ------------------------------------------------------------------------
158 
159  // ----- Timer --------------------------------------------------------
160  TStopwatch timer;
161  timer.Start();
162  // ------------------------------------------------------------------------
163 
164  // ---- Input files --------------------------------------------------------
165 
166  TString DigiFile = storePath+"/Lumi_digi_";
167  DigiFile += startEvent;
168  DigiFile += ".root";
169  TChain tdigiHits("pndsim");
170  tdigiHits.Add(DigiFile);
171 
172  TString recHit=storePath+"/Lumi_reco_";
173  recHit += startEvent;
174  recHit += ".root";
175  TChain tHits("pndsim");
176  tHits.Add(recHit);
177 
178  TString recHitmerged=storePath+"/Lumi_recoMerged_";
179  recHitmerged += startEvent;
180  recHitmerged += ".root";
181  TChain tHitsMerged("pndsim");
182  tHitsMerged.Add(recHitmerged);
183 
184  // TString recHitmerged=storePath+"/Lumi_reco_";
185  // recHitmerged += startEvent;
186  // recHitmerged += ".root";
187  // TChain tHitsMerged("pndsim");
188  // tHitsMerged.Add(recHitmerged);
189 
190 
191  TString trkCand = storePath+"/Lumi_TCand_";
192  trkCand += startEvent;
193  trkCand += ".root";
194  TChain tTrkCand("pndsim");
195  tTrkCand.Add(trkCand);
196 
197  TString recTrack = storePath+"/Lumi_Track_"; //NOTE: make sure kinematic filter was switched off, but filtering on hit base was done!
198  //TString recTrack = storePath+"/Lumi_TrackNotFiltered_";
199  recTrack += startEvent;
200  recTrack += ".root";
201  TChain tTrkRec("pndsim");
202  tTrkRec.Add(recTrack);
203 
204  TString simMC=storePath+"/Lumi_MC_";
205  simMC += startEvent;
206  simMC += ".root";
207  TChain tMC("pndsim");
208  tMC.Add(simMC);
209 
210  // ---------------------------------------------------------------------------------
211 
212  //--- Digitization info ------------------------------------------------------------
213  TClonesArray* fStripClusterArray;
214  fStripClusterArray = new TClonesArray("PndSdsClusterPixel");
215  tHits.SetBranchAddress("LMDPixelClusterCand",&fStripClusterArray);
216  TClonesArray* fStripDigiArray;
217  fStripDigiArray = new TClonesArray("PndSdsDigiPixel");
218  tdigiHits.SetBranchAddress("LMDPixelDigis",&fStripDigiArray);
219  //----------------------------------------------------------------------------------
220 
221  //--- Real Hits --------------------------------------------------------------------
222  TClonesArray* rechit_array=new TClonesArray("PndSdsMergedHit");
223  tHitsMerged.SetBranchAddress("LMDHitsMerged",&rechit_array); //Points for Tracks
224  // TClonesArray* rechit_array=new TClonesArray("PndSdsHit");
225  // tHitsMerged.SetBranchAddress("LMDHits",&rechit_array); //Points for Tracks
226  //----------------------------------------------------------------------------------
227 
228  //--- Track Candidate ---------------------------------------------------------------
229  TClonesArray* trkcand_array=new TClonesArray("PndTrackCand");
230  tTrkCand.SetBranchAddress("LMDTrackCand",&trkcand_array); //Points for Track Canidates
231  //-----------------------------------------------------------------------------------
232 
233  //--- Real tracks -------------------------------------------------------------------
234  // TClonesArray* rec_trk=new TClonesArray("PndLinTrack");
235  // tTrkRec.SetBranchAddress("LMDTrack",&rec_trk); //Tracks
236  TClonesArray* rec_trk=new TClonesArray("PndTrack");
237  tTrkRec.SetBranchAddress("LMDPndTrack",&rec_trk); //Tracks
238  //----------------------------------------------------------------------------------
239 
240  //--- MC info -----------------------------------------------------------------
241  TClonesArray* true_tracks=new TClonesArray("PndMCTrack");
242  tMC.SetBranchAddress("MCTrack",&true_tracks); //True Track to compare
243 
244  TClonesArray* true_points=new TClonesArray("PndSdsMCPoint");
245  tMC.SetBranchAddress("LMDPoint",&true_points); //True Points to compare
246  //----------------------------------------------------------------------------------
247 
248  //Histograms with residuals---------------------------------------------------------
249  //Create and fill histogramms
250  TObjArray *m_res_x = new TObjArray;
251  TObjArray *m_res_y = new TObjArray;
252  TObjArray *m_res_r = new TObjArray;
253 
254  TObjArray *m_res_mc_x = new TObjArray;
255  TObjArray *m_res_mc_y = new TObjArray;
256  TObjArray *m_res_mc_z = new TObjArray;
257 
258  double misal_scales[3];
259  const double scaleX = 0.05;
260  const double scaleY = 0.05;
261  const double scaleR = 0.05;
262 
263 
264 
265 
266  (scaleX == 0) ? misal_scales[0] = 1. : misal_scales[0] = 1.*scaleX;
267  (scaleY == 0) ? misal_scales[1] = 1. : misal_scales[1] = 1.*scaleY;
268  (scaleR == 0) ? misal_scales[2] = 1. : misal_scales[2] = 1.*scaleR;
269 
270  // TH2F *resx = new TH2F("resx","residuals vs. x; x, cm;#delta_{x}, cm",100,-15,15,5000,-2.5e-1,2.5e-1);
271  // TH2F *resy = new TH2F("resy","residuals vs. y; y, cm;#delta_{y}, cm",100,-15,15,5000,-2.5e-1,2.5e-1);
272 
273  TH2F *resx = new TH2F("resx","residuals vs. x; x, cm;#delta_{x}, cm",100,0,45,100,-5.0e-2,5.0e-2);
274  TH2F *resy = new TH2F("resy","residuals vs. y; y, cm;#delta_{y}, cm",100,-15,15,100,-5.0e-2,5.0e-2);
275 
276  TH2F *x_id = new TH2F("x_id","x vs. sensor id; id ;x, cm",400,0,400,1e3,0,45);
277  TH2F *y_id = new TH2F("y_id","y vs. sensor id; id ;y, cm",400,0,400,1e3,-15,15);
278  TH2F *z_id = new TH2F("z_id","z vs. sensor id; id ;z, cm",400,0,400,1e3,-10,60);
279  TNtuple *nhits = new TNtuple("nhits","hitsInfo","xrec:yrec:zrec:errxrec:erryrec:errzrec:xmc:ymc:zmc:half:plane:module:side");
280  // sens_id = mcpoint->GetSensorID(); // store the sensor id
281  // // calculate the plane and sensor on this plane
282  // lmddim->Get_sensor_by_id(sens_id, ihalf, iplane, imodule, iside, idie, isensor);
283 
284  TH2F *resxy = new TH2F("resxy","residuals y vs. residuals x; #delta_{x}, cm; #delta_{y}, cm",5000,-2.5e-1,2.5e-1,5000,-2.5e-1,2.5e-1);
285  TH1I *htrks = new TH1I("trks","number of trks/ev",2e1,0,2e1);
286  // for (unsigned int histID = 0; histID < 400; ++histID)
287  // {
288  // char histoResXName[100];
289  // char histoResYName[100];
290  // char histoResRName[100];
291 
292  // char histoResMCXName[100];
293  // char histoResMCYName[100];
294  // char histoResMCZName[100];
295  // sprintf(histoResXName, "residuals_x_%d",histID);
296  // TH1F* aHisto = new TH1F(histoResXName, histoResXName, 1000, -misal_scales[0], misal_scales[0]);
297  // m_res_x->Add(aHisto);
298 
299  // sprintf(histoResYName, "residuals_y_%d",histID);
300  // TH1F* aHisto2 = new TH1F(histoResYName, histoResYName, 1000, -misal_scales[1], misal_scales[1]);
301 
302  // m_res_y->Add(aHisto2);
303 
304  // sprintf(histoResRName, "residuals_r_%d",histID);
305  // TH1F* aHisto3 = new TH1F(histoResRName, histoResRName, 1000, -misal_scales[2], misal_scales[2]);
306  // m_res_r->Add(aHisto3);
307 
308  // sprintf(histoResMCXName, "residuals_mc_x_%d",histID);
309  // TH1F* aHistomc = new TH1F(histoResMCXName, histoResMCXName, 1000, -misal_scales[0], misal_scales[0]);
310  // m_res_mc_x->Add(aHistomc);
311 
312  // sprintf(histoResMCYName, "residuals_mc_y_%d",histID);
313  // TH1F* aHistomc2 = new TH1F(histoResMCYName, histoResMCYName, 1000, -misal_scales[1], misal_scales[1]);
314 
315  // m_res_mc_y->Add(aHistomc2);
316 
317  // sprintf(histoResMCZName, "residuals_mc_z_%d",histID);
318  // TH1F* aHistomc3 = new TH1F(histoResMCZName, histoResMCZName, 1000, -misal_scales[2], misal_scales[2]);
319  // m_res_mc_z->Add(aHistomc3);
320  // }
321  //----------------------------------------------------------------------------------
322 
323  //Load lumi geo params
324  PndLmdDim *lmddim = PndLmdDim::Instance();
325  // lmddim -> Read_transformation_matrices("matrices.txt", true);
326 
327  // TString mtx_perfect = "${VMCWORKDIR}/macro/lmd/matrices_perfect.txt";
328 
329  TString mtx_corr = storePath+"/matrices_corrected.txt";
330  TString mtx_perfect = storePath+"/trafo_matrices_lmd.dat";
331  // lmddim -> Read_transformation_matrices(mtx_perfect.Data(), false);
332  lmddim -> Read_transformation_matrices(mtx_perfect.Data(), false);
333  lmddim -> Read_transformation_matrices(mtx_corr.Data(), true);
334  //
335  // lmddim -> Read_transformation_matrices("/panda/pandaroot/macro/lmd/matrices_corrected.txt", true);
336  int icounter=0;
337  cout<<"In total there are "<<nEvents<<" events"<<endl;
338  int icounterALL=0;
339  for (Int_t j=0; j<nEvents; j++){
340  // Read REC tree -----------------------------------------------------------------
341  tTrkCand.GetEntry(j);
342  tTrkRec.GetEntry(j);
343  tHits.GetEntry(j);
344  tHitsMerged.GetEntry(j);
345  tdigiHits.GetEntry(j);
346  tMC.GetEntry(j);
348 
349  const int nRecTrks = rec_trk->GetEntriesFast();
350  icounterALL +=nRecTrks;
351  // if(nRecTrks>1) continue; //!!! TEST with 1 track/event only !!!
353  for (Int_t iN=0; iN<nRecTrks; iN++){
354  // PndLinTrack *trk_lin = (PndLinTrack*)rec_trk->At(iN);
355  // TVector3 startlintrk = trk_lin->GetStartVec();
356  // TVector3 dirlintrk = trk_lin->GetDirectionVec();
357  // Int_t candID = trk_lin->GetTCandID();
358  PndTrack *trkpnd = (PndTrack*)rec_trk->At(iN);
359  FairTrackParP fFittedTrkP = trkpnd->GetParamFirst();
360  TVector3 startlintrk(fFittedTrkP.GetX(),fFittedTrkP.GetY(),fFittedTrkP.GetZ());
361  TVector3 dirlintrk(fFittedTrkP.GetPx(),fFittedTrkP.GetPy(),fFittedTrkP.GetPz());
362  // double covMARS[6][6];
363  // fFittedTrkP.GetMARSCov(covMARS);
364  // TVector3 errMomRecLMD(sqrt(covMARS[0][0]),sqrt(covMARS[1][1]),sqrt(covMARS[2][2]));
365  // TVector3 errPosRecLMD(sqrt(covMARS[3][3]),sqrt(covMARS[4][4]),sqrt(covMARS[5][5]));
366  dirlintrk *= 1./dirlintrk.Mag();
367  int candID = trkpnd->GetRefIndex();
368  // cout<<"candID = "<<candID<<endl;
369  PndTrackCand *trkcand = (PndTrackCand*)trkcand_array->At(candID);
370  const int Ntrkcandhits= trkcand->GetNHits();
371  // cout<<"Ntrkcandhits = "<<Ntrkcandhits<<endl;
372  if(Ntrkcandhits<4) continue;
373  // if(Ntrkcandhits<3) continue; //!!! TEST with > 2 hits tracks only !!!
374  // if(Ntrkcandhits<2) continue; //!!! TEST
375  //if(Ntrkcandhits>4) continue; //!!! TEST with single hits only !!!
376  double phiMCgl;
377  //check if these hits are sutiable for sector aligment
378  int trkModules[Ntrkcandhits];
379  int trkHalfs[Ntrkcandhits];
380  int half1;
381  int module1;
382  for (Int_t iHit = 0; iHit < Ntrkcandhits; iHit++){
383  PndTrackCandHit candhit = (PndTrackCandHit)(trkcand->GetSortedHit(iHit));
384  Int_t hitID = candhit.GetHitId();
385  PndSdsHit* myHit = (PndSdsHit*)(rechit_array->At(hitID));
386  int sensorID = myHit->GetSensorID();
387  int ihalf, iplane, imodule, iside, idie, isensor;
388  // calculate the plane and sensor on this plane
389  lmddim->Get_sensor_by_id(sensorID, ihalf, iplane, imodule, iside, idie, isensor);
390  trkModules[iHit]=ihalf*5+imodule;
391  trkHalfs[iHit]=ihalf;
392  half1=ihalf;
393  module1=ihalf*5+imodule;
394  }
395  bool flagSector = true;
396  if(sectorPos<10){
397  for (Int_t iHit = 0; iHit < Ntrkcandhits; iHit++){
398  if(trkModules[iHit]!=sectorPos) flagSector=false;// TODO: check on diff trkModules[iHit] for alignment between sectors
399  }
400  }
401  else{
402  for (Int_t iHit = 0; iHit < Ntrkcandhits; iHit++){
403  if(trkHalfs[iHit]!=half1) flagSector=false;// TODO: check on diff trkModules[iHit] for alignment between sectors
404  if(trkModules[iHit]!=module1) flagSector=false;// TODO: check on diff trkModules[iHit] for alignment between sectors
405  }
406  }
407 
408  //(end) check if these hits are sutiable for sector aligment
409  if(!flagSector) continue;
410  icounter++;
411  if(iN==0) htrks->Fill(nRecTrks);//fill only if trk was accepted
412  for (Int_t iHit = 0; iHit < Ntrkcandhits; iHit++){
413  PndTrackCandHit candhit = (PndTrackCandHit)(trkcand->GetSortedHit(iHit));
414  Int_t hitID = candhit.GetHitId();
415  PndSdsHit* myHit = (PndSdsHit*)(rechit_array->At(hitID));
416  TVector3 HitPos = myHit->GetPosition();
417  int sensorID = myHit->GetSensorID();
418  int ihalf, iplane, imodule, iside, idie, isensor;
419  // calculate the plane and sensor on this plane
420  lmddim->Get_sensor_by_id(sensorID, ihalf, iplane, imodule, iside, idie, isensor);
421 
422  TVector3 HitPosLoc(lmddim->Transform_global_to_lmd_local(HitPos,false,true));
423  HitPosLoc = TVector3(lmddim->Transform_lmd_local_to_module_side(HitPosLoc,ihalf,0,imodule,0, false,true));
424  // TVector3 HitPosLoc = lmddim->Transform_global_to_sensor(HitPos, ihalf, iplane, imodule, iside, idie, isensor,false,true);
425  // HitPosLoc = lmddim->Transform_sensor_to_module_side(HitPosLoc,ihalf,0,imodule,0,idie, isensor,false,false);
426  // cout<<"ihalf, iplane, imodule, iside: "<<ihalf<<","<<iplane<<","<<imodule<<","<<iside<<endl;
427  TVector3 mcTop,mcTopOUT;
428  PndSdsClusterPixel* myCluster = (PndSdsClusterPixel*)(fStripClusterArray->At(myHit->GetClusterIndex()));
429  PndSdsDigiPixel* astripdigi = (PndSdsDigiPixel*)(fStripDigiArray->At(myCluster->GetDigiIndex(0)));
430  PndSdsMCPoint* MCPoint = (PndSdsMCPoint*)(true_points->At(astripdigi->GetIndex(0)));
431  // TVector3 MCtrue = MCPoint->GetPosition();
432  // MCtrue = lmddim->Transform_global_to_lmd_local(MCtrue,false,false);
433  // MCtrue = lmddim->Transform_lmd_local_to_module_side(MCtrue,ihalf,0,imodule,0, false,false);
434  // TVector3 recMCdiff = HitPosLoc - MCtrue;
435  // cout<<"rec - mc :"<<endl;
436  // recMCdiff.Print();
437 
438 
439  mcTop = MCPoint->GetPosition();
440  mcTopOUT = MCPoint->GetPositionOut();
441  // mcTop = TVector3(lmddim->Transform_global_to_lmd_local(mcTop, false, true));
442  // mcTopOUT = TVector3(lmddim->Transform_global_to_lmd_local(mcTopOUT, false, true));
443  TVector3 mcMid = mcTop+(mcTopOUT-mcTop)*0.5;
444 
445 
446 
447  // cout<<"ResRecMC: ";
448  // ResRecMC.Print();
449  // cout<<""<<endl;
450  double trkX = startlintrk.X()+dirlintrk.X()*(HitPos.Z()-startlintrk.Z());
451  double trkY = startlintrk.Y()+dirlintrk.Y()*(HitPos.Z()-startlintrk.Z());
452  double trkR = sqrt(trkX*trkX+trkY*trkY);
453  double hitR = sqrt(HitPos.X()*HitPos.X()+HitPos.Y()*HitPos.Y());
454 
455  // ((TH1F*)m_res_x->At(sensorID))->Fill(HitPos.X()-trkX);
456  // ((TH1F*)m_res_y->At(sensorID))->Fill(HitPos.Y()-trkY);
457  // ((TH1F*)m_res_r->At(sensorID))->Fill(hitR-trkR);
458  double resx_d = HitPos.X()-trkX;
459  double resy_d = HitPos.Y()-trkY;
460  resx->Fill(HitPos.X(),resx_d);
461  resy->Fill(HitPos.Y(),resy_d);
462  resxy->Fill(resx_d,resy_d);
463  // if(fabs(resx_d)>0.1 || fabs(resy_d)>0.1)
464  // cout<<"Track #"<<iN<<" from event #"<<j<<" has: resx="<<resx_d<<", resy="<<resy_d<<endl;
465 
466  x_id->Fill(sensorID,HitPos.X());
467  y_id->Fill(sensorID,HitPos.Y());
468  z_id->Fill(sensorID,HitPos.Z());
469 
470 
471  // int glModule = ihalf*4*5+(iplane)*5+imodule;// counted per one half
472  int glModule = iplane;// for alignment inside one sector
473  if(sectorPos>10) glModule = ihalf*4*5+(iplane)*5+imodule;// counted per one half
474  TMatrixD HitErr = myHit->GetCov();
475  // cout<<"HitErr:"<<endl;
476  // HitErr.Print();
477  TMatrixD HitErrLoc(lmddim->Transform_global_to_lmd_local(HitErr,true));
478  // cout<<"HitErrLoc:"<<endl;
479  // HitErrLoc.Print();
480  HitErrLoc = (lmddim->Transform_lmd_local_to_module_side(HitErrLoc,ihalf,0,imodule,0,true));
481  // cout<<"HitErrLoc:"<<endl;
482  // HitErrLoc.Print();
483  // TMatrixD HitErrLoc = lmddim->Transform_global_to_sensor(HitErr, ihalf, iplane, imodule, iside, idie, isensor,true);
484  // HitErrLoc = lmddim->Transform_sensor_to_module_side(HitErrLoc,ihalf,0,imodule,0,idie, isensor,false);
485 
486  // HitErrLoc = lmddim->Transform_lmd_local_to_module_side(HitErrLoc,ihalf,iplane,imodule,0,true);//TEST
487  TVector3 MCHitPosLoc(lmddim->Transform_global_to_lmd_local(mcMid,false,true));
488  MCHitPosLoc = TVector3(lmddim->Transform_lmd_local_to_module_side(MCHitPosLoc,ihalf,0,imodule,0, false,true));
489  nhits->Fill(HitPosLoc.X(),HitPosLoc.Y(),HitPosLoc.Z(),sqrt(HitErrLoc(0,0)),sqrt(HitErrLoc(1,1)),sqrt(HitErrLoc(2,2)),MCHitPosLoc.X(),MCHitPosLoc.Y(),MCHitPosLoc.Z(),ihalf,iplane,imodule,iside);
490  TVector3 ResRecMC = HitPosLoc - MCHitPosLoc;
491  // ((TH1F*)m_res_mc_x->At(sensorID))->Fill(ResRecMC.X());
492  // ((TH1F*)m_res_mc_y->At(sensorID))->Fill(ResRecMC.Y());
493  // ((TH1F*)m_res_mc_z->At(sensorID))->Fill(ResRecMC.Z());
494 
495  int endtrk=1;
496  if((Ntrkcandhits-iHit)>1) endtrk=0;
497  if(flagSector && output>0){
498  output<<HitPosLoc.X()<<" "<<sqrt(HitErrLoc(0,0))<<" ";
499  output<<HitPosLoc.Y()<<" "<<sqrt(HitErrLoc(1,1))<<" ";
500  // output<<HitPosLoc.Z()<<" "<<sqrt(HitErrLoc(2,2))<<" ";
501  output<<HitPosLoc.Z()<<" "<<sectorPos<<" "; //Z and sectorPos isn't used in Knossos
502  output<<glModule<<" "<<endtrk<<endl;
503  }
504  }
505  output<<" "<<endl;
506  }
507  }
509  output.close();
510  cout<<"This module has "<<icounter<<" trks from "<<icounterALL<<endl;
511  // for (unsigned int histID = 0; histID < 400; ++histID){
512  // ((TH1F*)m_res_x->At(histID))->Write();
513  // ((TH1F*)m_res_y->At(histID))->Write();
514  // ((TH1F*)m_res_r->At(histID))->Write();
515  // ((TH1F*)m_res_mc_x->At(histID))->Write();
516  // ((TH1F*)m_res_mc_y->At(histID))->Write();
517  // ((TH1F*)m_res_mc_z->At(histID))->Write();
518  // }
519 
520 
521 // TCanvas* c0 = new TCanvas("respl1x","residuals, pl#1",200,500,700,800);
522 // c0->SetFillColor(0);
523 // c0->SetBorderMode(0);
524 // c0->Divide(10,10);
525 // for (unsigned int histID = 0; histID < 100; ++histID){
526 // c0->cd(histID+1);
527 // ((TH1F*)m_res_x->At(histID))->GetXaxis()->SetTitle("(trk-hit)[x], cm");
528 // ((TH1F*)m_res_x->At(histID))->GetXaxis()->SetTitleSize(0.05);
529 // ((TH1F*)m_res_x->At(histID))->SetLineColor(1);
530 // ((TH1F*)m_res_x->At(histID))->Draw();
531 // }
532 // c0->Update();
533 // c0->Write();
534 
535 // TCanvas* c1 = new TCanvas("respl2x","residuals, pl#2",200,500,700,800);
536 // c1->SetFillColor(0);
537 // c1->SetBorderMode(0);
538 // c1->Divide(10,10);
539 // for (unsigned int histID = 100; histID < 200; ++histID){
540 // c1->cd(histID+1-100);
541 // ((TH1F*)m_res_x->At(histID))->GetXaxis()->SetTitle("(trk-hit)[x], cm");
542 // ((TH1F*)m_res_x->At(histID))->GetXaxis()->SetTitleSize(0.05);
543 // ((TH1F*)m_res_x->At(histID))->SetLineColor(1);
544 // ((TH1F*)m_res_x->At(histID))->Draw();
545 // }
546 // c1->Update();
547 // c1->Write();
548 
549 // TCanvas* c2 = new TCanvas("respl3x","residuals, pl#3",200,500,700,800);
550 // c2->SetFillColor(0);
551 // c2->SetBorderMode(0);
552 // c2->Divide(10,10);
553 // for (unsigned int histID = 200; histID < 300; ++histID){
554 // c2->cd(histID+1-200);
555 // ((TH1F*)m_res_x->At(histID))->GetXaxis()->SetTitle("(trk-hit)[x], cm");
556 // ((TH1F*)m_res_x->At(histID))->GetXaxis()->SetTitleSize(0.05);
557 // ((TH1F*)m_res_x->At(histID))->SetLineColor(1);
558 // ((TH1F*)m_res_x->At(histID))->Draw();
559 // }
560 // c2->Update();
561 // c2->Write();
562 
563 // TCanvas* c3 = new TCanvas("respl4x","residuals, pl#4",200,500,700,800);
564 // c3->SetFillColor(0);
565 // c3->SetBorderMode(0);
566 // c3->Divide(10,10);
567 // for (unsigned int histID = 300; histID < 400; ++histID){
568 // c3->cd(histID+1-300);
569 // ((TH1F*)m_res_x->At(histID))->GetXaxis()->SetTitle("(trk-hit)[x], cm");
570 // ((TH1F*)m_res_x->At(histID))->GetXaxis()->SetTitleSize(0.05);
571 // ((TH1F*)m_res_x->At(histID))->SetLineColor(1);
572 // ((TH1F*)m_res_x->At(histID))->Draw();
573 // }
574 // c3->Update();
575 // c3->Write();
576 
577 
578 // TCanvas* c01 = new TCanvas("respl1y","residuals, pl#1",200,500,700,800);
579 // c01->SetFillColor(0);
580 // c01->SetBorderMode(0);
581 // c01->Divide(10,10);
582 // for (unsigned int histID = 0; histID < 100; ++histID){
583 // c01->cd(histID+1);
584 // ((TH1F*)m_res_y->At(histID))->GetXaxis()->SetTitle("(trk-hit)[y], cm");
585 // ((TH1F*)m_res_y->At(histID))->GetXaxis()->SetTitleSize(0.05);
586 // ((TH1F*)m_res_y->At(histID))->SetLineColor(1);
587 // ((TH1F*)m_res_y->At(histID))->Draw();
588 // }
589 // c01->Update();
590 // c01->Write();
591 
592 // TCanvas* c11 = new TCanvas("respl2y","residuals, pl#2",200,500,700,800);
593 // c11->SetFillColor(0);
594 // c11->SetBorderMode(0);
595 // c11->Divide(10,10);
596 // int i11=1;
597 // for (unsigned int histID = 100; histID < 200; ++histID){
598 // c11->cd(i11);
599 // i11++;
600 // ((TH1F*)m_res_y->At(histID))->GetXaxis()->SetTitle("(trk-hit)[y], cm");
601 // ((TH1F*)m_res_y->At(histID))->GetXaxis()->SetTitleSize(0.05);
602 // ((TH1F*)m_res_y->At(histID))->SetLineColor(1);
603 // ((TH1F*)m_res_y->At(histID))->Draw();
604 // }
605 // c11->Update();
606 // c11->Write();
607 
608 // i11=1;
609 // TCanvas* c21 = new TCanvas("respl3y","residuals, pl#3",200,500,700,800);
610 // c21->SetFillColor(0);
611 // c21->SetBorderMode(0);
612 // c21->Divide(10,10);
613 // for (unsigned int histID = 200; histID < 300; ++histID){
614 // c21->cd(i11);
615 // i11++;
616 // ((TH1F*)m_res_y->At(histID))->GetXaxis()->SetTitle("(trk-hit)[y], cm");
617 // ((TH1F*)m_res_y->At(histID))->GetXaxis()->SetTitleSize(0.05);
618 // ((TH1F*)m_res_y->At(histID))->SetLineColor(1);
619 // ((TH1F*)m_res_y->At(histID))->Draw();
620 // }
621 // c21->Update();
622 // c21->Write();
623 
624 // i11=1;
625 // TCanvas* c31 = new TCanvas("respl4y","residuals, pl#4",200,500,700,800);
626 // c31->SetFillColor(0);
627 // c31->SetBorderMode(0);
628 // c31->Divide(10,10);
629 // for (unsigned int histID = 300; histID < 400; ++histID){
630 // c31->cd(i11);
631 // i11++;
632 // ((TH1F*)m_res_y->At(histID))->GetXaxis()->SetTitle("(trk-hit)[y], cm");
633 // ((TH1F*)m_res_y->At(histID))->GetXaxis()->SetTitleSize(0.05);
634 // ((TH1F*)m_res_y->At(histID))->SetLineColor(1);
635 // ((TH1F*)m_res_y->At(histID))->Draw();
636 // }
637 // c31->Update();
638 // c31->Write();
639 
640 // //---------------------------------------------
641 
642 // TCanvas* c0mc = new TCanvas("resmcpl1x","residuals, pl#1",200,500,700,800);
643 // c0mc->SetFillColor(0);
644 // c0mc->SetBorderMode(0);
645 // c0mc->Divide(10,10);
646 // for (unsigned int histID = 0; histID < 100; ++histID){
647 // c0mc->cd(histID+1);
648 // ((TH1F*)m_res_mc_x->At(histID))->GetXaxis()->SetTitle("(rec-mc)[x], cm");
649 // ((TH1F*)m_res_mc_x->At(histID))->GetXaxis()->SetTitleSize(0.05);
650 // ((TH1F*)m_res_mc_x->At(histID))->SetLineColor(1);
651 // ((TH1F*)m_res_mc_x->At(histID))->Draw();
652 // }
653 // c0mc->Update();
654 // c0mc->Write();
655 
656 // TCanvas* c1mc = new TCanvas("resmcpl2x","residuals, pl#2",200,500,700,800);
657 // c1mc->SetFillColor(0);
658 // c1mc->SetBorderMode(0);
659 // c1mc->Divide(10,10);
660 // for (unsigned int histID = 100; histID < 200; ++histID){
661 // c1mc->cd(histID+1-100);
662 // ((TH1F*)m_res_mc_x->At(histID))->GetXaxis()->SetTitle("(rec-mc)[x], cm");
663 // ((TH1F*)m_res_mc_x->At(histID))->GetXaxis()->SetTitleSize(0.05);
664 // ((TH1F*)m_res_mc_x->At(histID))->SetLineColor(1);
665 // ((TH1F*)m_res_mc_x->At(histID))->Draw();
666 // }
667 // c1mc->Update();
668 // c1mc->Write();
669 
670 // TCanvas* c2mc = new TCanvas("resmcpl3x","residuals, pl#3",200,500,700,800);
671 // c2mc->SetFillColor(0);
672 // c2mc->SetBorderMode(0);
673 // c2mc->Divide(10,10);
674 // for (unsigned int histID = 100; histID < 200; ++histID){
675 // c2mc->cd(histID+1-100);
676 // ((TH1F*)m_res_mc_x->At(histID))->GetXaxis()->SetTitle("(rec-mc)[x], cm");
677 // ((TH1F*)m_res_mc_x->At(histID))->GetXaxis()->SetTitleSize(0.05);
678 // ((TH1F*)m_res_mc_x->At(histID))->SetLineColor(1);
679 // ((TH1F*)m_res_mc_x->At(histID))->Draw();
680 // }
681 // c2mc->Update();
682 // c2mc->Write();
683 
684 // TCanvas* c3mc = new TCanvas("resmcpl4x","residuals, pl#4",200,500,700,800);
685 // c3mc->SetFillColor(0);
686 // c3mc->SetBorderMode(0);
687 // c3mc->Divide(10,10);
688 // for (unsigned int histID = 200; histID < 300; ++histID){
689 // c3mc->cd(histID+1-200);
690 // ((TH1F*)m_res_mc_x->At(histID))->GetXaxis()->SetTitle("(rec-mc)[x], cm");
691 // ((TH1F*)m_res_mc_x->At(histID))->GetXaxis()->SetTitleSize(0.05);
692 // ((TH1F*)m_res_mc_x->At(histID))->SetLineColor(1);
693 // ((TH1F*)m_res_mc_x->At(histID))->Draw();
694 // }
695 // c3mc->Update();
696 // c3mc->Write();
697 
698 
699 // TCanvas* c01mc = new TCanvas("resmcpl1y","residuals, pl#1",200,500,700,800);
700 // c01mc->SetFillColor(0);
701 // c01mc->SetBorderMode(0);
702 // c01mc->Divide(10,10);
703 // for (unsigned int histID = 0; histID < 100; ++histID){
704 // c01mc->cd(histID+1);
705 // ((TH1F*)m_res_mc_y->At(histID))->GetXaxis()->SetTitle("(rec-mc)[y], cm");
706 // ((TH1F*)m_res_mc_y->At(histID))->GetXaxis()->SetTitleSize(0.05);
707 // ((TH1F*)m_res_mc_y->At(histID))->SetLineColor(1);
708 // ((TH1F*)m_res_mc_y->At(histID))->Draw();
709 // }
710 // c01mc->Update();
711 // c01mc->Write();
712 
713 // TCanvas* c11mc = new TCanvas("resmcpl2y","residuals, pl#2",200,500,700,800);
714 // c11mc->SetFillColor(0);
715 // c11mc->SetBorderMode(0);
716 // c11mc->Divide(10,10);
717 // i11=1;
718 // for (unsigned int histID = 100; histID < 200; ++histID){
719 // c11mc->cd(i11);
720 // i11++;
721 // ((TH1F*)m_res_mc_y->At(histID))->GetXaxis()->SetTitle("(rec-mc)[y], cm");
722 // ((TH1F*)m_res_mc_y->At(histID))->GetXaxis()->SetTitleSize(0.05);
723 // ((TH1F*)m_res_mc_y->At(histID))->SetLineColor(1);
724 // ((TH1F*)m_res_mc_y->At(histID))->Draw();
725 // }
726 // c11mc->Update();
727 // c11mc->Write();
728 
729 // i11=1;
730 // TCanvas* c21mc = new TCanvas("resmcpl3y","residuals, pl#3",200,500,700,800);
731 // c21mc->SetFillColor(0);
732 // c21mc->SetBorderMode(0);
733 // c21mc->Divide(10,10);
734 // for (unsigned int histID = 200; histID < 300; ++histID){
735 // c21mc->cd(i11);
736 // i11++;
737 // ((TH1F*)m_res_mc_y->At(histID))->GetXaxis()->SetTitle("(rec-mc)[y], cm");
738 // ((TH1F*)m_res_mc_y->At(histID))->GetXaxis()->SetTitleSize(0.05);
739 // ((TH1F*)m_res_mc_y->At(histID))->SetLineColor(1);
740 // ((TH1F*)m_res_mc_y->At(histID))->Draw();
741 // }
742 // c21mc->Update();
743 // c21mc->Write();
744 
745 // i11=1;
746 // TCanvas* c31mc = new TCanvas("resmcpl4y","residuals, pl#4",200,500,700,800);
747 // c31mc->SetFillColor(0);
748 // c31mc->SetBorderMode(0);
749 // c31mc->Divide(10,10);
750 // for (unsigned int histID = 300; histID < 400; ++histID){
751 // c31mc->cd(i11);
752 // i11++;
753 // ((TH1F*)m_res_mc_y->At(histID))->GetXaxis()->SetTitle("(rec-mc)[y], cm");
754 // ((TH1F*)m_res_mc_y->At(histID))->GetXaxis()->SetTitleSize(0.05);
755 // ((TH1F*)m_res_mc_y->At(histID))->SetLineColor(1);
756 // ((TH1F*)m_res_mc_y->At(histID))->Draw();
757 // }
758 // c31mc->Update();
759 // c31mc->Write();
760 
761 
762 // TCanvas* c02mc = new TCanvas("resmcpl1z","residuals, pl#1",200,500,700,800);
763 // c02mc->SetFillColor(0);
764 // c02mc->SetBorderMode(0);
765 // c02mc->Divide(10,10);
766 // for (unsigned int histID = 0; histID < 100; ++histID){
767 // c02mc->cd(histID+1);
768 // ((TH1F*)m_res_mc_z->At(histID))->GetXaxis()->SetTitle("(rec-mc)[y], cm");
769 // ((TH1F*)m_res_mc_z->At(histID))->GetXaxis()->SetTitleSize(0.05);
770 // ((TH1F*)m_res_mc_z->At(histID))->SetLineColor(1);
771 // ((TH1F*)m_res_mc_z->At(histID))->Draw();
772 // }
773 // c02mc->Update();
774 // c02mc->Write();
775 
776 // TCanvas* c12mc = new TCanvas("resmcpl2z","residuals, pl#2",200,500,700,800);
777 // c12mc->SetFillColor(0);
778 // c12mc->SetBorderMode(0);
779 // c12mc->Divide(10,10);
780 // i11=1;
781 // for (unsigned int histID = 100; histID < 200; ++histID){
782 // c12mc->cd(i11);
783 // i11++;
784 // ((TH1F*)m_res_mc_z->At(histID))->GetXaxis()->SetTitle("(rec-mc)[y], cm");
785 // ((TH1F*)m_res_mc_z->At(histID))->GetXaxis()->SetTitleSize(0.05);
786 // ((TH1F*)m_res_mc_z->At(histID))->SetLineColor(1);
787 // ((TH1F*)m_res_mc_z->At(histID))->Draw();
788 // }
789 // c12mc->Update();
790 // c12mc->Write();
791 
792 // i11=1;
793 // TCanvas* c22mc = new TCanvas("resmcpl3z","residuals, pl#3",200,500,700,800);
794 // c22mc->SetFillColor(0);
795 // c22mc->SetBorderMode(0);
796 // c22mc->Divide(10,10);
797 // for (unsigned int histID = 200; histID < 300; ++histID){
798 // c22mc->cd(i11);
799 // i11++;
800 // ((TH1F*)m_res_mc_z->At(histID))->GetXaxis()->SetTitle("(rec-mc)[y], cm");
801 // ((TH1F*)m_res_mc_z->At(histID))->GetXaxis()->SetTitleSize(0.05);
802 // ((TH1F*)m_res_mc_z->At(histID))->SetLineColor(1);
803 // ((TH1F*)m_res_mc_z->At(histID))->Draw();
804 // }
805 // c22mc->Update();
806 // c22mc->Write();
807 
808 // i11=1;
809 // TCanvas* c32mc = new TCanvas("resmcpl4z","residuals, pl#4",200,500,700,800);
810 // c32mc->SetFillColor(0);
811 // c32mc->SetBorderMode(0);
812 // c32mc->Divide(10,10);
813 // for (unsigned int histID = 300; histID < 400; ++histID){
814 // c32mc->cd(i11);
815 // i11++;
816 // ((TH1F*)m_res_mc_z->At(histID))->GetXaxis()->SetTitle("(rec-mc)[y], cm");
817 // ((TH1F*)m_res_mc_z->At(histID))->GetXaxis()->SetTitleSize(0.05);
818 // ((TH1F*)m_res_mc_z->At(histID))->SetLineColor(1);
819 // ((TH1F*)m_res_mc_z->At(histID))->Draw();
820 // }
821 // c32mc->Update();
822 // c32mc->Write();
823  resx->Write();
824  resy->Write();
825  resxy->Write();
826  x_id->Write();
827  y_id->Write();
828  z_id->Write();
829  nhits->Write();
830  htrks->Write();
831  f->Close();
832 
833 }
main
int main(int argc, char **argv)
Definition: MvdMQ/run/runFileSampler.cxx:29
verboseLevel
int verboseLevel
Definition: Lars/runMvdSim.C:7
PndSdsHit::GetPosition
TVector3 GetPosition() const
Definition: PndSdsHit.h:93
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
PndSdsDigi::GetIndex
Int_t GetIndex(int i=0) const
Definition: PndSdsDigi.h:63
nhits
Int_t nhits
Definition: analysis_point_hit_fwendcap.C:256
PndLmdDim::Instance
static PndLmdDim * Instance()
Definition: PndLmdDim.cxx:249
PndTrackCand::GetSortedHit
PndTrackCandHit GetSortedHit(UInt_t i)
Definition: PndTrackCand.h:54
PndLmdDim::Transform_global_to_lmd_local
TVector3 Transform_global_to_lmd_local(const TVector3 &point, bool isvector=false, bool aligned=true)
Definition: PndLmdDim.cxx:3113
MCPoint
TFile * MCPoint
Definition: anaLmdDigi.C:25
startEvent
Int_t startEvent
Definition: emc/run_linkCreation.C:7
PndSdsMCPoint::GetPositionOut
TVector3 GetPositionOut() const
Definition: PndSdsMCPoint.h:91
storePath
TString storePath
Definition: runGeaneTestSimBox.cxx:50
output
FairParRootFileIo * output
Definition: sim_emc_apd.C:120
DigiFile
TString DigiFile
Definition: emc/run_linkCreation.C:14
PndSdsHit::GetCov
TMatrixD GetCov() const
Definition: PndSdsHit.h:98
PndSdsCluster::GetDigiIndex
Int_t GetDigiIndex(Int_t i) const
Definition: PndSdsCluster.h:40
PndLmdDim::Get_sensor_by_id
void Get_sensor_by_id(const int sensor_id, int &ihalf, int &iplane, int &imodule, int &iside, int &idie, int &isensor)
Definition: PndLmdDim.h:526
nEvents
Int_t nEvents
Definition: hit_dirc.C:11
timer
TStopwatch timer
Definition: hit_dirc.C:51
PndLmdDim::Transform_lmd_local_to_module_side
TVector3 Transform_lmd_local_to_module_side(const TVector3 &point, int ihalf, int iplane, int imodule, int iside, bool isvector=false, bool aligned=true)
Definition: PndLmdDim.cxx:3126
f
TFile * f
Definition: bump_analys.C:12
PndSdsMCPoint::GetPosition
TVector3 GetPosition() const
Definition: PndSdsMCPoint.h:90
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
PndTrack::GetRefIndex
Int_t GetRefIndex() const
Definition: PndTrack.h:36
PndSdsHit::GetSensorID
Int_t GetSensorID() const
Definition: PndSdsHit.h:90
PndTrackCandHit::GetHitId
Int_t GetHitId() const
Definition: PndTrackCandHit.h:55
PndTrack::GetParamFirst
FairTrackParP GetParamFirst()
Definition: PndTrack.h:49
TMatrixD
TMatrixT< double > TMatrixD
Definition: PndLmdDim.h:52
Generated on Wed Apr 3 2019 10:02:45 for FairRoot/PandaRoot by   doxygen 1.8.5