FEDRA emulsion software from the OPERA Collaboration
IO.C File Reference

Functions

void filltree (TTree *tree, EdbPVRec *al, int fillraw=0)
 
void getDataCouples (const char *file_name, const char *tree_name="couples", EdbPattern *pat, Bool_t MaxCHI2=kFALSE)
 
void getDataEdb (EdbRun *edbRun, EdbPVRec *pvol, TIndexCell *up, TIndexCell *down, float shrU, float shrD, float uoff[2], float doff[2], float plate[3])
 
void getDataFile (const char *file, EdbPatternsVolume *pvol, float zpat)
 
void getDataTree (const char *file, int npl, int entries[], float dz1, float dz2, float shr, EdbPatternsVolume *pvol)
 
void getDataTreeJ (const char *file, int npl, int entries[], float dz, EdbPatternsVolume *pvol)
 
void getPatternEdb (EdbRun *edbRun, EdbPattern *pat, TIndexCell *elist, float shr, float off[2], float dz)
 
TTree * inittree (char *tree_name="couples", char *file_name="aedb.root", char *mode="RECREATE")
 
void make_views_coord_map (EdbRun *edbRun, TIndexCell &up, TIndexCell &down)
 
void make_views_map (EdbRun *edbRun, TIndexCell &up, TIndexCell &down)
 

Function Documentation

◆ filltree()

void filltree ( TTree *  tree,
EdbPVRec al,
int  fillraw = 0 
)
456 {
457  EdbPatCouple *patc=0;
458 
459  int pid1,pid2;
460  EdbSegCouple *cp=0;
461  EdbSegP *s1=0;
462  EdbSegP *s2=0;
463  EdbSegP *s=0;
464 
465  tree->SetBranchAddress("pid1",&pid1);
466  tree->SetBranchAddress("pid2",&pid2);
467  tree->SetBranchAddress("cp" ,&cp);
468  tree->SetBranchAddress("s1." ,&s1);
469  tree->SetBranchAddress("s2." ,&s2);
470  tree->SetBranchAddress("s." ,&s );
471 
472  if(fillraw) {
473  // **** fill tree with raw segments ****
474  EdbPattern *pat=0;
475  int nic;
476  int nip=al->Npatterns();
477  for( int ip=0; ip<nip; ip++ ) {
478  pat = al->GetPattern(ip);
479  pid1 = pat->ID();
480  pid2 = -1;
481  nic=pat->N();
482  for( int ic=0; ic<nic; ic++ ) {
483  s1 = pat->GetSegment(ic);
484  tree->Fill();
485  }
486  }
487  }
488 
489  // **** fill tree with found couples ****
490  s = new EdbSegP();
491 
492  int nip=al->Ncouples();
493  for( int ip=0; ip<nip; ip++ ) {
494  patc = al->GetCouple(ip);
495  pid1 = patc->Pat1()->ID();
496  pid2 = patc->Pat2()->ID();
497 
498  int nic=patc->Ncouples();
499  for( int ic=0; ic<nic; ic++ ) {
500  // printf("%d %d\n",ip,ic);
501  cp = patc->GetSegCouple(ic);
502  s1 = patc->Pat1()->GetSegment(cp->ID1());
503  s2 = patc->Pat2()->GetSegment(cp->ID2());
504  //s = new EdbSegP(*s1);
505  //*s += *s2;
506  s->Set( ic, s1->X(), s1->Y(),
507  (s1->X()-s2->X())/(s1->Z()-s2->Z()),
508  (s1->Y()-s2->Y())/(s1->Z()-s2->Z()),
509  s1->W()+s2->W());
510  s->SetZ( s1->Z() );
511 
512 // s->Set( ic, (s1->X()+s2->X())/2., (s1->Y()+s2->Y())/2.,
513 // (s1->TX()+s2->TX())/2., (s1->TY()+s2->TY())/2.,
514 // s1->W()+s2->W());
515 
516 // s->SetZ( (s1->Z()+s2->Z())/2. );
517 
518  tree->Fill();
519  //delete s;
520  }
521  }
522 
523  tree->Write();
524  gDirectory->Purge();
525 }
EdbPatCouple * GetCouple(int i) const
Definition: EdbPVRec.h:190
Int_t Ncouples() const
Definition: EdbPVRec.h:188
Definition: EdbPVRec.h:27
EdbPattern * Pat2()
Definition: EdbPVRec.h:90
EdbPattern * Pat1()
Definition: EdbPVRec.h:89
EdbSegCouple * GetSegCouple(int i) const
Definition: EdbPVRec.h:85
int Ncouples() const
Definition: EdbPVRec.h:80
Definition: EdbPattern.h:280
int ID() const
Definition: EdbPattern.h:328
Int_t Npatterns() const
Definition: EdbPattern.h:380
EdbPattern * GetPattern(int id) const
Definition: EdbPattern.cxx:1887
Definition: EdbSegCouple.h:14
int ID2() const
Definition: EdbSegCouple.h:53
int ID1() const
Definition: EdbSegCouple.h:52
Definition: EdbSegP.h:18
Float_t X() const
Definition: EdbSegP.h:170
Float_t Z() const
Definition: EdbSegP.h:150
void SetZ(float z)
Definition: EdbSegP.h:122
Float_t Y() const
Definition: EdbSegP.h:171
Float_t W() const
Definition: EdbSegP.h:148
void Set(int id, float x, float y, float tx, float ty, float w, int flag)
Definition: EdbSegP.h:86
Int_t N() const
Definition: EdbPattern.h:89
EdbSegP * GetSegment(int i) const
Definition: EdbPattern.h:66
EdbSegP * s1
Definition: tlg2pattern.C:30
EdbSegP * s2
Definition: tlg2pattern.C:31
EdbSegP * s
Definition: tlg2pattern.C:32
EdbSegCouple * cp
Definition: tlg2pattern.C:29

◆ getDataCouples()

void getDataCouples ( const char *  file_name,
const char *  tree_name = "couples",
EdbPattern pat,
Bool_t  MaxCHI2 = kFALSE 
)
223 {
224  TTree *tree;
225  EdbSegCouple *cp = 0;
226  EdbSegP *s1 = 0;
227  EdbSegP *s2 = 0;
228  EdbSegP *s = 0;
229 
230  printf("get couples from file: %s\n",file_name);
231 
232  TFile *f = new TFile( file_name );
233  if (f) tree = (TTree*)f->Get(tree_name);
234 
235  int nentr = tree->GetEntries();
236 
237  //tree->SetBranchStatus("*" ,0 );
238  //tree->SetBranchStatus("s.*" ,1 );
239 
240  tree->SetBranchAddress("cp" , &cp );
241  //tree->SetBranchAddress("s1." , &s1 );
242  //tree->SetBranchAddress("s2." , &s2 );
243  tree->SetBranchAddress("s." , &s );
244 
245  int nseg = 0;
246 
247  for(int i=0; i<nentr; i++ ) {
248  tree->GetEntry(i);
249  if( !cut_seg_cp( s->X(),s->Y(),s->TX(),s->TY(),s->W()) ) continue;
250  if((cp->N1()>1) && MaxCHI2 ) continue;
251  if((cp->N2()>1) && MaxCHI2 ) continue;
252  s->SetZ(pat->Z());
253  pat->AddSegment( *s );
254  nseg++;
255  }
256  printf("%d (of %d) segments readed from file %s\n", nseg,nentr,file_name );
257  f->Close();
258 }
FILE * f
Definition: RecDispMC.C:150
int cut_seg_cp(float x0, float y0, float tx0, float ty0, float puls)
Definition: align.C:64
int N1() const
Definition: EdbSegCouple.h:54
int N2() const
Definition: EdbSegCouple.h:55
Float_t TX() const
Definition: EdbSegP.h:172
Float_t TY() const
Definition: EdbSegP.h:173
Float_t Z() const
Definition: EdbPattern.h:87
EdbSegP * AddSegment(int i, EdbSegP &s)
Definition: EdbPattern.cxx:71

◆ getDataEdb()

void getDataEdb ( EdbRun edbRun,
EdbPVRec pvol,
TIndexCell up,
TIndexCell down,
float  shrU,
float  shrD,
float  uoff[2],
float  doff[2],
float  plate[3] 
)
134 {
135  float base=plate[1];
136 
137  // float dzu = edbRun->GetHeader()->GetPlate()->GetUp() * shrU/2.;
138  float dzu = plate[2]/2.;
139  EdbPattern *patU = new EdbPattern(0.,0., base+dzu);
140  patU->SetID(1);
141 
142  getPatternEdb(edbRun,patU,up,shrU,uoff,dzu);
143  pvol->AddPattern(patU);
144 
145  // float dzd = edbRun->GetHeader()->GetPlate()->GetDown() * shrD/2.;
146  float dzd = plate[0]/2.;
147  EdbPattern *patD = new EdbPattern(0.,0., -dzd );
148  patD->SetID(2);
149 
150  getPatternEdb(edbRun,patD,down,shrD,doff,-dzd);
151  pvol->AddPattern(patD);
152 }
void getPatternEdb(EdbRun *edbRun, EdbPattern *pat, TIndexCell *elist, float shr, float off[2], float dz)
Definition: IO.C:155
void SetID(int id)
Definition: EdbPattern.h:318
void AddPattern(EdbPattern *pat)
Definition: EdbPattern.cxx:1867
Int_t plate
Definition: merge_Energy_SytematicSources_Electron.C:1

◆ getDataFile()

void getDataFile ( const char *  file,
EdbPatternsVolume pvol,
float  zpat 
)
262 {
263  EdbPattern *pat = new EdbPattern();
264 
265  FILE *ff = fopen( file, "r");
266  int id;
267  float ax,ay,x,y;
268  while( 5 == fscanf(ff,"%d %f %f %f %f", &id,&ax,&ay,&x,&y) ) {
269  pat->AddSegment( id, x, y, ax, ay );
270  }
271  printf("%d segments read from file %s\n", pat->N(),file);
272  fclose(ff);
273 
274  pat->SetZ(zpat);
275  pvol->AddPattern(pat);
276 
277  pvol->SetPatternsID();
278  pvol->FillTree(pat,pat,0);
279 
280 }
void SetPatternsID()
Definition: EdbPattern.cxx:1763
void SetZ(float z)
Definition: EdbPattern.h:41
fclose(pFile)
UInt_t id
Definition: tlg2pattern.C:118
TFile * file
Definition: write_pvr.C:3

◆ getDataTree()

void getDataTree ( const char *  file,
int  npl,
int  entries[],
float  dz1,
float  dz2,
float  shr,
EdbPatternsVolume pvol 
)
352 {
353  EdbPattern *pat = 0;
354  EdbSegP *segP = new EdbSegP();
355 
356  TFile *ff = new TFile(file);
357  TTree *tfxx = (TTree*)ff->Get("tfxx");
358 
359  NSFxxRawDB *rawdb=0;
360  NSSegRawDB *seg=0;
361 
362  tfxx->SetBranchStatus("*",1);
363 
364  tfxx->SetBranchAddress("fxx",&rawdb);
365 
366  int pl_pos=0, surf=0;
367  float zpl=0,ax=0;
368 
369  for(int ipl=0; ipl<npl; ipl++ ) {
370 
371  tfxx->GetEntry( entries[ipl] );
372 
373  if(rawdb->getPOS()==11) zpl=0;
374  if(rawdb->getPOS()==12) zpl=dz2;
375 
376  pat = new EdbPattern();
377  pat->SetZ(zpl);
378 
379  printf( "**** npl=%d, nseg=%d \n",npl, rawdb->N() );
380 
381  int nraw=rawdb->N();
382  for( int i=0; i<nraw; i++ ) {
383  seg = (NSSegRawDB*)(rawdb->getSegments()->At(i));
384 
385  //if( TMath::Abs(seg->getAY()) < gAYcut ) continue;
386  //if( seg->getAX()*10000 == -44. ) continue;
387  //if( seg->getAY()*10000 == 0. ) continue;
388 
389  if( seg->getAX() > .1) continue;
390  if( TMath::Abs(seg->getAY()) > .1 ) continue;
391  //if( seg->getPH()<75000 ) continue;
392 
393  //if(rawdb->getPOS()==12) ax=seg->getAX()-.05;
394  //else ax=seg->getAX();
395 
396 
397  int aoffx=-.01; // empiric angular offset
398  int aoffy=-.02; // empiric angular offset
399 
400  segP->Set( seg->getISG(),
401  seg->getX(),
402  seg->getY(),
403  seg->getAX()*shr-.01,
404  seg->getAY()*shr-.02 );
405  segP->SetZ(zpl);
406  segP->SetW(seg->getPH());
407  pat->AddSegment( *segP);
408  }
409 
410  pvol->AddPattern(pat);
411 
412  pvol->SetPatternsID();
413  pvol->FillTree(pat,pat,0);
414  }
415 
416  ff->Close();
417 }
void SetW(float w)
Definition: EdbSegP.h:129

◆ getDataTreeJ()

void getDataTreeJ ( const char *  file,
int  npl,
int  entries[],
float  dz,
EdbPatternsVolume pvol 
)
286 {
287  EdbPattern *pat = 0;
288  EdbSegP *segP = new EdbSegP();
289 
290  TFile *ff = new TFile(file);
291  TTree *tfxx = (TTree*)ff->Get("tfxx");
292 
293  NSFxxRawDB *rawdb=0;
294  NSSegRawDB *seg=0;
295 
296  tfxx->SetBranchStatus("*",1);
297 
298  tfxx->SetBranchAddress("fxx",&rawdb);
299 
300  float zpl=0,ax=0,shr=1.;
301  int pos;
302 
303  for(int ipl=0; ipl<npl; ipl++ ) {
304 
305  tfxx->GetEntry( entries[ipl] );
306 
307  pos = (int)(rawdb->getPOS()/10);
308 
309  zpl = (36-pos)*dz;
310  pat = new EdbPattern(0,0,zpl);
311 
312  printf( "**** %d (%d) nseg=%d run: %d event: %d plate = %d\n",
313  ipl,npl, rawdb->N(),rawdb->getRUN(),rawdb->getEVENT(),rawdb->getPOS() );
314 
315  int nraw=rawdb->N();
316  for( int i=0; i<nraw; i++ ) {
317  seg = (NSSegRawDB*)(rawdb->getSegments()->At(i));
318 
319  //if( TMath::Abs(seg->getAY()) < gAYcut ) continue;
320  if( seg->getAX()*10000 == -44. ) continue;
321  if( seg->getAY()*10000 == 0. ) continue;
322 
323  //if( seg->getAX()<.01) continue;
324  //if( TMath::Abs(seg->getAY()) > .1 ) continue;
325  //if( seg->getPH()<65000 ) continue;
326 
327  if(rawdb->getPOS()==12) ax=seg->getAX()-.05;
328  else ax=seg->getAX();
329 
330  segP->Set( seg->getISG(),
331  seg->getX(),
332  seg->getY(),
333  ax*shr,
334  seg->getAY()*shr );
335  segP->SetZ(zpl);
336  segP->SetW(seg->getPH());
337  pat->AddSegment( *segP);
338  }
339 
340  pvol->AddPattern(pat);
341 
342  pvol->SetPatternsID();
343  }
344 
345  ff->Close();
346 }
brick dz
Definition: RecDispMC.C:107

◆ getPatternEdb()

void getPatternEdb ( EdbRun edbRun,
EdbPattern pat,
TIndexCell elist,
float  shr,
float  off[2],
float  dz 
)
159 {
160  EdbSegP *segP = new EdbSegP();
161 
162  EdbView *view=edbRun->GetView();
163  EdbSegment *seg=0;
164 
165  float x,y,z,tx,ty;
166 
167  int nseg=0;
168  int nentr = edbRun->GetEntries();
169  //printf("nentr = %d\n",nentr);
170 
171  EdbViewHeader *head=0;
172  int entry;
173 
174  int nsegV;
175  int niu=elist->N();
176  for(int iu=0; iu<niu; iu++) {
177  entry = elist->At(iu)->Value();
178  view = edbRun->GetEntry(entry);
179 
180 
181  head = view->GetHeader();
182 
183  printf(" View %d : \t %d/%d ",
184  entry,head->GetNframesTop(),head->GetNframesBot() );
185  nseg=0;
186 
187  nsegV=view->Nsegments();
188  for(int j=0;j<nsegV;j++) {
189 
190  seg = view->GetSegment(j);
191 
192  if( !cut_seg(seg->GetX0(), seg->GetY0(),
193  seg->GetTx(), seg->GetTy(),
194  seg->GetPuls() ) ) continue;
195 
196  nseg++;
197  //segP->SetDZ( seg->GetDz() );
198 
199  tx = seg->GetTx()/shr + off[0];
200  ty = seg->GetTy()/shr + off[1];
201  x = (seg->GetX0()+view->GetXview()) + dz*tx;
202  y = (seg->GetY0()+view->GetYview()) + dz*ty;
203  segP->Set( seg->GetID(),x,y,tx,ty);
204  segP->SetZ( pat->Z() );
205  segP->SetW(seg->GetPuls());
206  segP->SetVid(entry,j);
207  pat->AddSegment( *segP);
208 
209  }
210 
211  printf("\t nseg = %d ( %d rejected ) \n", nseg, view->Nsegments()-nseg );
212 
213  }
214 
215 
216 }
TLegendEntry * entry
Definition: Canv_SYSTEMATICS_ALLCOMBINED__RMSEnergy__vs__Energy__ELECTRON.C:130
int GetEntries() const
Definition: EdbRun.h:135
EdbView * GetEntry(int entry, int ih=1, int icl=0, int iseg=1, int itr=0, int ifr=0)
Definition: EdbRun.cxx:489
EdbView * GetView() const
Definition: EdbRun.h:109
virtual Float_t GetX0() const
Definition: EdbSegment.h:37
virtual Float_t GetTx() const
Definition: EdbSegment.h:40
virtual Float_t GetY0() const
Definition: EdbSegment.h:38
virtual Float_t GetTy() const
Definition: EdbSegment.h:41
void SetVid(int vid, int sid)
Definition: EdbSegP.h:134
Definition: EdbSegment.h:61
Int_t GetPuls() const
Definition: EdbSegment.h:88
Int_t GetID() const
Definition: EdbSegment.h:90
Definition: EdbView.h:26
Int_t GetNframesTop() const
Definition: EdbView.h:120
Int_t GetNframesBot() const
Definition: EdbView.h:121
Definition: EdbView.h:134
EdbViewHeader * GetHeader() const
Definition: EdbView.h:163
Float_t GetXview() const
Definition: EdbView.h:193
Int_t Nsegments() const
Definition: EdbView.h:216
EdbSegment * GetSegment(int i) const
Definition: EdbView.h:219
Float_t GetYview() const
Definition: EdbView.h:194
Long_t Value() const
Definition: TIndexCell.h:79
TIndexCell const * At(Int_t narg, Int_t vind[]) const
Definition: TIndexCell.cpp:519
Int_t N() const
Definition: TIndexCell.cpp:344

◆ inittree()

TTree* inittree ( char *  tree_name = "couples",
char *  file_name = "aedb.root",
char *  mode = "RECREATE" 
)

423 {
424  TTree *tree;
425 
426  tree = (TTree*)gROOT->FindObject(tree_name);
427  if (!tree) {
428  TFile *f = new TFile(file_name,mode);
429  if (f) tree = (TTree*)f->Get(tree_name);
430  if(!tree) {
431 
432  tree = new TTree(tree_name,tree_name);
433 
434  int pid1,pid2;
435  EdbSegCouple *cp=0;
436  EdbSegP *s1=0;
437  EdbSegP *s2=0;
438  EdbSegP *s=0;
439 
440  tree->Branch("pid1",&pid1,"pid1/I");
441  tree->Branch("pid2",&pid2,"pid2/I");
442  tree->Branch("cp","EdbSegCouple",&cp,32000,99);
443  tree->Branch("s1.","EdbSegP",&s1,32000,99);
444  tree->Branch("s2.","EdbSegP",&s2,32000,99);
445  tree->Branch("s." ,"EdbSegP",&s,32000,99);
446 
447  }
448  }
449 
450  return tree;
451 }

◆ make_views_coord_map()

void make_views_coord_map ( EdbRun edbRun,
TIndexCell up,
TIndexCell down 
)
29 {
30  int nentr = edbRun->GetEntries();
31  printf("Make views coordinate map, nentr = %d\n",nentr);
32 
33 
34  TIndexCell upc;
35  TIndexCell downc;
36  Long_t v[3]; // x,y,entry
37  EdbViewHeader *head=0;
38 
39  Long_t xx=0, yy=0;
40  float cx = 1000., cy = 1000.; // 2x2 mm cells
41  float dx = 400. , dy = 400.; // 400 microns margins
42  float xv,yv;
43  int mx[9] = {0, 0, 0,-1, 1, -1,-1, 1, 1};
44  int my[9] = {0,-1, 1, 0, 0, -1, 1,-1, 1};
45 
46  for(int iv=0; iv<nentr; iv++ ) {
47 
48  head = edbRun->GetEntryHeader(iv);
49  yv = head->GetXview();
50  xv = head->GetYview();
51  xx = xv/cx;
52  yy = yv/cx;
53  v[0] = xx;
54  v[1] = yy;
55  v[2] = iv;
56 
57  if(head->GetNframesBot()==0) { // fill up views
58  upc.Add(3,v);
59  }
60  else if(head->GetNframesTop()==0) { // fill down views
61  downc.Add(3,v); // add center
62  for( int im=1; im<5; im++ ) { // add sides margins
63  v[0] = ( xv+dx*mx[im] ) / cx;
64  v[1] = ( yv+dy*my[im] ) / cy;
65  if( (v[0] != xx) || (v[1] != yy) )
66  downc.Add(3,v);
67  }
68  for( int im=5; im<9; im++ ) { // add angles margins
69  v[0] = ( xv+dx*mx[im] ) / cx;
70  v[1] = ( yv+dy*my[im] ) / cy;
71  if( (v[0] != xx) && (v[1] != yy) )
72  // if(!downc.Find(3,v))
73  downc.Add(3,v);
74  }
75 
76  }
77 
78  }
79 
80  upc.Sort();
81  downc.Sort();
82 
83  upc.SetName("x:y:entry");
84  downc.SetName("x:y:entry");
85 
86  upc.PrintStat();
87  downc.PrintStat();
88 
89  TIndexCell *clx=0;
90  TIndexCell *cly=0;
91 
92  int areac=0;
93  int nix,niy,nie;
94  nix=upc.N(1);
95  for(int ix=0; ix<nix; ix++) {
96  clx = upc.At(ix);
97  xx = clx->Value();
98  niy=clx->N(1);
99  for(int iy=0; iy<niy; iy++) {
100  cly = clx->At(iy);
101  yy = cly->Value();
102  areac++;
103 
104  nie = cly->N(1);
105  for(int ie=0; ie<nie; ie++) {
106  v[0]=areac;
107  v[1] = cly->At(ie)->Value();
108  up->Add(2,v);
109  }
110 
111  cly = downc.Find(xx)->Find(yy);
112  nie=cly->N(1);
113  for(int ie=0; ie<nie; ie++) {
114  v[0] = areac;
115  v[1] = cly->At(ie)->Value();
116  down->Add(2,v);
117  }
118  }
119  }
120 
121  up->Sort();
122  down->Sort();
123 
124  up->PrintStat();
125  down->PrintStat();
126 }
EdbViewHeader * GetEntryHeader(int entry) const
Definition: EdbRun.cxx:548
Float_t GetXview() const
Definition: EdbView.h:93
Float_t GetYview() const
Definition: EdbView.h:94
Definition: TIndexCell.h:19
void PrintStat() const
Definition: TIndexCell.cpp:264
void Sort(Int_t upto=kMaxInt)
Definition: TIndexCell.cpp:539
TIndexCell * Find(Int_t narg, Long_t varg[]) const
Definition: TIndexCell.cpp:613
Int_t Add(Int_t narg, Long_t varg[])
Definition: TIndexCell.cpp:602
void SetName(const char *varlist)
Definition: TIndexCell.cpp:549

◆ make_views_map()

void make_views_map ( EdbRun edbRun,
TIndexCell up,
TIndexCell down 
)
2 {
3  int nentr = edbRun->GetEntries();
4  printf("Make views entry map, nentr = %d\n",nentr);
5 
6  Long_t v[2]; // areaID,entry
7  EdbViewHeader *head=0;
8  for(int iv=0; iv<nentr; iv++ ) {
9 
10  head = edbRun->GetEntryHeader(iv);
11 
12  v[0]=head->GetAreaID();
13  v[1]=iv;
14 
15  if(head->GetNframesTop()==0) { // fill down views
16  down.Add(2,v);
17  }
18  else if(head->GetNframesBot()==0) { // fill up views
19  up.Add(2,v);
20  }
21  }
22 
23  //up->PrintStat();
24  //down->PrintStat();
25 }
Int_t GetAreaID() const
Definition: EdbView.h:91