StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StFgtLenTreeMaker.cxx
1 #include "StFgtLenTreeMaker.h"
2 #include "StRoot/StEvent/StFgtCollection.h"
3 #include "StRoot/StEvent/StFgtHitCollection.h"
4 #include "StRoot/StEvent/StFgtHit.h"
5 #include "StRoot/StFgtUtil/geometry/StFgtGeom.h"
6 #include "StRoot/StEvent/StEvent.h"
7 #include "StRoot/StEvent/StEventInfo.h"
8 #include <TH2D.h>
9 #include <TROOT.h>
10 #include <TStyle.h>
11 #include <TCanvas.h>
12 
13 
14 Int_t StFgtLenTreeMaker::Make()
15 {
16  StEvent* eventPtr = 0;
17  eventPtr = (StEvent*)GetInputDS("StEvent");
18  //(*outTxtFile) <<endl<<endl<<" ------new event: " << eventPtr->info()->id() << "----------------------" << " running nr: " << runningEvtNr << "------" << endl;
19 
20  Bool_t anycluster=false;
21  Int_t iTrk=0;
22 
23  for(int iD=0;iD<kFgtNumDiscs;iD++)
24  {
25  Ncl[iD]=0;
26  for(int iCl=0;iCl<20;iCl++)
27  {
28  cl_geoId[iD][iCl]=-1;
29  cl_seedType[iD][iCl]=-1;
30  cl_quad[iD][iCl]=-1;
31  cl_z[iD][iCl]=0.;cl_ez[iD][iCl]=0.;
32  cl_phi[iD][iCl]=0;cl_ephi[iD][iCl]=0;
33  cl_r[iD][iCl]=-1;cl_er[iD][iCl]=-1;
34  cl_charge[iD][iCl]=-1;cl_echarge[iD][iCl]=-1;
35  cl_numStrips[iD][iCl]=-1;
36  cl_tQuad[iD][iCl]=-1;
37  cl_tStrip[iD][iCl]=-1;
38  cl_layer[iD][iCl]=' ';
39  cl_key[iD][iCl]=-1;
40  cl_pointers[iD][iCl]=0;
41  maxadc[iD][iCl]=-1;seedadc[iD][iCl]=-1;
42  }
43  };
44 
45  Ntrk=0;
46  for(int iTr=0;iTr<4;iTr++)
47  {
48  tr_phi[iTr]=0.;tr_slope[iTr]=0.;tr_vtx[iTr]=0.;
49  tr_chi2[iTr]=0.;tr_ncluster[iTr]=0;
50  trkArray[iTr].phi=0.;trkArray[iTr].slope=0.;trkArray[iTr].vtx=0.;
51  trkArray[iTr].chi2=0.;trkArray[iTr].ncluster=0;
52  for(int iD=0;iD<6;iD++)
53  {
54  tr_iCl[iTr][iD]=-1;
55  trkArray[iTr].clArray[iD]=0;
56  };
57  }
58  Int_t ierr = StFgtQaMaker::Make();
59  cout <<"in plotter make " << endl;
60  for(int iD=0;iD<kFgtNumDiscs;iD++)
61  {
62  cout <<"trying to get clusters in disc " << iD << endl;
63  StFgtHitCollection* clusterCol=mFgtCollectionPtr->getHitCollection(iD);
64 
65  if(clusterCol)
66  {
67  anycluster=true;
68  Int_t iCl=0;
69  cout <<"got collection, looking at hits .. " <<endl;
70  const StSPtrVecFgtHit &hitVec=clusterCol->getHitVec();
71  StSPtrVecFgtHitConstIterator hitIter;
72  for(hitIter=hitVec.begin();hitIter != hitVec.end();hitIter++)
73  {
74  Int_t iq=(*hitIter)->getQuad();
75  Float_t phi=(*hitIter)->getPositionPhi();
76  Float_t ephi=(*hitIter)->getErrorPhi();
77  Float_t r=(*hitIter)->getPositionR();
78  Float_t er=(*hitIter)->getErrorR();
79  Float_t z=(*hitIter)->getPositionZ();
80  Float_t ez=(*hitIter)->getErrorZ();
81  Int_t geoId=(*hitIter)->getCentralStripGeoId();
82  Float_t charge=(*hitIter)->charge();
83  Float_t echarge=(*hitIter)->getChargeUncert();
84  Int_t numStrips=(*hitIter)->getStripWeightMap().size();
86  Int_t seedType=-99;
87  Bool_t containsSeed=false;
88  Int_t maxA=0;
89  Int_t seedA=0;
90  for(stripWeightMap_t::iterator it=(*hitIter)->getStripWeightMap().begin();it!=(*hitIter)->getStripWeightMap().end();it++)
91  {
92  if(it->first->getClusterSeedType()==kFgtSeedType1 || it->first->getClusterSeedType()==kFgtSeedType2 ||it->first->getClusterSeedType()==kFgtSeedType3) //require 1,2,3 strips
93  {
94  containsSeed=true;
95  if(seedA < it->first->getMaxAdc())seedA=it->first->getMaxAdc();
96  if(it->first->getClusterSeedType()==kFgtSeedType1)seedType=1;
97  if(it->first->getClusterSeedType()==kFgtSeedType2 && seedType!=1)seedType=2;
98  if(it->first->getClusterSeedType()==kFgtSeedType3 && seedType!=1 && seedType!=2)seedType=3;
99  }
100  if(maxA < it->first->getMaxAdc())maxA=it->first->getMaxAdc();
101  }
102  if(iCl<20){
103  maxadc[iD][iCl]=maxA;
104  seedadc[iD][iCl]=seedA;
105  };
106  if(containsSeed)
107  {
108  cout <<"cluster contains a seed " <<endl;
109  if((*hitIter)->getLayer()!='R')
110  {
111  continue;
112  }
113 
114  Short_t tDisc, tQuad,tStrip;
115  Char_t tLayer;
116  StFgtGeom::decodeGeoId(geoId,tDisc,tQuad,tLayer,tStrip);
117  if(iCl<20){
118  cout <<"filling cluster information " <<endl;
119  cl_geoId[iD][iCl]=geoId;
120  cl_seedType[iD][iCl]=seedType;
121  cl_quad[iD][iCl]=iq;
122  cl_z[iD][iCl]=z;cl_ez[iD][iCl]=ez;
123  cl_phi[iD][iCl]=phi;cl_ephi[iD][iCl]=ephi;
124  cl_r[iD][iCl]=r;cl_er[iD][iCl]=er;
125  cl_charge[iD][iCl]=charge;cl_echarge[iD][iCl]=echarge;
126  cl_numStrips[iD][iCl]=numStrips;
127  cl_tQuad[iD][iCl]=tQuad;
128  cl_tStrip[iD][iCl]=tStrip;
129  cl_layer[iD][iCl]=(*hitIter)->getLayer();
130  cl_key[iD][iCl]=(*hitIter)->getKey();
131  cl_pointers[iD][iCl]=(*hitIter);
132  }
133  iCl++;
134  }
135  }
136  Ncl[iD]=iCl;
137  }
138  }
139 
140 
141  for(Int_t iphi=0;iphi<2;iphi++)
142  {
143  htrk->Reset();
144  //Float_t minphi=;
145  //Float_t maxphi=;
146  Int_t trNcl=0;
147  Float_t maxA[6]={0.,0.,0.,0.,0.,0.};
148  Int_t kcl[6]={-1,-1,-1,-1,-1,-1};
149 
150  for(Int_t iD=0;iD<6;iD++)
151  {
152  Int_t Ntot=Ncl[iD];
153  if(Ntot>10)Ntot=10;
154  Float_t maxz=-1;
155  Float_t maxr=-1;
156  Float_t emaxr=-1;
157  for(Int_t iC=0;iC<Ntot;iC++)
158  {
159  //if(maxadc[iD][iC]>0. && cl_phi[iD][iC]>minphi && cl_phi[iD][iC]<maxphi)
160  if(maxadc[iD][iC]>maxA[iD] && cl_quad[iD][iC]==iphi)
161  {
162  //printf("%f %f \n",cl_z[iD][iC],cl_r[iD][iC]);
163  maxz=cl_z[iD][iC];
164  maxr=cl_r[iD][iC];emaxr=cl_er[iD][iC];
165  maxA[iD]=maxadc[iD][iC];
166  kcl[iD]=iC;
167  }
168  }
169  if(maxA[iD]>0.)
170  {
171  trNcl++;
172  Int_t zbin=htrk->FindBin(maxz);
173  htrk->SetBinContent(zbin,maxr);
174  Float_t errR=1.;
175  htrk->SetBinError(zbin,errR);
176  };
177  }
178  if(trNcl>2)
179  {
180  cout <<"Track found, fitting .. " <<endl;
181  f0->SetParameter(0,0);
182  f0->SetParameter(1,0);
183  htrk->SetStats(0);
184  htrk->Fit(f0,"Q");
185  htrk->SetMaximum(50);
186  trkArray[iTrk].slope=f0->GetParameter(1);
187  trkArray[iTrk].vtx=-f0->GetParameter(0)/f0->GetParameter(1);
188  trkArray[iTrk].chi2=f0->GetChisquare()/(trNcl-2);
189  trkArray[iTrk].ncluster=trNcl;
190  trkArray[iTrk].phi=iphi;
191  printf("slope=%f vtx=%f chi2=%f \n",trkArray[iTrk].slope,trkArray[iTrk].vtx,trkArray[iTrk].chi2);
192  tr_phi[iTrk]=trkArray[iTrk].phi;
193  tr_slope[iTrk]=trkArray[iTrk].slope;
194  tr_vtx[iTrk]=trkArray[iTrk].vtx;
195  tr_chi2[iTrk]=trkArray[iTrk].chi2;
196  tr_ncluster[iTrk]=trkArray[iTrk].ncluster;
197  for(Int_t iD=0;iD<6;iD++)
198  {
199  printf("%d \n",kcl[iD]);
200  tr_iCl[iTrk][iD]=kcl[iD];
201  if(kcl[iD]>-1)
202  {
203  trkArray[iTrk].clArray[iD]=cl_pointers[iD][kcl[iD]];
204  Int_t qq=trkArray[iTrk].clArray[iD]->getQuad();
205  Float_t pp=trkArray[iTrk].clArray[iD]->getPositionPhi();
206  Float_t rr=trkArray[iTrk].clArray[iD]->getPositionR();
207  printf("********##########*************** iTrk=%d ncluster=%d iD=%d quad=%d phi=%f r=%f",iTrk,trkArray[iTrk].ncluster,iD,qq,pp,rr);
208  }
209  };
210  iTrk++;
211  };
212  };
213  Ntrk=iTrk;
214  iEvt=runningEvtNr;
215  if(anycluster){
216  printf("*************** FILLING THE TREE %d **********************\n", iEvt);
217  tCl->Fill();
218  };
219 
220  FitFunc();
221  for(int iTr=0;iTr<Ntrk;iTr++)
222  {
223  for(int iD=0;iD<6;iD++)
224  {
225  if(trkArray[iTr].clArray[iD])
226  {
227  char hname[100];
228  StFgtHit* clTr=trkArray[iTr].clArray[iD];
229  for(stripWeightMap_t::iterator it=clTr->getStripWeightMap().begin();it!=clTr->getStripWeightMap().end();it++)
230  {
231  rdo=0;arm=apv=chn=stat=-1;
232  disk=quad=strip=-1.;layer=' ';
233  ordinate=lowerSpan=upperSpan=-1.;
234  Int_t dof=Ntimebin-4;
235 
236  it->first->getElecCoords( rdo, arm, apv, chn );
237  stat=0;//how to get status?
238  if(stat>0)continue;//stat=0 is good. Otherwise, skip.
239 
240  int geoId=it->first->getGeoId();
241  if (geoId>0){
242  StFgtGeom::decodeGeoId(geoId,disk,quad,layer,strip);
243  StFgtGeom::getPhysicalCoordinate(geoId,disk,quad,layer,ordinate,lowerSpan,upperSpan);
244  }
245  else{continue;};//strip was not mapped , we readout in the daq file some '0' form non-existing APVs.
246  //printf("%d %d %d %d %d %d \n",iEvt,rdo,arm,apv,chn,geoId);
247 
248  ped=it->first->getPed();
249  pedSig=it->first->getPedErr();
250  if(ped<=0.)continue;
251  for(Int_t is=0;is<7;is++){adc[is]=0.;};
252  for(Int_t is=0;is<Ntimebin;is++){
253  adc[is]=it->first->getAdc(is);
254  //adc[is]-=ped;//already subtracted
255  };
256  Bool_t pass=true;
257  //if(rdo==1 && arm==3 && apv==9)pass=false;
258  //if(rdo==1 && arm==4 && apv==21)pass=false;
259  //if(rdo==1 && arm==0 && apv>-1 && apv<10){}
260  //else if(rdo==2 && arm==0 && apv>-1 && apv<10){}
261  //else{pass=false;};
262 
263  if(pass){
264  chi2=-1.;tau=0.;t0=0.;beta=0.;offset=0.;errCode=0;
265  hh->Reset();
266  for(Int_t is=0;is<Ntimebin;is++){
267  hh->SetBinContent(is+1,adc[is]);
268  hh->SetBinError(is+1,pedSig);
269  }
270  sprintf(hname,"rdo%d_arm%d_apv%d_chn%d_%d",rdo,arm,apv,chn,iEvt);
271  hh->SetTitle(hname);
272  mmax=hh->GetMaximumBin()-1;
273  mmin=hh->GetMinimumBin()-1;
274  adcmax=hh->GetBinContent(mmax+1);
275  if(0){
276  for(Int_t is=0;is<Ntimebin;is++){
277  printf("%d ",it->first->getAdc(is));
278  };
279  printf("ped=%f \n",ped);
280  };
281  if(adcmax>fitThresh){
282  if(abs(mmax-mmin)==2 && mmin>0 && mmax>0 && mmin<Ntimebin-1 && mmax<Ntimebin-1){
283  Float_t middle1=(hh->GetBinContent(mmin)+hh->GetBinContent(mmin+2))/2.;
284  Float_t middle2=(hh->GetBinContent(mmax)+hh->GetBinContent(mmax+2))/2.;
285  if((middle1-hh->GetBinContent(mmin+1))/hh->GetBinError(mmin+1)>3. && (hh->GetBinContent(mmax+1)-middle2)/hh->GetBinError(mmax+1)>3.){errCode=1;}
286  }
287  Int_t highcnt=0;
288  for(Int_t is=0;is<Ntimebin;is++){
289  if(adc[is]>adcmax*0.95 && adcmax>20.*pedSig)highcnt++;
290  };
291  //if(highcnt>2){errCode=2;}
292 
293  if(!errCode){
294  InitFX();
295  hh->Fit(FX,"QI","",0.,Ntimebin);
296  chi2=FX->GetChisquare()/(Float_t)dof;
297  fmax=FX->GetMaximumX();
298  norm=FX->GetParameter(0);
299  tau=FX->GetParameter(1);
300  beta=FX->GetParameter(2);
301  offset=FX->GetParameter(3);
302  t0=FX->GetParameter(4);
303  };
304  tFgt->Fill();//save only fitted events
305  };
306  };
307  }
308  }
309  }
310  }
311  runningEvtNr++;
312  return ierr;
313 };
314 
315  StFgtLenTreeMaker::StFgtLenTreeMaker( const Char_t* name):runningEvtNr(0)
316 {
317  StFgtQaMaker( name, 0,0, "qName" );
318  fname="hFgt";
319  fitThresh=350;
320  Ntimebin=7;
321 };
322 
323 StFgtLenTreeMaker::~StFgtLenTreeMaker()
324 {
325  //delete histogram arrays
326 };
327 
328 
329 Int_t StFgtLenTreeMaker::Finish(){
330  //outTxtFile->close();
331  gStyle->SetPalette(1);
332  cout <<"cluster tree maker finish funciton " <<endl;
333  Int_t ierr = kStOk;
334 
335  tCl->Print();
336  fFgt->cd();
337  ierr=tCl->Write();
338 
339  tFgt->Print();
340  fFgt->cd();
341  tFgt->Write();
342 
343  fFgt->Close();
344  cout << "StFgtLenTreeMaker::Finish()" << endl;
345 
346  return ierr;
347 };
348 
349 Int_t StFgtLenTreeMaker::Init(){
350 
351  Int_t ierr = kStOk;
352  ierr=InitTree();
353 
354  htrk=new TH1F("htrk","Hits in FGT Discs;z (cm);r (cm)",300,-100,200);
355  f0=new TF1("f0","[0]+[1]*x",-500,500);
356  hh=new TH1F("hh","hh",Ntimebin,0,Ntimebin);
357 
358  return ierr;
359 };
360 
361 Int_t StFgtLenTreeMaker::InitTree(){
362 
363  Int_t ierr = kStOk;
364  if(fname.CompareTo("")==0){
365  LOG_ERROR << "No output file name given for the TTree" << endm;
366  ierr = kStErr;
367  }
368  else{
369  TString outName=fname+".tree.root";
370  fFgt = new TFile(outName,"recreate");
371  tCl = new TTree("tCl","TTree for FGT cluster analysis");
372  tCl->Branch("iEvt",&iEvt,"iEvt/I");
373  tCl->Branch("Ncl",Ncl,"Ncl[6]/I");
374  tCl->Branch("cl_geoId",cl_geoId,"cl_geoId[6][20]/I");
375  tCl->Branch("cl_seedType",cl_seedType,"cl_seedType[6][20]/I");
376  tCl->Branch("cl_quad",cl_quad,"cl_quad[6][20]/I");
377  tCl->Branch("cl_z",cl_z,"cl_z[6][20]/F");
378  tCl->Branch("cl_ez",cl_ez,"cl_ez[6][20]/F");
379  tCl->Branch("cl_phi",cl_phi,"cl_phi[6][20]/F");
380  tCl->Branch("cl_ephi",cl_ephi,"cl_ephi[6][20]/F");
381  tCl->Branch("cl_r",cl_r,"cl_r[6][20]/F");
382  tCl->Branch("cl_er",cl_er,"cl_er[6][20]/F");
383  tCl->Branch("cl_charge",cl_charge,"cl_charge[6][20]/F");
384  tCl->Branch("cl_echarge",cl_echarge,"cl_echarge[6][20]/F");
385  tCl->Branch("cl_numStrips",cl_numStrips,"cl_numStrips[6][20]/I");
386  tCl->Branch("cl_tStrip",cl_tStrip,"cl_tStrip[6][20]/I");
387  tCl->Branch("cl_layer",cl_layer,"cl_layer[6][20]/B");
388  tCl->Branch("cl_key",cl_key,"cl_key[6][20]/I");
389  tCl->Branch("maxadc",maxadc,"maxadc[6][20]/I");
390  tCl->Branch("seedadc",seedadc,"seedadc[6][20]/I");
391  tCl->Branch("Ntrk",&Ntrk,"Ntrk/I");
392  tCl->Branch("tr_phi",tr_phi,"tr_phi[4]/F");
393  tCl->Branch("tr_slope",tr_slope,"tr_slope[4]/F");
394  tCl->Branch("tr_vtx",tr_vtx,"tr_vtx[4]/F");
395  tCl->Branch("tr_chi2",tr_chi2,"tr_chi2[4]/F");
396  tCl->Branch("tr_ncluster",tr_ncluster,"tr_ncluster[4]/I");
397  tCl->Branch("tr_iCl",tr_iCl,"tr_iCl[4][6]/I");
398 
399  tFgt = new TTree("tFgt","TTree for FGT time shape analysis");
400 
401  tFgt->Branch("iEvt",&iEvt,"iEvt/I");
402  tFgt->Branch("rdo",&rdo,"rdo/I");
403  tFgt->Branch("arm",&arm,"arm/I");
404  tFgt->Branch("apv",&apv,"apv/I");
405  tFgt->Branch("chn",&chn,"chn/I");
406  tFgt->Branch("disk",&disk,"disk/S");
407  tFgt->Branch("quad",&quad,"quad/S");
408  tFgt->Branch("strip",&strip,"strip/S");
409  tFgt->Branch("stat",&stat,"stat/S");
410  tFgt->Branch("ordinate",&ordinate,"ordinate/D");
411  tFgt->Branch("lowerSpan",&lowerSpan,"lowerSpan/D");
412  tFgt->Branch("upperSpan",&upperSpan,"upperSpan/D");
413  tFgt->Branch("layer",&layer,"layer/B");
414  tFgt->Branch("adc",adc,"adc[7]/I");
415  tFgt->Branch("ped",&ped,"ped/D");
416  tFgt->Branch("pedSig",&pedSig,"pedSig/D");
417  tFgt->Branch("adcmax",&adcmax,"adcmax/I");
418  tFgt->Branch("mmin",&mmin,"mmin/I");
419  tFgt->Branch("mmax",&mmax,"mmax/I");
420  tFgt->Branch("chi2",&chi2,"chi2/F");
421  tFgt->Branch("fmax",&fmax,"fmax/F");
422  tFgt->Branch("norm",&norm,"norm/F");
423  tFgt->Branch("tau",&tau,"tau/F");
424  tFgt->Branch("t0",&t0,"t0/F");
425  tFgt->Branch("beta",&beta,"beta/F");
426  tFgt->Branch("offset",&offset,"offset/F");
427  tFgt->Branch("errCode",&errCode,"errCode/I");
428  };
429  return ierr;
430 };
431 
432 struct StFgtLenTreeMaker::MyFunc {
433  //par[0]=normalization
434  //par[1]=tau
435  //par[2]=exponent: fixed at 2
436  //par[3]=y-offset: fixed at 0
437  //par[4]=t-offset
438  MyFunc(TF1 * f): fFunc(f) {}
439  double Evaluate (double *x, double * par) const {
440  fFunc->SetParameter(0,par[0]);
441  fFunc->SetParameter(1,par[1]);
442  fFunc->SetParameter(2,par[2]);
443  fFunc->SetParameter(3,par[3]);
444  Float_t val=0;
445  if(x[0]-par[4]<0.){val=par[3];}
446  else{val=fFunc->Eval(x[0]-par[4]);};
447  return val;
448  }
449  TF1 * fFunc;
450 };
451 
452 
453 void StFgtLenTreeMaker::FitFunc()
454 {
455  fs = new TF1("fs","[0]*x**[2]*exp(-x/[1])+[3]",-100,100);
456  MyFunc * mf = new MyFunc(fs);
457  FX = new TF1("FX",mf,&MyFunc::Evaluate,-10.,10.,5);
458  FX->SetLineWidth(1);
459  FX->SetLineColor(2);
460  InitFX();
461 };
462 
463 void StFgtLenTreeMaker::InitFX()
464 {
465  FX->SetParameter(0,200.);
466  FX->SetParLimits(0,0.,100000.);
467  FX->SetParameter(1,1.);
468  FX->SetParLimits(1,0.1,10.);
469  FX->FixParameter(2,2.);
470  FX->SetParameter(4,0.);
471  FX->SetParLimits(4,-10.,17.);
472 };
473 
474 
475 ClassImp(StFgtLenTreeMaker);
StFgtQaMaker::Make
virtual Int_t Make()
Definition: StFgtQaMaker.cxx:98
kStErr
Definition: Stypes.h:44
kStOk
Definition: Stypes.h:41
Generated by   doxygen 1.8.5