HARPO  5.1.1
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
HarpoAnalysePrf.cxx
Go to the documentation of this file.
1 //
2 // File HarpoAnalysePrf.cxx
3 //
12 #include "HarpoAnalysePrf.h"
13 #include "HarpoConfig.h"
14 #include "HarpoDetSet.h"
15 #include "HarpoDebug.h"
16 #include "HarpoDccEvent.h"
17 #include "HarpoEvent.h"
18 #include "HarpoRecoEvent.h"
19 
20 #include "TFile.h"
21 #include "TStyle.h"
22 #include "TCanvas.h"
23 #include "TLatex.h"
24 #include "TGraph.h"
25 #include "TF1.h"
26 #include "TMath.h"
27 #include "TSystem.h"
28 
29 #include <cstdlib>
30 #include <cstring>
31 #include <cassert>
32 #include <fstream>
33 #include <iostream>
34 
35 ClassImp(HarpoAnalysePrf)
36 
37 void HarpoAnalysePrf::print()
38  {
39  unsigned int nd; // number of detectors
40  HarpoEventHeader *hdr = fEvt->GetHeader();
41 
42  assert(hdr != NULL);
43  hdr->print();
44 
45  for (nd = 0; nd < gkNDetectors; nd++) {
46  // if (fEvt->isdataExist(nd)) {
47  HarpoDccMap *plane = fEvt->GetDccMap(nd);
48  if (plane != NULL )plane->print();
49  }
50  }
51 
52  void HarpoAnalysePrf::process()
53  {
54  // Bool_t drawEvent = kFALSE;
55  nEvents++;
56  //Info("process","");
57 
58  // Process Reconstructed data (clusters, tracks, matched tracks)
59  HarpoRecoEvent* reco = fEvt->GetRecoEvent();
60  TClonesArray* clArray = reco->GetClustersArray();
61  //TClonesArray* trArray = reco->GetTracksArray();
62  // Int_t NtrX=reco->GetNtracksXEvt();
63  // Int_t NtrY=reco->GetNtracksYEvt();
64  // fHistNtracks->Fill(NtrX,NtrY);
65  // if(NtrX!=1) return;
66  // if(NtrY!=1) return;
67  Int_t nCl = clArray->GetEntries();
68  //Int_t nTr = trArray->GetEntries();
69  // for(Int_t itr = 0; itr<nCl; itr++){
70  // HarpoRecoTracks* track = (HarpoRecoTracks*)trArray->At(itr);
71  // Double_t id = track->GetNtrack();
72  // }
73  fNeventsCut++;
74  Double_t tmin[2], tmax[2];
75  tmin[0] = 512;
76  tmax[0] = 0;
77  tmin[1] = 512;
78  tmax[1] = 0;
79  for(Int_t icl = 0; icl<nCl; icl++){
80  HarpoRecoClusters* cluster = (HarpoRecoClusters*)clArray->At(icl);
81  Int_t type = cluster->GetType();
82  // if(type==0) continue;
83  Double_t qCl = cluster->GetQ();
84  if(qCl<fQclMin) continue;
85  Double_t xCl = cluster->GetMean();
86  Int_t index = cluster->GetIndex();
87  if(type == CCLUSTER && index>=NCHAN) continue;
88  Int_t plane = cluster->GetPlane();
89  Int_t width = cluster->GetWidth();
90  // Int_t iHist = Int_t((index-90)/20);
91  // if(iHist<0) continue;
92  // if(iHist>=16) continue;
93  if(width<3) continue;
94  if( cluster->GetQuality() != 0) continue;
95  if(type == TCLUSTER && index<tmin[plane]) tmin[plane] = index;
96  if(type == TCLUSTER && index>tmax[plane]) tmax[plane] = index;
97 
98  Double_t sig = cluster->GetSig();
99  //Info("process","sigma = %g",sig);
100  if(type == CCLUSTER)
101  fHistSigmaZVsZ[plane]->Fill(xCl,sig*sig);
102  else
103  fHistSigmaXVsZ[plane]->Fill(index,sig*sig);
104 
105  HarpoDccMap* m = fEvt->GetDccMap(plane);
106 
107  Double_t qCenter;
108  if(type == CCLUSTER)
109  qCenter = m->GetData(index,Int_t(xCl+0.5));
110  else
111  qCenter = m->GetData(Int_t(xCl+0.5),index);
112  Double_t qTot = qCenter;
113  Double_t xNorm = xCl - Int_t(xCl+0.5);
114  Double_t qNorm = qCenter/qCl;
115  //std::cout << plane << " " << type << std::endl;
116  if(type == CCLUSTER){
117  fHistTrf[plane]->Fill(xCl,xNorm,qCenter);
118  fHistTrfNorm[plane]->Fill(xCl,xNorm,qNorm);
119  // if(width>30) continue;
120  // if(width<5) continue;
121  fHistTrfCut[plane]->Fill(xCl,xNorm,qCenter);
122  fHistTrfNormCut[plane]->Fill(xCl,xNorm,qNorm);
123  // fNtupleTrf->Fill(j,xCl,q,qCl,index,width);
124  }
125  if(type == TCLUSTER){
126  fHistPrf[plane]->Fill(index,xNorm,qCenter);
127  fHistPrfNorm[plane]->Fill(index,xNorm,qNorm);
128  if(width<=15 && width>=2){
129  fHistPrfCut[plane]->Fill(index,xNorm,qCenter);
130  fHistPrfNormCut[plane]->Fill(index,xNorm,qNorm);
131  }
132  // fNtuplePrf->Fill(j,xCl,q,qCl,index,width);
133  }
134  for(Int_t sign = -1; sign<2; sign += 2){
135  Double_t qold = qCenter;
136  Bool_t up = kTRUE;
137  for(Int_t j=Int_t(xCl+0.5)+sign;TMath::Abs(j-xCl)<width;j += sign){ // Channels
138  if(j<0) continue;
139  if(type == CCLUSTER && j>=NADC) break;
140  if(type == TCLUSTER && j>=NCHAN) break;
141  Double_t q, noise;
142  if(type == CCLUSTER){
143  q = m->GetData(index,j);
144  noise = m->GetSigmaNoise(index);
145  }else{
146  q = m->GetData(j,index);
147  noise = m->GetSigmaNoise(j);
148  }
149  if(q<=-1000) continue;
150  if(q<qold - noise) up = kFALSE;
151  if(!up && q>qold + noise) break;
152  qTot += q;
153  xNorm = xCl - j;
154  qNorm = q/qCl;
155  //std::cout << plane << " " << type << std::endl;
156  if(type == CCLUSTER){
157  fHistTrf[plane]->Fill(xCl,xNorm,q);
158  fHistTrfNorm[plane]->Fill(xCl,xNorm,qNorm);
159  // if(width>30) continue;
160  // if(width<5) continue;
161  fHistTrfCut[plane]->Fill(xCl,xNorm,q);
162  fHistTrfNormCut[plane]->Fill(xCl,xNorm,qNorm);
163  // fNtupleTrf->Fill(j,xCl,q,qCl,index,width);
164  }
165  if(type == TCLUSTER){
166  fHistPrf[plane]->Fill(index,xNorm,q);
167  fHistPrfNorm[plane]->Fill(index,xNorm,qNorm);
168  if(width<=15 && width>=2){
169  fHistPrfCut[plane]->Fill(index,xNorm,q);
170  fHistPrfNormCut[plane]->Fill(index,xNorm,qNorm);
171  // fNtuplePrf->Fill(j,xCl,q,qCl,index,width);
172  }
173  }
174  }
175  }
176  //std::cout << qCl << " => " << qTot << std::endl;
177  }
178 
179  for(Int_t i = 0; i<2; i++){
180  fHistTmin[i]->Fill(tmin[i]);
181  fHistTmax[i]->Fill(tmax[i]);
182  }
183 
184  //std::cout << std::endl;
185  }
186 
187  void HarpoAnalysePrf::Init()
188  {
189 
190  fQclMin = 1000;
191  Double_t qClMin;
192  if ( !gHConfig->Lookup("Prf.QclMin",qClMin) )
193  Info("Init","Use default Minimum cluster size: %.0f",fQclMin);
194  else{
195  fQclMin = qClMin;
196  Info("Init","Minimum cluster size: %.0f",fQclMin);
197  }
198 
199  // Initialise histograms here
200  Int_t nBinsY1 = 100;
201  Double_t qmin1 = 0.01;
202  Double_t qmax1 = 2;
203  Double_t yBins1[nBinsY1+1];
204  Double_t alpha1 = TMath::Power(qmax1/qmin1,1./(nBinsY1+1));
205  Double_t qbin1 = qmin1;
206  for(Int_t bin = 0; bin <=nBinsY1; bin++){
207  yBins1[bin] = qbin1;
208  qbin1 *= alpha1;
209  yBins1[bin] = qmin1 + bin*(qmax1-qmin1)/(nBinsY1+1);
210  }
211  yBins1[nBinsY1] = qmax1;
212 
213  // Initialise histograms here
214  Int_t nBinsY12 = 200;
215  Double_t qmin12 = 0.001;
216  Double_t qmax12 = 1.;
217  Double_t yBins12[nBinsY12+1];
218  Double_t alpha12 = TMath::Power(qmax12/qmin12,1./(nBinsY12+1));
219  Double_t qbin12 = qmin12;
220  for(Int_t bin = 0; bin <=nBinsY12; bin++){
221  yBins12[bin] = qbin12;
222  qbin12 *= alpha12;
223  }
224  yBins12[nBinsY12] = qmax12;
225 
226  Int_t nBinsY2 = 100;
227  Double_t qmin2 = 1;
228  Double_t qmax2 = 4096;
229  Double_t yBins2[nBinsY2+1];
230  Double_t alpha2 = TMath::Power(qmax2/qmin2,1./(nBinsY2+1));
231  Double_t qbin2 = qmin2;
232  for(Int_t bin = 0; bin <=nBinsY2; bin++){
233  yBins2[bin] = qbin2;
234  qbin2 *= alpha2;
235  }
236 
237  fNeventsCut = 0;
238 
239  fHistNtracks = new TH2F("fHistNtracks",";Ntr_{X};Ntr_{Y}",10,0,10,10,0,10);
240  fHistVdrift2D = new TH2F("fHistVdrift2D",";Time [t.b.];Charge [ADC]",511,0,511,1000,0,20000);
241  fHistVdrift = new TH1F("fHistVdrift",";Time [t.b.];Charge [ADC]",511,0,511);
242  Int_t nbinsA = 128;
243  Double_t xminA = 0, xmaxA = 512;
244  Double_t binsA[nbinsA+1];
245  for(Int_t i = 0; i<nbinsA; i++)
246  binsA[i] = xminA + i*(xmaxA-xminA)/(nbinsA+1);
247  binsA[nbinsA] = xmaxA;
248 
249  Int_t nbinsB = 100;
250  Double_t xminB = -10, xmaxB = 10;
251  Double_t binsB[nbinsB+1];
252  for(Int_t i = 0; i<nbinsB; i++)
253  binsB[i] = xminB + i*(xmaxB-xminB)/(nbinsB+1);
254  binsB[nbinsB] = xmaxB;
255 
256  Int_t nbinsC = 400;
257  Double_t xminC = -80, xmaxC = 80;
258  Double_t binsC[nbinsC+1];
259  for(Int_t i = 0; i<nbinsC; i++)
260  binsC[i] = xminC + i*(xmaxC-xminC)/(nbinsC+1);
261  binsC[nbinsC] = xmaxC;
262 
263  fHistPrf[0] = new TH3F("fHistPrfX",";X_{pix}-X_{cl};Q_{pix}",nbinsA,binsA,nbinsB,binsB,nBinsY2,yBins2);
264  fHistPrf[1] = new TH3F("fHistPrfY",";X_{pix}-X_{cl};Q_{pix}",nbinsA,binsA,nbinsB,binsB,nBinsY2,yBins2);
265  fHistTrf[0] = new TH3F("fHistTrfX",";T_{pix}-T_{cl};Q_{pix}",nbinsA,binsA,nbinsC,binsC,nBinsY2,yBins2);
266  fHistTrf[1] = new TH3F("fHistTrfY",";T_{pix}-T_{cl};Q_{pix}",nbinsA,binsA,nbinsC,binsC,nBinsY2,yBins2);
267  fHistPrfNorm[0] = new TH3F("fHistPrfNormX",";X_{pix}-X_{cl};Q_{pix}/Q_{cl}",nbinsA,binsA,nbinsB,binsB,nBinsY1,yBins1);
268  fHistPrfNorm[1] = new TH3F("fHistPrfNormY",";X_{pix}-X_{cl};Q_{pix}/Q_{cl}",nbinsA,binsA,nbinsB,binsB,nBinsY1,yBins1);
269  fHistTrfNorm[0] = new TH3F("fHistTrfNormX",";T_{pix}-T_{cl};Q_{pix}/Q_{cl}",nbinsA,binsA,nbinsC,binsC,nBinsY12,yBins12);
270  fHistTrfNorm[1] = new TH3F("fHistTrfNormY",";T_{pix}-T_{cl};Q_{pix}/Q_{cl}",nbinsA,binsA,nbinsC,binsC,nBinsY12,yBins12);
271 
272  fHistPrfCut[0] = new TH3F("fHistPrfCutX",";X_{pix}-X_{cl};Q_{pix}",nbinsA,binsA,nbinsB,binsB,nBinsY2,yBins2);
273  fHistPrfCut[1] = new TH3F("fHistPrfCutY",";X_{pix}-X_{cl};Q_{pix}",nbinsA,binsA,nbinsB,binsB,nBinsY2,yBins2);
274  fHistTrfCut[0] = new TH3F("fHistTrfCutX",";T_{pix}-T_{cl};Q_{pix}",nbinsA,binsA,nbinsC,binsC,nBinsY2,yBins2);
275  fHistTrfCut[1] = new TH3F("fHistTrfCutY",";T_{pix}-T_{cl};Q_{pix}",nbinsA,binsA,nbinsC,binsC,nBinsY2,yBins2);
276  fHistPrfNormCut[0] = new TH3F("fHistPrfNormCutX",";X_{pix}-X_{cl};Q_{pix}/Q_{cl}",nbinsA,binsA,nbinsB,binsB,nBinsY1,yBins1);
277  fHistPrfNormCut[1] = new TH3F("fHistPrfNormCutY",";X_{pix}-X_{cl};Q_{pix}/Q_{cl}",nbinsA,binsA,nbinsB,binsB,nBinsY1,yBins1);
278  fHistTrfNormCut[0] = new TH3F("fHistTrfNormCutX",";T_{pix}-T_{cl};Q_{pix}/Q_{cl}",nbinsA,binsA,nbinsC,binsC,nBinsY12,yBins12);
279  fHistTrfNormCut[1] = new TH3F("fHistTrfNormCutY",";T_{pix}-T_{cl};Q_{pix}/Q_{cl}",nbinsA,binsA,nbinsC,binsC,nBinsY12,yBins12);
280 
281 
282  fHistSigmaXVsZ[0] = new TH2F("fHistSigmaXVsZX",";t [30ns t.b.]; #sigma_{cluster}",512,0,512,500,0,50);
283  fHistSigmaXVsZ[1] = new TH2F("fHistSigmaXVsZY",";t [30ns t.b.]; #sigma_{cluster}",512,0,512,500,0,50);
284  fHistSigmaZVsZ[0] = new TH2F("fHistSigmaZVsZX",";t [30ns t.b.]; #sigma_{cluster}",512,0,512,500,0,50);
285  fHistSigmaZVsZ[1] = new TH2F("fHistSigmaZVsZY",";t [30ns t.b.]; #sigma_{cluster}",512,0,512,500,0,50);
286 
287 
288  fHistTmin[0] = new TH1F("fHistTminX",";t_{0} [30ns t.b.];",512,0,512);
289  fHistTmin[1] = new TH1F("fHistTminY",";t_{0} [30ns t.b.];",512,0,512);
290  fHistTmax[0] = new TH1F("fHistTmaxX",";t_{0} [30ns t.b.];",512,0,512);
291  fHistTmax[1] = new TH1F("fHistTmaxY",";t_{0} [30ns t.b.];",512,0,512);
292 
293 
294  }
295 
296 void HarpoAnalysePrf::Save(char * /* mode */)
297  {
298 
299  /* TFile *hf = */ OpenHistFile("prf");
300 
301  fHistPrf[0]->Write();
302  fHistPrf[1]->Write();
303  fHistTrf[0]->Write();
304  fHistTrf[1]->Write();
305  fHistPrfNorm[0]->Write();
306  fHistPrfNorm[1]->Write();
307  fHistTrfNorm[0]->Write();
308  fHistTrfNorm[1]->Write();
309 
310  fHistPrfCut[0]->Write();
311  fHistPrfCut[1]->Write();
312  fHistTrfCut[0]->Write();
313  fHistTrfCut[1]->Write();
314  fHistPrfNormCut[0]->Write();
315  fHistPrfNormCut[1]->Write();
316  fHistTrfNormCut[0]->Write();
317  fHistTrfNormCut[1]->Write();
318 
319  fHistSigmaXVsZ[0]->Write();
320  fHistSigmaXVsZ[1]->Write();
321  fHistSigmaZVsZ[0]->Write();
322  fHistSigmaZVsZ[1]->Write();
323 
324  fHistTmin[0]->Write();
325  fHistTmin[1]->Write();
326  fHistTmax[0]->Write();
327  fHistTmax[1]->Write();
328 
329  // hf->Close();
330  // delete hf;
331  }
332 
HarpoAnalysePrf::fQclMin
Double_t fQclMin
Redefine empty default.
Definition: HarpoAnalysePrf.h:34
HarpoEvent::GetRecoEvent
HarpoRecoEvent * GetRecoEvent()
Definition: HarpoEvent.h:46
HarpoAnalysePrf::fHistPrf
TH3F * fHistPrf[2]
Definition: HarpoAnalysePrf.h:43
HarpoAnalysePrf::Save
void Save(char *mode=NULL)
Definition: HarpoAnalysePrf.cxx:296
NCHAN
#define NCHAN
Definition: HarpoDccMap.h:16
HarpoRecoClusters::GetMean
Double_t GetMean()
Definition: HarpoRecoEvent.h:39
HarpoAnalysePrf::fHistTmin
TH1F * fHistTmin[2]
Definition: HarpoAnalysePrf.h:56
HarpoDccMap::GetData
Double_t GetData(Int_t i, Int_t j)
Set/Get Data Cell.
Definition: HarpoDccMap.cxx:84
HarpoAnalysePrf::fHistSigmaXVsZ
TH2F * fHistSigmaXVsZ[2]
Definition: HarpoAnalysePrf.h:53
HarpoAnalysePrf::fHistSigmaZVsZ
TH2F * fHistSigmaZVsZ[2]
Definition: HarpoAnalysePrf.h:54
HarpoRecoEvent
Object containing the reconstruction information for one event (with array of HarpoRecoClusters Harpo...
Definition: HarpoRecoEvent.h:474
HarpoDccMap::print
virtual void print()
Definition: HarpoDccMap.cxx:124
HarpoEventHeader
FullEvent Header not scecific to the detectors The class is ....
Definition: HarpoEventHeader.h:27
HarpoAnalysePrf::fHistTrfNorm
TH3F * fHistTrfNorm[2]
Definition: HarpoAnalysePrf.h:46
HarpoAnalysePrf::fHistPrfCut
TH3F * fHistPrfCut[2]
Definition: HarpoAnalysePrf.h:48
HarpoEventHeader::print
virtual void print()
Definition: HarpoEventHeader.cxx:21
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
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
HarpoAnalysePrf::fHistPrfNorm
TH3F * fHistPrfNorm[2]
Definition: HarpoAnalysePrf.h:44
NADC
#define NADC
Definition: HarpoDccMap.h:18
TCLUSTER
#define TCLUSTER
Definition: HarpoRecoEvent.h:23
HarpoAnalysePrf
Analysis of space and time response of the readout.
Definition: HarpoAnalysePrf.h:21
HarpoConfigFile::Lookup
Bool_t Lookup(const char *path, Bool_t &val)
Lookup function for scalar values.
Definition: HarpoConfigFile.cxx:85
tmax
const ULong_t tmax
Definition: HarpoTimeStamp.cxx:18
HarpoAnalyse::fEvt
HarpoEvent * fEvt
Definition: HarpoAnalyse.h:70
HarpoAnalyse::nEvents
ULong_t nEvents
Definition: HarpoAnalyse.h:75
HarpoAnalysePrf::fHistTrf
TH3F * fHistTrf[2]
Definition: HarpoAnalysePrf.h:45
HarpoAnalysePrf::fHistTmax
TH1F * fHistTmax[2]
Definition: HarpoAnalysePrf.h:57
gkNDetectors
const ULong_t gkNDetectors
Definition: HarpoDet.h:14
CCLUSTER
#define CCLUSTER
Definition: HarpoRecoEvent.h:22
gHConfig
HarpoConfig * gHConfig
Definition: HarpoConfig.cxx:101
HarpoAnalysePrf::fHistTrfNormCut
TH3F * fHistTrfNormCut[2]
Definition: HarpoAnalysePrf.h:51
HarpoEvent::GetDccMap
HarpoDccMap * GetDccMap(Long_t plane=XDCC)
Definition: HarpoEvent.cxx:108
HarpoDccMap::GetSigmaNoise
Double_t GetSigmaNoise(Int_t i)
Definition: HarpoDccMap.cxx:68
HarpoRecoEvent::GetClustersArray
TClonesArray * GetClustersArray()
Definition: HarpoRecoEvent.h:594
HarpoAnalysePrf::fHistTrfCut
TH3F * fHistTrfCut[2]
Definition: HarpoAnalysePrf.h:50
HarpoAnalysePrf::fHistPrfNormCut
TH3F * fHistPrfNormCut[2]
Definition: HarpoAnalysePrf.h:49
Generated by   doxygen 1.8.6