00001
00002
00003
00004
00005
00006
00007
00009
00010 #include <cassert>
00011 #include <vector>
00012 #include "Algorithm/AlgConfig.h"
00013 #include "Candidate/CandContext.h"
00014 #include "Conventions/PlaneView.h"
00015 #include "MessageService/MsgService.h"
00016 #include "Validity/VldContext.h"
00017 #include "Calibrator/Calibrator.h"
00018 #include "CandSubShowerSR/AlgSubShowerSR.h"
00019 #include "CandSubShowerSR/CandSubShowerSRHandle.h"
00020 #include "CandSubShowerSR/ClusterType.h"
00021 #include "CandSubShowerSR/EMFluctuation.h"
00022 #include "CandFitShowerEM/BinFluctuationEM.h"
00023 #include "RecoBase/CandStripHandle.h"
00024 #include "TMath.h"
00025 #include "TMatrixD.h"
00026 #include "TVectorD.h"
00027 #include "TVector2.h"
00028 #include "TMatrixDSym.h"
00029 #include "TMatrixDSymEigen.h"
00030 #include "TCanvas.h"
00031 #include "TStyle.h"
00032
00033 #define D_EFF 6.0
00034 #define T_EFF 4.0
00035 #define EC_EFF 21.6297
00036 #define X0_EFF 4.1542
00037 #define RM_EFF 4.0717
00038 #define ST_WID 4.1
00039
00040 Double_t LongShwProfSamAng(double *,double *);
00041 Double_t TranShwProfSamAng(double *,double *);
00042 Double_t LongShwProfSamn(double *,double *);
00043 Double_t TranShwProfSamn(double *,double *);
00044 Bool_t SortStripByEnergy(Double_t *,Int_t,Int_t *);
00045 Double_t StpFluctuaionAng(double *,double *);
00046 Double_t Chi2(double *,double *);
00047
00048 ClassImp(AlgSubShowerSR)
00049
00050
00051 CVSID("$Id: AlgSubShowerSR.cxx,v 1.20 2006/08/10 11:35:30 cbs Exp $");
00052
00053
00054 AlgSubShowerSR::AlgSubShowerSR() :
00055 fPECut(1.5),fECut(0.1),fMaxDiffCut(1),fVtxDiffCut(2),
00056 fOut2RmPHCut(0.3),fTrkFraCut(0.8),fTrkLikeStpCut(1.3),
00057 fXtkFraCut(0.8),fMaxXTalkStpPH(2.0),fDoEMFit(1),fMIPperGeV(18.5)
00058 {
00059 flspl = new TH1F("flspl","flspl",500,-0.5,499.5);
00060 sflspl = new TH1F("sflspl","sflspl",500,-0.5,499.5);
00061 clusterlong = new TF1("clusterlong",LongShwProfSamAng,0,1,2);
00062 clustertran = new TF1("clustertran",TranShwProfSamAng,0,1,2);
00063 chi2func = new TF1("chi2func",Chi2,0,5000,1);
00064
00065 LongProf = new TH1F("longprof","Long Prof",1,0,1);
00066 LongHist = new TH1F("longhist","Long Hist",1,0,1);
00067 TranProf = new TH1F("tranprof","Tran Prof",1,0,1);
00068 TranHist = new TH1F("tranhist","Tran Hist",1,0,1);
00069 LongProfb = new TH1F("longprofb","Long Prof Both",1,0,1);
00070 LongHistb = new TH1F("longhistb","Long Hist Both",1,0,1);
00071 TranProfb = new TH1F("tranprofb","Tran Prof Both",1,0,1);
00072 TranHistb = new TH1F("tranhistb","Tran Hist Both",1,0,1);
00073 }
00074
00075
00076 AlgSubShowerSR::~AlgSubShowerSR()
00077 {
00078
00079 delete flspl;
00080 delete sflspl;
00081 delete clusterlong;
00082 delete clustertran;
00083 delete chi2func;
00084
00085 delete LongProf;
00086 delete LongHist;
00087 delete TranProf;
00088 delete TranHist;
00089 delete LongProfb;
00090 delete LongHistb;
00091 delete TranProfb;
00092 delete TranHistb;
00093
00094 LongPl.clear();
00095 LongPh.clear();
00096 LongFlu.clear();
00097 LongPhpe.clear();
00098 TranStp.clear();
00099 TranPh.clear();
00100 TranFlu.clear();
00101
00102 }
00103
00104
00105
00106 void AlgSubShowerSR::RunAlg(AlgConfig &ac, CandHandle &ch, CandContext &cx)
00107 {
00108
00109 MSG("SubShowerSR",Msg::kDebug) << "In AlgSubShowerSR::RunAlg" << endl;
00110 assert(ch.InheritsFrom("CandSubShowerSRHandle"));
00111 CandSubShowerSRHandle &cch = dynamic_cast<CandSubShowerSRHandle &>(ch);
00112 assert(cx.GetDataIn());
00113 assert(cx.GetDataIn()->InheritsFrom("TObjArray"));
00114 const TObjArray *tary = dynamic_cast<const TObjArray*>(cx.GetDataIn());
00115 if(tary->GetLast()<=-1) return;
00116 Int_t detector = 0;
00117 for (Int_t i=0; i<=tary->GetLast(); i++) {
00118 TObject *tobj = tary->At(i);
00119 assert(tobj->InheritsFrom("CandStripHandle"));
00120 CandStripHandle *striphandle = dynamic_cast<CandStripHandle*>(tobj);
00121 cch.AddDaughterLink(*striphandle);
00122 if(cch.GetPlaneView()==PlaneView::kUnknown)
00123 cch.SetPlaneView(striphandle->GetPlaneView());
00124 if (i==0) {
00125 const VldContext *vldcptr = striphandle->GetVldContext();
00126 assert(vldcptr);
00127 VldContext vldc = *vldcptr;
00128 if(vldc.GetDetector()==Detector::kFar) detector = 1;
00129 }
00130 }
00131
00132 fPECut = ac.GetDouble("PECut");
00133 fECut = ac.GetDouble("ECut");
00134 fMaxDiffCut = ac.GetDouble("MaxDiffCut");
00135 fVtxDiffCut = ac.GetDouble("VtxDiffCut");
00136 fOut2RmPHCut = ac.GetDouble("Out2RmPHCut");
00137 fTrkFraCut = ac.GetDouble("TrkFraCut");
00138 fTrkLikeStpCut = ac.GetDouble("TrkLikeStpCut");
00139 fXtkFraCut = ac.GetDouble("XtkFraCut");
00140 fMaxXTalkStpPH = ac.GetDouble("MaxXTalkStpPH");
00141 fDoEMFit = ac.GetInt("DoEMFit");
00142 fMIPperGeV = ac.GetDouble("MIPperGeV");
00143
00144 CalculateEnergyVertexAngle(cch,detector);
00145 SubShowerID(cch);
00146
00147 }
00148
00149
00150 void AlgSubShowerSR::CalculateEnergyVertexAngle(CandSubShowerSRHandle &cch,
00151 Int_t det)
00152 {
00153
00154 Int_t nstp = cch.GetNStrip();
00155 Int_t begpl = 500;
00156 Int_t endpl = -1;
00157 CandStripHandleItr stripItr(cch.GetDaughterIterator());
00158 while (CandStripHandle *stp = stripItr()){
00159 Int_t pln = stp->GetPlane();
00160 if(pln<begpl) begpl = pln;
00161 if(pln>endpl) endpl = pln;
00162 }
00163
00164 std::vector<Int_t> plane(nstp,0);
00165 std::vector<Int_t> strip(nstp,0);
00166 std::vector<Double_t> ph(nstp,0);
00167 std::vector<Double_t> ph_pe(nstp,0);
00168 std::vector<Double_t> z(nstp,0);
00169 std::vector<Double_t> tpos(nstp,0);
00170 std::vector<Double_t> time(nstp,0);
00171
00172 totw = 0.;
00173 Double_t totph = 0.;
00174 Int_t Ns0 = 0;
00175 Int_t icnt = 0;
00176
00177 maxStpPHinSS = 0;
00178 loweststp = 200;
00179 uppeststp = -200;
00180 Double_t lowestpl = begpl;
00181 Double_t uppestpl = endpl;
00182
00183 Calibrator& calibrator=Calibrator::Instance();
00184
00185 stripItr.Reset();
00186 while (CandStripHandle *stp = stripItr()) {
00187
00188 plane[icnt] = stp->GetPlane();
00189 strip[icnt] = stp->GetStrip();
00190 z[icnt] = stp->GetZPos();
00191 tpos[icnt] = stp->GetTPos();
00192 time[icnt] = stp->GetTime();
00193
00194 ph[icnt] = calibrator.GetMIP(stp->GetCharge(CalDigitType::kSigCorr));
00195 totw += ph[icnt];
00196
00197 ph_pe[icnt] = stp->GetCharge(CalDigitType::kPE);
00198 totph += ph_pe[icnt];
00199
00200 if(tpos[icnt]/T_EFF*100.>uppeststp) {
00201 uppeststp = tpos[icnt]/T_EFF*100.;
00202 uppestpl = plane[icnt];
00203 }
00204 if(tpos[icnt]/T_EFF*100.<loweststp) {
00205 loweststp = tpos[icnt]/T_EFF*100.;
00206 lowestpl = plane[icnt];
00207 }
00208 if(ph[icnt]>maxStpPHinSS) maxStpPHinSS = ph[icnt];
00209 if (ph_pe[icnt]<=fPECut) Ns0++;
00210 icnt++;
00211 }
00212
00213 maxdiff = 0.;
00214 vtxdiff = 0.;
00215 trkfra = 0;
00216 avgstpperpln = 0;
00217 xtkfra = 0;
00218 avgph = 0;
00219 out2Rmph = 0;
00220
00221 Int_t maxpl = -1;
00222 Double_t maxz = 0.;
00223 begz = 0.;
00224 endz = 30.;
00225 Int_t maxstp = -1;
00226 Double_t maxtpos = 0.;
00227 Int_t begstp = -1;
00228 Double_t mint = 999999999;
00229 Float_t normpos = 0;
00230
00231 Double_t Theta = 0;
00232 Theta_rad = 0;
00233 Double_t Theta1 = 0;
00234 Double_t Theta_rad1 = 0;
00235 Double_t avgdev = 0.;
00236 Double_t wavgdev = 0.;
00237 Double_t avgdev1 = 0.;
00238 Double_t wavgdev1 = 0.;
00239 Double_t avgdev2 = 0.;
00240 Double_t wavgdev2 = 0.;
00241 Double_t eval[2] = {0};
00242
00243 Eguess = totw*2./fMIPperGeV;
00244 if(Eguess==0) {
00245 MSG("SubShowerSR",Msg::kError) << "Eguess = 0" << endl;
00246 return;
00247 }
00248
00249 Double_t tmax = TMath::Log(Eguess*1000./EC_EFF)-0.858;
00250 Double_t tmaxpl = tmax*X0_EFF/D_EFF;
00251 flspl->Reset();
00252 sflspl->Reset();
00253
00254 for(Int_t i=0;i<nstp;i++){
00255 flspl->Fill(plane[i],ph[i]);
00256 }
00257 maxpl = Int_t(flspl->GetBinCenter(flspl->GetMaximumBin()));
00258 Int_t getz = 0;
00259 for(Int_t i=0;i<nstp;i++){
00260 if (plane[i]!=maxpl) sflspl->Fill(plane[i],ph[i]);
00261 else if (getz<1&&plane[i]==maxpl) {
00262 maxz = z[i];
00263 getz++;
00264 }
00265 }
00266 double secmaxpl = Int_t(sflspl->GetBinCenter(sflspl->GetMaximumBin()));
00267 maxdiff = TMath::Abs(maxpl-secmaxpl)/tmaxpl;
00268 vtxdiff = TMath::Abs(maxpl-begpl)/tmaxpl;
00269 MSG("SubShowerSR",Msg::kVerbose)
00270 << "vtxdiff " <<vtxdiff << " tmaxpl " << tmaxpl <<endl;
00271
00272 Int_t Ns = 0;
00273 Int_t Ns1 = 0;
00274 Double_t stpden = 0.;
00275 Double_t stpsp = 0.;
00276 double phplbeg = 0.;
00277 double phplend = 0.;
00278 Double_t avgbeg = 0.;
00279 Double_t posdif = 0.;
00280 double preavg = 0.;
00281 double prephpl = 0.;
00282 double prenstppl = 0.;
00283
00284 for (int pl = begpl;pl<=endpl;pl++){
00285 int nstppl = 0;
00286 double low = 200.;
00287 double high = -200.;
00288 double phpl = 0.;
00289 double avg = 0.;
00290 for(int i=0;i<nstp;i++){
00291 if(plane[i]==pl) {
00292 nstppl++;
00293 phpl += ph[i];
00294 avg +=tpos[i]*ph[i];
00295 if(tpos[i]/T_EFF*100.>high) high = tpos[i]/T_EFF*100.;
00296 if(tpos[i]/T_EFF*100.<low) low = tpos[i]/T_EFF*100.;
00297 }
00298 }
00299 if (nstppl>0) {
00300 Ns++;
00301 avgstpperpln += Double_t(nstppl);
00302 stpden += double(nstppl)*phpl/(high-low+1.);
00303 avg = avg/T_EFF*100./phpl;
00304 if (nstppl==1) Ns1++;
00305 if(pl>begpl && phplbeg!=0. && preavg!=0.) {
00306 double avgdif = ( (phpl + prephpl)*(avg - preavg) /
00307 (2.*double(nstppl + prenstppl) ) );
00308 stpsp += TMath::Abs(avgdif);
00309 if (avgdif>=0) posdif += 1;
00310 MSG("SubShowerSR",Msg::kVerbose)
00311 << "avg " << avg << " preavg " << preavg << " nstppl " << nstppl
00312 << " prenstppl " << prenstppl << " phpl " << phpl << " prephpl "
00313 << " avgdif " << avgdif << " posdif " << posdif
00314 << " stpsp " << stpsp << endl;
00315 }
00316 if(pl>=begpl&&phplbeg==0.) {
00317 phplbeg = phpl;
00318 avgbeg = avg;
00319 }
00320 if(pl<=endpl) phplend = phpl;
00321 preavg = avg;
00322 prephpl = phpl;
00323 prenstppl = nstppl;
00324 }
00325 }
00326 stpden = stpden/totw;
00327
00328 if (Ns>1) {
00329 stpsp = stpsp/(totw-(phplbeg+phplend)/2.);
00330 posdif = posdif/double(Ns-1);
00331 }
00332 else {
00333 stpsp = 100.;
00334 posdif = 0.5;
00335 }
00336 MSG("SubShowerSR",Msg::kVerbose) << "stpden " << stpden << " stpsp "
00337 << stpsp << " posdif " << posdif << endl;
00338
00339 if(endpl==begpl) {
00340 avgstpperpln = 999;
00341 trkfra = 0;
00342 }
00343 else {
00344 avgstpperpln /= (Double_t(endpl - begpl)/2. + 1.);
00345 trkfra = double(Ns1)/double(Ns);
00346 }
00347 xtkfra = double(Ns0)/double(nstp);
00348 avgph = totph/double(nstp);
00349 Float_t maxstpph = 0.;
00350 Float_t begstpph = 0.;
00351 Int_t nstpmaxpl = 0;
00352
00353 for(Int_t i=0;i<nstp;i++){
00354 if(plane[i]==maxpl&&ph[i]>maxstpph){
00355 maxstp = strip[i];
00356 maxstpph = ph[i];
00357 maxtpos = tpos[i];
00358 }
00359 if(plane[i]==maxpl){
00360 maxz = z[i];
00361 nstpmaxpl++;
00362 }
00363 if(plane[i]==begpl){
00364 if(ph[i]>begstpph){
00365 begz = z[i];
00366 begstp = strip[i];
00367 begstpph = ph[i];
00368 }
00369 if(time[i]<mint) mint = time[i];
00370 }
00371 Bool_t gotendz = false;
00372 if(!gotendz&&plane[i]==endpl){
00373 endz = z[i];
00374 gotendz = true;
00375 }
00376 }
00377
00378 Float_t ss[2][2] = {{0}};
00379 for (int m = 0;m<2;m++){
00380 for (int n = 0;n<2;n++){
00381 ss[m][n] = 0.;
00382 }
00383 }
00384
00385 Double_t cnt = 0;
00386 for(Int_t i=0;i<nstp;i++){
00387 Double_t dtpos = tpos[i] - maxtpos;
00388 Double_t dz = z[i] - maxz;
00389 ss[0][0]+=dz*dz*ph[i];
00390 ss[0][1]+=dtpos*dz*ph[i];
00391 ss[1][1]+=dtpos*dtpos*ph[i];
00392 cnt+=1;
00393 normpos += dtpos*dtpos+dz*dz;
00394 }
00395 ss[1][0]=ss[0][1];
00396
00397 if(totw*normpos!=0){
00398
00399 TMatrixDSym ms(2);
00400 for (int m = 0;m<2;m++){
00401 for (int n = 0;n<2;n++){
00402 ms[m][n] = 0.;
00403 }
00404 }
00405
00406 Double_t ang[2];
00407 for(Int_t i=0;i<2;i++) ang[i] = -999;
00408 for (Int_t i=0;i<2;i++){
00409 for (Int_t j=0;j<2;j++){
00410 ms(i,j)=ss[i][j]/totw/normpos;
00411 }
00412 }
00413
00414 if(ms.IsA()){
00415 TMatrixDSymEigen eigen(ms);
00416 TMatrixD EVect1 = eigen.GetEigenVectors();
00417 TVectorD EVal1 = eigen.GetEigenValues();
00418 TVector2 ev0(EVect1(0,0),EVect1(1,0));
00419 TVector2 ev1(EVect1(0,1),EVect1(1,1));
00420 ang[0] = ev0.Phi()*180/TMath::Pi();
00421 ang[1] = ev1.Phi()*180/TMath::Pi();
00422 eval[0] = EVal1(0);
00423 eval[1] = EVal1(1);
00424 }
00425 else {
00426 eval[0] = 0.;
00427 eval[1] = 0.;
00428 }
00429 for (Int_t i=0;i<2;i++){
00430 if (ang[i]>90&&ang[i]<270) ang[i]-=180;
00431 if (ang[i]>270&&ang[i]<=360) ang[i]-=360;
00432 if (ang[i]==90||ang[i]==270) ang[i]-=0.01;
00433 }
00434
00435 Theta = ang[0];
00436 Theta1 = ang[1];
00437 Theta_rad = (Theta/180)*TMath::Pi();
00438 Theta_rad1 = (Theta1/180)*TMath::Pi();
00439 double vs = maxtpos+tan(Theta_rad)*(begz-maxz);
00440 double vs1 = maxtpos+tan(Theta_rad1)*(begz-maxz);
00441 MSG("SubShowerSR",Msg::kVerbose) << "vs " << vs << " vs1 " << vs1
00442 << " avgbeg " << avgbeg*T_EFF/100.
00443 << " Theta " << Theta << endl;
00444 if(ss[0][0]==0. || (TMath::Abs(Theta)>TMath::Abs(Theta1) &&
00445 TMath::Abs(avgdev*T_EFF/100.-vs) > 2 *
00446 TMath::Abs(avgbeg*T_EFF/100.-vs1) &&
00447 (stpden<0.5 || (stpsp>0.5 &&
00448 (posdif>0. && posdif<1.)
00449 || Ns<=2)))){
00450 double swap = Theta;
00451 Theta = Theta1;
00452 Theta1 = swap;
00453 }
00454 Theta_rad = (Theta/180)*TMath::Pi();
00455 Theta_rad1 = (Theta1/180)*TMath::Pi();
00456 vs = maxtpos+tan(Theta_rad)*(begz-maxz);
00457 MSG("SubShowerSR",Msg::kVerbose)
00458 << "vs " << vs << " avgbeg " << avgbeg*T_EFF/100.
00459 << " Theta " << Theta << endl;
00460 if(vs>(uppeststp+(uppeststp-loweststp)*1.5/2.)*T_EFF/100.||
00461 vs<(loweststp-(uppeststp-loweststp)*1.5/2.)*T_EFF/100.
00462 ||vs<-4.||vs>4.) {
00463 double bigang = Theta;
00464 if(TMath::Abs(tan(Theta_rad))<TMath::Abs(tan(Theta_rad1))) bigang = Theta1;
00465 Theta = 0.;
00466 Theta1 = TMath::Abs(tan(bigang))/tan(bigang)*89.99;
00467 Theta_rad = (Theta/180)*TMath::Pi();
00468 Theta_rad1 = (Theta1/180)*TMath::Pi();
00469 }
00470 MSG("SubShowerSR",Msg::kVerbose)
00471 << "vs " << vs << " avgbeg " << avgbeg*T_EFF/100.
00472 << " Theta " << Theta << endl;
00473
00474 std::vector<Double_t> shwstpc(nstp,0);
00475 std::vector<Double_t> shwstpc1(nstp,0);
00476 std::vector<Double_t> shwstpc2(nstp,0);
00477 std::vector<Double_t> wshwstpc(nstp,0);
00478 for(Int_t i=0;i<nstp;i++) shwstpc[i] = -1.;
00479 for(Int_t i=0;i<nstp;i++) shwstpc1[i] = -1.;
00480 for(Int_t i=0;i<nstp;i++) shwstpc2[i] = -1.;
00481 for(Int_t i=0;i<nstp;i++) wshwstpc[i] = -1.;
00482
00483 avgdev = 0.;
00484 wavgdev = 0.;
00485
00486 Float_t offstp = 0;
00487 Float_t offstp2 = 0;
00488
00489 out2Rmph = 0;
00490 for(Int_t i=0;i<nstp;i++){
00491 shwstpc[i]=TMath::Abs(tpos[i]
00492 -(maxtpos+tan(Theta_rad)
00493 *(z[i]-maxz)))*TMath::Cos(Theta_rad);
00494 shwstpc2[i]=TMath::Abs(tpos[i]
00495 -(maxtpos+tan(Theta_rad1)
00496 *(z[i]-maxz)))*TMath::Cos(Theta_rad1);
00497 if (shwstpc[i]>2*RM_EFF/100.) out2Rmph+=ph[i];
00498 if (shwstpc[i]>1) offstp++;
00499 if (shwstpc2[i]>1) offstp2++;
00500 avgdev1 += shwstpc[i];
00501 avgdev2 += shwstpc2[i];
00502 wavgdev1 += shwstpc[i]*ph[i];
00503 wavgdev2 += shwstpc2[i]*ph[i];
00504 }
00505
00506 out2Rmph/=totw;
00507 avgdev1/=nstp;
00508 avgdev2/=nstp;
00509 wavgdev1/=totw;
00510 wavgdev2/=totw;
00511 avgdev = avgdev1;
00512 wavgdev = wavgdev1;
00513 offstp/=nstp;
00514 offstp2/=nstp;
00515 }
00516 double vs = maxtpos+tan(Theta_rad)*(begz-maxz);
00517
00518
00519
00520 if(vs<loweststp*T_EFF/100.-0.25 ||
00521 vs>uppeststp*T_EFF/100.+0.25) {
00522 vs = avgbeg*T_EFF/100.;
00523 Theta_rad = 0.;
00524 }
00525 if(cch.GetPlaneView()==PlaneView::kX ||
00526 cch.GetPlaneView()==PlaneView::kU) {
00527 cch.SetVtxU(vs);
00528 }
00529 if(cch.GetPlaneView()==PlaneView::kY ||
00530 cch.GetPlaneView()==PlaneView::kV) {
00531 cch.SetVtxV(vs);
00532 }
00533 cch.SetVtxZ(begz);
00534 cch.SetEndZ(endz);
00535 cch.SetVtxT(mint);
00536 cch.SetVtxPlane(begpl);
00537 cch.SetEndPlane(endpl);
00538 cch.SetEnergy(Eguess/2.);
00539 cch.SetSlope(tan(Theta_rad));
00540 cch.SetAvgDev(wavgdev);
00541
00542 LongPl.clear();
00543 LongPh.clear();
00544 LongFlu.clear();
00545 LongPhpe.clear();
00546 TranStp.clear();
00547 TranPh.clear();
00548 TranFlu.clear();
00549
00550 EMFluctuation EMFlu(Eguess);
00551 nn0l = 0;
00552 nn0t = 0;
00553
00554 if(cos(Theta_rad)!=0.){
00555 Bool_t xs = false;
00556 Int_t dpl = 2;
00557 if(det==1&&(begpl-249)*(endpl-249)<0) xs = true;
00558 for(Int_t j=begpl;j<=endpl;j+=dpl){
00559 dpl = 2;
00560 if(xs&&(begpl-249)*(j-249)<0) {
00561 dpl = 3;
00562 xs = false;
00563 }
00564 Double_t totphPl = 0.;
00565 Double_t totphPl_pe = 0.;
00566 Double_t plz = -1.;
00567 Double_t nstpPl = 0;
00568 for(Int_t i=0;i<nstp;i++){
00569 if(plane[i]==j && ph[i]>0) {
00570 totphPl+=ph[i];
00571 totphPl_pe+=ph_pe[i];
00572 plz = z[i];
00573 Double_t zt_vtx_corrected[2];
00574 zt_vtx_corrected[0] = z[i]-begz;
00575 zt_vtx_corrected[1] = tpos[i]-vs-(z[i]-begz)*tan(Theta_rad);
00576 nstpPl++;
00577 }
00578 }
00579 LongPl.push_back(plz-begz);
00580 LongPh.push_back(totphPl);
00581 LongPhpe.push_back(totphPl_pe);
00582 if(totphPl>0.) nn0l++;
00583 }
00584
00585 double z_adj = 0.;
00586 std::list<Double_t>::iterator PlIt = LongPl.begin();
00587 std::list<Double_t>::iterator PhIt = LongPh.begin();
00588 std::list<Double_t>::iterator PhpeIt = LongPhpe.begin();
00589 double t0 = (*PhIt);
00590 PhIt++;
00591 double t1 = (*PhIt);
00592 PhIt = LongPh.begin();
00593 if(vtxdiff>1.25) z_adj = Int_t((1-vtxdiff)*tmaxpl);
00594 if(t0>t1/2.&&t0<t1) z_adj = 1;
00595 else if(t0>=t1) {
00596 z_adj = 2;
00597 LongPl.push_front(0);
00598 LongPh.push_front(0.);
00599 LongPhpe.push_front(0.);
00600 }
00601
00602 Int_t inEM = 0;
00603 PlIt = LongPl.begin();
00604 while(PlIt!=LongPl.end()){
00605 if(z_adj!=0) {
00606 (*PlIt) += z_adj*D_EFF/100.;
00607 begz += z_adj*D_EFF/100.;
00608 endz += z_adj*D_EFF/100.;
00609 }
00610 double fluem = 0.;
00611 if ((*PlIt)>=0.) {
00612 fluem = EMFlu.CalcLongFluc((*PlIt)/cos(Theta_rad));
00613 inEM++;
00614 }
00615 double error = 0.;
00616 if ((*PhpeIt)<=0.001) error = fluem;
00617 else {
00618 MSG("SubShowerSR",Msg::kVerbose)
00619 << "fluem " << fluem << " *PhpeIt " << (*PhpeIt) << endl;
00620 error = TMath::Sqrt(TMath::Power(fluem,2) +
00621 TMath::Power(1./TMath::Sqrt(*PhpeIt),2));
00622 }
00623 if(inEM==1) error = TMath::Sqrt(TMath::Power(error,2)+1.);
00624 MSG("SubShowerSR",Msg::kVerbose)
00625 << (*PlIt) << " long ph pe " << (*PhpeIt)
00626 << " error " << error << " " << fluem << endl;
00627 LongFlu.push_back(error);
00628 PlIt++;
00629 PhpeIt++;
00630 }
00631 MSG("SubShowerSR",Msg::kVerbose)<<"z_adj "<<z_adj<<endl;
00632
00633 PlIt = LongPl.end();
00634 PlIt--;
00635 Double_t miss_z = (*PlIt);
00636 PlIt = LongPl.begin();
00637 t0 = *(LongPl.begin());
00638 while(miss_z-t0<2*tmax*X0_EFF/100.){
00639 miss_z+=2*D_EFF/100.;
00640 LongPl.push_back(miss_z);
00641 LongPh.push_back(0.);
00642 Double_t miss_flu = EMFlu.CalcLongFluc(miss_z/cos(Theta_rad));
00643 MSG("SubShowerSR",Msg::kVerbose)<<"miss_flu "<<miss_flu<<endl;
00644 LongFlu.push_back(miss_flu);
00645 }
00646 MSG("SubShowerSR",Msg::kVerbose)<< "0 loweststp " << loweststp
00647 << " uppeststp " << uppeststp << endl;
00648 Double_t Thetabd = atan(ST_WID/(Double_t(LongPl.size()) *
00649 D_EFF*2.))*180./TMath::Pi();
00650 Double_t dThetai = 1./180.*TMath::Pi();
00651 if(TMath::Abs(Thetabd)<5.) dThetai = TMath::Abs(Thetabd)/180.*TMath::Pi()/5.;
00652 Double_t maxph0 = 0.;
00653 Double_t maxvs0 = vs;
00654 Double_t maxtheta0 = Theta_rad;
00655 for(Int_t k=-6;k<=6;k++){
00656 for(Int_t j=-5;j<=5;j++){
00657 Double_t totph0 = 0.;
00658 for(Int_t i=0;i<nstp;i++){
00659 Double_t zt_vtx_corrected[2];
00660 zt_vtx_corrected[0] = z[i]-begz;
00661 zt_vtx_corrected[1] = tpos[i] - (vs+ST_WID/100.*k/4.) -
00662 (z[i]-begz) * tan(Theta_rad+j*dThetai);
00663 if(TMath::Abs(zt_vtx_corrected[1]-0.)<=T_EFF/2./100.&&ph[i]>0) {
00664 totph0+=ph[i];
00665 }
00666 }
00667 if(totph0>maxph0){
00668 maxph0 = totph0;
00669 maxvs0 = vs+ST_WID/100.*k/4.;
00670 maxtheta0 = Theta_rad+j*dThetai;
00671 MSG("SubShowerSR",Msg::kVerbose)
00672 <<"Theta_rad "<<Theta_rad<<" vs "<<vs
00673 <<" maxtheta0 "<<maxtheta0<<" maxvs0 "<<maxvs0
00674 <<" maxph0 "<<maxph0<<" j "<<j<<" k "<<k
00675 <<" dThetai "<<dThetai<<endl;
00676 }
00677 }
00678 }
00679 MSG("SubShowerSR",Msg::kDebug) <<"Theta_rad "<<Theta_rad
00680 <<" vs "<<vs
00681 <<" maxtheta0 "<<maxtheta0
00682 <<" maxvs0 "<<maxvs0<<endl;
00683 Theta_rad = maxtheta0;
00684 vs = maxvs0;
00685 for(Int_t j=int((loweststp-uppeststp)/2.)-1;
00686 j<=int((uppeststp-loweststp)/2.)+1;j++){
00687
00688 Double_t totphStp = 0.;
00689 Double_t totphStp_pe = 0.;
00690 Double_t stptpos = j*T_EFF/100.;
00691
00692 Double_t nstpTr = 0.;
00693 for(Int_t i=0;i<nstp;i++){
00694 Double_t zt_vtx_corrected[2];
00695 zt_vtx_corrected[0] = z[i]-begz;
00696 zt_vtx_corrected[1] = tpos[i]-vs-(z[i]-begz)*tan(Theta_rad);
00697
00698 if(TMath::Abs(zt_vtx_corrected[1]-stptpos)<=T_EFF/2./100.&&ph[i]>0) {
00699 totphStp+=ph[i];
00700 totphStp_pe+=ph_pe[i];
00701 nstpTr++;
00702 }
00703 }
00704 TranStp.push_back(stptpos);
00705 TranPh.push_back(totphStp);
00706 if(totphStp>0.) nn0t++;
00707 double fluem = EMFlu.CalcTranFluc(stptpos*cos(Theta_rad));
00708 double error = 0.;
00709 if (totphStp_pe<=0.001) error = fluem;
00710 else {
00711 MSG("SubShowerSR",Msg::kVerbose)<<"fluem "<<fluem
00712 <<" totphStp_pe "<<totphStp_pe<<endl;
00713 error = TMath::Sqrt(TMath::Power(fluem,2) +
00714 TMath::Power(1./TMath::Sqrt(totphStp_pe),2));
00715 MSG("SubShowerSR",Msg::kVerbose)<<"tran ph "<<totphStp_pe
00716 <<" error "<<error<<" "<<fluem<<endl;
00717 }
00718 TranFlu.push_back(error);
00719 }
00720
00721 std::list<Double_t>::iterator StpIt = TranStp.begin();
00722 Double_t miss_tn = *(TranStp.begin());
00723 StpIt = TranStp.end();
00724 StpIt--;
00725 Double_t miss_tp = *StpIt;
00726
00727 while(miss_tp-miss_tn<4*RM_EFF/100.){
00728 miss_tn-=T_EFF/100.;
00729 miss_tp+=T_EFF/100.;
00730 TranStp.push_front(miss_tn);
00731 TranStp.push_back(miss_tp);
00732 TranPh.push_front(0.);
00733 TranPh.push_back(0.);
00734 Double_t miss_flun = EMFlu.CalcLongFluc(miss_tn*cos(Theta_rad));
00735 Double_t miss_flup = EMFlu.CalcLongFluc(miss_tp*cos(Theta_rad));
00736 MSG("SubShowerSR",Msg::kVerbose)<<"miss_flun "<<miss_flun
00737 <<" miss_flup "<<miss_flup<<endl;
00738 TranFlu.push_front(miss_flun);
00739 TranFlu.push_back(miss_flup);
00740 loweststp--;
00741 uppeststp++;
00742 }
00743 }
00744 }
00745
00746
00747
00748 void AlgSubShowerSR::SubShowerID(CandSubShowerSRHandle &cch){
00749
00750 ClusterType::ClusterType_t id = ClusterType::kUnknown;
00751
00752 if((xtkfra>fXtkFraCut || avgph<2*fPECut) && maxStpPHinSS<fMaxXTalkStpPH )
00753 id = ClusterType::kXTalk;
00754 else if(trkfra>fTrkFraCut || avgstpperpln<fTrkLikeStpCut)
00755 id = ClusterType::kTrkLike;
00756 else if(cch.GetEnergy()<fECut)
00757 id = ClusterType::kHadLike;
00758 else if(maxdiff>fMaxDiffCut)
00759 id = ClusterType::kHadLike;
00760 else if(vtxdiff>fVtxDiffCut)
00761 id = ClusterType::kHadLike;
00762 else if(out2Rmph>fOut2RmPHCut)
00763 id = ClusterType::kHadLike;
00764 else id = ClusterType::kEMLike;
00765
00766 cch.SetClusterID(id);
00767 MSG("SubShowerSR",Msg::kVerbose)<<"id "<<id<<" E "<<cch.GetEnergy()<<endl;
00768
00769 bothfitprob = 0.;
00770 Bool_t makeplots = false;
00771 MSG("SubShowerSR",Msg::kDebug)
00772 << "nn0l " << nn0l << " LongPl.size() " << LongPl.size()
00773 << " nn0t " << nn0t << " TranStp.size() " << TranStp.size() << endl;
00774
00775
00776 if(fDoEMFit &&
00777 (id==ClusterType::kEMLike ||
00778 (id==ClusterType::kHadLike && cch.GetEnergy()>=fECut)) &&
00779 LongPl.size()-nn0l <= 3 && TranStp.size()-nn0t <= 3 ){
00780
00781 MSG("SubShowerSR",Msg::kDebug)<<"in fitting "<<endl;
00782
00783 std::list<Double_t>::iterator PlIt = LongPl.begin();
00784 std::list<Double_t>::iterator PhIt = LongPh.begin();
00785 std::list<Double_t>::iterator PluIt = LongFlu.begin();
00786 std::list<Double_t>::iterator StpIt = TranStp.begin();
00787 std::list<Double_t>::iterator StphIt = TranPh.begin();
00788 std::list<Double_t>::iterator StpluIt = TranFlu.begin();
00789
00790 MSG("SubShowerSR",Msg::kVerbose)
00791 << "funct boundary " << ((loweststp-uppeststp)/2.-1)*T_EFF-T_EFF/2.
00792 << " " << ((uppeststp-loweststp)/2.+1)*T_EFF+T_EFF/2. << endl;
00793
00794 Double_t thetarange = Theta_rad*0.05;
00795 if(thetarange<10./180.*TMath::Pi()) thetarange = 10./180.*TMath::Pi();
00796 Double_t erange = Eguess*0.05;
00797 if(erange<0.5) erange = 0.5;
00798
00799 clusterlong->SetRange(0.-D_EFF/2.,(endz-begz+(*PlIt))*100+D_EFF);
00800 clusterlong->SetParameter(0,Eguess);
00801 clusterlong->SetParLimits(0,Eguess-erange,Eguess+erange);
00802 clusterlong->SetParameter(1,Theta_rad);
00803 clusterlong->SetParLimits(1,Theta_rad-thetarange,Theta_rad+thetarange);
00804
00805 clustertran->SetRange(((loweststp-uppeststp)/2.-1)*T_EFF,
00806 ((uppeststp-loweststp)/2.+1)*T_EFF);
00807 clustertran->SetParameter(0,Eguess);
00808 clustertran->SetParLimits(0,Eguess-erange,Eguess+erange);
00809 clustertran->SetParameter(1,Theta_rad);
00810 clustertran->SetParLimits(1,Theta_rad-thetarange,Theta_rad+thetarange);
00811
00812 LongProf->Reset();
00813 LongProf->SetBins(int((endz-begz+(*PlIt))*100+2*D_EFF),
00814 0-D_EFF,(endz-begz+(*PlIt))*100+D_EFF);
00815 LongHist->Reset();
00816 LongHist->SetBins(int((endz-begz+(*PlIt))*100+2*D_EFF),
00817 0-D_EFF,(endz-begz+(*PlIt))*100+D_EFF);
00818 TranProf->Reset();
00819 TranProf->SetBins(int((uppeststp-loweststp+2)*T_EFF),
00820 ((loweststp-uppeststp)/2.-1)*T_EFF,
00821 ((uppeststp-loweststp)/2.+1)*T_EFF);
00822 TranHist->Reset();
00823 TranHist->SetBins(int((uppeststp-loweststp+2)*T_EFF),
00824 ((loweststp-uppeststp)/2.-1)*T_EFF,
00825 ((uppeststp-loweststp)/2.+1)*T_EFF);
00826 LongProfb->Reset();
00827 LongProfb->SetBins(int((endz-begz+(*PlIt))*100+2*D_EFF),
00828 0-D_EFF,(endz-begz+(*PlIt))*100+D_EFF);
00829 LongHistb->Reset();
00830 LongHistb->SetBins(int((endz-begz+(*PlIt))*100+2*D_EFF),
00831 0-D_EFF,(endz-begz+(*PlIt))*100+D_EFF);
00832 TranProfb->Reset();
00833 TranProfb->SetBins(int((uppeststp-loweststp+2)*T_EFF),
00834 ((loweststp-uppeststp)/2.-1)*T_EFF,
00835 ((uppeststp-loweststp)/2.+1)*T_EFF);
00836 TranHistb->Reset();
00837 TranHistb->SetBins(int((uppeststp-loweststp+2)*T_EFF),
00838 ((loweststp-uppeststp)/2.-1)*T_EFF,
00839 ((uppeststp-loweststp)/2.+1)*T_EFF);
00840
00841 MSG("SubShowerSR",Msg::kVerbose) << "loweststp " << loweststp
00842 << " uppeststp " << uppeststp << endl;
00843
00844 Double_t Minchi2Both = 100000;
00845 Double_t BestEBoth = Eguess;
00846 Double_t BestThetaBoth = Theta_rad;
00847 double bestnormlb = 1;
00848 double bestnormtb = 1;
00849 Bool_t fittedboth = false;
00850 Double_t dE = 0.1;
00851 Double_t Eini = Eguess-erange;
00852 Double_t Ein = Eini;
00853
00854 MSG("SubShowerSR",Msg::kDebug)
00855 << "done initial stuff " <<Eguess << " " << Ein
00856 << " " << Theta_rad << " " << thetarange << endl;
00857
00858 while(Ein>=Eini&&Ein<=Eguess+erange){
00859 MSG("SubShowerSR",Msg::kVerbose)<<"Ein "<<Ein<<endl;
00860 clusterlong->SetParameter(0,Ein);
00861 Double_t dTheta = 2.5/180.*TMath::Pi();
00862 Double_t Thetaini = Theta_rad-thetarange;
00863 Double_t Thetain = Thetaini;
00864 while(Thetain>=Thetaini&&Thetain<=Theta_rad+thetarange){
00865 MSG("SubShowerSR",Msg::kVerbose)<<"Thetain "<<Thetain<<endl;
00866 if(cos(Thetain)!=0.){
00867 clusterlong->SetParameter(1,Thetain);
00868 Double_t norml = 0.0;
00869 std::vector<Double_t> funcl(LongPl.size(),0);
00870 std::vector<Double_t> errl(LongPl.size(),0);
00871 Int_t nfpl = 0;
00872 UInt_t j = 0;
00873 Int_t nlong = 0;
00874 PhIt = LongPh.begin();
00875 while(PhIt!=LongPh.end()){
00876 if((*PhIt)>0.) nlong++;
00877 PhIt++;
00878 }
00879 PhIt = LongPh.begin();
00880 PlIt = LongPl.begin();
00881 while(PlIt!=LongPl.end()){
00882 funcl[j] = 0.;
00883 if((*PlIt)>=0.) funcl[j] = clusterlong->Eval((*PlIt)*100.);
00884 MSG("SubShowerSR",Msg::kVerbose)<<"funcl[j] "<<funcl[j]<<endl;
00885 norml += funcl[j];
00886 if (funcl[j]>0.||(*PlIt)<0.) nfpl++;
00887 PlIt++;
00888 j++;
00889 }
00890 MSG("SubShowerSR",Msg::kVerbose)<<"nfpl "<<nfpl<<endl;
00891 Double_t chi2Long = 0.;
00892 PlIt = LongPl.begin();
00893 PhIt = LongPh.begin();
00894 PluIt = LongFlu.begin();
00895 j = 0;
00896 MSG("SubShowerSR",Msg::kVerbose)
00897 << "nfpl " << nfpl << " norml " << norml << " nlong " << nlong
00898 << " Ein " << Ein << " Thetain " << Thetain << endl;
00899
00900 if(double(nfpl)/double(LongPl.size())>0.5 &&
00901 norml>1e-3 && nlong>1){
00902
00903 while(PlIt!=LongPl.end()){
00904 Double_t data = (*PhIt)/totw*norml;
00905 Double_t err = funcl[j]*(*PluIt);
00906 if ((*PluIt)==0.0) err = 0.0001;
00907 UInt_t sub = j;
00908 while(funcl[sub]<0.05*norml&&j==0) {
00909 sub++;
00910 err = (funcl[sub]+funcl[j])/2.*(*PluIt);
00911 MSG("SubShowerSR",Msg::kVerbose)<<"err "<<err<<" sub "<<endl;
00912 if(err>0.0001) break;
00913 else if(sub==LongPl.size()-1||(sub-j)>=3){
00914 err = 0.0001;
00915 break;
00916 }
00917 }
00918 sub = j;
00919 while(funcl[sub]<0.05*norml&&j>0) {
00920 sub--;
00921 err = (funcl[sub]+funcl[j])/2.*(*PluIt);
00922 MSG("SubShowerSR",Msg::kVerbose)<<"err "<<err<<" sub "<<endl;
00923 if(err>0.0001) break;
00924 else if(sub==0||(j-sub)>=3){
00925 err = 0.0001;
00926 break;
00927 }
00928 }
00929
00930 MSG("SubShowerSR",Msg::kDebug)
00931 << " LongPl.size() " << LongPl.size()
00932 << " Eguess " <<Eguess << " norml " << norml
00933 << " Theta " << Theta_rad << " j " << j
00934 << " LongPl.at(j) " << (*PlIt)
00935 << " err " << funcl[j] << " " << funcl[j-1]
00936 << " " << funcl[j+1] << " " << (*PluIt)
00937 << " err " << err << " data - funcl " << data-funcl[j] << endl;
00938
00939 if(err<0.0001) err = 0.0001;
00940 MSG("SubShowerSR",Msg::kVerbose)
00941 <<"err "<<err<<" funcl[j] "<<funcl[j]<<endl;
00942
00943 if(j==0)
00944 errl[j] = TMath::Sqrt(TMath::Power(err,2) +
00945 TMath::Power((funcl[j+1] -
00946 funcl[j])/2.,2));
00947 else if(j==LongPl.size()-1)
00948 errl[j] = TMath::Sqrt(TMath::Power(err,2) +
00949 TMath::Power((funcl[j] -
00950 funcl[j-1])/2.,2));
00951 else
00952 errl[j] = TMath::Sqrt(TMath::Power(err,2) +
00953 TMath::Power((funcl[j+1] -
00954 funcl[j-1])/4.,2));
00955
00956 MSG("SubShowerSR",Msg::kVerbose)
00957 << "data-funcl[j] " << data-funcl[j]
00958 << " errl[j] " << errl[j] << endl;
00959
00960 chi2Long += ( TMath::Power(data-funcl[j],2) /
00961 TMath::Power(errl[j],2) );
00962
00963 MSG("SubShowerSR",Msg::kVerbose)
00964 << "chi2long " << (*PlIt) << " " << data
00965 << " " << funcl[j] << " " << (*PluIt) << " " << err
00966 << " " << TMath::Power(data-funcl[j],2)/TMath::Power(err,2)
00967 << " " << chi2Long << endl;
00968
00969 PlIt++;
00970 PhIt++;
00971 PluIt++;
00972 j++;
00973 }
00974
00975 MSG("SubShowerSR",Msg::kVerbose)
00976 << "long fit " << chi2Long/(LongPl.size()-1)
00977 << "norml " << LongPl.size() << " " << norml << endl;
00978 }
00979 else chi2Long = 100000;
00980
00981 Double_t normt = 0.0;
00982 std::vector<Double_t> funct(TranStp.size(),0);
00983 std::vector<Double_t> errt(TranStp.size(),0);
00984 Int_t nfpt = 0;
00985 j = 0;
00986 Int_t ntran = 0;
00987 StphIt = TranPh.begin();
00988 while(StphIt!=TranPh.end()){
00989 if(*StphIt>0.) ntran++;
00990 StphIt++;
00991 }
00992 StphIt = TranPh.begin();
00993 StpIt = TranStp.begin();
00994 while(StpIt!=TranStp.end()){
00995 funct[j] = clustertran->Eval(*StpIt*100.);
00996 normt += funct[j];
00997 if(funct[j]>0.) nfpt++;
00998 StpIt++;
00999 j++;
01000 }
01001
01002 Double_t chi2Tran = 0.;
01003 MSG("SubShowerSR",Msg::kVerbose)<<"nfpt "<<nfpt<<" "<<normt<<endl;
01004 if(double(nfpt)/double(TranStp.size())>0.5 &&
01005 normt>1e-3 && ntran>1){
01006
01007 StpIt = TranStp.begin();
01008 StphIt = TranPh.begin();
01009 StpluIt = TranFlu.begin();
01010 j = 0;
01011
01012 while(StpIt!=TranStp.end()){
01013 Double_t data = *StphIt/totw*normt;
01014 Double_t err = funct[j]*(*StpluIt);
01015 if (*StpluIt==0.0) err = 0.0001;
01016 UInt_t sub = j;
01017 while(funct[sub]<0.05*normt&&*StpIt<=0) {
01018 sub++;
01019 err = (funct[sub]+funct[j])/2.*(*StpluIt);
01020 if(err>0.0001) break;
01021 else if(sub==TranStp.size()-1||(sub-j)>=3){
01022 err = 0.0001;
01023 break;
01024 }
01025 }
01026 sub = j;
01027 while(funct[sub]<0.05*normt && *StpIt>0) {
01028 sub--;
01029 err = (funct[sub]+funct[j])/2.*(*StpluIt);
01030 if(err>0.0001) break;
01031 else if(sub==0||(j-sub)>=3){
01032 err = 0.0001;
01033 break;
01034 }
01035 }
01036 MSG("SubShowerSR",Msg::kDebug)
01037 << "funct " << j << " " << *StpIt << " " << *StpluIt
01038 << " err " << " " << funct[j] << " " << funct[j-1]
01039 << " " << funct[j+1] << " err " << err
01040 << " data - funcl " << data-funcl[j] << endl;
01041
01042 if(err<0.0001) err = 0.0001;
01043 if(j==0)
01044 errt[j] = TMath::Sqrt(TMath::Power(err,2) +
01045 TMath::Power((funct[j+1] -
01046 funct[j])/2.,2));
01047 else if(j==TranStp.size()-1)
01048 errt[j] = TMath::Sqrt(TMath::Power(err,2) +
01049 TMath::Power((funct[j] -
01050 funct[j-1])/2.,2));
01051 else
01052 errt[j] = TMath::Sqrt(TMath::Power(err,2) +
01053 TMath::Power((funct[j+1] -
01054 funct[j-1])/4.,2));
01055
01056 MSG("SubShowerSR",Msg::kVerbose) << "data-funcl[j] " <<
01057 data-funcl[j] << " errt[j] " << errt[j] << endl;
01058
01059 chi2Tran += ( TMath::Power(data-funct[j],2) /
01060 (TMath::Power(errt[j],2)) );
01061
01062 MSG("SubShowerSR",Msg::kVerbose) << "chi2tran " <<
01063 TMath::Power(data-funcl[j],2)/TMath::Power(err,2)
01064 << " " << chi2Tran << endl;
01065
01066 StpIt++;
01067 StphIt++;
01068 StpluIt++;
01069 j++;
01070 }
01071
01072 MSG("SubShowerSR",Msg::kVerbose)
01073 << "tran fit " << chi2Tran/(TranStp.size()-1)
01074 << "normt " << TranStp.size() << " " << normt << endl;
01075 }
01076 else chi2Tran = 100000;
01077
01078 bothfitprob = 0.;
01079 Double_t chi2Both = 0.;
01080
01081 MSG("SubShowerSR",Msg::kVerbose) << "chi2Tran " <<chi2Tran
01082 << " chi2Long " << chi2Long << endl;
01083
01084 chi2Both = chi2Tran+chi2Long;
01085 if((double(nfpl)/double(LongPl.size())>0.5&&norml>1e-3&&nlong>1)&&
01086 (double(nfpt)/double(TranStp.size())>0.5&&normt>1e-3&&ntran>1)&&
01087 chi2Both<=Minchi2Both){
01088 Minchi2Both = chi2Both;
01089 BestEBoth = Ein;
01090 BestThetaBoth = Thetain;
01091 bestnormtb = normt;
01092 bestnormlb = norml;
01093 fittedboth = true;
01094 if(makeplots){
01095 PlIt = LongPl.begin();
01096 PhIt = LongPh.begin();
01097 PluIt = LongFlu.begin();
01098 j = 0;
01099 while(PlIt!=LongPl.end()){
01100 double data = (*PhIt)/totw*norml;
01101 Int_t bin = LongProfb->FindBin((*PlIt)*100.);
01102 LongProfb->SetBinContent(bin,funcl[j]);
01103 LongProfb->SetBinError(bin,errl[j]);
01104 LongHistb->SetBinContent(bin,data);
01105 LongHistb->SetBinError(bin,0.0);
01106 PlIt++;
01107 PhIt++;
01108 PluIt++;
01109 j++;
01110 }
01111 StpIt = TranStp.begin();
01112 StphIt = TranPh.begin();
01113 StpluIt = TranFlu.begin();
01114 j = 0;
01115
01116 while(StpIt!=TranStp.end()){
01117 double data = *StphIt/totw*normt;
01118 Int_t bin = TranProfb->FindBin(*StpIt*100.);
01119 TranProfb->SetBinContent(bin,funct[j]);
01120 TranProfb->SetBinError(bin,errt[j]);
01121 TranHistb->SetBinContent(bin,data);
01122 TranHistb->SetBinError(bin,0.);
01123 StpIt++;
01124 StphIt++;
01125 StpluIt++;
01126 j++;
01127 }
01128 }
01129 }
01130 }
01131 Thetain+=dTheta;
01132 }
01133 Ein+=dE;
01134 }
01135
01136 if(fittedboth){
01137 if(LongPl.size()>1&&TranStp.size()>1){
01138 Int_t ndfLong = LongPl.size();
01139 Int_t ndfTran = TranStp.size();
01140 Int_t ndfBoth = ndfLong + ndfTran;
01141 chi2func->SetParameter(0,ndfBoth);
01142 if(fittedboth && Minchi2Both/ndfBoth<100.)
01143 bothfitprob = 1.0 - chi2func->Integral(0.0,Minchi2Both);
01144 else bothfitprob = 0.0;
01145
01146 MSG("SubShowerSR",Msg::kDebug)
01147 << "Best Fit: Both- " << BestEBoth << " GeV "
01148 << BestThetaBoth << " rad with chi2/ndof "
01149 << Minchi2Both/ndfBoth << " ndf " << ndfBoth
01150 << " of probability " << bothfitprob
01151 << " comparing to " << Eguess << " GeV "
01152 << Theta_rad << " rad" << endl;
01153
01154 }
01155
01156 if(makeplots){
01157 gStyle->SetOptStat(0);
01158 char name[256];
01159 TCanvas *can = new TCanvas();
01160 can->cd();
01161 can->Divide(2,1);
01162 can->cd(1);
01163 LongHistb->SetMaximum(bestnormlb);
01164 LongHistb->Draw("EP");
01165 LongProfb->SetMarkerColor(4);
01166 LongProfb->Draw("EPsame");
01167 can->cd(2);
01168 TranHistb->SetMaximum(bestnormtb);
01169 TranHistb->Draw("EP");
01170 TranProfb->SetMarkerColor(4);
01171 TranProfb->Draw("EPsame");
01172 if(bestnormlb+bestnormtb>0){
01173 sprintf(name,"Fit%f_%fb.eps",Eguess,Theta_rad);
01174 can->Print(name);
01175 }
01176 delete can;
01177 }
01178 }
01179 }
01180 else {
01181 bothfitprob = 0.;
01182 }
01183 MSG("SubShowerSR",Msg::kVerbose) << "long " << LongPl.size()
01184 << " tran " <<TranStp.size()
01185 << " id " << id
01186 << " eng " << cch.GetEnergy()
01187 << " nstrip " << cch.GetNStrip()
01188 << " Prob " << bothfitprob << endl;
01189 cch.SetProbEM(Float_t(bothfitprob));
01190 }
01191
01192
01193 void AlgSubShowerSR::Trace(const char * ) const
01194 {
01195 }
01196
01197
01198 Double_t LongShwProfSamAng(double *x,double *par){
01199 Double_t plz = x[0];
01200 Double_t E = par[0]*1000;
01201 if(E<=0.) return 0;
01202 Double_t tmax = log(E/EC_EFF)-0.858;
01203 Double_t theta = par[1];
01204 Double_t theta_rt = TMath::ATan(TMath::Tan(theta)*(T_EFF*X0_EFF)/(D_EFF*RM_EFF));
01205 Double_t delta_z = 0.1*D_EFF;
01206 if(plz/X0_EFF/cos(theta_rt)>2*tmax) delta_z = 0.2*D_EFF;
01207 Double_t z = (plz-D_EFF/2.)/X0_EFF;
01208 Double_t dE = 0.;
01209 for(Int_t i=0;i<=int(D_EFF/delta_z);i++){
01210 Double_t t[1] = {0.};
01211 if (cos(theta_rt)!=0.) t[0] = z/cos(theta_rt);
01212 if(t[0]>0.01*delta_z/X0_EFF)
01213 dE += LongShwProfSamn(t,par)*cos(theta)*delta_z;
01214 MSG("SubShowerSR",Msg::kVerbose)<<"z "<<z<<" t[0] "<<t[0]<<" dE "<<dE<<endl;
01215 z+=delta_z/X0_EFF;
01216 }
01217 return dE;
01218 }
01219
01220
01221 Double_t TranShwProfSamAng(double *x,double *par){
01222
01223 Double_t stptpos = x[0];
01224 Double_t E = par[0]*1000;
01225 if(E<=0) return 0;
01226 Double_t tmax = log(E/EC_EFF)-0.858;
01227 MSG("SubShowerSR",Msg::kVerbose)<<"E "<<E<<" tmax "<<tmax<<endl;
01228 Double_t theta = par[1];
01229 Double_t theta_rt = TMath::ATan(TMath::Tan(theta)*(T_EFF*X0_EFF)/(D_EFF*RM_EFF));
01230 Double_t delta_z = 0.25*D_EFF;
01231 Double_t delta_tpos = 0.1*T_EFF;
01232 if(TMath::Abs(stptpos/RM_EFF*cos(theta_rt))>1.) delta_tpos = 0.2*T_EFF;
01233 Double_t z = 0;
01234 Double_t tpos = (stptpos-T_EFF/2.)/RM_EFF;
01235 Double_t dE = 0.;
01236 Double_t cutoff = 1e-3;
01237 for(Int_t i=0;i<=int(T_EFF/delta_tpos);i++){
01238 Double_t rt[2] = {0.,0.};
01239 rt[1] = tpos*cos(theta_rt);
01240 Double_t allplE = 0.;
01241 Bool_t keepgoing = true;
01242 z = 0.;
01243 MSG("SubShowerSR",Msg::kVerbose)<<"rt[1] "<<rt[1]<<" "<<theta_rt<<endl;
01244 while(keepgoing){
01245 if (cos(theta_rt)!=0.) rt[0] = z/cos(theta_rt)+tpos*sin(theta_rt);
01246 else {
01247 rt[0] = tpos*sin(theta_rt);
01248 keepgoing = false;
01249 }
01250 MSG("SubShowerSR",Msg::kVerbose)<<"rt[0] "<<rt[0]<<endl;
01251 Double_t fracE = 0.;
01252
01253 fracE = TranShwProfSamn(rt,par);
01254 allplE += fracE*cos(theta)*delta_tpos*delta_z;
01255 if(fracE<cutoff||z<0.||z>1.5*tmax) keepgoing = false;
01256 if(z>1.25*tmax) delta_z = 0.5*D_EFF;
01257 z+=delta_z/X0_EFF;
01258 }
01259 dE += allplE;
01260 tpos+=delta_tpos/RM_EFF;
01261 }
01262 return dE;
01263 }
01264
01265
01266 Double_t LongShwProfSamn(double *x,double *par){
01267
01268 Double_t t = x[0];
01269 Double_t E = par[0]*1000.;
01270 if(E<=0.||t<=0.) return 0;
01271 Double_t Z = 24.9814;
01272 Double_t Ec = 20.7941;
01273 Double_t Fs = 0.692366;
01274 Double_t ehat = 0.8752;
01275
01276 Double_t lny = TMath::Log(E/Ec);
01277 Double_t alpha = 0.21 + (0.492+2.38/Z)*lny;
01278 Double_t T = lny - 0.858;
01279
01280 T = T - 0.59/Fs - 0.53*(1-ehat);
01281 alpha = alpha - 0.444/Fs;
01282
01283 Double_t belta = (alpha-1)/T;
01284 if(belta<=0.||alpha<=0) return 0.;
01285 Double_t f = TMath::Power((belta*t),(alpha-1))
01286 *belta*TMath::Exp(-(belta*t))/TMath::Gamma(alpha);
01287 MSG("SubShowerSR",Msg::kVerbose)<<"belta "<<belta<<" t "<<t<<" alpha "<<alpha<<" "<<TMath::Gamma(alpha)<<" "<<f<<endl;
01288 return f;
01289
01290 }
01291
01292
01293 Double_t TranShwProfSamn(double *x, double *par){
01294
01295 Double_t t = x[0];
01296 Double_t r = TMath::Abs(x[1]);
01297 Double_t E = par[0]*1000.;
01298 if(E<=0) return 0;
01299 MSG("SubShowerSR",Msg::kVerbose)<<"rte "<<t<<" "<<r<<" "<<E<<endl;
01300 Double_t Z = 24.9814;
01301 Double_t Ec = 21.6297;
01302 Double_t Fs = 0.692366;
01303 Double_t ehat = 0.8752;
01304
01305 Double_t lny = log(E/Ec);
01306 Double_t T = lny - 0.858;
01307 T = T - 0.59/Fs - 0.53*(1-ehat);
01308 if(T<=0||t/T>100.) return 0;
01309 MSG("SubShowerSR",Msg::kVerbose)<<"T "<<T<<endl;
01310 Double_t tau = t/T;
01311
01312 Double_t z1 = 0.0251+0.00319*log(E);
01313 Double_t z2 = 0.1162+(-0.000381*Z);
01314
01315 Double_t k1 = 0.659+(-0.00309*Z);
01316 Double_t k2 = 0.645;
01317 Double_t k3 = -2.59;
01318 Double_t k4 = 0.3585+0.0421*log(E);
01319
01320 Double_t p1 = 2.632+(-0.00094*Z);
01321 Double_t p2 = 0.401+0.00187*Z;
01322 Double_t p3 = 1.313+(-0.0686*log(E));
01323 if(tau-k2<-1.) return 0;
01324 MSG("SubShowerSR",Msg::kVerbose)<<"k "<<k1<<" "<<k2<<" "<<k3<<" "<<k4<<endl;
01325 Double_t Rc = z1+z2*tau;
01326 Double_t Rt = k1*(exp(k3*(tau-k2))+exp(k4*(tau-k2)));
01327 Double_t p = p1*exp((p2-tau)/p3-exp((p2-tau)/p3));
01328 MSG("SubShowerSR",Msg::kVerbose)<<"Rt0 "<<Rt<<endl;
01329 MSG("SubShowerSR",Msg::kVerbose)<<ehat<<" "<<Fs<<" "<<tau<<endl;
01330 Rc = Rc - 0.0203*(1.-ehat)+(0.0397/Fs)*exp(-tau);
01331 MSG("SubShowerSR",Msg::kVerbose)<<"Rc "<<Rc<<endl;
01332 Rt = Rt - 0.14*(1-ehat)-(0.495/Fs)*exp(-tau);
01333 MSG("SubShowerSR",Msg::kVerbose)<<"Rt "<<Rt<<endl;
01334 p = p + (1-ehat)*(0.348-(0.642/Fs)*exp(-TMath::Power(tau-1,2)));
01335 MSG("SubShowerSR",Msg::kVerbose)<<"p "<<p<<endl;
01336 Double_t fc = 0;
01337 Double_t ft = 0;
01338 if (r==0.) r = 0.01;
01339 if(r!=0) {
01340 fc = 2*r*Rc*Rc/TMath::Power(r*r+Rc*Rc,2);
01341 ft = 2*r*Rt*Rt/TMath::Power(r*r+Rt*Rt,2);
01342 }
01343
01344 Double_t f = p*fc+(1-p)*ft;
01345 return f;
01346 }
01347
01348
01349 Bool_t SortStripByEnergy(Double_t *eng0,Int_t nstp,Int_t *sort){
01350 std::vector<Double_t> eng(nstp,0);
01351 for (Int_t i=0; i<=nstp-1; i++) {
01352 sort[i] = -1;
01353 eng[i] = 0.;
01354 }
01355 for (Int_t i=0; i<=nstp-1; i++) {
01356 Int_t j = 0;
01357 while(eng0[i]<=eng[j]){
01358 j++;
01359 }
01360 MSG("SubShowerSR",Msg::kVerbose)<< "i " << i << " " << eng0[i]
01361 << " j " << j << " " << eng[j] << endl;
01362 if(eng[j]>0){
01363 Int_t k = nstp-j;
01364 while(k>0){
01365 if(eng[j+k-1]>0){
01366 eng[j+k] = eng[j+k-1];
01367 sort[j+k] = sort[j+k-1];
01368 }
01369 k--;
01370 }
01371 }
01372 eng[j] = eng0[i];
01373 sort[j] = i;
01374 }
01375 return true;
01376 }
01377
01378
01379 Double_t StpFluctuaionAng(double *x,double *par){
01380 Double_t plz = x[0];
01381 Double_t z = plz;
01382 Double_t stptpos = x[1];
01383 Double_t tpos = stptpos;
01384 Double_t E = par[0]*1000;
01385 if(E<=0) return 0;
01386 Double_t theta = par[1];
01387
01388 BinFluctuationEM fluEM = BinFluctuationEM(par[0]);
01389 Double_t rt[2] = {0.,0.};
01390 rt[1] = TMath::Abs(tpos*cos(theta)*100./T_EFF);
01391 if (cos(theta)!=0.) rt[0] = (z/cos(theta)+tpos*sin(theta))*100./D_EFF;
01392 else {
01393 rt[0] = tpos*sin(theta)*100./D_EFF;
01394 }
01395 if(rt[0]<0.) rt[0] = 0.;
01396 MSG("SubShowerSR",Msg::kVerbose)<<"t r "<<plz<<" "<<stptpos<<" "<<theta<<" "<<rt[0]<<" "<<rt[1]<<endl;
01397 return fluEM.CalcFluctuation(rt[0],rt[1]);
01398 }
01399
01400
01401 Double_t Chi2(double *x,double *par){
01402
01403 return TMath::Exp(-x[0]/2.)*TMath::Power(x[0],par[0]/2. - 1)/(TMath::Power(2.,par[0]/2.)*TMath::Gamma(par[0]/2.));
01404
01405 }
