FairRoot/PandaRoot
PndSciTAnaIdeal.cxx
Go to the documentation of this file.
1 // -------------------------------------------------------------------------
2 // ----- PndSciTAnaIdeal source file -----
3 // ----- Created 10/01/06 by A.Sanchez -----
4 // -------------------------------------------------------------------------
5 
6 #include <cmath>
7 
8 #include "TClonesArray.h"
9 #include "TGeoManager.h"
10 #include "FairRootManager.h"
11 #include "FairGeoNode.h"
12 #include "PndSciTAnaIdeal.h"
13 #include "PndSciTPoint.h"
14 
15 #include "../tracking/riemannfit/PndRiemannHit.h"
16 
17 #include "RhoBase/RhoCandidate.h"
18 #include "PndPidCandidate.h"
19 #include "RhoTools/TEventShape.h"
20 #include "RhoBase/RhoCandList.h"
21 
22 #include "RhoBase/RhoFactory.h"
23 //#include "RhoBase/TRho.h"
24 
25 #include "FairRunAna.h"
26 #include "FairRuntimeDb.h"
27 #include "PndMCTrack.h"
28 //Root includes
29 #include "TROOT.h"
30 #include "TString.h"
31 
32 // ----- Default constructor -------------------------------------------
33 PndSciTAnaIdeal::PndSciTAnaIdeal() :
34  FairTask("Ideal PndSciT Hit Producer"),evt(0),friemann(false),fLHcut(false),fScFtime(1.), fstt(false) {
35  fFileName = "";
36 
37  fBranchName = "SciTPoint";
38 
39 
40  histo = new TH2F("h1","h1",300,-2,2,300,0,10);
41  histpi = new TH2F("h","h",300,0,1.5,300,0,10);
42  histk = new TH2F("hk","hk",200,0,2,200,0,1.5);
43  histp = new TH2F("hp","hp",300,0,1.5,300,0,10);
44  //histo = new TH2F("h1","h1",3000,0,0.02,3000,0,0.02);
45  map1 = new TH2F("h2","h2",3000,0,30,3000,0,1);
46  map2 = new TH2F("h3","h3",100,-4,4,100,-4,4);
47  map3 = new TH2F("h4","h4",100,-4,4,100,-4,4);
48  map4 = new TH2F("h5","h5",100,-4,4,100,-4,4);
49 
50  Dpi1 = new TH1F("hpi1","hpi1",300,0,2.);
51  Dpi2 = new TH1F("hpi2","hpi2",300,0,40.);
52  Dp2 = new TH1F("hp2","hp2",300,0,2);
53  Dpix = new TH1F("hpix","hpix",300,0,2.);
54 }
55 // -------------------------------------------------------------------------
56 
57 
58 
59 // ----- Destructor ----------------------------------------------------
60 PndSciTAnaIdeal::~PndSciTAnaIdeal()
61 { delete(histo);
62 delete fread;
63 delete fGeoH;
64 
65 }
66 
67 void PndSciTAnaIdeal::SetParContainers()
68 {
69  // Get Base Container
70  // FairRunAna* ana = FairRunAna::Instance();
71 // FairRuntimeDb* rtdb=ana->GetRuntimeDb();
72 //
73 
74 }
75 
76 InitStatus PndSciTAnaIdeal::ReInit()
77 {InitStatus stat=kERROR;
78  return stat;
79 }
80 
81 
82 // ----- Public method Init --------------------------------------------
83 InitStatus PndSciTAnaIdeal::Init()
84 {
85  // Get RootManager
86  FairRootManager* ioman = FairRootManager::Instance();
87 
88  if ( ! ioman )
89  {
90  std::cout << "-E- PndSciTAnaIdeal::Init: "
91  << "RootManager not instantiated!" << std::endl;
92  return kFATAL;
93  }
94 
95  // Get input array
96  fPointArray = (TClonesArray*) ioman->GetObject(fBranchName);
97 
98  if ( ! fPointArray )
99  {
100  std::cout << "-W- PndSciTAnaIdeal::Init: "
101  << "No SciTPoint array!" << std::endl;
102  return kERROR;
103  }
104 
105 
106 
107 
108  fmcArray=(TClonesArray*) ioman->GetObject("MCTrack");
109  if ( ! fmcArray )
110  {
111  std::cout << "-W- PndSciTAnaIdeal::Init: "
112  << "No SciTPoint array!" << std::endl;
113  return kERROR;
114  }
115  // Create and register output array
116 
117 
118 
119  std::cout << "-I- PndSciTAnaIdeal: Intialisation successfull" << std::endl;
120  return kSUCCESS;
121 }
122 // -------------------------------------------------------------------------
123 
124 
125 // ----- Public method Exec --------------------------------------------
126 void PndSciTAnaIdeal::Exec(Option_t*)
127 {
128 
129  ++evt;
130  fTrigArray->Delete();
131  fHitOutputArray->Delete();
132  if (fMicroCandidates->GetEntriesFast() != 0) fMicroCandidates->Delete();
133 
134 
135  std::cout << "**** PndSciTAnaIdeal::Exec *****" << evt<<std::endl;
136 
137 
138  TRandom *r = new TRandom; Double_t gam;
139 
140 
141  TVector3 vecs,pos,vecsV,vecsc,posV,momV;
142  int mcpdg = -1;
143  int scpdg = -1;
144  int Idtof,Idsci;
145  Int_t event;
146  event=-1;
147 
148  Double_t mult,timeSciT,timeFib,time,Ssci,Stof,length,bet;
149  Double_t R;
150 
151  Double_t sigt,sigL,sigt1,sigL1,x,y,z,sigx,sigy,sigz;
152 
153  timeSciT=0.;timeFib=0.;x=0.;gam=0.;
154  Stof=0.;Ssci=0.;R =0.;sigx=sigt=sigL=sigt1=sigL1=0.;
155 
156  Int_t iStrip = 0;
157 
158  std::vector<unsigned int> tpchits;
159  std::vector<unsigned int> scifhits;
160  std::vector<unsigned int> tofhits;
161  std::vector<Double_t> reltime;
162 
163 
164  Int_t ntrack = fTrackCandArray->GetEntriesFast();
165  Int_t nhit = fSTPointArray->GetEntriesFast();
166  Int_t nhyp = fHPointArray->GetEntriesFast();
167 
168 
169  unsigned int detId, sfdetId ,fhitId,lhitId,
170  n,det,hit,tpdetid,tphitid;
171  Int_t MotherId;
172 
173  detId = fhitId = lhitId =n =sfdetId =det=hit=-1;
174  Int_t kaon;
175 
176  kaon=0;
177  PndSciTMassTrig* trig = new PndSciTMassTrig();
178  Double_t Phi;
179  Phi=0.;
180 
181  Double_t Lk,Lpi,Lp,Lmu,Le;
182  Double_t tk,tp,tpi,te,tmu;
183 
184 
185 
186 
187 
188  for (Int_t iPoint = 0; iPoint < ntrack; iPoint++)
189  {
190  GFTrackCand* track = (GFTrackCand*) fTrackCandArray->At(iPoint);
191  // std::cout << " Ideal Hit Producer -Point-: " << nPoints << std::endl;
192  if (track==0) continue;
193 
194  PndSciTPoint* Hit;
195 
196 
197  Hit = (PndSciTPoint*)(fPointArray->At(lhitId));
198 
199 
200 
201  if(Hit==0){
202  //std::cout<<" to ze "<<std::endl;
203  continue;
204  }
205 
206 
207 
208 
209  PndMCTrack* mctruth=(PndMCTrack*)fmcArray->At(Hit->GetTrackID());
210  if(mctruth==0)continue;
211  //AddAndExpand(point->GetTrackID(),1,iPoint);
212 
213 
214  MotherId= mctruth->GetMotherID();
215 
216  //if(MotherId!=-1)continue;
217 
218 
219  mcpdg= mctruth->GetPdgCode();
220 
221  double q;
222  double charg;
223  double pfq;
224 
225  if(mcpdg>5000){
226 
227  q = GetChargeIon(mcpdg);
228  //std::cout<<" particle "<<q<<std::endl;
229  }
230  else{ q=TDatabasePDG::Instance()->GetParticle(mcpdg)->Charge()/3.;
231  //std::cout<<" particle "<<q<<std::endl;
232  }
233  fQ=q;
234  TVector3 pos1;
235  sfHit->Position(pos1);
236  TVector3 pos2;
237  Hit->Position(pos2);
238  PndRiemannTrack* RT;
239 
240  double x0;
241  double y0;
242  Double_t S;
243  S=0;
244  Double_t sinth;
245 
246  if(friemann){
247 
248  RT = GetRiemannTrack(track);
249  TVectorT<double> dVec = RT->orig();
250  x0 = dVec[0];
251  y0 = dVec[1];
252 
253 
254  Double_t param[3]={x0,y0,RT->r()};
255 
256  //tpcHit->Position(pos1);
257 
258  Phi = GetDPhi(param,pos1,pos2);
259  sinth = RT->dip(); //sintheta=cosLambda
260  Double_t tanL = sqrt(1.-sinth*sinth)/sinth;
261  //(RT->sign()>0 ? pfq = -1.: pfq = 1.);//tanL
262 
263 
264  S=RT->r()*Phi/sinth;
265 
266  TVector2 centre(x0, y0);
267  Double_t fi0 = centre.Phi(); // 0 - 2Pi
268  Double_t d0 = fQ*(TMath::Sqrt(x0*x0 + y0*y0) - param[2]);
269  TVector2 pivot(-x0, -y0);
270  TVector2 pnt2(pos2.X() - x0,pos2.Y() - y0);
271  Double_t phi2 = pnt2.DeltaPhi(pivot);
272  Double_t pT = fQ*0.3*RT->r()*0.01;
273 
274  fvecp.SetXYZ(-pT*sin(phi2+fi0),pT*cos(phi2+fi0),pT*tanL);
275 
276  }
277 
278  // std::cout << "Riemann Fit r: " << RT.r() << " offset: "
279  // << x<< " " << y << " tandip "
280  // << RT.m()<<" z0 "<<RT.t()<<" phi "<<fabs(S)<<std::endl;
281  // std::cout << " Pt "<<0.3*RT.r()*0.01<<" P "<<(0.3*RT.r()*0.01)/sinth<<std::endl;
282  // std::cout << " sinth "<<sinth<<" sinth v "<<momV.Perp()/momV.Mag()<<std::endl;
283 
284  momV = mctruth->GetMomentum();
285  //Hit->Momentum(momV);
286 
287  if(friemann)x = (0.3*RT->r()*0.01)/sinth; //helix riemanntrack
288  else x = momV.Mag();//y = momV.y();z = momV.z();
289 
290  if(x>2.)continue;
291 
292  R = fabs(momV.Perp()/0.3*q);
293  //if(2*R<0.4335)continue;
294 
295  sigx = r->Gaus(0,0.01* momV.Mag());
296 
297 
298  x += sigx;
299  sigt=sigt1=0.;
300  te = tk = tp = tmu = tpi = 0;
301  Le = Lp = Lmu = Lpi = Lk = 0;
302 
303  event = Hit->GetEventID();
304  //sigt = 0.08;//100ps
305  if(fScFtime==1)sigt = 0.;
306  else if(fScFtime>=0.08) sigt = gRandom->Gaus(0,0.08);
307  timeSciT = Hit->GetTime();
308  timeSciT +=sigt;
309 
310  Stof = Hit->GetLength();
311  //sigL = 0.01*Stof;//1%
312  sigL = r->Gaus(0,0.01*Stof);
313 
314  Stof +=sigL;
315 
316  //kaon selector
317  //sigma =200ps/80ps;->sigkaon=0.150GeV,meanKaon=0.489Gev
318  //Cut on low mass range=0.355 Gev
319 
320  // **** first fiber ring *****
321 
322  Ssci = sfHit->GetLength();
323  //sigL1 = 0.01*Ssci;//1%
324  sigL1 = r->Gaus(0,0.01*Ssci);
325  Ssci +=sigL1;
326  // double
327  timeFib= sfHit->GetTime();
328  // //sigt1 = 0.200;//0.450ns//1ns
329  if(fScFtime==1)sigt1 = 0.;
330  else sigt1 = gRandom->Gaus(0,fScFtime);
331  timeFib +=sigt1;
332 
333  time=fabs(timeSciT-timeFib);
334  //time=fabs(timeSciT);
335 
336  length = fabs(Stof-Ssci);
337 
338  //std::cout<<"length "<<length*momV.Perp()/momV.Mag()<<std::endl;
339 
340  bet = (length*0.01)/(0.3*time);//c=0.3 m/ns;bet=v/c
341  //bet = (fabs(S)*0.01)/(0.3*time);
342  double m;
343  if(friemann)m = x*sqrt(((0.09*time*time)/(fabs(S)*0.01*fabs(S)*0.01))-1.);
344  else m = x*sqrt(((0.09*time*time)/(length*0.01*length*0.01))-1.);
345 
346  //if(friemann)std::cout<<" charge "<<fQ<<" prefitter "
347  //<<double(TMath::Sign(1.,pfq))<<" dip "<<RT->sign()<<std::endl;
348 
349  //if(mcpdg==321)++kaon;
350  //std::cout<<" pdg "<<mcpdg<<std::endl;
351  if(q>0)charg=1.;
352  if(q<0)charg=-1.;
353  Double_t qp = charg * x;
354  TLorentzVector lv;
355 
356  if(friemann)lv.SetXYZT(fvecp.x(),fvecp.y(), fvecp.z(),m);
357  else lv.SetXYZT(momV.x(),momV.y(), momV.z(),m);
358 
359  tk = ((length*0.01)/0.3)*sqrt(((0.490*0.490)/(x*x))+1.);
360  tpi = ((length*0.01)/0.3)*sqrt(((0.139*0.139)/(x*x))+1.);
361  tp = ((length*0.01)/0.3)*sqrt(((0.938*0.938)/(x*x))+1.);
362  tmu = ((length*0.01)/0.3)*sqrt(((0.105*0.105)/(x*x))+1.);
363  te = ((length*0.01)/0.3)*sqrt(((0.0005*0.0005)/(x*x))+1.);
364 
365  Lk = TMath::Gaus(tk,time,fScFtime)/(TMath::Gaus(tk,time,fScFtime)+
366  TMath::Gaus(tpi,time,fScFtime)+
367  TMath::Gaus(tp,time,fScFtime)+
368  TMath::Gaus(tmu,time,fScFtime)+
369  TMath::Gaus(te,time,fScFtime));
370 
371  Lp = TMath::Gaus(tp,time,fScFtime)/(TMath::Gaus(tk,time,fScFtime)+
372  TMath::Gaus(tpi,time,fScFtime)+
373  TMath::Gaus(tp,time,fScFtime)+
374  TMath::Gaus(tmu,time,fScFtime)+
375  TMath::Gaus(te,time,fScFtime));
376 
377  Lpi = TMath::Gaus(tpi,time,fScFtime)/(TMath::Gaus(tk,time,fScFtime)+
378  TMath::Gaus(tpi,time,fScFtime)+
379  TMath::Gaus(tp,time,fScFtime)+
380  TMath::Gaus(tmu,time,fScFtime)+
381  TMath::Gaus(te,time,fScFtime));
382 
383  Lmu = TMath::Gaus(tmu,time,fScFtime)/(TMath::Gaus(tk,time,fScFtime)+
384  TMath::Gaus(tpi,time,fScFtime)+
385  TMath::Gaus(tp,time,fScFtime)+
386  TMath::Gaus(tmu,time,fScFtime)+
387  TMath::Gaus(te,time,fScFtime));
388 
389  Le = TMath::Gaus(te,time,fScFtime)/(TMath::Gaus(tk,time,fScFtime)+
390  TMath::Gaus(tpi,time,fScFtime)+
391  TMath::Gaus(tp,time,fScFtime)+
392  TMath::Gaus(tmu,time,fScFtime)+
393  TMath::Gaus(te,time,fScFtime));
394 
395  //std::cout << " Pt v "<<momV.Perp()<<" P v "<<momV.Mag()<<std::endl;
396  //if (bet<1.&&bet>0.&&time<100&&time>1&&length<2000&&length>0&&q!=0)
397  //if (bet<1.&&bet>0.&&q!=0)
398  //{
399  if (bet==1.){
400  gam=99999;
401  }else gam= 1./sqrt(1.-(bet*bet));
402 
403  if(fLHcut){ if(Lk>0.4)++kaon;
404  }else{
405  if(mcpdg==321)++kaon;
406  }
407 
408  Dpi2->Fill(time);
409  histo->Fill(charg*x,m/x);// 1/(gam*bet));
410  if(q<0)histp->Fill(x,m/x);// 1/(gam*bet));
411  if(q>0)histpi->Fill(x,m/x);// 1/(gam*bet));
412 
413 
414  //if(q>0)
415  histk->Fill(x,m*m);
416  //if(mcpdg==321||mcpdg==211)
417  if(q<0)Dpi1->Fill(m);//x/(gam*bet));
418  if(q>0)Dp2->Fill(m);//(x/(gam*bet)));
419 
420  // Int_t count = fTrigArray->GetEntriesFast();
421  // trig->SetEvt(event);
422  // trig = new ((*fTrigArray)[count])PndSciTMassTrig(Hit->GetTrackID(),iPoint,bet,charg*x,time,length);
423  //}
424 
425  /*
426  Int_t count = fTrigArray->GetEntriesFast();
427 
428  trig = new ((*fTrigArray)[count])PndSciTMassTrig(Hit->GetTrackID(),iPoint,
429  bet,qp ,time,length);
430  trig->SetEvt(event);
431  trig->SetP(momV);*/
432 
433  Int_t micsize = fMicroCandidates->GetEntriesFast();
434  PndPidCandidate *micro= new ((*fMicroCandidates)[micsize]) PndPidCandidate((Int_t)fQ,
435  pos1,lv);
436  micro->SetFirstHit(pos1);
437  micro->SetLastHit(pos2);
438  if(friemann)micro->SetTrackLength(fabs(S));
439  else micro->SetTrackLength(length);
440  micro->SetMcIndex(Hit->GetTrackID());
441  micro->SetTrackIndex(iPoint);
442  micro->SetSciTStopTime(time);
443  micro->SetSciTM2(m*m);
444  micro->SetPionPidLH(Lpi);
445  micro->SetKaonPidLH(Lk);
446  micro->SetProtonPidLH(Lp);
447  micro->SetMuonPidLH(Lmu);
448  micro->SetElectronPidLH(Le);
449 
450 
451 
452  }
453 
454  Int_t count = fTrigArray->GetEntriesFast();
455  new ((*fTrigArray)[count])PndSciTMassTrig(event,kaon);
456 
457 
458  std::cout<<"evt "<<evt<<" k "<<kaon<<std::endl;
459 
460  std::cout<<" entries of PndSciTMassTrig "<<fTrigArray->GetEntriesFast()<<std::endl;
461  std::cout<<" entries of PndSciTMicro "<<fMicroCandidates->GetEntriesFast()<<std::endl;
462 
463  // *** Appl. kaon multiplicity Trigger to Stopped Xi- events
464  // *** cheking the efficiency of tof system.
465  // ** trigger on Silicon tracker
466 
467  PndHypPoint* stg = (PndHypPoint*)(fHPointArray->At(0));
468  if(stg!=0){
469 
470  if((stg->GetEventID()==event)&&kaon>0){
471  std::cout<<" tof trig -> si "<<event<<" "<<stg->GetEventID()<<std::endl;
472 
473  for(Int_t iHypHit=0;iHypHit<nhyp;++iHypHit)
474  {
475  fCurrentHypPoint =(PndHypPoint*)fHPointArray->At(iHypHit);
476 
477  int size = fHitOutputArray->GetEntriesFast();
478  InitTransMat();
479  fCurrentHypPoint->PositionIn(posV);
480  smearLocal(posV);
481  TVector3 dposLocal(0.01,0.01,0.05);
482 
483  new ((*fHitOutputArray)[size]) PndHypHit(fCurrentHypPoint->GetVolumeID(),
484  (fCurrentHypPoint->GetDetName()).Data(),
485  posV,dposLocal,
486  iHypHit,
487  fCurrentHypPoint->GetEnergyLoss(),1);
488 
489  }
490  }
491  }
492 
493  if(fFileName!="") SetSciTEventCorr(event,kaon);
494 
495  // Event summary
496 
497  fsciFPointArray->Delete();
498  fSTPointArray->Delete();
499  fHPointArray->Delete();
500  fPointArray->Delete();
501  fmcArray->Delete();
502  fTrackCandArray->Delete();
503  fTpcArray->Delete();
504 
505 
506 
507 }
508 // -------------------------------------------------------------------------
509 void PndSciTAnaIdeal::Finish() {
510  if(fFileName!=""){
511 
512  fFile->Write();
513  fFile->Close();
514 
515  delete fEvt;
516  delete fFile;
517  }
518 
519 
520 
521  cout<<" -I PndHyp::FinishRun():closing and deleting fFile fEvt "<<endl;
522  }
523 void PndSciTAnaIdeal::SetSciTEventCorr(Int_t ev, Int_t mult)
524 {
525  TVector3 vecs,pos;
526  Int_t nhit = fSTPointArray->GetEntriesFast();
527  PndHypPoint* stg = (PndHypPoint*)(fSTPointArray->At(0));
528 
529  if(stg!=0){
530  if(ev==stg->GetEventID()&&(mult>0))
531  {
532  for (Int_t ihit = 0; ihit < nhit; ihit++)
533  {
534  PndHypPoint* Hit = (PndHypPoint*)(fSTPointArray->At(ihit));
535  if(Hit==0){
536  //std::cout<<" to ze "<<std::endl;
537  continue;
538  }
539  std::cout<<" it "<<Hit->GetEventID()<<" ev "<<ev<<endl;
540  //if(mult==0) continue;
541 
542  fEvt->Clear();
543  Int_t cnt = 0;
544 
545  fread->GetFragment(ev);
546 
547  TLorentzVector P,V;
548  Int_t ntr = fread->GetNtr();
549  Int_t pid,AI,ZI,LI;
550  Double_t mf;
551  TVector3 pv;
552  Hit->PositionOut(pos);
553  Hit->MomentumIn(vecs);
554 
555  Dpix->Fill(vecs.Mag());
556 
557  for (Int_t i=0;i<ntr;i++){
558 
559  fread->GetData(i, pv,pid,mf);
560  Int_t tt;
561 
562  Double_t Energy=0.;
563  Energy = sqrt(pv.Mag2()+(mf*mf));
564 
565  //cout<<" track "<<i<<" "<<pid<<" "<<mf<<endl;
566  //cout<<" track "<<pv.x()<<" "<<pv.y()<<" "<<pv.z()<<endl;
567 
568  P.SetPxPyPzE(pv.x(),pv.y(),pv.z(),Energy);
569  V.SetX( pos.x()); V.SetY(pos.y()); V.SetZ(pos.z());
570 
571 
572  fread->GetAZH(pid,AI,ZI,LI);
573  //cout<<" track "<<i<<" "<<pid<<" A "<<AI<<" Z "<<ZI<<" L "<<LI<<endl;
574 
575 
576  THParticle fparticle(pid,1,ev,0,mf,AI,ZI,LI,P,V);
577  new((*fEvt)[cnt++]) THParticle(fparticle);
578  }
579 
580  activeCnt = cnt;
581  t->Fill();
582 
583 
584 
585  }
586  }
587  }
588 }
589 
590 void PndSciTAnaIdeal::InitTransMat()
591 {
592  gGeoManager->cd(fGeoH->GetPath( fCurrentHypPoint->GetDetName() ).Data());
593  fCurrentTransMat = gGeoManager->GetCurrentMatrix();
594 
595 
596 }
597 
598 
599 void PndSciTAnaIdeal::smear(TVector3& pos)
600 {
602 
603  double sigx=gRandom->Gaus(0,0.01);
604  double sigy=gRandom->Gaus(0,0.01);
605  double sigz=gRandom->Gaus(0,0.005);
606 
607  double x = pos.x() + sigx;
608  double y = pos.y() + sigy;
609  double z = pos.z() + sigz;
610 
611 
612  pos.SetXYZ(x,y,z);
613 
614  return;
615 }
616 
617 
618 void PndSciTAnaIdeal::smearLocal(TVector3& pos)
619 {
621 
622  Double_t posLab[3], posSens[3];
623 
624  posLab[0]=pos.x(); posLab[1]=pos.y(); posLab[2]=pos.z();
625  fCurrentTransMat->MasterToLocal(posLab,posSens);
626 
627  pos.SetXYZ(posSens[0],posSens[1],posSens[2]);
628 
629  smear(pos); // apply a gaussian
630 
631  posSens[0]=pos.x(); posSens[1]=pos.y(); posSens[2]=pos.z();
632  fCurrentTransMat->LocalToMaster(posSens,posLab);
633  pos.SetXYZ(posLab[0],posLab[1],posLab[2]);
634 
635  return;
636 }
637 PndRiemannTrack* PndSciTAnaIdeal::GetRiemannTrack(GFTrackCand* cand)
638 {
639  PndRiemannTrack* result= new PndRiemannTrack();
640  //result->SetVerbose(fVerbose);
641  for (int i = 0; i < cand->getNHits(); i++){
642  unsigned int detId, hitId;
643  cand->getHit(i, detId, hitId);
644  TVector3 pos(0.,0.,0.);
645  TVector3 dpos(0.,0.,0.);
646 
647  FairMCPoint* myHit;
648 
649  if (detId == 1){
650  myHit = (FairMCPoint*)fsciFPointArray->At(hitId);
651  }
652  else if (detId == 2){
653  myHit = (FairMCPoint*)fTpcArray->At(hitId);
654  }
655  else if (detId == 3){
656  myHit = (FairMCPoint*)fPointArray->At(hitId);
657  }
658 
659  if (myHit != 0){
660 
661  PndRiemannHit hit;
662  hit.setXYZ(myHit->GetX(),myHit->GetY(),myHit->GetZ());
663  hit.setDXYZ(0.01,0.01,0.05);
664  result->addHit(hit);
665  }
666  }
667 
668  result->refit();
669  result->szFit();
670 
671  return result;
672 }
673 
674 Double_t PndSciTAnaIdeal::GetDPhi(Double_t* par,TVector3 &v1,TVector3 &v2)
675 {
676  Double_t DPhi;
677  DPhi=0.;
678  Double_t r,Q,x0,y0;
679  x0=par[0];
680  y0=par[1];
681  r=par[2];
682  (fQ>0 ? Q=1.: Q=-1.);
683 
684  TVector2 centre(x0, y0);
685  Double_t fi0 = centre.Phi(); // 0 - 2Pi
686  Double_t d0 = Q*(TMath::Sqrt(x0*x0 + y0*y0) - r);
687  TVector3 pivot(-x0, -y0,0.);
688 
689 
690 
691  TVector3 pnt1(v1.X() - x0,v1.Y() - y0,0.);
692  Double_t phi1;// = pnt1.DeltaPhi(pivot);
693  TVector3 pnt2(v2.X() - x0, v2.Y() - y0,0.);
694  Double_t phi2 ;//= pnt2.DeltaPhi(pivot);
695  Double_t cosPhi1 = (pnt1*pivot)/(pnt1.Mag()*pivot.Mag());
696  Double_t cosPhi2 = (pnt2*pivot)/(pnt2.Mag()*pivot.Mag());
697 
698  //std::cout<<" phi1 "<<cosPhi1<<" "<<cosPhi2<<std::endl;
699 
700  if(cosPhi1 >=0.) phi1=acos(cosPhi1);
701  else phi1 = TMath::Pi() + acos(cosPhi1);
702  if(cosPhi2 >=0.) phi2=acos(cosPhi2);
703  else phi2 = TMath::Pi()+acos(cosPhi2);
704  //pnt2.DeltaPhi(pnt1);//phi2-phi1;
705  if(phi2-phi1>0) DPhi= (phi2-phi1);
706  if(phi2-phi1<0) DPhi = (phi2-phi1);
707 
708  //std::cout<<" phi1 "<<acos(cosPhi1)<<" "<<acos(cosPhi2)<<std::endl;
709 
710  return DPhi;
711 
712 }
713 
714 // ___________________________________________________________________________
715 
716 
717 
718 Int_t PndSciTAnaIdeal::GetChargeIon(Int_t ion)
719 {
720  Int_t A,Z,L;
721 
722  if(ion>1000000000&&(ion<1010000000))
723  { ion -= 1000000000;
724  Z = ion/10000;
725  ion -= 10000*Z;
726  A = ion/10;
727  //cout<<" ion charge "<<Z<<endl;
728 
729  return Z;
730 
731  }
732 if((ion>1010000000||ion>1020000000))
733  {
734  ion -= 1000000000;
735  L = ion/10000000;
736  ion -= 10000000*L;
737  Z = ion/10000;
738  ion -= 10000*Z;
739  A = ion/10;
740  //cout<<L<<" hypernuclei charge "<<Z<<endl;
741 
742  return Z;
743 
744  }
745 }
746 void PndSciTAnaIdeal::WriteHistograms(){
747  TFile* file = FairRootManager::Instance()->GetOutFile();
748  file->cd();
749  file->mkdir("tof_trig");
750  file->cd("tof_trig");
751 
752  histo->Write();
753  delete histo;
754  histo=NULL;
755 
756  histpi->Write();
757  delete histpi;
758  histpi=NULL;
759 
760  histk->Write();
761  delete histk;
762  histk=NULL;
763 
764  histp->Write();
765  delete histp;
766  histp=NULL;
767 
768  map1->Write();
769  delete map1;
770  map1=NULL;
771  map2->Write();
772  delete map2;
773  map2=NULL;
774  map3->Write();
775  delete map3;
776  map3=NULL;
777  map4->Write();
778  delete map4;
779  map4=NULL;
780 
781  Dpi1->Write();
782  delete Dpi1;
783  Dpi1 = NULL;
784 
785  Dpi2->Write();
786  delete Dpi2;
787  Dpi2 = NULL;
788 
789  Dp2->Write();
790  delete Dp2;
791  Dp2 = NULL;
792 
793  Dpix->Write();
794  delete Dpix;
795  Dpix = NULL;
796  }
797 
798 
799 ClassImp(PndSciTAnaIdeal)
pos
TVector3 pos
Definition: GammaSpectraAnalysis.C:43
acos
friend F32vec4 acos(const F32vec4 &a)
Definition: P4_F32vec4.h:113
PndPidCandidate::SetTrackIndex
void SetTrackIndex(int idx)
Definition: PndPidCandidate.h:326
x0
Double_t x0
Definition: checkhelixhit.C:70
nhit
int nhit
Definition: anal_hit_digi_cluster_fwendcap.C:48
PndHypPoint::MomentumIn
void MomentumIn(TVector3 &mom)
Definition: PndHypPoint.h:126
PndSciTAnaIdeal::fCurrentHypPoint
PndHypPoint * fCurrentHypPoint
Definition: PndSciTAnaIdeal.h:86
cos
friend F32vec4 cos(const F32vec4 &a)
Definition: P4_F32vec4.h:112
PndSciTAnaIdeal::fTrackCandArray
TClonesArray * fTrackCandArray
Definition: PndSciTAnaIdeal.h:83
PndSciTAnaIdeal::fGeoH
PndHypGeoHandling * fGeoH
Definition: PndSciTAnaIdeal.h:68
PndSciTAnaIdeal::smearLocal
void smearLocal(TVector3 &pos)
Definition: PndSciTAnaIdeal.cxx:618
PndSciTAnaIdeal::Exec
virtual void Exec(Option_t *opt)
Definition: PndSciTAnaIdeal.cxx:126
PndSciTAnaIdeal::fread
HypStatDecay * fread
Definition: PndSciTAnaIdeal.h:119
r
double r
Definition: RiemannTest.C:14
i
Int_t i
Definition: run_full.C:25
m
__m128 m
Definition: P4_F32vec4.h:28
file
TFile * file
Definition: drc/old/draw_geom.C:12
PndSciTAnaIdeal::fEvt
TClonesArray * fEvt
Definition: PndSciTAnaIdeal.h:115
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
track
PndRiemannTrack track
Definition: RiemannTest.C:33
HypStatDecay::GetFragment
void GetFragment(Int_t primGen)
Definition: HypStatDecay.cxx:102
CAMath::Sqrt
static T Sqrt(const T &x)
Definition: PndCAMath.h:37
PndMCTrack::GetPdgCode
Int_t GetPdgCode() const
Definition: PndMCTrack.h:73
sin
friend F32vec4 sin(const F32vec4 &a)
Definition: P4_F32vec4.h:111
ev
int ev
Definition: GammaSpectraAnalysis.C:44
PndSciTAnaIdeal::SetParContainers
virtual void SetParContainers()
Definition: PndSciTAnaIdeal.cxx:67
GFTrackCand::getNHits
unsigned int getNHits() const
Definition: GFTrackCand.h:113
n
int n
Definition: HypStatDecay.cxx:21
pid
int pid()
Definition: detectors/disc/pid.C:1
par
Double_t par[3]
Definition: anaLmdCluster.C:146
PndPidCandidate::SetMcIndex
void SetMcIndex(int idx)
Definition: PndPidCandidate.h:323
PndMCTrack::GetMomentum
TVector3 GetMomentum() const
Definition: PndMCTrack.h:78
PndSciTAnaIdeal::InitTransMat
void InitTransMat()
Input file name.
Definition: PndSciTAnaIdeal.cxx:590
evt
int evt
Definition: checkhelixhit.C:36
PndSciTAnaIdeal::fPointArray
TClonesArray * fPointArray
Definition: PndSciTAnaIdeal.h:78
PndHypPoint::GetVolumeID
Int_t GetVolumeID() const
Definition: PndHypPoint.h:91
PndSciTAnaIdeal::fTpcArray
TClonesArray * fTpcArray
Definition: PndSciTAnaIdeal.h:87
gGeoManager
TGeoManager * gGeoManager
Definition: createRootGeoFile.C:97
PndPidCandidate::SetFirstHit
void SetFirstHit(TVector3 &pos)
Definition: PndPidCandidate.h:310
PndRiemannHit::setDXYZ
void setDXYZ(double dx, double dy, double dz)
Definition: PndRiemannHit.cxx:70
PndSciTAnaIdeal::fsciFPointArray
TClonesArray * fsciFPointArray
Definition: PndSciTAnaIdeal.h:82
d0
Double_t d0
Definition: checkhelixhit.C:59
vecs
TVector3 vecs
Definition: GammaSpectraAnalysis.C:43
PndSciTAnaIdeal::GetRiemannTrack
PndRiemannTrack * GetRiemannTrack(GFTrackCand *cand)
Definition: PndSciTAnaIdeal.cxx:637
PndSciTAnaIdeal::fmcArray
TClonesArray * fmcArray
Definition: PndSciTAnaIdeal.h:98
PndHypPoint::PositionIn
void PositionIn(TVector3 &pos)
Definition: PndHypPoint.h:125
PndHypGeoHandling::GetPath
TString GetPath(TString id)
for a given ID the path is returned
Definition: PndHypGeoHandling.cxx:76
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
GFTrackCand::getHit
void getHit(unsigned int i, unsigned int &detId, unsigned int &hitId) const
Get detector ID and cluster index (hitId) for hit number i.
Definition: GFTrackCand.h:84
PndSciTAnaIdeal::GetDPhi
Double_t GetDPhi(Double_t *par, TVector3 &v1, TVector3 &v2)
Definition: PndSciTAnaIdeal.cxx:674
y0
Double_t y0
Definition: checkhelixhit.C:71
PndRiemannHit::setXYZ
void setXYZ(double x, double y, double z)
Definition: PndRiemannHit.cxx:57
PndSciTAnaIdeal::fCurrentTransMat
TGeoHMatrix * fCurrentTransMat
Definition: PndSciTAnaIdeal.h:85
PndSciTAnaIdeal::fTrigArray
TClonesArray * fTrigArray
Definition: PndSciTAnaIdeal.h:81
PndRiemannTrack::refit
void refit(bool withErrorCalc=true)
Definition: PndRiemannTrack.cxx:172
GFTrackCand
Track candidate – a list of cluster indices.
Definition: GFTrackCand.h:55
PndSciTAnaIdeal::fSTPointArray
TClonesArray * fSTPointArray
Definition: PndSciTAnaIdeal.h:80
PndSciTAnaIdeal::fHPointArray
TClonesArray * fHPointArray
Definition: PndSciTAnaIdeal.h:79
z
Double_t z
Definition: createRootFscGeometryFile.C:146
PndSciTAnaIdeal::fFileName
const Char_t * fFileName
Gives Access to the Statistical decay products.
Definition: PndSciTAnaIdeal.h:120
PndSciTAnaIdeal::Finish
virtual void Finish()
Definition: PndSciTAnaIdeal.cxx:509
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
gam
CandList gam
Definition: softtrigger_kin5.C:144
mcpdg
int mcpdg
Definition: GammaSpectraAnalysis.C:44
PndRiemannTrack::r
double r() const
Definition: PndRiemannTrack.cxx:416
HypStatDecay::GetNtr
Int_t GetNtr() const
Definition: HypStatDecay.h:57
PndHypPoint::PositionOut
void PositionOut(TVector3 &pos)
Definition: PndHypPoint.h:128
hit
int hit(Int_t nEvents=0, TString inFile="sim.root", TString parFile="par.root", TString inDigi="digi.root", TString outFile="hit.root", Int_t timeBased=0)
Definition: hit.C:1
PndSciTAnaIdeal::GetChargeIon
Int_t GetChargeIon(Int_t ion)
Definition: PndSciTAnaIdeal.cxx:718
P
GeV c P
Definition: ReadMesDecCorrRes_PiT.C:42
PndSciTAnaIdeal::fMicroCandidates
TClonesArray * fMicroCandidates
Definition: PndSciTAnaIdeal.h:88
x
Double_t x
Definition: createRootFscGeometryFile.C:159
HypStatDecay::GetAZH
void GetAZH(Int_t ion, Int_t &AI, Int_t &ZI, Int_t &L)
Definition: HypStatDecay.cxx:224
PndHypPoint::GetEventID
Int_t GetEventID() const
Definition: PndHypPoint.h:90
HypStatDecay::GetData
void GetData(Int_t tr, TVector3 &p, Int_t &pid, Double_t &mass)
search fragment in event
Definition: HypStatDecay.cxx:207
PndSciTAnaIdeal::fHitOutputArray
TClonesArray * fHitOutputArray
Definition: PndSciTAnaIdeal.h:84
PndSciTAnaIdeal::smear
void smear(TVector3 &pos, TVector3 &dpos)
v1
TVector3 v1
Definition: bump_analys.C:40
cnt
Int_t cnt
Definition: hist-t7.C:106
PndRiemannTrack::orig
TVectorD orig() const
Definition: PndRiemannTrack.cxx:405
PndHypPoint::GetDetName
TString GetDetName() const
Definition: PndHypPoint.h:119
ntrack
Int_t ntrack
Definition: reco_apd.C:21
ClassImp
ClassImp(PndAnaContFact)
Z
double Z
Definition: anaLmdDigi.C:68
v2
TVector3 v2
Definition: bump_analys.C:40
PndSciTAnaIdeal::ReInit
virtual InitStatus ReInit()
Definition: PndSciTAnaIdeal.cxx:76
PndRiemannTrack::szFit
void szFit(bool withErrorCalc=true)
Definition: PndRiemannTrack.cxx:441
PndRiemannTrack::addHit
void addHit(PndRiemannHit &hit)
Definition: PndRiemannTrack.cxx:102
PndPidCandidate::SetLastHit
void SetLastHit(TVector3 &pos)
Definition: PndPidCandidate.h:315
count
int count
Definition: GammaSpectraAnalysis.C:48
PndSciTAnaIdeal::SetSciTEventCorr
void SetSciTEventCorr(Int_t ev, Int_t mult)
Definition: PndSciTAnaIdeal.cxx:523
y
Double_t y
Definition: createRootFscGeometryFile.C:159
PndSciTAnaIdeal::Init
virtual InitStatus Init()
Definition: PndSciTAnaIdeal.cxx:83
PndMCTrack::GetMotherID
Int_t GetMotherID() const
Definition: PndMCTrack.h:74
Pi
Double_t Pi
Definition: PndGeoHypGeCrystal.cxx:10
MotherId
Int_t MotherId
Definition: GammaSpectraAnalysis.C:39
R
Double_t R
Definition: checkhelixhit.C:61
mult
Double_t mult
Definition: GammaSpectraAnalysis.C:45
Generated on Sun Aug 18 2019 03:01:43 for FairRoot/PandaRoot by   doxygen 1.8.8