HARPO  5.1.1
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
HarpoAnalysePattern.cxx
Go to the documentation of this file.
1 //
2 // File HarpoAnalysePattern.cxx
3 //
11 #include "HarpoAnalysePattern.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 #include "HarpoSimEvent.h"
20 #include "TpcSimMCEvent.h"
21 #include "MakeNiceHisto.h"
22 
23 #include "TFile.h"
24 #include "TStyle.h"
25 #include "TCanvas.h"
26 #include "TRootEmbeddedCanvas.h"
27 #include "TLatex.h"
28 #include "TGraph.h"
29 #include "TF1.h"
30 #include "TMath.h"
31 #include "TArrow.h"
32 #include "TGLabel.h"
33 #include "TSystem.h"
34 
35 #include <cstdlib>
36 #include <cstring>
37 #include <cassert>
38 #include <fstream>
39 #include <iostream>
40 
41 ClassImp(HarpoAnalysePattern)
42 
43 void HarpoAnalysePattern::print()
44  {
45  unsigned int nd; // number of detectors
46  HarpoEventHeader *hdr = fEvt->GetHeader();
47 
48  assert(hdr != NULL);
49  hdr->print();
50 
51  for (nd = 0; nd < gkNDetectors; nd++) {
52  // if (fEvt->isdataExist(nd)) {
53  HarpoDccMap *plane = fEvt->GetDccMap(nd);
54  if (plane != NULL )plane->print();
55  }
56  }
57 
58 void HarpoAnalysePattern::process()
59 {
60  // Bool_t drawEvent = kFALSE;
61  nEvents++;
62 
63  // Process RAW data
64  for(Int_t ndet=0;ndet<2;ndet++) { // Detectors 0 and 1 = X and Y planes
65  if (!gHDetSet->isExist(ndet)) continue;
66  HarpoDccMap* m = fEvt->GetDccMap(ndet);
67  if ( m == NULL ) continue;
68  for(Int_t iV = 0; iV<kMaxNvertex; iV++){
69  fHistRoute[ndet][iV]->Reset();
70  fHistRoute2[ndet][iV]->Reset();
71  fHistShortest[ndet][iV]->Reset();
72  for(Int_t j=0;j<NCHAN;j++){ // Channels
73  //Double_t qch = 0;
74  for(Int_t i=0;i<NADC;i++){ // Time bins
75  Double_t q = m->GetData(j,i);
76  if(q <= -1000) continue;
77  fHistShortest[ndet][iV]->SetBinContent(i/fRebinZ+1,j/fRebinX+1,1000000);
78  }
79  }
80  }
81 
82  Int_t lengthmax = 0, i1max = -1, j1max = -1, i2max = -1, j2max = -1;
83  for(Int_t j1=0;j1<NCHAN/fRebinX;j1++){ // Channels
84  for(Int_t i1=0;i1<NADC/fRebinZ;i1++){ // Time bins
85  if(fHistShortest[ndet][0]->GetBinContent(i1+1,j1+1)<=0) continue;
86  TH2I* h = (TH2I*)fHistShortest[ndet][0]->Clone("hShortest");
87  GetShortest(h,1,i1+1,j1+1);
88  // Int_t lengthmax = 0, imax = -1, jmax = -1;
89  for(Int_t j2=0;j2<NCHAN/fRebinX;j2++){ // Channels
90  for(Int_t i2=0;i2<NADC/fRebinZ;i2++){ // Time bins
91  Int_t l = h->GetBinContent(i2+1,j2+1);
92  if(l <= 0) continue;
93  if(l >10000) continue;
94  if(l>lengthmax){
95  lengthmax = l;
96  i1max = i1+1;
97  j1max = j1+1;
98  i2max = i2+1;
99  j2max = j2+1;
100  }
101  }
102  }
103  }
104  }
105  fHistRoute2[ndet][3]->SetBinContent(i1max,j1max,10000);
106  fHistRoute2[ndet][3]->SetBinContent(i2max,j2max,10000);
107 
108  }
109 
110  //return;
111 
112 
113  if(fEvt->GetRecoEvent()){ // Reconstructed data (clusters, tracks, matched tracks)
114  HarpoRecoEvent* reco = fEvt->GetRecoEvent();
115  TClonesArray* vArray = reco->GetVertexArray();
116  Int_t nV = vArray->GetEntries();
117  Int_t nVp[2];
118  nVp[0] = 0;
119  nVp[1] = 0;
120  for(Int_t iV = 0; iV<nV; iV++){
121  HarpoRecoVertex* vertex = (HarpoRecoVertex*)vArray->At(iV);
122  Int_t plane = vertex->GetPlane();
123  Double_t xV = vertex->GetXvertex();
124  Double_t zV = vertex->GetZvertex();
125  Int_t iv = nVp[plane];
126  nVp[plane]++;
127  Info("process","%d %d",Int_t(zV/fRebinZ)+1,Int_t(xV/fRebinX)+1);
128  GetShortest(fHistShortest[plane][iv],1,Int_t(zV/fRebinZ)+1,Int_t(xV/fRebinX)+1);
129  Int_t lengthmax = 0, imax = -1, jmax = -1;
130  for(Int_t j=0;j<NCHAN/fRebinX;j++){ // Channels
131  for(Int_t i=0;i<NADC/fRebinZ;i++){ // Time bins
132  Int_t l = fHistShortest[plane][iv]->GetBinContent(i+1,j+1);
133  if(l <= 0) continue;
134  if(l >10000) continue;
135  if(l>lengthmax){
136  lengthmax = l;
137  imax = i+1;
138  jmax = j+1;
139  }
140  }
141  }
142  if(lengthmax<=1) continue;
143  GetRoute(fHistRoute[plane][iv],fHistShortest[plane][iv],lengthmax+1,imax,jmax);
144  Int_t lengthmax2 = 0, imax2 = -1, jmax2 = -1;
145  for(Int_t j=0;j<NCHAN/fRebinX;j++){ // Channels
146  for(Int_t i=0;i<NADC/fRebinZ;i++){ // Time bins
147  if(fHistRoute[plane][iv]->GetBinContent(i+1,j+1)>0) continue;
148  Int_t l = fHistShortest[plane][iv]->GetBinContent(i+1,j+1);
149  if(l <= 0) continue;
150  if(l >10000) continue;
151  if(l>lengthmax2){
152  lengthmax2 = l;
153  imax2 = i+1;
154  jmax2 = j+1;
155  }
156  }
157  }
158  if(lengthmax2<=2) continue;
159  GetRoute(fHistRoute2[plane][iv],fHistShortest[plane][iv],lengthmax2+1,imax2,jmax2);
160  }
161  }
162 
163 
164 }
165 
166 
167 
168 void HarpoAnalysePattern::GetShortest(TH2I* h, Int_t length, Int_t i, Int_t j)
169 {
170 
171  // if(i<0) return;
172  // if(j<0) return;
173  // if(i>kNx) return;
174  // if(i<kNz) return;
175 
176  if(h->GetBinContent(i,j)>length) h->SetBinContent(i,j,length);
177  else return;
178  GetShortest(h,length+1,i-1,j);
179  GetShortest(h,length+1,i,j-1);
180  GetShortest(h,length+1,i+1,j);
181  GetShortest(h,length+1,i,j+1);
182 
183 }
184 
185 
186 
187 void HarpoAnalysePattern::GetRoute(TH2I* hRoute, TH2I* h, Int_t length, Int_t i, Int_t j)
188 {
189 
190  Info("GetRoute","%d %d %d / %d",i,j,Int_t(h->GetBinContent(i,j)),length);
191  if(h->GetBinContent(i,j)!=length-1) return;
192  if(hRoute->GetBinContent(i,j)>0) return;
193  hRoute->SetBinContent(i,j,100000);
194  GetRoute(hRoute,h,length-1,i-1,j);
195  GetRoute(hRoute,h,length-1,i,j-1);
196  GetRoute(hRoute,h,length-1,i+1,j);
197  GetRoute(hRoute,h,length-1,i,j+1);
198 }
199 
200 
201 
202 
203 void HarpoAnalysePattern::Init()
204 {
205 
206  fChooseRebinX = 0;
207 
208  fRebinX = 6;//5;
209  fRebinZ = 12;//10;
210 
211  Long64_t rebinX = 10;
212  if( ! gHConfig->Lookup("ClusteringBlocs.rebinX",rebinX ))
213  Info("Constructor","Use default RebinX %d",fRebinX);
214  else
215  fRebinX = rebinX;
216 
217  Long64_t rebinZ = 10;
218  if( ! gHConfig->Lookup("ClusteringBlocs.rebinZ",rebinZ ))
219  Info("Constructor","Use default RebinZ %d",fRebinZ);
220  else
221  fRebinZ = rebinZ;
222 
223  for(Int_t i = 0; i<kMaxNvertex; i++){
224  fHistShortest[XPLANE][i] = new TH2I("fHistShortestX","",NADC/fRebinZ,0,NADC,NCHAN/fRebinX,0,NCHAN);
225  fHistShortest[YPLANE][i] = new TH2I("fHistShortestY","",NADC/fRebinZ,0,NADC,NCHAN/fRebinX,0,NCHAN);
226  fHistRoute[XPLANE][i] = new TH2I("fHistRouteX","",NADC/fRebinZ,0,NADC,NCHAN/fRebinX,0,NCHAN);
227  fHistRoute[YPLANE][i] = new TH2I("fHistRouteY","",NADC/fRebinZ,0,NADC,NCHAN/fRebinX,0,NCHAN);
228  fHistRoute2[XPLANE][i] = new TH2I("fHistRoute2X","",NADC/fRebinZ,0,NADC,NCHAN/fRebinX,0,NCHAN);
229  fHistRoute2[YPLANE][i] = new TH2I("fHistRoute2Y","",NADC/fRebinZ,0,NADC,NCHAN/fRebinX,0,NCHAN);
230  }
231 
232 }
233 
234 
235 
236 void HarpoAnalysePattern::DisplayAnalysis(TRootEmbeddedCanvas* ecTab, TGListBox* /* infobox */ )
237 {
238  // in hrecomonitorgui: fills the canvas on the "Analysis" tab
239  TCanvas* c = ecTab->GetCanvas();
240  c->GetPad(1)->Delete();
241  c->GetPad(2)->Delete();
242  c->Divide(2);
243 
244 
245  HarpoRecoEvent* reco = fEvt->GetRecoEvent();
246  TClonesArray* vArray = reco->GetVertexArray();
247  Int_t nV = vArray->GetEntries();
248  TClonesArray* clArray = reco->GetClustersArray();
249  Int_t nCl = clArray->GetEntries();
250 
251  for(Int_t plane = 0; plane<2;plane++){
252  TVirtualPad* cMap = c->GetPad(plane+1);
253  cMap->Divide(2,2);
254  for(Int_t i = 0; i<4; i++){
255  MakeNiceHisto(fHistShortest[plane][i],cMap->GetPad(i+1),"col");
256  fHistRoute[plane][i]->SetLineColor(kRed);
257  fHistRoute[plane][i]->Draw("cont3 same");
258  fHistRoute2[plane][i]->SetLineColor(kBlue);
259  fHistRoute2[plane][i]->Draw("cont3 same");
260  }
261 
262  for(Int_t iV = 0; iV<nV; iV++){
263  HarpoRecoVertex* vertex = (HarpoRecoVertex*)vArray->At(iV);
264  if(plane!=vertex->GetPlane()) continue;
265  Double_t xV = vertex->GetXvertex();
266  Double_t zV = vertex->GetZvertex();
267  Double_t nx = vertex->GetPxVertex();
268  Double_t nz = vertex->GetPzVertex();
269  for(Int_t k = 0; k<5; k++){
270  Double_t X = 0, Z = 0;
271  for(Int_t icl1 = 0; icl1<nCl; icl1++){
272  HarpoRecoClusters* cluster1 = (HarpoRecoClusters*)clArray->At(icl1);
273  if(plane!=cluster1->GetPlane()) continue;
274  Double_t xV1 = cluster1->GetX();
275  Double_t zV1 = cluster1->GetZ();
276  // for(Int_t icl2 = 0; icl2<nCl; icl2++){
277  // HarpoRecoClusters* cluster2 = (HarpoRecoClusters*)clArray->At(icl2);
278  // if(plane!=cluster2->GetPlane()) continue;
279  // Double_t xV2 = cluster2->GetX()-xV;
280  // Double_t zV2 = cluster2->GetZ()-zV;
281  if(nx*(xV1-xV) + nz*(zV1-zV)>0){
282  X += xV1 - xV;
283  Z += zV1 - zV;
284  }else{
285  X += xV - xV1;
286  Z += zV - zV1;
287  }
288  //}
289  }
290 
291  cout << "(X,Z) = (" << X << ", " << Z << ")" << endl;
292  Double_t norm = TMath::Sqrt(X*X+Z*Z);
293  nx = X/norm;
294  nz = Z/norm;
295  }
296 
297  cout << "n = (" << nx << ", " << nz << ")" << endl;
298  cMap->cd(1);
299  TArrow* arrow = new TArrow();
300  arrow->SetLineColor(kGray);
301  arrow->SetFillColor(kGray);
302  arrow->SetArrowSize(0.01);
303  arrow->DrawArrow(zV,xV,zV+30*nz,xV+30*nx);
304  }
305 
306  }
307 
308  c->Update();
309 }
310 
311 
312 void HarpoAnalysePattern::ConfigFrame(TGMainFrame* fMain, Int_t id)
313 {
314  // in hrecomonitorgui: Creates a popup window for analysis configuration
315  // You must define SetConfig() properly
316 
317  UInt_t xsize = 200;
318  UInt_t ysize2 = 20;
319  UInt_t ysize = 10*ysize2;
320  TGTransientFrame* main = new TGTransientFrame(gClient->GetRoot(), fMain, xsize, ysize);
321  main->Connect("CloseWindow()", "HarpoRecoMonitorGui", main, "CloseWindow()");
322  main->DontCallClose(); // to avoid double deletions.
323 
324  // use hierarchical cleaning
325  main->SetCleanup(kDeepCleanup);
326 
327  TGVerticalFrame* frame = new TGVerticalFrame(main,xsize,ysize,kVerticalFrame);
328 
329  // Object layout options
330  TGLayoutHints* fLayout1 = new TGLayoutHints(kLHintsLeft ,5,5,5,5);
331  TGLayoutHints* fLayout2 = new TGLayoutHints(kLHintsRight ,5,5,5,5);
332  TGLayoutHints* fLayout3 = new TGLayoutHints(kLHintsTop | kLHintsExpandX ,5,5,5,5);
333  TGLayoutHints* fLayout4 = new TGLayoutHints(kLHintsExpandX | kLHintsExpandY,5,5,5,5);
334 
335  //________ DO NOT MODIFY ABOVE _____________________________
336 
337 
338 
339 
340  // Title of the analysis
341  TGLabel* fAnalysisLabel = new TGLabel(frame, "Template analysis");
342 
343  // Create a line for choosing value of parameter
344  TGHorizontalFrame* rebinXFrame = new TGHorizontalFrame(frame,xsize,ysize2,kHorizontalFrame);
345  TGLabel* rebinXLabel = new TGLabel(rebinXFrame, "fRebinX");
346  fChooseRebinX = new TGNumberEntry(rebinXFrame);
347  fChooseRebinX->SetNumber(fRebinX);
348  TGHorizontalFrame* rebinZFrame = new TGHorizontalFrame(frame,xsize,ysize2,kHorizontalFrame);
349  TGLabel* rebinZLabel = new TGLabel(rebinZFrame, "fRebinZ");
350  fChooseRebinZ = new TGNumberEntry(rebinZFrame);
351  fChooseRebinZ->SetNumber(fRebinZ);
352 
353  main->AddFrame(frame,fLayout4);
354  frame->AddFrame(fAnalysisLabel,fLayout3);
355  frame->AddFrame(rebinXFrame,fLayout3);
356  rebinXFrame->AddFrame(rebinXLabel,fLayout1);
357  rebinXFrame->AddFrame(fChooseRebinX,fLayout2);
358  frame->AddFrame(rebinZFrame,fLayout3);
359  rebinZFrame->AddFrame(rebinZLabel,fLayout1);
360  rebinZFrame->AddFrame(fChooseRebinZ,fLayout2);
361 
362 
363 
364  //________ DO NOT MODIFY BELOW _____________________________
365  // Button to validate configuration
366  TGTextButton* setConf = new TGTextButton(frame,"Save Config",id);
367  setConf->Associate(fMain);
368 
369  frame->AddFrame(setConf,fLayout3);
370 
371  main->MapSubwindows();
372  main->MapWindow();
373  main->Resize();
374  return;
375 }
376 
377 
378 
379 void HarpoAnalysePattern::SetConfig()
380 {
381  // Update the configuration according to the values in the popup window
382 
383  if(!fChooseRebinX) return;
384 
385  fRebinX = fChooseRebinX->GetNumber();
386  fRebinZ = fChooseRebinZ->GetNumber();
387  for(Int_t i = 0; i<kMaxNvertex; i++){
388  fHistShortest[XPLANE][i]->Delete();
389  fHistShortest[YPLANE][i]->Delete();
390  fHistRoute[XPLANE][i]->Delete();
391  fHistRoute[YPLANE][i]->Delete();
392  fHistRoute2[XPLANE][i]->Delete();
393  fHistRoute2[YPLANE][i]->Delete();
394  fHistShortest[XPLANE][i] = new TH2I("fHistShortestX","",NADC/fRebinZ,0,NADC,NCHAN/fRebinX,0,NCHAN);
395  fHistShortest[YPLANE][i] = new TH2I("fHistShortestY","",NADC/fRebinZ,0,NADC,NCHAN/fRebinX,0,NCHAN);
396  fHistRoute[XPLANE][i] = new TH2I("fHistRouteX","",NADC/fRebinZ,0,NADC,NCHAN/fRebinX,0,NCHAN);
397  fHistRoute[YPLANE][i] = new TH2I("fHistRouteY","",NADC/fRebinZ,0,NADC,NCHAN/fRebinX,0,NCHAN);
398  fHistRoute2[XPLANE][i] = new TH2I("fHistRoute2X","",NADC/fRebinZ,0,NADC,NCHAN/fRebinX,0,NCHAN);
399  fHistRoute2[YPLANE][i] = new TH2I("fHistRoute2Y","",NADC/fRebinZ,0,NADC,NCHAN/fRebinX,0,NCHAN);
400  }
401 }
402 
403 
404 void HarpoAnalysePattern::Save(char * /* mode */ )
405  {
406 
407  OpenHistFile("pattern");
408 
409  // Save histograms here
410 
411  }
HarpoEvent::GetRecoEvent
HarpoRecoEvent * GetRecoEvent()
Definition: HarpoEvent.h:46
XPLANE
#define XPLANE
Definition: HarpoConfig.h:24
HarpoAnalysePattern::kMaxNvertex
static const Int_t kMaxNvertex
Definition: HarpoAnalysePattern.h:45
NCHAN
#define NCHAN
Definition: HarpoDccMap.h:16
HarpoAnalysePattern::fHistRoute2
TH2I * fHistRoute2[2][kMaxNvertex]
Definition: HarpoAnalysePattern.h:49
HarpoRecoVertex::GetZvertex
Double_t GetZvertex()
Definition: HarpoRecoEvent.h:288
HarpoDccMap::GetData
Double_t GetData(Int_t i, Int_t j)
Set/Get Data Cell.
Definition: HarpoDccMap.cxx:84
main
int main(int argc, char **argv)
Definition: drawAnalyse.cxx:25
HarpoDetSet::isExist
Bool_t isExist(ULong_t det)
Detecror date exist //! Number of Real Detectors.
Definition: HarpoDetSet.h:33
HarpoRecoEvent
Object containing the reconstruction information for one event (with array of HarpoRecoClusters Harpo...
Definition: HarpoRecoEvent.h:474
HarpoRecoVertex
2D vertex object, containing position, angle and associated track numbers, and quality info ...
Definition: HarpoRecoEvent.h:263
HarpoDccMap::print
virtual void print()
Definition: HarpoDccMap.cxx:124
HarpoAnalysePattern::GetRoute
void GetRoute(TH2I *hRoute, TH2I *h, Int_t length, Int_t i, Int_t j)
Definition: HarpoAnalysePattern.cxx:187
HarpoEventHeader
FullEvent Header not scecific to the detectors The class is ....
Definition: HarpoEventHeader.h:27
HarpoRecoVertex::GetPzVertex
Double_t GetPzVertex()
Definition: HarpoRecoEvent.h:290
HarpoAnalysePattern::GetShortest
void GetShortest(TH2I *h, Int_t length, Int_t i, Int_t j)
Redefine empty default.
Definition: HarpoAnalysePattern.cxx:168
HarpoAnalysePattern::fChooseRebinZ
TGNumberEntry * fChooseRebinZ
Definition: HarpoAnalysePattern.h:52
HarpoRecoEvent::GetVertexArray
TClonesArray * GetVertexArray()
Definition: HarpoRecoEvent.h:606
HarpoEventHeader::print
virtual void print()
Definition: HarpoEventHeader.cxx:21
HarpoAnalysePattern::fHistShortest
TH2I * fHistShortest[2][kMaxNvertex]
Definition: HarpoAnalysePattern.h:47
HarpoRecoClusters
Cluster object, containing position, charge and quality information.
Definition: HarpoRecoEvent.h:28
HarpoAnalysePattern::DisplayAnalysis
void DisplayAnalysis(TRootEmbeddedCanvas *ecTab, TGListBox *infobox)
Definition: HarpoAnalysePattern.cxx:236
HarpoAnalyse::OpenHistFile
TFile * OpenHistFile(const char *ananame)
Definition: HarpoAnalyse.cxx:150
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
MakeNiceHisto
TVirtualPad * MakeNiceHisto(TH1 *hist, TVirtualPad *c1, const char *opt, Bool_t copy)
Definition: MakeNiceHisto.cxx:221
HarpoAnalysePattern::Save
void Save(char *mode=NULL)
Definition: HarpoAnalysePattern.cxx:404
NADC
#define NADC
Definition: HarpoDccMap.h:18
YPLANE
#define YPLANE
Definition: HarpoConfig.h:25
HarpoConfigFile::Lookup
Bool_t Lookup(const char *path, Bool_t &val)
Lookup function for scalar values.
Definition: HarpoConfigFile.cxx:85
HarpoAnalyse::fEvt
HarpoEvent * fEvt
Definition: HarpoAnalyse.h:70
HarpoAnalysePattern::fChooseRebinX
TGNumberEntry * fChooseRebinX
Definition: HarpoAnalysePattern.h:51
HarpoAnalyse::nEvents
ULong_t nEvents
Definition: HarpoAnalyse.h:75
HarpoAnalysePattern
Dummy analysis to run as test and example. Give basic histograms of the data.
Definition: HarpoAnalysePattern.h:22
HarpoRecoVertex::GetPxVertex
Double_t GetPxVertex()
Definition: HarpoRecoEvent.h:289
gkNDetectors
const ULong_t gkNDetectors
Definition: HarpoDet.h:14
gHConfig
HarpoConfig * gHConfig
Definition: HarpoConfig.cxx:101
HarpoEvent::GetDccMap
HarpoDccMap * GetDccMap(Long_t plane=XDCC)
Definition: HarpoEvent.cxx:108
HarpoAnalysePattern::fHistRoute
TH2I * fHistRoute[2][kMaxNvertex]
Definition: HarpoAnalysePattern.h:48
HarpoRecoVertex::GetXvertex
Double_t GetXvertex()
Definition: HarpoRecoEvent.h:287
HarpoRecoEvent::GetClustersArray
TClonesArray * GetClustersArray()
Definition: HarpoRecoEvent.h:594
gHDetSet
R__EXTERN HarpoDetSet * gHDetSet
Definition: HarpoDetSet.h:71
HarpoAnalysePattern::ConfigFrame
void ConfigFrame(TGMainFrame *fMain, Int_t id)
Definition: HarpoAnalysePattern.cxx:312
Generated by   doxygen 1.8.6