#include <AlgSubShowerSR.h>
Inheritance diagram for AlgSubShowerSR:
Public Member Functions | |
| AlgSubShowerSR () | |
| virtual | ~AlgSubShowerSR () |
| virtual void | RunAlg (AlgConfig &ac, CandHandle &ch, CandContext &cx) |
| virtual void | Trace (const char *c) const |
Private Member Functions | |
| void | CalculateEnergyVertexAngle (CandSubShowerSRHandle &cch, Int_t det) |
| void | SubShowerID (CandSubShowerSRHandle &cch) |
Private Attributes | |
| Double_t | fPECut |
| Double_t | fECut |
| Double_t | fMaxDiffCut |
| Double_t | fVtxDiffCut |
| Double_t | fOut2RmPHCut |
| Double_t | fTrkFraCut |
| Double_t | fTrkLikeStpCut |
| Double_t | fXtkFraCut |
| Bool_t | fDoEMFit |
| Double_t | fMIPperGeV |
| Double_t | maxStpPHinSS |
| Double_t | maxdiff |
| Double_t | vtxdiff |
| Double_t | trkfra |
| Double_t | avgstpperpln |
| Double_t | xtkfra |
| Double_t | avgph |
| Double_t | out2Rmph |
| Double_t | bothfitprob |
| Double_t | stpdismean |
| Double_t | stpdisrms |
| std::list< Double_t > | LongPl |
| std::list< Double_t > | LongPh |
| std::list< Double_t > | LongFlu |
| std::list< Double_t > | LongPhpe |
| std::list< Double_t > | TranStp |
| std::list< Double_t > | TranPh |
| std::list< Double_t > | TranFlu |
| Double_t | loweststp |
| Double_t | uppeststp |
| Double_t | begz |
| Double_t | endz |
| Double_t | Theta_rad |
| Double_t | totw |
| Double_t | Eguess |
| Int_t | nn0l |
| Int_t | nn0t |
|
|
Definition at line 59 of file AlgSubShowerSR.cxx. 00059 : 00060 fPECut(1.5),fECut(0.1),fMaxDiffCut(1),fVtxDiffCut(2), 00061 fOut2RmPHCut(0.3),fTrkFraCut(0.8),fTrkLikeStpCut(1.3),fXtkFraCut(0.8), 00062 fDoEMFit(1),fMIPperGeV(18.5) 00063 { 00064 }
|
|
|
Definition at line 67 of file AlgSubShowerSR.cxx. References LongFlu, LongPh, LongPhpe, LongPl, TranFlu, TranPh, and TranStp. 00068 {
00069 LongPl.clear();
00070 LongPh.clear();
00071 LongFlu.clear();
00072 LongPhpe.clear();
00073 TranStp.clear();
00074 TranPh.clear();
00075 TranFlu.clear();
00076 }
|
|
||||||||||||
|
Definition at line 123 of file AlgSubShowerSR.cxx. References avgph, avgstpperpln, begz, EMFluctuation::CalcLongFluc(), EMFluctuation::CalcTranFluc(), D_EFF, EC_EFF, Eguess, endz, CandRecoHandle::GetBegPlane(), CandStripHandle::GetCharge(), CandHandle::GetDaughterIterator(), CandRecoHandle::GetEndPlane(), Calibrator::GetMIP(), CandRecoHandle::GetNStrip(), CandStripHandle::GetPlane(), CandSubShowerSRHandle::GetPlaneView(), CandStripHandle::GetStrip(), CandStripHandle::GetTime(), CandStripHandle::GetTPos(), CandStripHandle::GetZPos(), Calibrator::Instance(), LongFlu, LongPh, LongPhpe, LongPl, loweststp, maxdiff, maxStpPHinSS, Munits::ms, MSG, nn0l, nn0t, out2Rmph, RM_EFF, CandSubShowerSRHandle::SetAvgDev(), CandRecoHandle::SetEndZ(), CandSubShowerSRHandle::SetEnergy(), CandSubShowerSRHandle::SetSlope(), CandRecoHandle::SetVtxPlane(), CandRecoHandle::SetVtxT(), CandRecoHandle::SetVtxU(), CandRecoHandle::SetVtxV(), CandRecoHandle::SetVtxZ(), ST_WID, T_EFF, Theta_rad, totw, TranFlu, TranPh, TranStp, trkfra, uppeststp, vtxdiff, X0_EFF, and xtkfra. Referenced by RunAlg(). 00125 {
00126
00127 Int_t nstp = cch.GetNStrip();
00128 Int_t begpl = cch.GetBegPlane();
00129 Int_t endpl = cch.GetEndPlane();
00130 Int_t *plane = new Int_t[nstp];
00131 Int_t *strip = new Int_t[nstp];
00132 Double_t *ph = new Double_t[nstp];
00133 Double_t *ph_pe = new Double_t[nstp];
00134 Double_t *z = new Double_t[nstp];
00135 Double_t *tpos = new Double_t[nstp];
00136 Double_t *time = new Double_t[nstp];
00137
00138 totw = 0.;
00139 Double_t totph = 0.;
00140 Int_t Ns0 = 0;
00141 Int_t icnt = 0;
00142
00143 maxStpPHinSS = 0;
00144 loweststp = 200;
00145 uppeststp = -200;
00146 Double_t lowestpl = begpl;
00147 Double_t uppestpl = endpl;
00148 Calibrator& calibrator=Calibrator::Instance();
00149 CandStripHandleItr stripItr(cch.GetDaughterIterator());
00150 while (CandStripHandle *stp = stripItr()) {
00151 plane[icnt] = stp->GetPlane();
00152 strip[icnt] = stp->GetStrip();
00153 z[icnt] = stp->GetZPos();
00154 tpos[icnt] = stp->GetTPos();
00155 time[icnt] = stp->GetTime();
00156
00157 ph[icnt] = calibrator.GetMIP(stp->GetCharge(CalDigitType::kSigCorr));
00158 totw += ph[icnt];
00159 ph_pe[icnt] = stp->GetCharge(CalDigitType::kPE);
00160 totph += ph_pe[icnt];
00161
00162 if(tpos[icnt]/T_EFF*100.>uppeststp) {
00163 uppeststp = tpos[icnt]/T_EFF*100.;
00164 uppestpl = plane[icnt];
00165 }
00166 if(tpos[icnt]/T_EFF*100.<loweststp) {
00167 loweststp = tpos[icnt]/T_EFF*100.;
00168 lowestpl = plane[icnt];
00169 }
00170 if(ph[icnt]>maxStpPHinSS) maxStpPHinSS = ph[icnt];
00171 if (ph_pe[icnt]<=fPECut) Ns0++;
00172 icnt++;
00173 }
00174
00175 maxdiff = 0.;
00176 vtxdiff = 0.;
00177 trkfra = 0;
00178 avgstpperpln = 0;
00179 xtkfra = 0;
00180 avgph = 0;
00181 out2Rmph = 0;
00182
00183 Int_t maxpl = -1;
00184 Double_t maxz = 0.;
00185 begz = 0.;
00186 endz = 30.;
00187 Int_t maxstp = -1;
00188 Double_t maxtpos = 0.;
00189 Int_t begstp = -1;
00190 Double_t mint = 999999999;
00191 Float_t normpos = 0;
00192
00193 Double_t Theta = 0;
00194 Theta_rad = 0;
00195 Double_t Theta1 = 0;
00196 Double_t Theta_rad1 = 0;
00197 Double_t avgdev = 0.;
00198 Double_t wavgdev = 0.;
00199 Double_t avgdev1 = 0.;
00200 Double_t wavgdev1 = 0.;
00201 Double_t avgdev2 = 0.;
00202 Double_t wavgdev2 = 0.;
00203 Double_t eval[2] = {0};
00204
00205 Eguess = totw*2./fMIPperGeV;
00206 if(Eguess==0) {
00207 MSG("SubShowerSR",Msg::kError) << "Eguess = 0" << endl;
00208 return;
00209 }
00210
00211 Double_t tmax = log(Eguess*1000./EC_EFF)-0.858;
00212 Double_t tmaxpl = tmax*X0_EFF/D_EFF;
00213 TH1F *flspl = new TH1F("flspl","flspl",500,-0.5,499.5);
00214 TH1F *sflspl = new TH1F("sflspl","sflspl",500,-0.5,499.5);
00215
00216 for(Int_t i=0;i<nstp;i++){
00217 flspl->Fill(plane[i],ph[i]);
00218 }
00219 maxpl = Int_t(flspl->GetBinCenter(flspl->GetMaximumBin()));
00220 Int_t getz = 0;
00221 for(Int_t i=0;i<nstp;i++){
00222 if (plane[i]!=maxpl) sflspl->Fill(plane[i],ph[i]);
00223 else if (getz<1&&plane[i]==maxpl) {
00224 maxz = z[i];
00225 getz++;
00226 }
00227 }
00228 double secmaxpl = Int_t(sflspl->GetBinCenter(sflspl->GetMaximumBin()));
00229 maxdiff = TMath::Abs(maxpl-secmaxpl)/tmaxpl;
00230 vtxdiff = TMath::Abs(maxpl-begpl)/tmaxpl;
00231 MSG("SubShowerSR",Msg::kVerbose)
00232 << "vtxdiff " <<vtxdiff << " tmaxpl " << tmaxpl <<endl;
00233
00234 delete flspl;
00235 delete sflspl;
00236 flspl = NULL;
00237 sflspl = NULL;
00238
00239 Int_t Ns = 0;
00240 Int_t Ns1 = 0;
00241 Double_t stpden = 0.;
00242 Double_t stpsp = 0.;
00243 double phplbeg = 0.;
00244 double phplend = 0.;
00245 Double_t avgbeg = 0.;
00246 Double_t posdif = 0.;
00247 double preavg = 0.;
00248 double prephpl = 0.;
00249 double prenstppl = 0.;
00250 for (int pl = begpl;pl<=endpl;pl++){
00251 int nstppl = 0;
00252 double low = 200.;
00253 double high = -200.;
00254 double phpl = 0.;
00255 double avg = 0.;
00256 for(int i=0;i<nstp;i++){
00257 if(plane[i]==pl) {
00258 nstppl++;
00259 phpl += ph[i];
00260 avg +=tpos[i]*ph[i];
00261 if(tpos[i]/T_EFF*100.>high) high = tpos[i]/T_EFF*100.;
00262 if(tpos[i]/T_EFF*100.<low) low = tpos[i]/T_EFF*100.;
00263 }
00264 }
00265 if (nstppl>0) {
00266 Ns++;
00267 avgstpperpln += Double_t(nstppl);
00268 stpden += double(nstppl)*phpl/(high-low+1.);
00269 avg = avg/T_EFF*100./phpl;
00270 if (nstppl==1) Ns1++;
00271 if(pl>begpl && phplbeg!=0. && preavg!=0.) {
00272 double avgdif = ( (phpl + prephpl)*(avg - preavg) /
00273 (2.*double(nstppl + prenstppl) ) );
00274 stpsp += TMath::Abs(avgdif);
00275 if (avgdif>=0) posdif += 1;
00276 MSG("SubShowerSR",Msg::kVerbose)
00277 << "avg " << avg << " preavg " << preavg << " nstppl " << nstppl
00278 << " prenstppl " << prenstppl << " phpl " << phpl << " prephpl "
00279 << " avgdif " << avgdif << " posdif " << posdif
00280 << " stpsp " << stpsp << endl;
00281 }
00282 if(pl>=begpl&&phplbeg==0.) {
00283 phplbeg = phpl;
00284 avgbeg = avg;
00285 }
00286 if(pl<=endpl) phplend = phpl;
00287 preavg = avg;
00288 prephpl = phpl;
00289 prenstppl = nstppl;
00290 }
00291 }
00292 stpden = stpden/totw;
00293
00294 if (Ns>1) {
00295 stpsp = stpsp/(totw-(phplbeg+phplend)/2.);
00296 posdif = posdif/double(Ns-1);
00297 }
00298 else {
00299 stpsp = 100.;
00300 posdif = 0.5;
00301 }
00302 MSG("SubShowerSR",Msg::kVerbose) << "stpden " << stpden << " stpsp "
00303 << stpsp << " posdif " << posdif << endl;
00304
00305 if(endpl==begpl) {
00306 avgstpperpln = 999;
00307 trkfra = 0;
00308 }
00309 else {
00310 avgstpperpln /= (Double_t(endpl - begpl)/2. + 1.);
00311 trkfra = double(Ns1)/double(Ns);
00312 }
00313 xtkfra = double(Ns0)/double(nstp);
00314 avgph = totph/double(nstp);
00315 Float_t maxstpph = 0.;
00316 Float_t begstpph = 0.;
00317 Int_t nstpmaxpl = 0;
00318
00319 for(Int_t i=0;i<nstp;i++){
00320 if(plane[i]==maxpl&&ph[i]>maxstpph){
00321 maxstp = strip[i];
00322 maxstpph = ph[i];
00323 maxtpos = tpos[i];
00324 }
00325 if(plane[i]==maxpl){
00326 maxz = z[i];
00327 nstpmaxpl++;
00328 }
00329 if(plane[i]==begpl){
00330 if(ph[i]>begstpph){
00331 begz = z[i];
00332 begstp = strip[i];
00333 begstpph = ph[i];
00334 }
00335 if(time[i]<mint) mint = time[i];
00336 }
00337 Bool_t gotendz = false;
00338 if(!gotendz&&plane[i]==endpl){
00339 endz = z[i];
00340 gotendz = true;
00341 }
00342 }
00343
00344 Float_t ss[2][2] = {{0}};
00345 for (int m = 0;m<2;m++){
00346 for (int n = 0;n<2;n++){
00347 ss[m][n] = 0.;
00348 }
00349 }
00350
00351 Double_t cnt = 0;
00352 for(Int_t i=0;i<nstp;i++){
00353 Double_t dtpos = tpos[i] - maxtpos;
00354 Double_t dz = z[i] - maxz;
00355 ss[0][0]+=dz*dz*ph[i];
00356 ss[0][1]+=dtpos*dz*ph[i];
00357 ss[1][1]+=dtpos*dtpos*ph[i];
00358 cnt+=1;
00359 normpos += dtpos*dtpos+dz*dz;
00360 }
00361 ss[1][0]=ss[0][1];
00362
00363 if(totw*normpos!=0){
00364
00365 TMatrixDSym ms(2);
00366 for (int m = 0;m<2;m++){
00367 for (int n = 0;n<2;n++){
00368 ms[m][n] = 0.;
00369 }
00370 }
00371 Double_t *ang = new Double_t[2];
00372
00373 for(Int_t i=0;i<2;i++) ang[i] = -999;
00374 for (Int_t i=0;i<2;i++){
00375 for (Int_t j=0;j<2;j++){
00376 ms(i,j)=ss[i][j]/totw/normpos;
00377 }
00378 }
00379
00380 if(ms.IsA()){
00381 TMatrixDSymEigen eigen(ms);
00382 TMatrixD EVect1 = eigen.GetEigenVectors();
00383 TVectorD EVal1 = eigen.GetEigenValues();
00384 TVector2 *ev0 = new TVector2(EVect1(0,0),EVect1(1,0));
00385 TVector2 *ev1 = new TVector2(EVect1(0,1),EVect1(1,1));
00386 ang[0] = ev0->Phi()*180/TMath::Pi();
00387 ang[1] = ev1->Phi()*180/TMath::Pi();
00388 eval[0] = EVal1(0);
00389 eval[1] = EVal1(1);
00390 delete ev0;
00391 delete ev1;
00392 }
00393 else {
00394 eval[0] = 0.;
00395 eval[1] = 0.;
00396 }
00397 for (Int_t i=0;i<2;i++){
00398 if (ang[i]>90&&ang[i]<270) ang[i]-=180;
00399 if (ang[i]>270&&ang[i]<=360) ang[i]-=360;
00400 if (ang[i]==90||ang[i]==270) ang[i]-=0.01;
00401 }
00402
00403 Theta = ang[0];
00404 Theta1 = ang[1];
00405 Theta_rad = (Theta/180)*TMath::Pi();
00406 Theta_rad1 = (Theta1/180)*TMath::Pi();
00407 double vs = maxtpos+tan(Theta_rad)*(begz-maxz);
00408 double vs1 = maxtpos+tan(Theta_rad1)*(begz-maxz);
00409 MSG("SubShowerSR",Msg::kVerbose) << "vs " << vs << " vs1 " << vs1
00410 << " avgbeg " << avgbeg*T_EFF/100.
00411 << " Theta " << Theta << endl;
00412 if(ss[0][0]==0. || (TMath::Abs(Theta)>TMath::Abs(Theta1) &&
00413 TMath::Abs(avgdev*T_EFF/100.-vs) > 2 *
00414 TMath::Abs(avgbeg*T_EFF/100.-vs1) &&
00415 (stpden<0.5 || (stpsp>0.5 &&
00416 (posdif>0. && posdif<1.)
00417 || Ns<=2)))){
00418 double swap = Theta;
00419 Theta = Theta1;
00420 Theta1 = swap;
00421 }
00422 Theta_rad = (Theta/180)*TMath::Pi();
00423 Theta_rad1 = (Theta1/180)*TMath::Pi();
00424 vs = maxtpos+tan(Theta_rad)*(begz-maxz);
00425 MSG("SubShowerSR",Msg::kVerbose)
00426 << "vs " << vs << " avgbeg " << avgbeg*T_EFF/100.
00427 << " Theta " << Theta << endl;
00428 if(vs>(uppeststp+(uppeststp-loweststp)*1.5/2.)*T_EFF/100.||
00429 vs<(loweststp-(uppeststp-loweststp)*1.5/2.)*T_EFF/100.
00430 ||vs<-4.||vs>4.) {
00431 double bigang = Theta;
00432 if(TMath::Abs(tan(Theta_rad))<TMath::Abs(tan(Theta_rad1))) bigang = Theta1;
00433 Theta = 0.;
00434 Theta1 = TMath::Abs(tan(bigang))/tan(bigang)*89.99;
00435 Theta_rad = (Theta/180)*TMath::Pi();
00436 Theta_rad1 = (Theta1/180)*TMath::Pi();
00437 }
00438 MSG("SubShowerSR",Msg::kVerbose)
00439 << "vs " << vs << " avgbeg " << avgbeg*T_EFF/100.
00440 << " Theta " << Theta << endl;
00441 Double_t *shwstpc = new Double_t[nstp];
00442 Double_t *shwstpc1 = new Double_t[nstp];
00443 Double_t *shwstpc2 = new Double_t[nstp];
00444 Double_t *wshwstpc = new Double_t[nstp];
00445 for(Int_t i=0;i<nstp;i++) shwstpc[i] = -1.;
00446 for(Int_t i=0;i<nstp;i++) shwstpc1[i] = -1.;
00447 for(Int_t i=0;i<nstp;i++) shwstpc2[i] = -1.;
00448 for(Int_t i=0;i<nstp;i++) wshwstpc[i] = -1.;
00449
00450 avgdev = 0.;
00451 wavgdev = 0.;
00452
00453 Float_t offstp = 0;
00454 Float_t offstp2 = 0;
00455
00456 out2Rmph = 0;
00457 for(Int_t i=0;i<nstp;i++){
00458 shwstpc[i]=TMath::Abs(tpos[i]
00459 -(maxtpos+tan(Theta_rad)
00460 *(z[i]-maxz)))*TMath::Cos(Theta_rad);
00461 shwstpc2[i]=TMath::Abs(tpos[i]
00462 -(maxtpos+tan(Theta_rad1)
00463 *(z[i]-maxz)))*TMath::Cos(Theta_rad1);
00464 if (shwstpc[i]>2*RM_EFF/100.) out2Rmph+=ph[i];
00465 if (shwstpc[i]>1) offstp++;
00466 if (shwstpc2[i]>1) offstp2++;
00467 avgdev1 += shwstpc[i];
00468 avgdev2 += shwstpc2[i];
00469 wavgdev1 += shwstpc[i]*ph[i];
00470 wavgdev2 += shwstpc2[i]*ph[i];
00471 }
00472
00473 out2Rmph/=totw;
00474 avgdev1/=nstp;
00475 avgdev2/=nstp;
00476 wavgdev1/=totw;
00477 wavgdev2/=totw;
00478 avgdev = avgdev1;
00479 wavgdev = wavgdev1;
00480 offstp/=nstp;
00481 offstp2/=nstp;
00482 delete [] ang;
00483 delete [] shwstpc;
00484 delete [] shwstpc1;
00485 delete [] shwstpc2;
00486 delete [] wshwstpc;
00487 }
00488 double vs = maxtpos+tan(Theta_rad)*(begz-maxz);
00489
00490 if(vs<-4.||vs>4.) {
00491 vs = avgbeg*T_EFF/100.;
00492 Theta_rad = 0.;
00493 }
00494 if(cch.GetPlaneView()==PlaneView::kX ||
00495 cch.GetPlaneView()==PlaneView::kU) {
00496 cch.SetVtxU(vs);
00497 }
00498 if(cch.GetPlaneView()==PlaneView::kY ||
00499 cch.GetPlaneView()==PlaneView::kV) {
00500 cch.SetVtxV(vs);
00501 }
00502 cch.SetVtxZ(begz);
00503 cch.SetEndZ(endz);
00504 cch.SetVtxT(mint);
00505 cch.SetVtxPlane(begpl);
00506 cch.SetEnergy(Eguess/2.);
00507 cch.SetSlope(tan(Theta_rad));
00508 cch.SetAvgDev(wavgdev);
00509
00510 LongPl.clear();
00511 LongPh.clear();
00512 LongFlu.clear();
00513 LongPhpe.clear();
00514 TranStp.clear();
00515 TranPh.clear();
00516 TranFlu.clear();
00517
00518 EMFluctuation EMFlu(Eguess);
00519 nn0l = 0;
00520 nn0t = 0;
00521
00522 if(cos(Theta_rad)!=0.){
00523 Bool_t xs = false;
00524 Int_t dpl = 2;
00525 if(det==1&&(begpl-249)*(endpl-249)<0) xs = true;
00526 for(Int_t j=begpl;j<=endpl;j+=dpl){
00527 dpl = 2;
00528 if(xs&&(begpl-249)*(j-249)<0) {
00529 dpl = 3;
00530 xs = false;
00531 }
00532 Double_t totphPl = 0.;
00533 Double_t totphPl_pe = 0.;
00534 Double_t plz = -1.;
00535 Double_t nstpPl = 0;
00536 for(Int_t i=0;i<nstp;i++){
00537 if(plane[i]==j && ph[i]>0) {
00538 totphPl+=ph[i];
00539 totphPl_pe+=ph_pe[i];
00540 plz = z[i];
00541 Double_t zt_vtx_corrected[2];
00542 zt_vtx_corrected[0] = z[i]-begz;
00543 zt_vtx_corrected[1] = tpos[i]-vs-(z[i]-begz)*tan(Theta_rad);
00544 nstpPl++;
00545 }
00546 }
00547 LongPl.push_back(plz-begz);
00548 LongPh.push_back(totphPl);
00549 LongPhpe.push_back(totphPl_pe);
00550 if(totphPl>0.) nn0l++;
00551 }
00552
00553 double z_adj = 0.;
00554 std::list<Double_t>::iterator PlIt = LongPl.begin();
00555 std::list<Double_t>::iterator PhIt = LongPh.begin();
00556 std::list<Double_t>::iterator PhpeIt = LongPhpe.begin();
00557 double t0 = (*PhIt);
00558 PhIt++;
00559 double t1 = (*PhIt);
00560 PhIt = LongPh.begin();
00561 if(vtxdiff>1.25) z_adj = Int_t((1-vtxdiff)*tmaxpl);
00562 if(t0>t1/2.&&t0<t1) z_adj = 1;
00563 else if(t0>=t1) {
00564 z_adj = 2;
00565 LongPl.push_front(0);
00566 LongPh.push_front(0.);
00567 LongPhpe.push_front(0.);
00568 }
00569
00570 Int_t inEM = 0;
00571 PlIt = LongPl.begin();
00572 while(PlIt!=LongPl.end()){
00573 if(z_adj!=0) {
00574 (*PlIt) += z_adj*D_EFF/100.;
00575 begz += z_adj*D_EFF/100.;
00576 endz += z_adj*D_EFF/100.;
00577 }
00578 double fluem = 0.;
00579 if ((*PlIt)>=0.) {
00580 fluem = EMFlu.CalcLongFluc((*PlIt)/cos(Theta_rad));
00581 inEM++;
00582 }
00583 double error = 0.;
00584 if ((*PhpeIt)<=0.001) error = fluem;
00585 else {
00586 MSG("SubShowerSR",Msg::kVerbose)
00587 << "fluem " << fluem << " *PhpeIt " << (*PhpeIt) << endl;
00588 error = TMath::Sqrt(TMath::Power(fluem,2) +
00589 TMath::Power(1./TMath::Sqrt(*PhpeIt),2));
00590 }
00591 if(inEM==1) error = TMath::Sqrt(TMath::Power(error,2)+1.);
00592 MSG("SubShowerSR",Msg::kVerbose)
00593 << (*PlIt) << " long ph pe " << (*PhpeIt)
00594 << " error " << error << " " << fluem << endl;
00595 LongFlu.push_back(error);
00596 PlIt++;
00597 PhpeIt++;
00598 }
00599 MSG("SubShowerSR",Msg::kVerbose)<<"z_adj "<<z_adj<<endl;
00600
00601 PlIt = LongPl.end();
00602 PlIt--;
00603 Double_t miss_z = (*PlIt);
00604 PlIt = LongPl.begin();
00605 t0 = *(LongPl.begin());
00606 while(miss_z-t0<2*tmax*X0_EFF/100.){
00607 miss_z+=2*D_EFF/100.;
00608 LongPl.push_back(miss_z);
00609 LongPh.push_back(0.);
00610 Double_t miss_flu = EMFlu.CalcLongFluc(miss_z/cos(Theta_rad));
00611 MSG("SubShowerSR",Msg::kVerbose)<<"miss_flu "<<miss_flu<<endl;
00612 LongFlu.push_back(miss_flu);
00613 }
00614 MSG("SubShowerSR",Msg::kVerbose)<< "0 loweststp " << loweststp
00615 << " uppeststp " << uppeststp << endl;
00616 Double_t Thetabd = atan(ST_WID/(Double_t(LongPl.size()) *
00617 D_EFF*2.))*180./TMath::Pi();
00618 Double_t dThetai = 1./180.*TMath::Pi();
00619 if(TMath::Abs(Thetabd)<5.) dThetai = TMath::Abs(Thetabd)/180.*TMath::Pi()/5.;
00620 Double_t maxph0 = 0.;
00621 Double_t maxvs0 = vs;
00622 Double_t maxtheta0 = Theta_rad;
00623 for(Int_t k=-6;k<=6;k++){
00624 for(Int_t j=-5;j<=5;j++){
00625 Double_t totph0 = 0.;
00626 for(Int_t i=0;i<nstp;i++){
00627 Double_t zt_vtx_corrected[2];
00628 zt_vtx_corrected[0] = z[i]-begz;
00629 zt_vtx_corrected[1] = tpos[i] - (vs+ST_WID/100.*k/4.) -
00630 (z[i]-begz) * tan(Theta_rad+j*dThetai);
00631 if(TMath::Abs(zt_vtx_corrected[1]-0.)<=T_EFF/2./100.&&ph[i]>0) {
00632 totph0+=ph[i];
00633 }
00634 }
00635 if(totph0>maxph0){
00636 maxph0 = totph0;
00637 maxvs0 = vs+ST_WID/100.*k/4.;
00638 maxtheta0 = Theta_rad+j*dThetai;
00639 MSG("SubShowerSR",Msg::kVerbose)
00640 <<"Theta_rad "<<Theta_rad<<" vs "<<vs
00641 <<" maxtheta0 "<<maxtheta0<<" maxvs0 "<<maxvs0
00642 <<" maxph0 "<<maxph0<<" j "<<j<<" k "<<k
00643 <<" dThetai "<<dThetai<<endl;
00644 }
00645 }
00646 }
00647 MSG("SubShowerSR",Msg::kDebug) <<"Theta_rad "<<Theta_rad
00648 <<" vs "<<vs
00649 <<" maxtheta0 "<<maxtheta0
00650 <<" maxvs0 "<<maxvs0<<endl;
00651 Theta_rad = maxtheta0;
00652 vs = maxvs0;
00653 for(Int_t j=int((loweststp-uppeststp)/2.)-1;
00654 j<=int((uppeststp-loweststp)/2.)+1;j++){
00655
00656 Double_t totphStp = 0.;
00657 Double_t totphStp_pe = 0.;
00658 Double_t stptpos = j*T_EFF/100.;
00659
00660 Double_t nstpTr = 0.;
00661 for(Int_t i=0;i<nstp;i++){
00662 Double_t zt_vtx_corrected[2];
00663 zt_vtx_corrected[0] = z[i]-begz;
00664 zt_vtx_corrected[1] = tpos[i]-vs-(z[i]-begz)*tan(Theta_rad);
00665
00666 if(TMath::Abs(zt_vtx_corrected[1]-stptpos)<=T_EFF/2./100.&&ph[i]>0) {
00667 totphStp+=ph[i];
00668 totphStp_pe+=ph_pe[i];
00669 nstpTr++;
00670 }
00671 }
00672 TranStp.push_back(stptpos);
00673 TranPh.push_back(totphStp);
00674 if(totphStp>0.) nn0t++;
00675 double fluem = EMFlu.CalcTranFluc(stptpos*cos(Theta_rad));
00676 double error = 0.;
00677 if (totphStp_pe<=0.001) error = fluem;
00678 else {
00679 MSG("SubShowerSR",Msg::kVerbose)<<"fluem "<<fluem
00680 <<" totphStp_pe "<<totphStp_pe<<endl;
00681 error = TMath::Sqrt(TMath::Power(fluem,2) +
00682 TMath::Power(1./TMath::Sqrt(totphStp_pe),2));
00683 MSG("SubShowerSR",Msg::kVerbose)<<"tran ph "<<totphStp_pe
00684 <<" error "<<error<<" "<<fluem<<endl;
00685 }
00686 TranFlu.push_back(error);
00687 }
00688
00689 std::list<Double_t>::iterator StpIt = TranStp.begin();
00690 Double_t miss_tn = *(TranStp.begin());
00691 StpIt = TranStp.end();
00692 StpIt--;
00693 Double_t miss_tp = *StpIt;
00694
00695 while(miss_tp-miss_tn<4*RM_EFF/100.){
00696 miss_tn-=T_EFF/100.;
00697 miss_tp+=T_EFF/100.;
00698 TranStp.push_front(miss_tn);
00699 TranStp.push_back(miss_tp);
00700 TranPh.push_front(0.);
00701 TranPh.push_back(0.);
00702 Double_t miss_flun = EMFlu.CalcLongFluc(miss_tn*cos(Theta_rad));
00703 Double_t miss_flup = EMFlu.CalcLongFluc(miss_tp*cos(Theta_rad));
00704 MSG("SubShowerSR",Msg::kVerbose)<<"miss_flun "<<miss_flun
00705 <<" miss_flup "<<miss_flup<<endl;
00706 TranFlu.push_front(miss_flun);
00707 TranFlu.push_back(miss_flup);
00708 loweststp--;
00709 uppeststp++;
00710 }
00711 }
00712
00713 delete [] plane;
00714 delete [] strip;
00715 delete [] ph;
00716 delete [] ph_pe;
00717 delete [] z;
00718 delete [] tpos;
00719 delete [] time;
00720 }
|
|
||||||||||||||||
|
Implements AlgBase. Definition at line 80 of file AlgSubShowerSR.cxx. References CandHandle::AddDaughterLink(), CalculateEnergyVertexAngle(), fDoEMFit, fECut, fMaxDiffCut, fMIPperGeV, fOut2RmPHCut, fPECut, fTrkFraCut, fTrkLikeStpCut, fVtxDiffCut, fXtkFraCut, CandContext::GetDataIn(), VldContext::GetDetector(), Registry::GetDouble(), Registry::GetInt(), CandStripHandle::GetPlaneView(), CandSubShowerSRHandle::GetPlaneView(), CandHandle::GetVldContext(), MSG, CandSubShowerSRHandle::SetPlaneView(), and SubShowerID(). 00081 {
00082
00083 MSG("SubShowerSR",Msg::kDebug) << "In AlgSubShowerSR::RunAlg" << endl;
00084 assert(ch.InheritsFrom("CandSubShowerSRHandle"));
00085 CandSubShowerSRHandle &cch = dynamic_cast<CandSubShowerSRHandle &>(ch);
00086 assert(cx.GetDataIn());
00087 assert(cx.GetDataIn()->InheritsFrom("TObjArray"));
00088 const TObjArray *tary = dynamic_cast<const TObjArray*>(cx.GetDataIn());
00089 if(tary->GetLast()<=-1) return;
00090 Int_t detector = 0;
00091 for (Int_t i=0; i<=tary->GetLast(); i++) {
00092 TObject *tobj = tary->At(i);
00093 assert(tobj->InheritsFrom("CandStripHandle"));
00094 CandStripHandle *striphandle = dynamic_cast<CandStripHandle*>(tobj);
00095 cch.AddDaughterLink(*striphandle);
00096 if(cch.GetPlaneView()==PlaneView::kUnknown)
00097 cch.SetPlaneView(striphandle->GetPlaneView());
00098 if (i==0) {
00099 const VldContext *vldcptr = striphandle->GetVldContext();
00100 assert(vldcptr);
00101 VldContext vldc = *vldcptr;
00102 if(vldc.GetDetector()==DetectorType::kFar) detector = 1;
00103 }
00104 }
00105
00106 fPECut = ac.GetDouble("PECut");
00107 fECut = ac.GetDouble("ECut");
00108 fMaxDiffCut = ac.GetDouble("MaxDiffCut");
00109 fVtxDiffCut = ac.GetDouble("VtxDiffCut");
00110 fOut2RmPHCut = ac.GetDouble("Out2RmPHCut");
00111 fTrkFraCut = ac.GetDouble("TrkFraCut");
00112 fTrkLikeStpCut = ac.GetDouble("TrkLikeStpCut");
00113 fXtkFraCut = ac.GetDouble("XtkFraCut");
00114 fDoEMFit = ac.GetInt("DoEMFit");
00115 fMIPperGeV = ac.GetDouble("MIPperGeV");
00116
00117 CalculateEnergyVertexAngle(cch,detector);
00118 SubShowerID(cch);
00119
00120 }
|
|
|
Definition at line 724 of file AlgSubShowerSR.cxx. References avgph, avgstpperpln, begz, bothfitprob, Chi2(), D_EFF, Eguess, endz, fDoEMFit, fTrkFraCut, fXtkFraCut, CandSubShowerSRHandle::GetEnergy(), CandRecoHandle::GetNStrip(), LongFlu, LongPh, LongPl, LongShwProfSamAng(), loweststp, maxdiff, maxStpPHinSS, MSG, name, nn0l, nn0t, out2Rmph, CandHandle::Print(), CandSubShowerSRHandle::SetClusterID(), CandSubShowerSRHandle::SetProbEM(), sub(), T_EFF, Theta_rad, totw, TranFlu, TranPh, TranShwProfSamAng(), TranStp, trkfra, uppeststp, vtxdiff, and xtkfra. Referenced by RunAlg(). 00724 {
00725
00726 ClusterType::ClusterType_t id = ClusterType::kUnknown;
00727
00728 if((xtkfra>fXtkFraCut || avgph<2*fPECut) && maxStpPHinSS<2.0 )
00729 id = ClusterType::kXTalk;
00730 else if(trkfra>fTrkFraCut || avgstpperpln<fTrkLikeStpCut)
00731 id = ClusterType::kTrkLike;
00732 else if(cch.GetEnergy()<fECut)
00733 id = ClusterType::kHadLike;
00734 else if(maxdiff>fMaxDiffCut)
00735 id = ClusterType::kHadLike;
00736 else if(vtxdiff>fVtxDiffCut)
00737 id = ClusterType::kHadLike;
00738 else if(out2Rmph>fOut2RmPHCut)
00739 id = ClusterType::kHadLike;
00740 else id = ClusterType::kEMLike;
00741
00742 cch.SetClusterID(id);
00743 MSG("SubShowerSR",Msg::kVerbose)<<"id "<<id<<" E "<<cch.GetEnergy()<<endl;
00744
00745 bothfitprob = 0.;
00746 Bool_t makeplots = false;
00747 MSG("SubShowerSR",Msg::kDebug)
00748 << "nn0l " << nn0l << " LongPl.size() " << LongPl.size()
00749 << " nn0t " << nn0t << " TranStp.size() " << TranStp.size() << endl;
00750
00751 //checks before trying EM fit:
00752 if(fDoEMFit &&
00753 (id==ClusterType::kEMLike ||
00754 (id==ClusterType::kHadLike && cch.GetEnergy()>=fECut)) &&
00755 LongPl.size()-nn0l <= 3 && TranStp.size()-nn0t <= 3 ){
00756
00757 MSG("SubShowerSR",Msg::kDebug)<<"in fitting "<<endl;
00758
00759 std::list<Double_t>::iterator PlIt = LongPl.begin();
00760 std::list<Double_t>::iterator PhIt = LongPh.begin();
00761 std::list<Double_t>::iterator PluIt = LongFlu.begin();
00762 std::list<Double_t>::iterator StpIt = TranStp.begin();
00763 std::list<Double_t>::iterator StphIt = TranPh.begin();
00764 std::list<Double_t>::iterator StpluIt = TranFlu.begin();
00765
00766 TF1 *clusterlong = NULL;
00767 if(gROOT->FindObject("clusterlong")){
00768 clusterlong = (TF1*) gROOT->FindObject("clusterlong");
00769 }
00770 else {
00771 clusterlong = new TF1("clusterlong",
00772 LongShwProfSamAng,
00773 0-D_EFF,(endz-begz+(*PlIt))*100+D_EFF,2);
00774 }
00775 TF1 *clustertran = NULL;
00776 if(gROOT->FindObject("clustertran")){
00777 clustertran = (TF1*) gROOT->FindObject("clustertran");
00778 }
00779 else {
00780 clustertran = new TF1("clustertran",
00781 TranShwProfSamAng,
00782 ((loweststp-uppeststp)/2.-1)*T_EFF,
00783 ((uppeststp-loweststp)/2.+1)*T_EFF,2);
00784 }
00785
00786 MSG("SubShowerSR",Msg::kVerbose)
00787 << "funct boundary " << ((loweststp-uppeststp)/2.-1)*T_EFF-T_EFF/2.
00788 << " " << ((uppeststp-loweststp)/2.+1)*T_EFF+T_EFF/2. << endl;
00789
00790 Double_t thetarange = Theta_rad*0.05;
00791 if(thetarange<10./180.*TMath::Pi()) thetarange = 10./180.*TMath::Pi();
00792 Double_t erange = Eguess*0.05;
00793 if(erange<0.5) erange = 0.5;
00794
00795 clusterlong->SetRange(0.-D_EFF/2.,(endz-begz+(*PlIt))*100+D_EFF);
00796 clusterlong->SetParameter(0,Eguess);
00797 clusterlong->SetParLimits(0,Eguess-erange,Eguess+erange);
00798 clusterlong->SetParameter(1,Theta_rad);
00799 clusterlong->SetParLimits(1,Theta_rad-thetarange,Theta_rad+thetarange);
00800
00801 clustertran->SetRange(((loweststp-uppeststp)/2.-1)*T_EFF,
00802 ((uppeststp-loweststp)/2.+1)*T_EFF);
00803 clustertran->SetParameter(0,Eguess);
00804 clustertran->SetParLimits(0,Eguess-erange,Eguess+erange);
00805 clustertran->SetParameter(1,Theta_rad);
00806 clustertran->SetParLimits(1,Theta_rad-thetarange,Theta_rad+thetarange);
00807
00808 TH1F *LongProf = new TH1F("longprof","Long Prof",
00809 int((endz-begz+(*PlIt))*100+2*D_EFF),
00810 0-D_EFF,(endz-begz+(*PlIt))*100+D_EFF);
00811 TH1F *LongHist = new TH1F("longhist","Long Hist",
00812 int((endz-begz+(*PlIt))*100+2*D_EFF),
00813 0-D_EFF,(endz-begz+(*PlIt))*100+D_EFF);
00814 TH1F *TranProf = new TH1F("tranprof","Tran Prof",
00815 int((uppeststp-loweststp+2)*T_EFF),
00816 ((loweststp-uppeststp)/2.-1)*T_EFF,
00817 ((uppeststp-loweststp)/2.+1)*T_EFF);
00818 TH1F *TranHist = new TH1F("tranhist","Tran Hist",
00819 int((uppeststp-loweststp+2)*T_EFF),
00820 ((loweststp-uppeststp)/2.-1)*T_EFF,
00821 ((uppeststp-loweststp)/2.+1)*T_EFF);
00822 TH1F *LongProfb = new TH1F("longprofb","Long Prof Both",
00823 int((endz-begz+(*PlIt))*100+2*D_EFF),
00824 0-D_EFF,(endz-begz+(*PlIt))*100+D_EFF);
00825 TH1F *LongHistb = new TH1F("longhistb","Long Hist Both",
00826 int((endz-begz+(*PlIt))*100+2*D_EFF),
00827 0-D_EFF,(endz-begz+(*PlIt))*100+D_EFF);
00828 TH1F *TranProfb = new TH1F("tranprofb","Tran Prof Both",
00829 int((uppeststp-loweststp+2)*T_EFF),
00830 ((loweststp-uppeststp)/2.-1)*T_EFF,
00831 ((uppeststp-loweststp)/2.+1)*T_EFF);
00832 TH1F *TranHistb = new TH1F("tranhistb","Tran Hist Both",
00833 int((uppeststp-loweststp+2)*T_EFF),
00834 ((loweststp-uppeststp)/2.-1)*T_EFF,
00835 ((uppeststp-loweststp)/2.+1)*T_EFF);
00836 MSG("SubShowerSR",Msg::kVerbose) << "loweststp " << loweststp
00837 << " uppeststp " << uppeststp << endl;
00838
00839 Double_t Minchi2Both = 100000;
00840 Double_t BestEBoth = Eguess;
00841 Double_t BestThetaBoth = Theta_rad;
00842 double bestnormlb = 1;
00843 double bestnormtb = 1;
00844 Bool_t fittedboth = false;
00845 Double_t dE = 0.1;
00846 Double_t Eini = Eguess-erange;
00847 Double_t Ein = Eini;
00848
00849 MSG("SubShowerSR",Msg::kDebug)
00850 << "done initial stuff " <<Eguess << " " << Ein
00851 << " " << Theta_rad << " " << thetarange << endl;
00852
00853 while(Ein>=Eini&&Ein<=Eguess+erange){
00854 MSG("SubShowerSR",Msg::kVerbose)<<"Ein "<<Ein<<endl;
00855 clusterlong->SetParameter(0,Ein);
00856 Double_t dTheta = 2.5/180.*TMath::Pi();
00857 Double_t Thetaini = Theta_rad-thetarange;
00858 Double_t Thetain = Thetaini;
00859 while(Thetain>=Thetaini&&Thetain<=Theta_rad+thetarange){
00860 MSG("SubShowerSR",Msg::kVerbose)<<"Thetain "<<Thetain<<endl;
00861 if(cos(Thetain)!=0.){
00862 clusterlong->SetParameter(1,Thetain);
00863 Double_t norml = 0.0;
00864 Double_t *funcl = new Double_t[LongPl.size()];
00865 Double_t *errl = new Double_t[LongPl.size()];
00866 Int_t nfpl = 0;
00867 UInt_t j = 0;
00868 Int_t nlong = 0;
00869 PhIt = LongPh.begin();
00870 while(PhIt!=LongPh.end()){
00871 if((*PhIt)>0.) nlong++;
00872 PhIt++;
00873 }
00874 PhIt = LongPh.begin();
00875 PlIt = LongPl.begin();
00876 while(PlIt!=LongPl.end()){
00877 funcl[j] = 0.;
00878 if((*PlIt)>=0.) funcl[j] = clusterlong->Eval((*PlIt)*100.);
00879 MSG("SubShowerSR",Msg::kVerbose)<<"funcl[j] "<<funcl[j]<<endl;
00880 norml += funcl[j];
00881 if (funcl[j]>0.||(*PlIt)<0.) nfpl++;
00882 PlIt++;
00883 j++;
00884 }
00885 MSG("SubShowerSR",Msg::kVerbose)<<"nfpl "<<nfpl<<endl;
00886 Double_t chi2Long = 0.;
00887 PlIt = LongPl.begin();
00888 PhIt = LongPh.begin();
00889 PluIt = LongFlu.begin();
00890 j = 0;
00891 MSG("SubShowerSR",Msg::kVerbose)
00892 << "nfpl " << nfpl << " norml " << norml << " nlong " << nlong
00893 << " Ein " << Ein << " Thetain " << Thetain << endl;
00894
00895 if(double(nfpl)/double(LongPl.size())>0.5 &&
00896 norml>1e-3 && nlong>1){
00897
00898 while(PlIt!=LongPl.end()){
00899 Double_t data = (*PhIt)/totw*norml;
00900 Double_t err = funcl[j]*(*PluIt);
00901 if ((*PluIt)==0.0) err = 0.0001;
00902 UInt_t sub = j;
00903 while(funcl[sub]<0.05*norml&&j==0) {
00904 sub++;
00905 err = (funcl[sub]+funcl[j])/2.*(*PluIt);
00906 MSG("SubShowerSR",Msg::kVerbose)<<"err "<<err<<" sub "<<endl;
00907 if(err>0.0001) break;
00908 else if(sub==LongPl.size()-1||(sub-j)>=3){
00909 err = 0.0001;
00910 break;
00911 }
00912 }
00913 sub = j;
00914 while(funcl[sub]<0.05*norml&&j>0) {
00915 sub--;
00916 err = (funcl[sub]+funcl[j])/2.*(*PluIt);
00917 MSG("SubShowerSR",Msg::kVerbose)<<"err "<<err<<" sub "<<endl;
00918 if(err>0.0001) break;
00919 else if(sub==0||(j-sub)>=3){
00920 err = 0.0001;
00921 break;
00922 }
00923 }
00924
00925 MSG("SubShowerSR",Msg::kDebug)
00926 << " LongPl.size() " << LongPl.size()
00927 << " Eguess " <<Eguess << " norml " << norml
00928 << " Theta " << Theta_rad << " j " << j
00929 << " LongPl.at(j) " << (*PlIt)
00930 << " err " << funcl[j] << " " << funcl[j-1]
00931 << " " << funcl[j+1] << " " << (*PluIt)
00932 << " err " << err << " data - funcl " << data-funcl[j] << endl;
00933
00934 if(err<0.0001) err = 0.0001;
00935 MSG("SubShowerSR",Msg::kVerbose)
00936 <<"err "<<err<<" funcl[j] "<<funcl[j]<<endl;
00937
00938 if(j==0)
00939 errl[j] = TMath::Sqrt(TMath::Power(err,2) +
00940 TMath::Power((funcl[j+1] -
00941 funcl[j])/2.,2));
00942 else if(j==LongPl.size()-1)
00943 errl[j] = TMath::Sqrt(TMath::Power(err,2) +
00944 TMath::Power((funcl[j] -
00945 funcl[j-1])/2.,2));
00946 else
00947 errl[j] = TMath::Sqrt(TMath::Power(err,2) +
00948 TMath::Power((funcl[j+1] -
00949 funcl[j-1])/4.,2));
00950
00951 MSG("SubShowerSR",Msg::kVerbose)
00952 << "data-funcl[j] " << data-funcl[j]
00953 << " errl[j] " << errl[j] << endl;
00954
00955 chi2Long += ( TMath::Power(data-funcl[j],2) /
00956 TMath::Power(errl[j],2) );
00957
00958 MSG("SubShowerSR",Msg::kVerbose)
00959 << "chi2long " << (*PlIt) << " " << data
00960 << " " << funcl[j] << " " << (*PluIt) << " " << err
00961 << " " << TMath::Power(data-funcl[j],2)/TMath::Power(err,2)
00962 << " " << chi2Long << endl;
00963
00964 PlIt++;
00965 PhIt++;
00966 PluIt++;
00967 j++;
00968 }
00969
00970 MSG("SubShowerSR",Msg::kVerbose)
00971 << "long fit " << chi2Long/(LongPl.size()-1)
00972 << "norml " << LongPl.size() << " " << norml << endl;
00973 }
00974 else chi2Long = 100000;
00975
00976 Double_t normt = 0.0;
00977 Double_t *funct = new Double_t[TranStp.size()];
00978 Double_t *errt = new Double_t[TranStp.size()];
00979 Int_t nfpt = 0;
00980 j = 0;
00981 Int_t ntran = 0;
00982 StphIt = TranPh.begin();
00983 while(StphIt!=TranPh.end()){
00984 if(*StphIt>0.) ntran++;
00985 StphIt++;
00986 }
00987 StphIt = TranPh.begin();
00988 StpIt = TranStp.begin();
00989 while(StpIt!=TranStp.end()){
00990 funct[j] = clustertran->Eval(*StpIt*100.);
00991 normt += funct[j];
00992 if(funct[j]>0.) nfpt++;
00993 StpIt++;
00994 j++;
00995 }
00996
00997 Double_t chi2Tran = 0.;
00998 MSG("SubShowerSR",Msg::kVerbose)<<"nfpt "<<nfpt<<" "<<normt<<endl;
00999 if(double(nfpt)/double(TranStp.size())>0.5 &&
01000 normt>1e-3 && ntran>1){
01001
01002 StpIt = TranStp.begin();
01003 StphIt = TranPh.begin();
01004 StpluIt = TranFlu.begin();
01005 j = 0;
01006
01007 while(StpIt!=TranStp.end()){
01008 Double_t data = *StphIt/totw*normt;
01009 Double_t err = funct[j]*(*StpluIt);
01010 if (*StpluIt==0.0) err = 0.0001;
01011 UInt_t sub = j;
01012 while(funct[sub]<0.05*normt&&*StpIt<=0) {
01013 sub++;
01014 err = (funct[sub]+funct[j])/2.*(*StpluIt);
01015 if(err>0.0001) break;
01016 else if(sub==TranStp.size()-1||(sub-j)>=3){
01017 err = 0.0001;
01018 break;
01019 }
01020 }
01021 sub = j;
01022 while(funct[sub]<0.05*normt && *StpIt>0) {
01023 sub--;
01024 err = (funct[sub]+funct[j])/2.*(*StpluIt);
01025 if(err>0.0001) break;
01026 else if(sub==0||(j-sub)>=3){
01027 err = 0.0001;
01028 break;
01029 }
01030 }
01031 MSG("SubShowerSR",Msg::kDebug)
01032 << "funct " << j << " " << *StpIt << " " << *StpluIt
01033 << " err " << " " << funct[j] << " " << funct[j-1]
01034 << " " << funct[j+1] << " err " << err
01035 << " data - funcl " << data-funcl[j] << endl;
01036
01037 if(err<0.0001) err = 0.0001;
01038 if(j==0)
01039 errt[j] = TMath::Sqrt(TMath::Power(err,2) +
01040 TMath::Power((funct[j+1] -
01041 funct[j])/2.,2));
01042 else if(j==TranStp.size()-1)
01043 errt[j] = TMath::Sqrt(TMath::Power(err,2) +
01044 TMath::Power((funct[j] -
01045 funct[j-1])/2.,2));
01046 else
01047 errt[j] = TMath::Sqrt(TMath::Power(err,2) +
01048 TMath::Power((funct[j+1] -
01049 funct[j-1])/4.,2));
01050
01051 MSG("SubShowerSR",Msg::kVerbose) << "data-funcl[j] " <<
01052 data-funcl[j] << " errt[j] " << errt[j] << endl;
01053
01054 chi2Tran += ( TMath::Power(data-funct[j],2) /
01055 (TMath::Power(errt[j],2)) );
01056
01057 MSG("SubShowerSR",Msg::kVerbose) << "chi2tran " <<
01058 TMath::Power(data-funcl[j],2)/TMath::Power(err,2)
01059 << " " << chi2Tran << endl;
01060
01061 StpIt++;
01062 StphIt++;
01063 StpluIt++;
01064 j++;
01065 }
01066
01067 MSG("SubShowerSR",Msg::kVerbose)
01068 << "tran fit " << chi2Tran/(TranStp.size()-1)
01069 << "normt " << TranStp.size() << " " << normt << endl;
01070 }
01071 else chi2Tran = 100000;
01072
01073 bothfitprob = 0.;
01074 Double_t chi2Both = 0.;
01075
01076 MSG("SubShowerSR",Msg::kVerbose) << "chi2Tran " <<chi2Tran
01077 << " chi2Long " << chi2Long << endl;
01078
01079 chi2Both = chi2Tran+chi2Long;
01080 if((double(nfpl)/double(LongPl.size())>0.5&&norml>1e-3&&nlong>1)&&
01081 (double(nfpt)/double(TranStp.size())>0.5&&normt>1e-3&&ntran>1)&&
01082 chi2Both<=Minchi2Both){
01083 Minchi2Both = chi2Both;
01084 BestEBoth = Ein;
01085 BestThetaBoth = Thetain;
01086 bestnormtb = normt;
01087 bestnormlb = norml;
01088 fittedboth = true;
01089 if(makeplots){
01090 PlIt = LongPl.begin();
01091 PhIt = LongPh.begin();
01092 PluIt = LongFlu.begin();
01093 j = 0;
01094 while(PlIt!=LongPl.end()){
01095 double data = (*PhIt)/totw*norml;
01096 Int_t bin = LongProfb->FindBin((*PlIt)*100.);
01097 LongProfb->SetBinContent(bin,funcl[j]);
01098 LongProfb->SetBinError(bin,errl[j]);
01099 LongHistb->SetBinContent(bin,data);
01100 LongHistb->SetBinError(bin,0.0);
01101 PlIt++;
01102 PhIt++;
01103 PluIt++;
01104 j++;
01105 }
01106 StpIt = TranStp.begin();
01107 StphIt = TranPh.begin();
01108 StpluIt = TranFlu.begin();
01109 j = 0;
01110
01111 while(StpIt!=TranStp.end()){
01112 double data = *StphIt/totw*normt;
01113 Int_t bin = TranProfb->FindBin(*StpIt*100.);
01114 TranProfb->SetBinContent(bin,funct[j]);
01115 TranProfb->SetBinError(bin,errt[j]);
01116 TranHistb->SetBinContent(bin,data);
01117 TranHistb->SetBinError(bin,0.);
01118 StpIt++;
01119 StphIt++;
01120 StpluIt++;
01121 j++;
01122 }
01123 }
01124 }
01125 delete [] funcl;
01126 delete [] errl;
01127 delete [] funct;
01128 delete [] errt;
01129 }
01130 Thetain+=dTheta;
01131 }
01132 Ein+=dE;
01133 }
01134
01135 if(fittedboth){
01136 if(LongPl.size()>1&&TranStp.size()>1){
01137 TF1 *form = new TF1("form",Chi2,0,5000,1);
01138 Int_t ndfLong = LongPl.size();
01139 Int_t ndfTran = TranStp.size();
01140 Int_t ndfBoth = ndfLong+ndfTran;
01141 form->SetParameter(0,ndfBoth);
01142 if(fittedboth&&Minchi2Both/ndfBoth<100.)
01143 bothfitprob = 1.-form->Integral(0.0,Minchi2Both);
01144 else bothfitprob = 0.0;
01145 delete form;
01146
01147 MSG("SubShowerSR",Msg::kDebug)
01148 << "Best Fit: Both- " << BestEBoth << " GeV "
01149 << BestThetaBoth << " rad with chi2/ndof "
01150 << Minchi2Both/ndfBoth << " ndf " << ndfBoth
01151 << " of probability " << bothfitprob
01152 << " comparing to " << Eguess << " GeV "
01153 << Theta_rad << " rad" << endl;
01154 }
01155 if(makeplots){
01156 gStyle->SetOptStat(0);
01157 char name[256];
01158 TCanvas *can = new TCanvas();
01159 can->cd();
01160 can->Divide(2,1);
01161 can->cd(1);
01162 LongHistb->SetMaximum(bestnormlb);
01163 LongHistb->Draw("EP");
01164 LongProfb->SetMarkerColor(4);
01165 LongProfb->Draw("EPsame");
01166 can->cd(2);
01167 TranHistb->SetMaximum(bestnormtb);
01168 TranHistb->Draw("EP");
01169 TranProfb->SetMarkerColor(4);
01170 TranProfb->Draw("EPsame");
01171 if(bestnormlb+bestnormtb>0){
01172 sprintf(name,"Fit%f_%fb.eps",Eguess,Theta_rad);
01173 can->Print(name);
01174 }
01175 delete can;
01176 }
01177 }
01178 delete LongHist;
01179 LongHist = NULL;
01180 delete LongProf;
01181 LongProf = NULL;
01182 delete TranHist;
01183 TranHist = NULL;
01184 delete TranProf;
01185 TranProf = NULL;
01186 delete LongHistb;
01187 LongHistb = NULL;
01188 delete LongProfb;
01189 LongProfb = NULL;
01190 delete TranHistb;
01191 TranHistb = NULL;
01192 delete TranProfb;
01193 TranProfb = NULL;
01194 delete clusterlong;
01195 clusterlong = NULL;
01196 delete clustertran;
01197 clustertran = NULL;
01198 }
01199 else {
01200 bothfitprob = 0.;
01201 }
01202 MSG("SubShowerSR",Msg::kVerbose) << "long " << LongPl.size()
01203 << " tran " <<TranStp.size()
01204 << " id " << id
01205 << " eng " << cch.GetEnergy()
01206 << " nstrip " << cch.GetNStrip()
01207 << " Prob " << bothfitprob << endl;
01208 cch.SetProbEM(Float_t(bothfitprob));
01209 }
|
|
|
Reimplemented from AlgBase. Definition at line 1212 of file AlgSubShowerSR.cxx. 01213 {
01214 }
|
|
|
Definition at line 50 of file AlgSubShowerSR.h. Referenced by CalculateEnergyVertexAngle(), and SubShowerID(). |
|
|
Definition at line 48 of file AlgSubShowerSR.h. Referenced by CalculateEnergyVertexAngle(), and SubShowerID(). |
|
|
Definition at line 66 of file AlgSubShowerSR.h. Referenced by CalculateEnergyVertexAngle(), and SubShowerID(). |
|
|
Definition at line 52 of file AlgSubShowerSR.h. Referenced by SubShowerID(). |
|
|
Definition at line 70 of file AlgSubShowerSR.h. Referenced by CalculateEnergyVertexAngle(), and SubShowerID(). |
|
|
Definition at line 67 of file AlgSubShowerSR.h. Referenced by CalculateEnergyVertexAngle(), and SubShowerID(). |
|
|
Definition at line 41 of file AlgSubShowerSR.h. Referenced by RunAlg(), and SubShowerID(). |
|
|
Definition at line 34 of file AlgSubShowerSR.h. Referenced by RunAlg(). |
|
|
Definition at line 35 of file AlgSubShowerSR.h. Referenced by RunAlg(). |
|
|
Definition at line 42 of file AlgSubShowerSR.h. Referenced by RunAlg(). |
|
|
Definition at line 37 of file AlgSubShowerSR.h. Referenced by RunAlg(). |
|
|
Definition at line 33 of file AlgSubShowerSR.h. Referenced by RunAlg(). |
|
|
Definition at line 38 of file AlgSubShowerSR.h. Referenced by RunAlg(), and SubShowerID(). |
|
|
Definition at line 39 of file AlgSubShowerSR.h. Referenced by RunAlg(). |
|
|
Definition at line 36 of file AlgSubShowerSR.h. Referenced by RunAlg(). |
|
|
Definition at line 40 of file AlgSubShowerSR.h. Referenced by RunAlg(), and SubShowerID(). |
|
|
Definition at line 58 of file AlgSubShowerSR.h. Referenced by CalculateEnergyVertexAngle(), SubShowerID(), and ~AlgSubShowerSR(). |
|
|
Definition at line 57 of file AlgSubShowerSR.h. Referenced by CalculateEnergyVertexAngle(), SubShowerID(), and ~AlgSubShowerSR(). |
|
|
Definition at line 59 of file AlgSubShowerSR.h. Referenced by CalculateEnergyVertexAngle(), and ~AlgSubShowerSR(). |
|
|
Definition at line 56 of file AlgSubShowerSR.h. Referenced by CalculateEnergyVertexAngle(), SubShowerID(), and ~AlgSubShowerSR(). |
|
|
Definition at line 64 of file AlgSubShowerSR.h. Referenced by CalculateEnergyVertexAngle(), and SubShowerID(). |
|
|
Definition at line 45 of file AlgSubShowerSR.h. Referenced by CalculateEnergyVertexAngle(), and SubShowerID(). |
|
|
Definition at line 44 of file AlgSubShowerSR.h. Referenced by CalculateEnergyVertexAngle(), and SubShowerID(). |
|
|
Definition at line 71 of file AlgSubShowerSR.h. Referenced by CalculateEnergyVertexAngle(), and SubShowerID(). |
|
|
Definition at line 72 of file AlgSubShowerSR.h. Referenced by CalculateEnergyVertexAngle(), and SubShowerID(). |
|
|
Definition at line 51 of file AlgSubShowerSR.h. Referenced by CalculateEnergyVertexAngle(), and SubShowerID(). |
|
|
Definition at line 53 of file AlgSubShowerSR.h. |
|
|
Definition at line 54 of file AlgSubShowerSR.h. |
|
|
Definition at line 68 of file AlgSubShowerSR.h. Referenced by CalculateEnergyVertexAngle(), and SubShowerID(). |
|
|
Definition at line 69 of file AlgSubShowerSR.h. Referenced by CalculateEnergyVertexAngle(), and SubShowerID(). |
|
|
Definition at line 62 of file AlgSubShowerSR.h. Referenced by CalculateEnergyVertexAngle(), SubShowerID(), and ~AlgSubShowerSR(). |
|
|
Definition at line 61 of file AlgSubShowerSR.h. Referenced by CalculateEnergyVertexAngle(), SubShowerID(), and ~AlgSubShowerSR(). |
|
|
Definition at line 60 of file AlgSubShowerSR.h. Referenced by CalculateEnergyVertexAngle(), SubShowerID(), and ~AlgSubShowerSR(). |
|
|
Definition at line 47 of file AlgSubShowerSR.h. Referenced by CalculateEnergyVertexAngle(), and SubShowerID(). |
|
|
Definition at line 65 of file AlgSubShowerSR.h. Referenced by CalculateEnergyVertexAngle(), and SubShowerID(). |
|
|
Definition at line 46 of file AlgSubShowerSR.h. Referenced by CalculateEnergyVertexAngle(), and SubShowerID(). |
|
|
Definition at line 49 of file AlgSubShowerSR.h. Referenced by CalculateEnergyVertexAngle(), and SubShowerID(). |
1.3.9.1