HARPO  5.1.1
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
HarpoAnalyseBertrand.cxx
Go to the documentation of this file.
1 //
2 // File HarpoAnalyseBertrand.cxx
3 //
10 #include "HarpoAnalyseBertrand.h"
11 #include "HarpoConfig.h"
12 #include "HarpoDetSet.h"
13 #include "HarpoDebug.h"
14 #include "HarpoDccEvent.h"
15 #include "Pmm2Event.h"
16 #include "HarpoEvent.h"
17 #include "HarpoRecoEvent.h"
18 #include "HarpoSimEvent.h"
19 #include "TpcSimIonisationTrack.h"
20 #include "MakeNiceHisto.h"
21 
22 #include "TFile.h"
23 #include "TStyle.h"
24 #include "TCanvas.h"
25 #include "TRootEmbeddedCanvas.h"
26 #include "TGListBox.h"
27 #include "TLatex.h"
28 #include "TGraph.h"
29 #include "TF1.h"
30 #include "TTree.h"
31 #include "TRandom.h"
32 #include "TMath.h"
33 #include "TGLabel.h"
34 #include "TSystem.h"
35 
36 #include <vector>
37 #include <cstdlib>
38 #include <cstring>
39 #include <cassert>
40 #include <fstream>
41 #include <iostream>
42 
43 ClassImp(HarpoAnalyseBertrand)
44 
45 
46 void HarpoAnalyseBertrand::print()
47  {
48  unsigned int nd; // number of detectors
49  HarpoEventHeader *hdr = fEvt->GetHeader();
50 
51  assert(hdr != NULL);
52  hdr->print();
53 
54  for (nd = 0; nd < gkNDetectors; nd++) {
55  // if (fEvt->isdataExist(nd)) {
56  HarpoDccMap *plane = fEvt->GetDccMap(nd);
57  if (plane != NULL )plane->print();
58  }
59  }
60 
62 {
63 
64 
65 
66 
67 /* Raw Charge pour X & Y*/
68 
69  // Bool_t drawEvent = kFALSE;
70  nEvents++;
71 
72  // // Process RAW data
73  // for(Int_t ndet=0;ndet<2;ndet++) { // Detectors 0 and 1 = X and Y planes
74  // if (!gHDetSet->isExist(ndet)) continue;
75  // HarpoDccMap* m = fEvt->GetDccMap(ndet);
76  // if ( m == NULL ) continue;
77  // for(Int_t j=0;j<NCHAN;j++){ // Channels
78  // Double_t qch = 0;
79  // for(Int_t i=0;i<NADC;i++){ // Time bins
80  // Double_t q = m->GetData(j,i);
81  // if(q <= -1000) continue;
82  // hQraw[ndet]->Fill(q);
83  // }
84  // }
85  // }
86 
87 
88 
89 
90 // /* PPM2 Charge par Chan */
91 
92 // if(gHDetSet->isExist(PMM2)) { // PMm2 (trigger) data
93 // Pmm2Event *anaEvt = static_cast<Pmm2Event *>(fEvt->GetDetEvent(PMM2));
94 // Int_t triggerType = anaEvt->GetTriggerType(); // Trigger line (with mesh, without, traversing tracks, ...)
95 
96 // int imes;
97 // int esize = anaEvt->GetHeader()->eventSize; // number of PMs hit
98 // Pmm2MesVect * pm = anaEvt->GetMesurements(); // data for each included PM
99 // for(imes=0;imes<esize;imes++) {
100 // int ch = pm->at(imes).getChNum(); // channel
101 // int q = pm->at(imes).getCharge(); // charge
102 // hPmm2->Fill(ch,q);
103 // }
104 // }
105 
106 
107 
108 
109 
110 /* Charge par Cluster */
111 
112  if(fEvt->GetRecoEvent()){ // Reconstructed data (clusters, tracks, matched tracks)
113  HarpoRecoEvent* reco = fEvt->GetRecoEvent();
114  // Clusters in the event (both X and Y directions)
115  TClonesArray* clArray = reco->GetClustersArray();
116  Int_t nCl = clArray->GetEntries();
117  // Tracks in the event (both X and Y directions)
118  // TClonesArray* trArray = reco->GetTracksArray();
119  // Int_t nTr = trArray->GetEntries();
120 
121  qmyn = 0;
122  Double_t qmync = 0;
123  Int_t nTmp =0;
124 
125  for(Int_t icl = 0; icl<nCl; icl++){
126  HarpoRecoClusters* cluster = (HarpoRecoClusters*)clArray->At(icl);
127  Double_t q = cluster->GetQ();
128  Int_t type = cluster->GetType();
129  Int_t plane = cluster->GetPlane();
130  if(type==0) continue;
131  hQ[plane]->Fill(q);
132  if(plane==1) continue;
133  qmyn+=q;
134  nTmp++;
135  }
136 
137  for(Int_t icl = 0; icl<nCl; icl++){
138  HarpoRecoClusters* cluster = (HarpoRecoClusters*)clArray->At(icl);
139  Double_t q = cluster->GetQ();
140  Int_t type = cluster->GetType();
141  Int_t plane = cluster->GetPlane();
142  if(type==0 || plane==1) continue;
143  qmync+=((q-(qmyn/nTmp))*(q-(qmyn/nTmp)));
144  }
145  qsig=TMath::Sqrt((qmync/(nTmp-1)));
146  hQxy->Add(hQ[0], hQ[1], 1, 1);
147 
148 
149  }
150 
151 
152 /* Angle */
153 
154  // for(Int_t itr = 0; itr<nTr; itr++){
155  // HarpoRecoTracks* track = (HarpoRecoTracks*)trArray->At(itr);
156  // Double_t id = track->GetNtrack();
157  // Double_t angle = track->GetAngleTrack();
158  // hAngle->Fill(angle);
159  // }
160  // }
161 
162 
163 
164 
165 /* Charge par Angle */
166 
167  // if(fEvt->GetRecoEvent()){ // Reconstructed data (clusters, tracks, matched tracks)
168  // HarpoRecoEvent* reco = fEvt->GetRecoEvent();
169  // TClonesArray* clArray = reco->GetClustersArray();
170  // Int_t nCl = clArray->GetEntries();
171  // TClonesArray* trArray = reco->GetTracksArray();
172  // Int_t nTr = trArray->GetEntries();
173  // for(Int_t jtr = 0; jtr<nTr; jtr++){
174  // HarpoRecoTracks* track = (HarpoRecoTracks*)trArray->At(jtr);
175  // Int_t trpl = track->GetPlane();
176  // Double_t trid = track->GetNtrack();
177  // for(Int_t jcl = 0; jcl<nCl; jcl++){
178  // HarpoRecoClusters* cluster = (HarpoRecoClusters*)clArray->At(jcl);
179  // Int_t clpl = cluster->GetPlane();
180  // Double_t cltr = cluster->GetIdClusterTrack();
181  // Int_t type = cluster->GetType();
182 
183  // if(clpl==trpl && cltr==trid && type==1){
184  // Double_t angle = track->GetAngleTrack();
185  // Double_t q = cluster->GetQ();
186  // hAQ[clpl]->Fill(angle,q);
187  // }
188  // }
189  // }
190  // }
191 
192 
193 
194 
195 /* r et P(r) pour MCEvent */
196 
197  if (gHDetSet->isExist(SIMDET)){
199  TpcSimMCEvent* mcEvent = simEvt->GetMCEvent();
200  Int_t nMCtr = mcEvent->GetNtracks();
201  for(Int_t i = 0; i<nMCtr; i++){
202  TpcSimMCTrack* track = mcEvent->GetTrack(i);
203 
204  //Double_t r[3] = {track->GetX0(),track->GetY0(),-track->GetZ0());
205  //Double_t Pr[3] = {track->GetPx(),track->GetPy(),-track->GetPz());
206  Double_t xx = track->GetX0();
207  Double_t yy = track->GetY0();
208  Double_t zz = track->GetZ0();
209  // Double_t px = track->GetPx() ;
210  // Double_t py = track->GetPy();
211  // Double_t pz = track->GetPz();
212  // Double_t R = sqrt(zz+yy+xx);
213  // Double_t P = sqrt(pz*pz+py*py+px*px);
214  hxyz->Fill(xx,yy,-zz);
215  // hPxyz->Fill(px,py,-pz);
216  // hPx->Fill(px,xx);
217  }
218  }
219 
220 
221 
222 
223 
224  /* Ne- pour SimEvent */
225 
226  if (gHDetSet->isExist(SIMDET)){
228  Int_t nIonTr = simEvt->GetNtracks();
229  for(Int_t i = 0; i<nIonTr; i++){
230  TpcSimIonisationTrack* tr = simEvt->GetIonisationTrack(i);;
231  Int_t nPoints = tr->GetNpoints();
232 
233  emyn = 0;
234  Double_t emync = 0;
235 
236  for(Int_t j = 0; j<nPoints; j++){
237  TpcSimIonisationPoint* point = tr->GetPoint(j);
238  Double_t Ne = point->GetNe();
239  hNeSim->Fill(Ne);
240  emyn+=Ne;
241  }
242 
243  for(Int_t k=0; k<nPoints; k++){
244  TpcSimIonisationPoint* point = tr->GetPoint(k);
245  Double_t Ne = point->GetNe();
246  emync+=(Ne-(emyn/nPoints))*(Ne-(emyn/nPoints));
247  }
248  esig = TMath::Sqrt((emync/(nPoints-1)));
249 
250  }
251  }
252 
253 
254 
255 
256 
257  /*Dist X_recocluster - X_sim */
258  if (gHDetSet->isExist(SIMDET)){
260  Int_t nIonTr = simEvt->GetNtracks();
261 
262  HarpoRecoEvent* reco = fEvt->GetRecoEvent();
263  TClonesArray* clArray = reco->GetClustersArray();
264  Int_t nCl = clArray->GetEntries();
265 
266  for(Int_t zz = 0; zz<nCl; zz++){
267  HarpoRecoClusters* cluster = (HarpoRecoClusters*)clArray->At(zz);
268  Int_t zcl = cluster->GetIndex();
269  // Double_t q = cluster->GetQ();
270  Int_t type = cluster->GetType();
271  Int_t plane = cluster->GetPlane();
272  Double_t cx = cluster->GetMean();
273  cx += 1.2;
274  zcl = 24.5 + 0.920749*zcl;
275 
276  // zcl = 29.4122 + 0.877193*zcl;
277  if(type==1 && plane==0){
278  TEST->Fill(zcl,cx);
279  }
280  }
281 
282 
283  for(Int_t i = 0; i<nIonTr; i++){
284  TpcSimIonisationTrack* tr = simEvt->GetIonisationTrack(i);;
285  Int_t nPoints = tr->GetNpoints();
286 
287  for(Int_t n = 0; n<nPoints; n++){
288  TpcSimIonisationPoint* point = tr->GetPoint(n);
289  Double_t x = point->GetX();
290  x+=151;
291  Double_t y = point->GetY();
292  y+=151;
293  Double_t z = point->GetZ();
294  z+=89;
295 
296  TEST2->Fill(z,x);
297 
298  for(Int_t clz=0; clz<nCl; clz++){
299  HarpoRecoClusters* cluster = (HarpoRecoClusters*)clArray->At(clz);
300  Int_t zcl = cluster->GetIndex();
301  zcl = 24.5 + 0.920749*zcl;
302  Double_t q = cluster->GetQ();
303  Int_t type = cluster->GetType();
304  Int_t plane = cluster->GetPlane();
305  if(type==1 && plane==0 && TMath::Abs(zcl-z)<=0.25){
306  Double_t cx = cluster->GetMean();
307  // cx+=1.5;
308  cx += 1.2;
309  hclx->Fill(cx-x);
310  hclxa->Fill(cx-x,q);
311  }
312  }
313  }
314  }
315  }
316 
317 
318 
319 
320 
321  Double_t data[7];
322 
323 /*Dist X(puis Y) _reco - X(puis Y)_sim */
324  if (gHDetSet->isExist(SIMDET)){
326  Int_t nIonTr = simEvt->GetNtracks();
327  HarpoDccMap* m1 = fEvt->GetDccMap(0);
328  // HarpoDccMap* m2 = fEvt->GetDccMap(1);
329  for(Int_t i = 0; i<nIonTr; i++){
330  TpcSimIonisationTrack* tr = simEvt->GetIonisationTrack(i);;
331  Int_t nPoints = tr->GetNpoints();
332 
333  TpcSimMCEvent* mcEvent = simEvt->GetMCEvent();
334  TpcSimMCTrack* track = mcEvent->GetTrack(i);
335  Double_t px = track->GetPx();
336  Double_t py = track->GetPy();
337  Double_t pz = track->GetPz();
338  Double_t tx = px/pz;
339  Double_t ty = py/pz;
340  Ang[0] = tx;
341  Ang[1] = ty;
342 
343  for(Int_t n = 0; n<nPoints; n++){
344  TpcSimIonisationPoint* point = tr->GetPoint(n);
345  Double_t x = point->GetX();
346  x+=151;
347  Double_t y = point->GetY();
348  y+=151;
349  Double_t z = point->GetZ();
350  z+=80;
351  for(Int_t j=0;j<NCHAN;j++){ // Channels
352  Double_t qx = m1->GetData(j,z);
353  //Double_t qy = m2->GetData(j,z);
354  if (qx>=300 && TMath::Abs(tx)<=0.79){
355  hdxq->Fill(j-x,qx);
356  hdx->Fill(j-x);
357  hdxa->Fill(j-x,TMath::Cos(tx));
358  }
359  }
360 
361  for(Int_t k=0; k<NADC;k++){ // Time bins
362  Double_t qz1 = m1->GetData(x,k);
363  if (qz1>=300){
364  hdzq->Fill(k-z,qz1);
365  hdz->Fill(k-z);
366  hdza->Fill(k-z,TMath::Abs(TMath::Sin(tx)));
367  }
368  }
369 
370  // Cut Hist to better Mean/RMS
371  // Double_t cc = ((hxmxq->GetMaximum())/(3.5));
372  // for(Int_t n=5; n<395; n++){
373  // Double_t c = hxmxq->GetBinContent(n);
374  // if(c>=cc){
375  // hcut->SetBinContent(n,c);
376  // }
377  // }
378  // H -> GET FIRST/LAST BIN ABOVE (XX)....
379 
380  //Computation of RMS/Mean by hand, in a selected range [-5,5] : Delta X
381  Double_t L=0;
382  Double_t I=0;
383  Double_t II=0;
384 
385  for(Int_t p=200;p<=401;p++){
386  Double_t b = hdxq->GetBinContent(p);
387  L+= b;
388  I+= 0.05*(b*(p-300.5));
389  }
390  for(Int_t p=200;p<=401;p++){
391  Double_t b = hdxq->GetBinContent(p);
392  II+=b*((0.05*(p-300.5))-(I/L))*((0.05*(p-300.5))-(I/L));
393  }
394  Double_t si = TMath::Sqrt(II/L);
395 
396  //Computation of RMS/Mean by hand, in a selected range [-12,12] : Delta Z
397  Double_t zL=0;
398  Double_t zI=0;
399  Double_t zII=0;
400 
401  for(Int_t p=100;p<=201;p++){
402  Double_t b = hdzq->GetBinContent(p);
403  zL+= b;
404  zI+= 0.1*(b*(p-150.5));
405  }
406  for(Int_t p=100;p<=201;p++){
407  Double_t b = hdzq->GetBinContent(p);
408  zII+=b*((0.1*(p-150.5))-(zI/zL))*((0.1*(p-150.5))-(zI/zL));
409  }
410  Double_t zsi = TMath::Sqrt(zII/zL);
411 
412  // hxmxq->GetXaxis()->SetRangeUser(-7,7);
413  // hymyq->GetXaxis()->SetRangeUser(-7,7);
414  // hzmzq->GetXaxis()->SetRangeUser(-12,12);
415  // hcut->GetXaxis()->SetRangeUser(-5,5);
416  data[0]=z;
417  data[1]=(z+(zI/zL));
418  data[2]=zsi;
419  data[3]= x;
420  data[4]= (x + (I/L));
421  data[5]= si;
422  data[6]= tx;
423  if (data[3]<1 || data[0]<1) continue;
424  Ntp->Fill(data);
425  }
426  }
427  }
428 
429 
431 /* Angle RECO - Angle SIM */
432 
433 
434  if (gHDetSet->isExist(SIMDET)){
436  Int_t nIonTr = simEvt->GetNtracks();
437  HarpoRecoEvent* reco = fEvt->GetRecoEvent();
438  TClonesArray* trArray = reco->GetTracksArray();
439  Int_t nTr = trArray->GetEntries();
440 
441  for(Int_t i = 0; i<nIonTr; i++){
442  TpcSimMCEvent* mcEvent = simEvt->GetMCEvent();
443  TpcSimMCTrack* mtrack = mcEvent->GetTrack(i);
444  Double_t px = mtrack->GetPx();
445  Double_t py = mtrack->GetPy();
446  Double_t pz = mtrack->GetPz();
447  Double_t tx = (px/pz);
448  Double_t ty = (py/pz);
449 
450  if(nIonTr ==0 ) continue;
451 
452  for(Int_t itr = 0; itr<nTr; itr++){
453  HarpoRecoTracks* track = (HarpoRecoTracks*)trArray->At(itr);
454  Int_t pla = track->GetPlane();
455  Double_t angle = track->GetAngleTrack();
456  if(pla==0){
457  hresx->Fill(angle-tx);
458  hresxa->Fill(angle-tx,TMath::Abs(TMath::Cos(tx)));
459  }
460  if(pla==1){
461  hresy->Fill(angle-ty);
462  hresya->Fill(angle-ty,TMath::Abs(TMath::Cos(ty)));
463  }
464  }
465  }
466  }
467 
468 
469 
470 
471 
472 
473 
474 /* 3D Track pour SimEvent */
475 
476  // if (gHDetSet->isExist(SIMDET)){
477  // HarpoSimEvent *simEvt = (HarpoSimEvent *)(fEvt->GetDetEvent(SIMDET));
478  // Int_t nIonTr = simEvt->GetNtracks();
479  // for(Int_t i = 0; i<nIonTr; i++){
480  // TpcSimIonisationTrack* tr = simEvt->GetIonisationTrack(i);;
481  // Int_t nPoints = tr->GetNpoints();
482  // for(Int_t j = 0; j<nPoints; j++){
483  // TpcSimIonisationPoint* point = tr->GetPoint(j);
484  // Double_t x = point->GetX();
485  // x = (x+151);
486  // Double_t y = point->GetY();
487  // y = (y+151);
488  // Double_t z = point->GetZ();
489  // z = (z+80);
490  // h3->Fill(z,x,y);
491  // }
492  // }
493  // }
494 
495 
496 
497 
498  /* Ne- sommés par pas de dz=1mm pour SimEvent */
499 
500  if (gHDetSet->isExist(SIMDET)){
502  Int_t nIonTr = simEvt->GetNtracks();
503  for(Int_t i = 0; i<nIonTr; i++){
504  TpcSimIonisationTrack* tr = simEvt->GetIonisationTrack(i);;
505  Int_t nPoints = tr->GetNpoints();
506  if (nPoints==0) continue;
507 
508  TpcSimIonisationPoint* pnti = tr->GetPoint(0);
509  Double_t zi = pnti->GetZ();
510  TpcSimIonisationPoint* pntz = tr->GetPoint(nPoints-1);
511  Double_t zf = pntz->GetZ();
512 
513 
514  // TH1F* htest = new TH1F("htest","",600,-100,500);
515  // // Search End of ITrack
516  // if (nPoints==0) continue;
517  // for(Int_t k = 0; k<nPoints; k++){
518  // TpcSimIonisationPoint* ptk = tr->GetPoint(k);
519  // Double_t zk = ptk->GetZ();
520  // Double_t nek = ptk->GetNe();
521  // htest->Fill(zk);
522  // }
523 
524  // htest->Delete();
525 
526  // if(binf <2) break;
527  // Double_t zf = binf;
528 
529  for(Int_t az = zi; az<=zf; az++){
530  Double_t Nez = 0;
531  for(Int_t j = 0; j<nPoints; j++){
532  TpcSimIonisationPoint* point = tr->GetPoint(j);
533  Int_t z = point->GetZ();
534  if(z<az || z>=(az+1)) continue;
535  Double_t nz = point->GetNe();
536  Nez = Nez + nz;
537  }
538  hNez->Fill(Nez);
539  }
540  }
541  }
542 
543  // hSig->Fill(esig,qsig);
544  hMean->Fill(emyn,qmyn);
545  hRappMean->Fill(qmyn/emyn);
546  // hRappSig->Fill(qsig/esig);
547 
548 }
549 
550 
552 {
553 
554  fChooseNbins = 0;
555 
556  // Initialise histograms here
557  fNbins = 120;
558  Double_t fQmin = 10, fQmax = 1e7;
559 
560  Long64_t Nbins;
561  if( ! gHConfig->Lookup("template.nbins",Nbins ))
562  Info("Constructor","Use default fnbinr %d",fNbins);
563  else
564  fNbins = Nbins;
565 
566  Double_t qmin;
567  if( ! gHConfig->Lookup("template.qmin",qmin ))
568  Info("Constructor","Use default Qmin %.3g",fQmin);
569  else
570  fQmin = qmin;
571 
572  Double_t qmax;
573  if( ! gHConfig->Lookup("template.qmax",qmax ))
574  Info("Constructor","Use default Qmax %.3g",fQmax);
575  else
576  fQmax = qmax;
577 
578 
579  //hQ[XPLANE] = HistLog("hQx","",fNbins,fQmin,fQmax);
580  //hQ[YPLANE] = HistLog("hQy","",fNbins,fQmin,fQmax);
581  hQ[XPLANE] = new TH1F("Q_x","",500,0,2000);
582  hQ[YPLANE] = new TH1F("Q_y","",500,0,2000);
583  hQxy = new TH1F("hQ_xy","",500,0,2000);
584  // hQraw[XPLANE] = new TH1F("Qraw_X","",500,-500,5500);
585  // hQraw[YPLANE] = new TH1F("Qraw_Y","",500,-500,5500);
586  // hPmm2 = new TH2F("hPmm2","",16,0,16,128,0,4096);
587 
588  hNeSim = new TH1F("NeSim","",35,0,35);
589  hNez = new TH1F("Nez","",100,0,100);
590 
591  // hSig = new TH2F("SigmaQNe","",100,0,5000,10000,0,100000000);
592  hMean = new TH2F("MeanQNe","",100,0,20000,1000,0,10000000);
593  hRappMean = new TH1F("RapportMeanQNe","",50,150,350);
594  // hRappSig = new TH1F("RapportSigQNe","",300,0,800);
595 
596  // hAngle = new TH1F("Angle","",20,-5,5);
597  // hAQ[XPLANE] = new TH1F("Ang(q)_x","",50,-5,5);
598  // hAQ[YPLANE] = new TH1F("Ang(q)_y","",50,-5,5);
599 
600  hxyz = new TH3F("hxyz","",500,-200,200,500,-200,200,500,-200,200);
601  // hPx = new TH2F("hP_x","",500,-1,1,500,-200,200);
602  // hPxyz = new TH3F("hPxyz","",500,-1,1,500,-1,1,500,-1,1);
603 
604  // h3 = new TH3F("h3D","",512,0,512,288,0,288,288,0,288);
605 
606 
607  // hymy = new TH1F("hy-y","",400,-50,50);
608  //hymyq = new TH1F("hdeltayq","",400,-50,50);
609  //hymya = new TH2F("hdeltaya","",200,-1,1,400,-40,40);
610 
611  hdx = new TH1F("hdx","",600,-15,15);
612  hdxa = new TH1F("hdxa","",600,-15,15);
613  hdxq = new TH1F("hdxq","",600,-15,15);
614 
615  hdz = new TH1F("hdz","",100,-20,20);
616  hdza = new TH1F("hdza","",100,-20,20);
617  hdzq = new TH1F("hdzq","",100,-20,20);
618 
619  // hcut = new TH1F("hcut","",400,-50,50);
620  // hdeltaq = new TH2F("hdeltaq","",400,-20,20,400,00,3000);
621 
622  hclx = new TH1F("hclx","",300,-15,15);
623  hclxa = new TH1F("hclxa","",300,-15,15);
624 
625  hresx = new TH1F("hresx","",250,-0.25,0.25);
626  hresy = new TH1F("hresy","",250,-0.25,0.25);
627  hresxa = new TH1F("hresxa","",125,-0.25,0.25);
628  hresya = new TH1F("hresya","",125,-0.25,0.25);
629 
630  TEST = new TH2F("TEST","",512,0,512,302,0,302);
631  TEST2 = new TH2F("TEST2","",512,0,512,302,0,302);
632 
633 
634  Ntp = new TNtupleD("Ntp","","zs:zm:sigz:xs:xm:sigx:alpha");
635  Ntp->SetDirectory(0);
636 
637 }
638 
639 void HarpoAnalyseBertrand::Save(char * /* mode */)
640  {
641 
642  OpenHistFile("bertrand");
643 
644  // Save histograms here
645  // hQ[XPLANE]->Write();
646  // hQ[YPLANE]->Write();
647  hQxy->Write();
648  // hQraw[XPLANE]->Write();
649  // hQraw[YPLANE]->Write();
650 
651  hNeSim->Write();
652  hNez->Write();
653 
654  // hSig->Write();
655  hMean->Write();
656  hRappMean->Write();
657  // hRappSig->Write();
658 
659  // hAngle->Write();
660  // hPmm2->Write();
661  // hAQ[XPLANE]->Write();
662  // hAQ[YPLANE]->Write();
663 
664  //hPx->Write();
665  hxyz->Write();
666  //hPxyz->Write();
667 
668  // hymyq->Write();
669  // hymy->Write()
670  // hymya->Write();
671 
672  hdx->Write();
673  hdxa->Write();
674  hdxq->Write();
675 
676  hdz->Write();
677  hdza->Write();
678  hdzq->Write();
679 
680  // hcut->Write();
681  // hdeltaq->Write();
682 
683  hclx->Write();
684  hclxa->Write();
685 
686  hresx->Write();
687  hresy->Write();
688  hresxa->Write();
689  hresya->Write();
690 
691  TEST->Write();
692  TEST2->Write();
693 
694  Ntp->Write();
695  // h3->Write();
696 }
697 
698 
699 
700 
701 void HarpoAnalyseBertrand::DisplayAnalysis(TRootEmbeddedCanvas* ecTab, TGListBox* infobox)
702 {
703  // in hrecomonitorgui: fills the canvas on the "Analysis" tab
704  TCanvas* c = ecTab->GetCanvas();
705  c->GetPad(1)->Delete();
706  c->GetPad(2)->Delete();
707  c->Divide(2);
708 
709  // c->GetPad(1)->Divide(2);
710  // c->GetPad(2)->Divide(2);
711  //c->GetPad(2)->Divide(2,2);
712 
713  // c->GetPad(1)->SetPad(0., 0., 0.40, 1);
714  // c->GetPad(2)->SetPad(0.41, 0, 1, 1 );
715 
716 
717 
718  infobox->NewEntry("");
719  infobox->NewEntry("Left:");
720 
721 
722 
723 
724  TH2F* ha = new TH2F("ha","",512,0,512,288,0,288);
725  TH2F* hb = new TH2F("hb","",512,0,512,288,0,288);
726  //
727  TH2F* h1 = new TH2F("h1","",512,0,512,288,0,288);
728  TH2F* h2 = new TH2F("h2","",512,0,512,288,0,288);
729  //
730  TH2F* hX = new TH2F("hX","",512,0,512,288,0,288);
731  TH2F* hY = new TH2F("hY","",512,0,512,288,0,288);
732 
733  // TH2F* hBB = new TH2F("hBB","",512,0,512,288,0,288);
734 
735  //
736  // // TH3F* h3 = new TH3F("h3D","",512,0,512,288,0,288,288,0,288);
737  // // TH3F* h3bin = new TH3F("h3Dbin","",512,0,512,288,0,288,288,0,288);
738  // // TH3F* hreco = new TH3F("h3Dreco","",512,0,512,288,0,288,288,0,288);
739  // //
740 
741 
742 
743 
744 /* DccMap et remplissage (X puis Y) */
745 
746  HarpoDccMap* m1 = fEvt->GetDccMap(0);
747  for(Int_t i=1;i<NADC;i++){ //Time bins
748  for(Int_t j=0;j<NCHAN;j++){ // Channels
749  Double_t q = m1->GetData(j,i);
750  ha->SetBinContent(i+1,j+1,q);
751  }
752  }
753  ha->SetMinimum(0);
754 
755  HarpoDccMap* m2 = fEvt->GetDccMap(1);
756  for(Int_t i=1;i<NADC;i++){ //Time bins
757  for(Int_t j=0;j<NCHAN;j++){ // Channels
758  Double_t q = m2->GetData(j,i);
759  hb->SetBinContent(i+1,j+1,q);
760  }
761  }
762  hb->SetMinimum(0);
763 
764 
765 
766  // // // for(Int_t clz=0; clz<nCl; clz++){
767  // // // HarpoRecoClusters* cluster = (HarpoRecoClusters*)clArray->At(clz);
768  // // // Int_t zcl = cluster->GetIndex();
769  // // // Double_t q = cluster->GetQ();
770  // // // q+=100;
771  // // // Int_t type = cluster->GetType();
772  // // // Int_t plane = cluster->GetPlane();
773  // // // Double_t cx = cluster->GetMean();
774  // // // Double_t w = cluster->GetWidth();
775  // // // infobox->NewEntry(Form("NCL = %d",nCl));
776  // // // if(plane==0 && type==1){
777  // // // infobox->NewEntry(Form("I %d ; M %g",zcl,cx));
778  // // // }
779  // // // hBB->SetBinContent(zcl+1,cx+1,q);
780 
781  // // // }
782  // // // hBB->SetMinimum(0);
783 
784 
785 
786 /* Sim Track avec IonPoint (XZ,YZ,XYZ) */
787 
788  if (gHDetSet->isExist(SIMDET)){
790  Int_t nIonTr = simEvt->GetNtracks();
791  for(Int_t i = 0; i<nIonTr; i++){
792  TpcSimIonisationTrack* tr = simEvt->GetIonisationTrack(i);;
793  Int_t nPoints = tr->GetNpoints();
794  for(Int_t j = 0; j<nPoints; j++){
795  TpcSimIonisationPoint* point = tr->GetPoint(j);
796  Double_t x = point->GetX();
797  x+= 151;
798  Double_t y = point->GetY();
799  y+= 151;
800  Double_t z = point->GetZ();
801  z+= 89;
802  h1->SetBinContent(z,x,5000);
803  h2->SetBinContent(z,y,5000);
804  // h3bin->SetBinContent(z,x,y,5000);
805  // h3->Fill(z,x,y);
806  }
807  }
808  h1->SetMinimum(0);
809  h2->SetMinimum(0);
810  }
811 
812 
813  //
814  //MakeNiceHisto(h3,c->GetPad(2)->GetPad(1),"");
815  //
816  /*Dist X(puis Y) _reco - X(puis Y)_sim */
817  // if (gHDetSet->isExist(SIMDET)){
818  // HarpoSimEvent *simEvt = (HarpoSimEvent *)(fEvt->GetDetEvent(SIMDET));
819  // Int_t nIonTr = simEvt->GetNtracks();
820  // for(Int_t i = 0; i<nIonTr; i++){
821  // TpcSimIonisationTrack* tr = simEvt->GetIonisationTrack(i);;
822  // Int_t nPoints = tr->GetNpoints();
823  // for(Int_t n = 0; n<nPoints; n++){
824  // TpcSimIonisationPoint* point = tr->GetPoint(n);
825  // Double_t x = point->GetX();
826  // x = (x+151);
827  // Double_t y = point->GetY();
828  // y = (y+151);
829  // Double_t z = point->GetZ();
830  // z = (z+80);
831  // for(Int_t j=0;j<NCHAN;j++){ // Channels
832  // Double_t qx = m1->GetData(j,z);
833  // Double_t qy = m2->GetData(j,z);
834  // if (qx>=100){
835  // hxmx->Fill(j-y);
836  // }
837  // if (qy>=100){
838  // hymy->Fill(j-x);
839  // }
840  // }
841  // }
842  // }
843  // }
844 
845 
846 
847 
848  /* Histo 3D : Bins reconstruits s'ils sont dans l'entourage d'un bin de l'IonTrack */
849  //
850  // HarpoDccMap* m1 = fEvt->GetDccMap(0);
851  // HarpoDccMap* m2 = fEvt->GetDccMap(1);
852  // for(Int_t i=1;i<NADC;i++){ //Time bins
853  // for(Int_t j=3;j<NCHAN;j++){ // Channels
854  // for(Int_t k=3;k<NCHAN;k++){
855  // Double_t q2 = m2->GetData(k,i);
856  // Double_t q1 = m1->GetData(j,i);
857  // Double_t q3 = 0;
858  // Double_t zero = 0;
859  // for(Int_t p=-2; p<3;p++){
860  // q3+=h3bin->GetBinContent(i+p,k,j);
861  // q3+=h3bin->GetBinContent(i,k+p,j);
862  // q3+=h3bin->GetBinContent(i,k,j+p);
863  // q3+=h3bin->GetBinContent(i+p,k+p,j);
864  // q3+=h3bin->GetBinContent(i+p,k,j+p);
865  // q3+=h3bin->GetBinContent(i+p,k+p,j-p);
866  // q3+=h3bin->GetBinContent(i+p,k+p,j+p);
867  // }
868  // if ( q3<15000) continue;
869  // Double_t q4 = (TMath::Max(zero,m1->GetData(j,i))+TMath::Max(zero,m2->GetData(k,i)));
870  // hreco->SetBinContent(i+1,k+1,j+1,q4);
871  // }
872  // }
873  // }
874  //MakeNiceHisto(hreco,c->GetPad(2)->GetPad(2),"col2");
875 // //
876 
877 
878 
879 
880 /* Track X (puis Y) Sim ET Dcc */
881 
882  for(Int_t i = 0; i <=512 ; i++ ){
883  for(Int_t j = 0; j <=288 ; j++ ){
884  Double_t zero = 0;
885 
886  Double_t na = TMath::Max(zero,ha->GetBinContent(i,j));
887  Double_t nb = TMath::Max(zero,hb->GetBinContent(i,j));
888 
889  Double_t n1 = h1->GetBinContent(i,j);
890  Double_t n2 = h2->GetBinContent(i,j);
891 
892  hX->SetBinContent(i,j,100*(n1+na));
893  hY->SetBinContent(i,j,100*(n2+nb));
894  }
895  }
896  hX->SetMinimum(0);
897  hY->SetMinimum(0);
898 
899  // // MakeNiceHisto(hxmx,c->GetPad(2)->GetPad(1),"");
900  // // MakeNiceHisto(hymy,c->GetPad(2)->GetPad(2),"");
901 
902 
903  // infobox->NewEntry("-L: Track X(Dcc) & Track Y(IonTr) [2D]");
904  // infobox->NewEntry("-R: Track Y(Dcc) & Track X(IonTr) [2D]");
905  hX->GetXaxis()->SetTitle("Time bin");
906  hX->GetYaxis()->SetTitle("Spatial Channel");
907  hY->GetXaxis()->SetTitle("Time bin");
908  hY->GetYaxis()->SetTitle("Spatial Channel");
909  MakeNiceHisto(hX,c->GetPad(1),"colz");
910  MakeNiceHisto(hY,c->GetPad(2),"colz");
911 
912  hX->Delete();
913  hY->Delete();
914  ha->Delete();
915  hb->Delete();
916  h1->Delete();
917  h2->Delete();
918  //hBB->Delete();
919 
920  // hxmx->Delete();
921  // hymy->Delete();
922 
923  // h3->Delete();
924  // hreco->Delete();
925  // h3bin->Delete();
926 
927 
928 
929  infobox->NewEntry("Right");
930 
931  // infobox->NewEntry("Lfar: x0(Px) [2D]");
932  // hPx->GetYaxis()->SetTitle("x0");
933  // hPx->GetXaxis()->SetTitle("Proj.x P norm");
934  // MakeNiceHisto(hPx,c->GetPad(1)->GetPad(1),"box");
935 
936  // infobox->NewEntry("Lmid: Angle Distribution [1D]");
937  // hAngle->GetYaxis()->SetTitle("Tracks (X & Y)");
938  // hAngle->GetXaxis()->SetTitle("Angle [rad]");
939  // MakeNiceHisto(hAngle,c->GetPad(1)->GetPad(2),"");
940 
941  // infobox->NewEntry("Lfar: Sum of Cluster Charge Distribution [1D]");
942  // hQxy->GetYaxis()->SetTitle("Entries");
943  // hQxy->GetXaxis()->SetTitle("Q_x + Q_y [ADC]");
944  // MakeNiceHisto(hQxy,c->GetPad(1)->GetPad(1),"box");
945 
946  // infobox->NewEntry("Lmid: Angle Distribution [1D]");
947  // hQ[YPLANE]->GetYaxis()->SetTitle("Entries");
948  // hQ[YPLANE]->GetXaxis()->SetTitle("Q_y [ADC]");
949  // MakeNiceHisto(hQ[YPLANE],c->GetPad(1)->GetPad(2),"");
950 
951  // infobox->NewEntry("-TL: Cluster Charge Distribution [1D]");
952  // hQ[XPLANE]->GetYaxis()->SetTitle("Entries");
953  // hQ[XPLANE]->GetXaxis()->SetTitle("Q_x [ADC]");
954  // MakeNiceHisto(hQ[XPLANE],c->GetPad(2)->GetPad(1),"");
955 
956  // infobox->NewEntry("-TR: Sim e- Distribution for dz=1mm by IonTrack [1D]");
957  // hNez->GetYaxis()->SetTitle("Entries");
958  // hNez->GetXaxis()->SetTitle("e- Simulated @ z");
959  // MakeNiceHisto(hNez,c->GetPad(2)->GetPad(2),"");
960 
961  // infobox->NewEntry("Rmid: e- Simulated Distribution on L loop by Ionisation Track [1D]");
962  // hNeL->GetYaxis()->SetTitle("Entries");
963  // hNeL->GetXaxis()->SetTitle("e- Simulated @ L");
964  // MakeNiceHisto(hNeL,c->GetPad(2)->GetPad(1),"");
965 
966  // infobox->NewEntry("-BL: Ratio (Sum Q)/(Sum Ne) Distribution [1D]");
967  // hRappMean->GetYaxis()->SetTitle("Entries");
968  //hRappMean->GetXaxis()->SetTitle("Q/Ne");
969  // MakeNiceHisto(hRappMean,c->GetPad(2)->GetPad(3),"");
970 
971  // infobox->NewEntry("Rfar: e- Simulated Distribution on each Ionisation Point [1D]");
972  // hNeSim->GetYaxis()->SetTitle("Entries");
973  // hNeSim->GetXaxis()->SetTitle("e- Simulated");
974  // MakeNiceHisto(hNeSim,c->GetPad(2)->GetPad(2),"");
975 
976  // infobox->NewEntry("-BR: Angle by Charge (X Track) [2D]");
977  // hAQ[0]->GetYaxis()->SetTitle("Charge [ADC]");
978  // hAQ[0]->GetXaxis()->SetTitle("Angle [rad]");
979  // MakeNiceHisto(hAQ[0],c->GetPad(2)->GetPad(2),"col");
980 
981  infobox->NewEntry("");
982  infobox->NewEntry(Form("Event # %lu",nEvents));
983  infobox->NewEntry(Form("Rapport Q/Ne- : %g",(qmyn/emyn)));
984  infobox->NewEntry(Form("Angle 1 : %g ; Angle 2 : %g ",(Ang[0]),(Ang[1])));
985  infobox->NewEntry(Form("NADC : %u ; NCHAN : %u ",NADC,NCHAN));
986 
987  c->Update();
988 
989 }
990 void HarpoAnalyseBertrand::ConfigFrame(TGMainFrame* fMain, Int_t id)
991 {
992  // in hrecomonitorgui: Creates a popup window for analysis configuration
993  // You must define SetConfig() properly
994 
995  UInt_t xsize = 200;
996  UInt_t ysize2 = 20;
997  UInt_t ysize = 10*ysize2;
998  TGTransientFrame* main = new TGTransientFrame(gClient->GetRoot(), fMain, xsize, ysize);
999  main->Connect("CloseWindow()", "HarpoRecoMonitorGui", main, "CloseWindow()");
1000  main->DontCallClose(); // to avoid double deletions.
1001 
1002  // use hierarchical cleaning
1003  main->SetCleanup(kDeepCleanup);
1004 
1005  TGVerticalFrame* frame = new TGVerticalFrame(main,xsize,ysize,kVerticalFrame);
1006 
1007  // Object layout options
1008  TGLayoutHints* fLayout1 = new TGLayoutHints(kLHintsLeft ,5,5,5,5);
1009  TGLayoutHints* fLayout2 = new TGLayoutHints(kLHintsRight ,5,5,5,5);
1010  TGLayoutHints* fLayout3 = new TGLayoutHints(kLHintsTop | kLHintsExpandX ,5,5,5,5);
1011  TGLayoutHints* fLayout4 = new TGLayoutHints(kLHintsExpandX | kLHintsExpandY,5,5,5,5);
1012 
1013  //________ DO NOT MODIFY ABOVE _____________________________
1014 
1015 
1016 
1017 
1018  // Title of the analysis
1019  TGLabel* fAnalysisLabel = new TGLabel(frame, "Template analysis");
1020 
1021  // Create a line for choosing value of parameter
1022  TGHorizontalFrame* nBinsFrame = new TGHorizontalFrame(frame,xsize,ysize2,kHorizontalFrame);
1023  TGLabel* nBinsLabel = new TGLabel(nBinsFrame, "fNbins");
1024  fChooseNbins = new TGNumberEntry(nBinsFrame);
1025  fChooseNbins->SetNumber(fNbins);
1026 
1027  main->AddFrame(frame,fLayout4);
1028  frame->AddFrame(fAnalysisLabel,fLayout3);
1029  frame->AddFrame(nBinsFrame,fLayout3);
1030  nBinsFrame->AddFrame(nBinsLabel,fLayout1);
1031  nBinsFrame->AddFrame(fChooseNbins,fLayout2);
1032 
1033 
1034 
1035  //________ DO NOT MODIFY BELOW _____________________________
1036  // Button to validate configuration
1037  TGTextButton* setConf = new TGTextButton(frame,"Save Config",id);
1038  setConf->Associate(fMain);
1039 
1040  frame->AddFrame(setConf,fLayout3);
1041 
1042  main->MapSubwindows();
1043  main->MapWindow();
1044  main->Resize();
1045  return;
1046 }
1047 
1048 
1049 
1051 {
1052  // Update the configuration according to the values in the popup window
1053  if(!fChooseNbins) return;
1054  fNbins = fChooseNbins->GetNumber();
1055 }
void Save(char *mode=NULL)
TpcSimIonisationTrack * GetIonisationTrack(Int_t iTr)
HarpoRecoEvent * GetRecoEvent()
Definition: HarpoEvent.h:46
#define XPLANE
Definition: HarpoConfig.h:24
Double_t GetY0()
Definition: TpcSimMCEvent.h:24
#define NCHAN
Definition: HarpoDccMap.h:16
Double_t GetMean()
Double_t GetPz()
Definition: TpcSimMCEvent.h:32
Int_t GetNtracks()
Definition: HarpoSimEvent.h:54
Double_t GetData(Int_t i, Int_t j)
Set/Get Data Cell.
Definition: HarpoDccMap.cxx:84
int main(int argc, char **argv)
Definition: drawAnalyse.cxx:25
Double_t GetZ0()
Definition: TpcSimMCEvent.h:25
Bool_t isExist(ULong_t det)
Detecror date exist //! Number of Real Detectors.
Definition: HarpoDetSet.h:33
Double_t GetAngleTrack()
Track object, containing position, angle, charge and quality information.
Object containing the reconstruction information for one event (with array of HarpoRecoClusters Harpo...
virtual void print()
Double_t GetPy()
Definition: TpcSimMCEvent.h:31
TpcSimIonisationPoint * GetPoint(Int_t i)
FullEvent Header not scecific to the detectors The class is ....
A class store HARPO row DCC event data and header. End provide access metods to the row data...
Definition: HarpoSimEvent.h:24
virtual void print()
TGNumberEntry * fChooseNbins
Cluster object, containing position, charge and quality information.
TFile * OpenHistFile(const char *ananame)
static int type
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
void DisplayAnalysis(TRootEmbeddedCanvas *ecTab, TGListBox *infobox)
Redefine empty default.
void ConfigFrame(TGMainFrame *fMain, Int_t id)
HarpoDetEvent * GetDetEvent(Long_t plane=XDCC)
Definition: HarpoEvent.cxx:93
Data from Keller temperuture and pressure sensors.
Definition: HarpoDet.h:22
TVirtualPad * MakeNiceHisto(TH1 *hist, TVirtualPad *c1, const char *opt, Bool_t copy)
Int_t GetNtracks()
TpcSimMCTrack * GetTrack(Int_t i)
Definition: TpcSimMCEvent.h:91
#define NADC
Definition: HarpoDccMap.h:18
Double_t GetX0()
Definition: TpcSimMCEvent.h:23
#define YPLANE
Definition: HarpoConfig.h:25
TpcSimMCEvent * GetMCEvent()
Definition: HarpoSimEvent.h:57
Bool_t Lookup(const char *path, Bool_t &val)
Lookup function for scalar values.
HarpoEvent * fEvt
Definition: HarpoAnalyse.h:70
ULong_t nEvents
Definition: HarpoAnalyse.h:75
const ULong_t gkNDetectors
Definition: HarpoDet.h:14
HarpoConfig * gHConfig
TClonesArray * GetTracksArray()
HarpoDccMap * GetDccMap(Long_t plane=XDCC)
Definition: HarpoEvent.cxx:108
Double_t GetPx()
Definition: TpcSimMCEvent.h:30
TClonesArray * GetClustersArray()
R__EXTERN HarpoDetSet * gHDetSet
Definition: HarpoDetSet.h:71