HARPO  5.1.1
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
HarpoAnalyseResolution.cxx
Go to the documentation of this file.
1 //
2 // File HarpoAnalyseResolution.cxx
3 //
11 #include "HarpoAnalyseResolution.h"
12 #include "HarpoConfig.h"
13 #include "HarpoDetSet.h"
14 #include "HarpoDebug.h"
15 #include "HarpoDccEvent.h"
16 #include "Pmm2Event.h"
17 #include "HarpoEvent.h"
18 #include "HarpoRecoEvent.h"
19 
20 #include "HarpoSimEvent.h"
21 #include "TpcSimMCEvent.h"
22 #include "MakeNiceHisto.h"
23 
24 #include "TFile.h"
25 #include "TStyle.h"
26 #include "TCanvas.h"
27 #include "TLatex.h"
28 #include "TGraph.h"
29 #include "TF1.h"
30 #include "TMath.h"
31 #include "TSystem.h"
32 #include "TApplication.h"
33 
34 #include <cstdlib>
35 #include <cstring>
36 #include <cassert>
37 #include <fstream>
38 #include <iostream>
39 
40 ClassImp(HarpoAnalyseResolution);
41 
42 void HarpoAnalyseResolution::print()
43  {
44  unsigned int nd; // number of detectors
45  HarpoEventHeader *hdr = fEvt->GetHeader();
46 
47  assert(hdr != NULL);
48  hdr->print();
49 
50  for (nd = 0; nd < gkNDetectors; nd++) {
51  // if (fEvt->isdataExist(nd)) {
52  HarpoDccMap *plane = fEvt->GetDccMap(nd);
53  if (plane != NULL )plane->print();
54  }
55  }
56 
57 void HarpoAnalyseResolution::process()
58 {
59  // Bool_t drawEvent = kFALSE;
60  nEvents++;
61 
62 
63  if(nEvents%1000 == 0)
64  Info("process","Event %ld", nEvents);
65 
66  //???????????????????????????????
67  if(!fEvt->GetRecoEvent()) return;
68  // Reconstructed data (clusters, tracks, matched tracks)
69  HarpoRecoEvent* reco = fEvt->GetRecoEvent();
70 
71  TClonesArray* clArray = reco->GetClustersArray();
72  Int_t nCl = clArray->GetEntries();
73  TClonesArray* trArray = reco->GetTracksArray();
74  // Int_t NtrX=reco->GetNtracksXEvt();
75  // Int_t NtrY=reco->GetNtracksYEvt();
76 
77 
78  // if(reco->GetNclXEvt()<20)
79  // return;
80  // if(reco->GetQtotXEvt()/reco->GetNclXEvt()<2000)
81  // return;
82 
83 
84  if (gHDetSet->isExist(SIMDET)){
85 
86  HarpoSimEvent *simEvt = (HarpoSimEvent *)(fEvt->GetDetEvent(SIMDET));
87  TpcSimMCEvent *mcEvt = simEvt->GetMCEvent();
88  Int_t nTr = mcEvt->GetNtracks();
89  for(Int_t i = 0; i<nTr; i++){
90  TpcSimMCTrack* tr = mcEvt->GetTrack(i);
91  Double_t pX = tr->GetPx();
92  Double_t pY = tr->GetPy();
93  Double_t pZ = tr->GetPz();
94  Double_t x0 = tr->GetX0();
95  Double_t y0 = tr->GetY0();
96  // Double_t z0 = tr->GetZ0();
97 
98 
99  hAngleSimX->Fill(-TMath::ATan(pY/pZ));
100  Double_t thetaSim[2];
101  thetaSim[0] = -TMath::ATan(pY/pZ);
102  thetaSim[1] = -TMath::ATan(pX/pZ);
103  Double_t rhoSim[2];
104  rhoSim[0] = -((511./2. - 61.3)*TMath::Tan(thetaSim[0]) - y0 - 5.5)*TMath::Cos(thetaSim[0]);
105  rhoSim[1] = -((511./2. - 61.3)*TMath::Tan(thetaSim[1]) - x0 - 5.8)*TMath::Cos(thetaSim[1]);
106 
107  for(Int_t icl = 0; icl<nCl; icl++){
108  HarpoRecoClusters* cluster = (HarpoRecoClusters*)clArray->At(icl);
109  Int_t plane = cluster->GetPlane();
110  //if(TMath::Sin(theta[plane])<=0) continue;
111  Int_t type = cluster->GetType();
112  Double_t xCl = -1.0, zCl = -1.0 , delta = -1.0;
113  Double_t qCl = cluster->GetQ();
114  if(type == CCLUSTER){
115  zCl = cluster->GetMean();
116  xCl = cluster->GetIndex();
117  // delta = zCl - z0Sim[plane] - xCl/TMath::Tan(thetaSim[plane]) + 511/2;
118  delta = zCl - (rhoSim[plane]/TMath::Sin(thetaSim[plane]) - (xCl - 144.)/TMath::Tan(thetaSim[plane]) + 511/2);
119  }
120  if(type == TCLUSTER){
121  xCl = cluster->GetMean();
122  zCl = cluster->GetIndex();
123  //delta = xCl - x0Sim[plane] - zCl*TMath::Tan(thetaSim[plane]) + 144;
124  delta = xCl - (rhoSim[plane]/TMath::Cos(thetaSim[plane]) - (zCl - 511/2)*TMath::Tan(thetaSim[plane]) + 144);
125  }
126  if(type == BLOCCLUSTER){
127  xCl = cluster->GetMean();
128  zCl = cluster->GetIndex();
129  //delta = xCl - x0Sim[plane] - zCl*TMath::Tan(thetaSim[plane]) + 144;
130  delta = xCl - (rhoSim[plane]/TMath::Cos(thetaSim[plane]) - (zCl - 511/2)*TMath::Tan(thetaSim[plane]) + 144);
131  }
132  // cluster->SetXfit(0,cluster->GetMean() - delta);
133  if(cluster->GetQuality() != 0) continue;
134  if(cluster->GetWidth() < 3) continue;
135  Bool_t separators = (TMath::Abs(xCl - 94) <= 7) || (TMath::Abs(xCl - 194) <= 7);
136  Bool_t GemSectors =
137  (TMath::Abs(xCl - 28) <= 6) ||
138  (TMath::Abs(xCl - 61) <= 6) ||
139  (TMath::Abs(xCl - 127) <= 6) ||
140  (TMath::Abs(xCl - 160) <= 6) ||
141  (TMath::Abs(xCl - 226) <= 6) ||
142  (TMath::Abs(xCl - 259) <= 6);
143  hDeltaSim[type][plane]->Fill(delta);
144 
145 
146  Bool_t spaceCut = !separators && !(plane==0 && GemSectors) && xCl>25 && xCl<263;
147  Bool_t angleCut = TMath::Cos(thetaSim[plane])>0.5;
148  if(type == CCLUSTER)
149  angleCut = TMath::Cos(thetaSim[plane])>0.5 && TMath::Cos(thetaSim[plane])<0.8;
150 
151  if(spaceCut)
152  hDeltaSimVsA[type][plane]->Fill(TMath::Cos(thetaSim[plane]),delta);
153  //hDeltaSimVsA[type][plane]->Fill(theta[plane],delta);
154  if(angleCut){
155  if(spaceCut){
156  hDeltaSimVsT[type][plane]->Fill(zCl,delta);
157  hDeltaSimVsQ[type][plane]->Fill(qCl,delta);
158  }
159  hDeltaSimVsX[type][plane]->Fill(xCl,delta);
160  hDeltaSimVsXT[type][plane]->Fill(xCl,zCl,delta);
161  }
162  }
163  hAngleSimY->Fill(-TMath::ATan(pX/pZ));
164 
165  }
166  }
167 
168  // if(NtrX!=1) return;
169  // if(NtrY!=1) return;
170  const Int_t nTr = trArray->GetEntries();
171  if(gHarpoDebug>1)
172  Info("process","nTr = %d",nTr);
173  if(nTr != 2) return;
174 
175  Double_t angle[2], xstart[2], zstart[2];//, theta[2];
176  // Double_t angleCut = 100;//0.5;
177  for(Int_t itr = 0; itr<nTr; itr++){
178  // HarpoRecoTracks* track = (HarpoRecoTracks*)trArray->At(itr);
179  HarpoRecoTracks* track = (HarpoRecoTracks*)trArray->At(itr);
180  Int_t plane = track->GetPlane();
181  angle[plane] = track->GetAngleTrack();
182  // xstart[plane] = track->GetXstart();
183  xstart[plane] = track->GetX0();
184  zstart[plane] = 0;//track->GetZstart();
185  //theta[plane] = track->GetMtheta();
186 
187  }
188 
189  if(gHarpoDebug>1)
190  Info("process","angles = %g %g",angle[0],angle[1]);
191  // if(angle[0]==0) return;
192  // if(angle[1]==0) return;
193 
194 
195  // hAngleX->Fill(theta[0]);
196  // hAngleY->Fill(theta[1]);
197  hAngleX->Fill(angle[0]);
198  hAngleY->Fill(angle[1]);
199 
200 
201  for(Int_t icl = 0; icl<nCl; icl++){
202  HarpoRecoClusters* cluster = (HarpoRecoClusters*)clArray->At(icl);
203  Int_t plane = cluster->GetPlane();
204  Int_t type = cluster->GetType();
205  Double_t xCl = -1.0, zCl = -1.0, delta = -1.0;
206  Double_t qCl = cluster->GetQ();
207  if(type == CCLUSTER){
208  zCl = cluster->GetMean();
209  xCl = cluster->GetIndex();
210  delta = zCl - zstart[plane] - (xCl - xstart[plane])/angle[plane];
211  }
212  if(type == TCLUSTER){
213  xCl = cluster->GetMean();
214  zCl = cluster->GetIndex();
215  delta = xCl - xstart[plane] - (zCl - zstart[plane])*angle[plane];
216  }
217  if(type == BLOCCLUSTER){
218  xCl = cluster->GetX();
219  zCl = cluster->GetZ();
220  delta = (xCl - xstart[plane] - (zCl - zstart[plane])*angle[plane])/TMath::Sqrt(1+angle[plane]*angle[plane]);
221  }
222  if(gHarpoDebug>1)
223  Info("process","%g %g %g",xCl,zCl,delta);
224  if(cluster->GetQuality() != 0) continue;
225  if(cluster->GetWidth() < 3) continue;
226  // Double_t deltaZ = zCl - zstart[plane] - (xCl - xstart[plane])/angle[plane];
227  // Double_t deltaX = xCl - xstart[plane] - (zCl - zstart[plane])*angle[plane];
228  Bool_t separators = (TMath::Abs(xCl - 94) <= 7) || (TMath::Abs(xCl - 194) <= 7);
229  Bool_t GemSectors =
230  (TMath::Abs(xCl - 28) <= 6) ||
231  (TMath::Abs(xCl - 61) <= 6) ||
232  (TMath::Abs(xCl - 127) <= 6) ||
233  (TMath::Abs(xCl - 160) <= 6) ||
234  (TMath::Abs(xCl - 226) <= 6) ||
235  (TMath::Abs(xCl - 259) <= 6);
236  // if(type == 0){
237  // hDeltaZ[plane]->Fill(deltaZ);
238  // if(!separators && !(plane==0 && GemSectors))
239  // hDeltaZvsA[plane]->Fill(TMath::Cos(theta[plane]),deltaZ);
240  // if(TMath::Cos(theta[plane])<0.2){
241  // if(!separators && !(plane==0 && GemSectors))
242  // hDeltaZvsT[plane]->Fill(zCl,deltaZ);
243  // hDeltaZvsX[plane]->Fill(xCl,deltaZ);
244  // }
245  // }else{
246  // hDeltaX[plane]->Fill(deltaX);
247  // if(!separators && !(plane==0 && GemSectors))
248  // hDeltaXvsA[plane]->Fill(TMath::Cos(theta[plane]),deltaX);
249  // if(TMath::Cos(theta[plane])<0.2){
250  // if(!separators && !(plane==0 && GemSectors))
251  // hDeltaXvsT[plane]->Fill(zCl,deltaX);
252  // hDeltaXvsX[plane]->Fill(xCl,deltaX);
253  // }
254  // }
255  hDelta[type][plane]->Fill(delta);
256 
257 
258  Bool_t spaceCut = !separators && !(plane==0 && GemSectors) && xCl>25 && xCl<263;
259  // Bool_t angleCut = TMath::Cos(theta[plane])>0.5;
260  Bool_t angleCut = TMath::Abs(angle[plane])<1;
261  if(type == CCLUSTER)
262  angleCut = TMath::Abs(angle[plane])>0.1 && TMath::Abs(angle[plane])<1;
263  // if(type == BLOCCLUSTER)
264  // angleCut = 1;
265  // angleCut = TMath::Cos(theta[plane])>0.5 && TMath::Cos(theta[plane])<0.8;
266 
267  if(spaceCut)
268  hDeltaVsA[type][plane]->Fill(angle[plane],delta);
269  // hDeltaVsA[type][plane]->Fill(TMath::Cos(theta[plane]),delta);
270  //hDeltaVsA[type][plane]->Fill(theta[plane],delta);
271  if(angleCut){
272  if(spaceCut){
273  hDeltaVsT[type][plane]->Fill(zCl,delta);
274  hDeltaVsQ[type][plane]->Fill(qCl,delta);
275  }
276  hDeltaVsX[type][plane]->Fill(xCl,delta);
277  hDeltaVsXT[type][plane]->Fill(xCl,zCl,delta);
278  }
279  }
280 
281 
282  if (gHDetSet->isExist(SIMDET)){
283 
284  HarpoSimEvent *simEvt = (HarpoSimEvent *)(fEvt->GetDetEvent(SIMDET));
285  TpcSimMCEvent *mcEvt = simEvt->GetMCEvent();
286  Int_t nTr = mcEvt->GetNtracks();
287  for(Int_t i = 0; i<nTr; i++){
288  TpcSimMCTrack* tr = mcEvt->GetTrack(i);
289  Double_t pX = tr->GetPx();
290  Double_t pY = tr->GetPy();
291  Double_t pZ = tr->GetPz();
292  // Double_t x0 = tr->GetX0();
293  // Double_t y0 = tr->GetY0();
294  // Double_t z0 = tr->GetZ0();
295 
296  Double_t thetaSim[2];
297  thetaSim[0] = -TMath::ATan(pY/pZ);
298  thetaSim[1] = -TMath::ATan(pX/pZ);
299  // Double_t x0Sim[2];
300  // x0Sim[0] = y0 + 5.5 + 61.3*TMath::Tan(thetaSim[0]);
301  // x0Sim[1] = x0 + 5.8 + 61.3*TMath::Tan(thetaSim[1]);
302  // Double_t z0Sim[2];
303  // z0Sim[0] = (y0 + 149.5)/TMath::Tan(thetaSim[0]) + 61.3;
304  // z0Sim[1] = (x0 + 149.8)/TMath::Tan(thetaSim[1]) + 61.3;
305  // Double_t rhoSim[2];
306  // rhoSim[0] = -((511./2. - 61.3)*TMath::Tan(thetaSim[0]) - y0 - 5.5)*TMath::Cos(thetaSim[0]);
307  // rhoSim[1] = -((511./2. - 61.3)*TMath::Tan(thetaSim[1]) - x0 - 5.8)*TMath::Cos(thetaSim[1]);
308 
309  // hResAngleX->Fill(theta[0] - thetaSim[0]);
310  // hResAngleY->Fill(theta[1] - thetaSim[1]);
311  hResAngleX->Fill(angle[0] - thetaSim[0]);
312  hResAngleY->Fill(angle[1] - thetaSim[1]);
313 
314  hCompAngleX->Fill(TMath::Tan(angle[0]), pY/pZ);
315  hCompAngleY->Fill(TMath::Tan(angle[1]), pX/pZ);
316 
317  if(pZ)
318  hCompLength->Fill(TMath::Abs(TMath::Cos(angle[0])*TMath::Cos(angle[1])),pZ);
319 
320 
321  hAngleSimCutX->Fill(thetaSim[0]);
322  hAngleSimCutY->Fill(thetaSim[0]);
323 
324  // hResRhoX->Fill(rho[0]/TMath::Cos(theta[0]) + (511./2. - 61.3)*TMath::Tan(theta[0]) - y0 - 5.5);
325  // hResRhoY->Fill(rho[1]/TMath::Cos(theta[1]) + (511./2. - 61.3)*TMath::Tan(theta[1]) - x0 - 5.8);
326  // hResRhoX->Fill(rho[0]/TMath::Cos(angle[0]) + (511./2. - 61.3)*TMath::Tan(angle[0]) - y0 - 5.5);
327  // hResRhoY->Fill(rho[1]/TMath::Cos(angle[1]) + (511./2. - 61.3)*TMath::Tan(angle[1]) - x0 - 5.8);
328 
329  // hRhoX->Fill(rho[0]/TMath::Cos(theta[0]) - y0 - z0*TMath::Tan(thetaX),TMath::Tan(theta[0]));
330  // hRhoY->Fill(rho[1]/TMath::Cos(theta[1]) + (511./2. - 60)*TMath::Tan(theta[1]) - x0,TMath::Tan(theta[1]));
331 
332  for(Int_t i = 0; i<200; i++){
333  // hRhoX->Fill(rho[0]/TMath::Cos(theta[0]) + (511./2. - 70 + i*0.1)*TMath::Tan(theta[0]) - y0,i*0.1);
334  // hRhoY->Fill(rho[1]/TMath::Cos(theta[1]) + (511./2. - 70 + i*0.1)*TMath::Tan(theta[1]) - x0,i*0.1);
335  // hRhoX->Fill(rho[0]/TMath::Cos(angle[0]) + (511./2. - 70 + i*0.1)*TMath::Tan(angle[0]) - y0,i*0.1);
336  // hRhoY->Fill(rho[1]/TMath::Cos(angle[1]) + (511./2. - 70 + i*0.1)*TMath::Tan(angle[1]) - x0,i*0.1);
337  }
338  // hRhoX->Fill(rho[1]/TMath::Cos(theta[1]) - x0,TMath::Tan(theta[1]));
339 
340  // for(Int_t icl = 0; icl<nCl; icl++){
341  // HarpoRecoClusters* cluster = (HarpoRecoClusters*)clArray->At(icl);
342  // Int_t plane = cluster->GetPlane();
343  // //if(TMath::Sin(theta[plane])<=0) continue;
344  // Int_t type = cluster->GetType();
345  // Double_t xCl, zCl, delta;
346  // Double_t qCl = cluster->GetQ();
347  // if(type == 0){
348  // zCl = cluster->GetMean();
349  // xCl = cluster->GetIndex();
350  // // delta = zCl - z0Sim[plane] - xCl/TMath::Tan(thetaSim[plane]) + 511/2;
351  // delta = zCl - (rhoSim[plane]/TMath::Sin(thetaSim[plane]) - (xCl - 144.)/TMath::Tan(thetaSim[plane]) + 511/2);
352  // }else{
353  // xCl = cluster->GetMean();
354  // zCl = cluster->GetIndex();
355  // //delta = xCl - x0Sim[plane] - zCl*TMath::Tan(thetaSim[plane]) + 144;
356  // delta = xCl - (rhoSim[plane]/TMath::Cos(thetaSim[plane]) - (zCl - 511/2)*TMath::Tan(thetaSim[plane]) + 144);
357  // }
358  // cluster->SetFit(0,cluster->GetMean() - delta);
359  // if(cluster->GetQuality() != 0) continue;
360  // if(cluster->GetWidth() < 3) continue;
361  // Bool_t separators = (TMath::Abs(xCl - 94) <= 7) || (TMath::Abs(xCl - 194) <= 7);
362  // Bool_t GemSectors =
363  // (TMath::Abs(xCl - 28) <= 6) ||
364  // (TMath::Abs(xCl - 61) <= 6) ||
365  // (TMath::Abs(xCl - 127) <= 6) ||
366  // (TMath::Abs(xCl - 160) <= 6) ||
367  // (TMath::Abs(xCl - 226) <= 6) ||
368  // (TMath::Abs(xCl - 259) <= 6);
369  // hDeltaSim[type][plane]->Fill(delta);
370 
371 
372  // Bool_t spaceCut = !separators && !(plane==0 && GemSectors) && xCl>25 && xCl<263;
373  // Bool_t angleCut = TMath::Cos(theta[plane])>0.5;
374  // if(type == 0)
375  // angleCut = TMath::Cos(theta[plane])>0.5 && TMath::Cos(theta[plane])<0.8;
376 
377  // if(spaceCut)
378  // hDeltaSimVsA[type][plane]->Fill(TMath::Cos(theta[plane]),delta);
379  // //hDeltaSimVsA[type][plane]->Fill(theta[plane],delta);
380  // if(angleCut){
381  // if(spaceCut){
382  // hDeltaSimVsT[type][plane]->Fill(zCl,delta);
383  // hDeltaSimVsQ[type][plane]->Fill(qCl,delta);
384  // }
385  // hDeltaSimVsX[type][plane]->Fill(xCl,delta);
386  // hDeltaSimVsXT[type][plane]->Fill(xCl,zCl,delta);
387  // }
388  // }
389 
390  }
391 
392 
393 
394 
395  // Int_t nIonTr = simEvt->GetNtracks();
396  // for(Int_t i = 0; i<nIonTr; i++){
397  // TpcSimIonisationTrack* tr = simEvt->GetIonisationTrack(i);
398  // Int_t nPoints = tr->GetNpoints();
399  // for(Int_t j = 0; j<nPoints; j++){
400  // TpcSimIonisationPoint* point = tr->GetPoint(j);
401  // }
402  // }
403  }
404 
405 
406 }
407 
408 void HarpoAnalyseResolution::Init()
409 {
410 
411  hAngle = new TH1F("hAngle","map",200,-10,10);
412  hResAngleX = new TH1F("hResAngleX",";#sigma(#theta_{X})",1000,-TMath::Pi()/10,TMath::Pi()/10);
413  hResAngleY = new TH1F("hResAngleY",";#sigma(#theta_{Y})",1000,-TMath::Pi()/10,TMath::Pi()/10);
414  hAngleX = new TH1F("hAngleX",";#sigma(#theta_{X})",200,-TMath::Pi(),TMath::Pi());
415  hAngleY = new TH1F("hAngleY",";#sigma(#theta_{Y})",200,-TMath::Pi(),TMath::Pi());
416  hAngleSimX = new TH1F("hAngleSimX",";#sigma(#theta_{X})",200,-TMath::Pi(),TMath::Pi());
417  hAngleSimY = new TH1F("hAngleSimY",";#sigma(#theta_{Y})",200,-TMath::Pi(),TMath::Pi());
418  hAngleSimCutX = new TH1F("hAngleSimCutX",";#sigma(#theta_{X})",200,-TMath::Pi(),TMath::Pi());
419  hAngleSimCutY = new TH1F("hAngleSimCutY",";#sigma(#theta_{Y})",200,-TMath::Pi(),TMath::Pi());
420 
421 
422  hResRhoX = new TH1F("hResRhoX",";#sigma(#rho_{Y})",200,-10,10);
423  hResRhoY = new TH1F("hResRhoY",";#sigma(#rho_{Y})",200,-10,10);
424  hRhoX = new TH2F("hRhoX",";#rho_{X};X",2000,-20,20,200,0,20);//,-TMath::Pi(),TMath::Pi());
425  hRhoY = new TH2F("hRhoY",";#rho_{Y};Y",2000,-20,20,200,0,20);//200,-TMath::Pi(),TMath::Pi());
426 
427  hCompAngleX = new TH2F("hCompAngleX",";tan(#alpha_{X});p_{Y}/p_{Z}",2000,-20,20,2000,-20,20);//200,-TMath::Pi(),TMath::Pi());
428  hCompAngleY = new TH2F("hCompAngleY",";tan(#alpha_{Y});p_{X}/p_{Z}",2000,-20,20,2000,-20,20);//200,-TMath::Pi(),TMath::Pi());
429 
430  hCompLength = new TH2F("hCompLength",";cos(#alpha_{X})#timescos(#alpha_{Y});1/p_{Z}",100,0,1,100,0,1);//200,-TMath::Pi(),TMath::Pi());
431 
432 
433  Double_t max = 200;
434  // hDeltaX[0] = new TH1F("hDeltaXX",";#delta_{X}",1000,-max,max);
435  // hDeltaX[1] = new TH1F("hDeltaXY",";#delta_{X}",1000,-max,max);
436  // hDeltaZ[0] = new TH1F("hDeltaZX",";#delta_{Z}",1000,-max,max);
437  // hDeltaZ[1] = new TH1F("hDeltaZY",";#delta_{Z}",1000,-max,max);
438  // hDeltaXvsA[0] = new TH2F("hDeltaXvsAX",";#delta_{X}",100,0,1,1000,-max,max);
439  // hDeltaXvsA[1] = new TH2F("hDeltaXvsAY",";#delta_{X}",100,0,1,1000,-max,max);
440  // hDeltaZvsA[0] = new TH2F("hDeltaZvsAX",";#delta_{Z}",100,0,1,1000,-max,max);
441  // hDeltaZvsA[1] = new TH2F("hDeltaZvsAY",";#delta_{Z}",100,0,1,1000,-max,max);
442 
443  // hDeltaXvsT[0] = new TH2F("hDeltaXvsTX",";Z;#delta_{X}",512,0,512,1000,-max,max);
444  // hDeltaXvsT[1] = new TH2F("hDeltaXvsTY",";Z;#delta_{X}",512,0,512,1000,-max,max);
445  // hDeltaXvsX[0] = new TH2F("hDeltaXvsXX",";Z;#delta_{X}",288,0,288,1000,-max,max);
446  // hDeltaXvsX[1] = new TH2F("hDeltaXvsXY",";Z;#delta_{X}",288,0,288,1000,-max,max);
447 
448 
449  const char* type[3] = {"Z","X",""};
450  const char* plane[2] = {"X","Y"};
451  for(Int_t i = 0; i<3; i++){
452  for(Int_t j = 0; j<2; j++){
453  max = 20;
454  hDelta[i][j] = new TH1F(Form("hDelta%s%s",type[i],plane[j]),Form(";#delta_{%s}",type[i]),1000,-max,max);
455  // hDeltaVsA[i][j] = new TH2F(Form("hDelta%svsA%s",type[i],plane[j]),Form(";Angle;#delta_{%s}",type[i]),100,0,1,1000,-max,max);
456  hDeltaVsA[i][j] = new TH2F(Form("hDelta%svsA%s",type[i],plane[j]),Form(";Angle;#delta_{%s}",type[i]),100,-TMath::Pi(),TMath::Pi(),1000,-max,max);
457  hDeltaVsT[i][j] = new TH2F(Form("hDelta%svsT%s",type[i],plane[j]),Form(";T_{drift};#delta_{%s}",type[i]),512,0,512,1000,-max,max);
458  hDeltaVsQ[i][j] = new TH2F(Form("hDelta%svsQ%s",type[i],plane[j]),Form(";Q_{cl};#delta_{%s}",type[i]),1000,0,10000,1000,-max,max);
459  hDeltaVsX[i][j] = new TH2F(Form("hDelta%svsX%s",type[i],plane[j]),Form(";Channel;#delta_{%s}",type[i]),288,0,288,1000,-max,max);
460  hDeltaVsXT[i][j] = new TH3F(Form("hDelta%svsXT%s",type[i],plane[j]),Form(";Channel;T_{drift};#delta_{%s}",type[i]),288,0,288,512,0,512,200,-max,max);
461 
462  max = 20;
463  hDeltaSim[i][j] = new TH1F(Form("hDeltaSim%s%s",type[i],plane[j]),Form(";#delta_{%s}",type[i]),1000,-max,max);
464  // hDeltaSimVsA[i][j] = new TH2F(Form("hDeltaSim%svsA%s",type[i],plane[j]),Form(";Angle;#delta_{%s}",type[i]),100,0,1,1000,-max,max);
465  hDeltaSimVsA[i][j] = new TH2F(Form("hDeltaSim%svsA%s",type[i],plane[j]),Form(";Angle;#delta_{%s}",type[i]),100,-TMath::Pi(),TMath::Pi(),1000,-max,max);
466  hDeltaSimVsT[i][j] = new TH2F(Form("hDeltaSim%svsT%s",type[i],plane[j]),Form(";T_{drift};#delta_{%s}",type[i]),512,0,512,1000,-max,max);
467  hDeltaSimVsQ[i][j] = new TH2F(Form("hDeltaSim%svsQ%s",type[i],plane[j]),Form(";Q_{cl};#delta_{%s}",type[i]),1000,0,10000,1000,-max,max);
468  hDeltaSimVsX[i][j] = new TH2F(Form("hDeltaSim%svsX%s",type[i],plane[j]),Form(";Channel;#delta_{%s}",type[i]),288,0,288,1000,-max,max);
469  hDeltaSimVsXT[i][j] = new TH3F(Form("hDeltaSim%svsXT%s",type[i],plane[j]),Form(";Channel;T_{drift};#delta_{%s}",type[i]),288,0,288,512,0,512,200,-max,max);
470  }
471  }
472  // const Int_t nbinsQcl = 500;
473  // Double_t xminQcl = 10;
474  // Double_t xmaxQcl = 1e5;
475  // Double_t logxminQcl = TMath::Log(xminQcl);
476  // Double_t logxmaxQcl = TMath::Log(xmaxQcl);
477  // Double_t binwidthQcl = (logxmaxQcl-logxminQcl)/nbinsQcl;
478  // Double_t xbinsQcl[nbinsQcl+1];
479  // xbinsQcl[0] = xminQcl;
480  // for (Int_t i=1;i<=nbinsQcl;i++)
481  // xbinsQcl[i] = xminQcl + TMath::Exp(logxminQcl+i*binwidthQcl);
482  // hQcl[0] = new TH1F("hQclX",";Q_{cl}",nbinsQcl,xbinsQcl);
483  // hQcl[1] = new TH1F("hQclY",";Q_{cl}",nbinsQcl,xbinsQcl);
484 
485  }
486 
487 void HarpoAnalyseResolution::Save(char * /* mode */ )
488  {
489 
490 
491 
492  // TApplication* theApp = new TApplication("App", 0, 0);
493  // gStyle->SetOptStat(0);
494 
495  // TCanvas* c1 = new TCanvas();
496  // c1->Divide(2);
497  // MakeNiceHisto(hAngleSimX,c1->GetPad(1));
498  // MakeNiceHisto(hAngleSimY,c1->GetPad(2));
499  // TCanvas* c = new TCanvas();
500  // c->Divide(2);
501  // MakeNiceHisto(hResRhoX,c1->GetPad(1));
502  // MakeNiceHisto(hResRhoY,c1->GetPad(2));
503  // // MakeNiceHisto(hRhoX,c->GetPad(1),"colz");
504  // // MakeNiceHisto(hRhoY,c->GetPad(2),"colz");
505  // theApp->Run();
506 
507 
508 
509  /* TFile* hf = */ OpenHistFile("resolution");
510 
511  // Save histograms here
512  for(Int_t i = 0; i<3; i++){
513  for(Int_t j = 0; j<2; j++){
514  hDelta[i][j]->Write();
515  hDeltaVsA[i][j]->Write();
516  hDeltaVsX[i][j]->Write();
517  hDeltaVsT[i][j]->Write();
518  hDeltaVsQ[i][j]->Write();
519  hDeltaVsXT[i][j]->Write();
520 
521  hDeltaSim[i][j]->Write();
522  hDeltaSimVsA[i][j]->Write();
523  hDeltaSimVsX[i][j]->Write();
524  hDeltaSimVsT[i][j]->Write();
525  hDeltaSimVsQ[i][j]->Write();
526  hDeltaSimVsXT[i][j]->Write();
527  }
528  }
529 
530  hAngle->Write();
531  hResAngleX->Write();
532  hResAngleY->Write();
533  hAngleX->Write();
534  hAngleY->Write();
535  hAngleSimX->Write();
536  hAngleSimY->Write();
537  hAngleSimCutX->Write();
538  hAngleSimCutY->Write();
539  hCompAngleX->Write();
540  hCompAngleY->Write();
541  hCompLength->Write();
542 
543  hResRhoX->Write();
544  hResRhoY->Write();
545 
546  }
547 
HarpoEvent::GetRecoEvent
HarpoRecoEvent * GetRecoEvent()
Definition: HarpoEvent.h:46
TpcSimMCTrack::GetY0
Double_t GetY0()
Definition: TpcSimMCEvent.h:24
HarpoRecoClusters::GetMean
Double_t GetMean()
Definition: HarpoRecoEvent.h:39
TpcSimMCTrack::GetPz
Double_t GetPz()
Definition: TpcSimMCEvent.h:32
HarpoDetSet::isExist
Bool_t isExist(ULong_t det)
Detecror date exist //! Number of Real Detectors.
Definition: HarpoDetSet.h:33
HarpoRecoTracks::GetAngleTrack
Double_t GetAngleTrack()
Definition: HarpoRecoEvent.h:136
HarpoRecoTracks
Track object, containing position, angle, charge and quality information.
Definition: HarpoRecoEvent.h:125
HarpoRecoEvent
Object containing the reconstruction information for one event (with array of HarpoRecoClusters Harpo...
Definition: HarpoRecoEvent.h:474
BLOCCLUSTER
#define BLOCCLUSTER
Definition: HarpoRecoEvent.h:24
HarpoDccMap::print
virtual void print()
Definition: HarpoDccMap.cxx:124
TpcSimMCTrack::GetPy
Double_t GetPy()
Definition: TpcSimMCEvent.h:31
HarpoEventHeader
FullEvent Header not scecific to the detectors The class is ....
Definition: HarpoEventHeader.h:27
HarpoSimEvent
A class store HARPO row DCC event data and header. End provide access metods to the row data...
Definition: HarpoSimEvent.h:24
HarpoEventHeader::print
virtual void print()
Definition: HarpoEventHeader.cxx:21
HarpoRecoTracks::GetX0
Double_t GetX0()
Definition: HarpoRecoEvent.h:137
HarpoRecoClusters
Cluster object, containing position, charge and quality information.
Definition: HarpoRecoEvent.h:28
HarpoAnalyse::OpenHistFile
TFile * OpenHistFile(const char *ananame)
Definition: HarpoAnalyse.cxx:150
type
static int type
Definition: HarpoConfig.cxx:106
HarpoAnalyseResolution
Dummy analysis to run as test and example. Give basic histograms of the data.
Definition: HarpoAnalyseResolution.h:20
HarpoDccMap
unpacked dcc data The class contains the data map for DCC or Feminos The data is stored as a 2D TMatr...
Definition: HarpoDccMap.h:29
HarpoEvent::GetDetEvent
HarpoDetEvent * GetDetEvent(Long_t plane=XDCC)
Definition: HarpoEvent.cxx:93
SIMDET
Data from Keller temperuture and pressure sensors.
Definition: HarpoDet.h:22
TpcSimMCEvent::GetNtracks
Int_t GetNtracks()
Definition: TpcSimMCEvent.h:100
TpcSimMCEvent::GetTrack
TpcSimMCTrack * GetTrack(Int_t i)
Definition: TpcSimMCEvent.h:91
TpcSimMCTrack::GetX0
Double_t GetX0()
Definition: TpcSimMCEvent.h:23
TCLUSTER
#define TCLUSTER
Definition: HarpoRecoEvent.h:23
gHarpoDebug
Long64_t gHarpoDebug
Definition: HarpoDebug.cxx:9
HarpoSimEvent::GetMCEvent
TpcSimMCEvent * GetMCEvent()
Definition: HarpoSimEvent.h:57
HarpoEvent::GetHeader
HarpoEventHeader * GetHeader()
Definition: HarpoEvent.cxx:80
HarpoAnalyse::fEvt
HarpoEvent * fEvt
Definition: HarpoAnalyse.h:70
HarpoAnalyse::nEvents
ULong_t nEvents
Definition: HarpoAnalyse.h:75
HarpoAnalyseResolution::print
void print()
Overloaded method which do all job.
Definition: HarpoAnalyseResolution.cxx:42
HarpoAnalyseResolution::hDelta
TH1F * hDelta[3][2]
Redefine empty default.
Definition: HarpoAnalyseResolution.h:40
gkNDetectors
const ULong_t gkNDetectors
Definition: HarpoDet.h:14
CCLUSTER
#define CCLUSTER
Definition: HarpoRecoEvent.h:22
HarpoRecoEvent::GetTracksArray
TClonesArray * GetTracksArray()
Definition: HarpoRecoEvent.h:597
HarpoEvent::GetDccMap
HarpoDccMap * GetDccMap(Long_t plane=XDCC)
Definition: HarpoEvent.cxx:108
TpcSimMCTrack::GetPx
Double_t GetPx()
Definition: TpcSimMCEvent.h:30
HarpoRecoEvent::GetClustersArray
TClonesArray * GetClustersArray()
Definition: HarpoRecoEvent.h:594
gHDetSet
R__EXTERN HarpoDetSet * gHDetSet
Definition: HarpoDetSet.h:71
Generated by   doxygen 1.8.6