loon: AngClusterFitAna.cxx Source File

00001 
00002 
00003 
00004 
00005 
00006 
00007 
00008 
00009 
00010 
00011 
00012 
00013 
00014 #include <iostream>
00015 #include <algorithm>
00016 #include <TCanvas.h>
00017 #include <TRotation.h>
00018 #include <TFile.h>
00019 #include <TH2F.h>
00020 #include <TH1F.h>
00021 #include <TH3F.h>
00022 #include <TPolyMarker3D.h>
00023 #include <TPolyLine3D.h>
00024 #include "StandardNtuple/NtpStRecord.h"
00025 #include "CandNtupleSR/NtpSRRecord.h"
00026 #include "CandNtupleSR/NtpSREvent.h"
00027 #include "CandNtupleSR/NtpSRTrack.h"
00028 #include "CandNtupleSR/NtpSRStrip.h"
00029 #include "MessageService/MsgService.h"
00030 #include "NueAna/AngClusterFitAna.h"
00031 #include "AnalysisNtuples/ANtpDefaultValue.h"
00032 
00033 CVSID("$Id: AngClusterFitAna.cxx,v 1.21 2007/03/01 16:38:49 rhatcher Exp $");
00034 
00035 using std::cout;
00036 using std::endl;
00037 
00038 ClassImp(AngClusterFitAna)
00039 
00040 const Float_t kDetectorPitch=0.05949; //Munits::meters
00041 
00042 const Float_t kStripWidth=0.041; //Munits::meters
00043 
00044 const Float_t kRadLengthToPlane=1.54; //X0->planes conversion factor
00045 
00046 const Float_t kMolRadDet=0.0391; //Characteristic Moliere Radius 
00047                                  //of the Minos Detectors
00048 
00049 const Float_t kMeuToGeV=1.0;  //Conversion factor from MEU to GeV
00050 
00051 AngClusterFitAna::AngClusterFitAna(AngClusterFit &acf) :
00052     fLongHitHisto(0),
00053     fLongHitHistoRad(0),
00054     fTransHitHisto(0),    
00055     fTransHitHistoRad(0),
00056     fScatterAxis(0),
00057     fLongFit(0),
00058     fTransFit(0),
00059     fAngClusterFit(acf)
00060 {
00061         
00062 }
00063 
00064 AngClusterFitAna::~AngClusterFitAna()
00065 {
00066     //Prevent memory leaks
00067     if(fLongHitHisto!=0){
00068         fLongHitHisto->Delete();
00069         fLongHitHisto=0;
00070     }
00071     
00072     if(fLongHitHistoRad !=0){
00073         fLongHitHistoRad->Delete();
00074     }
00075     if(fTransHitHisto !=0){
00076         fTransHitHisto->Delete();
00077         fTransHitHisto=0;
00078     }
00079     
00080     if(fTransHitHistoRad !=0){
00081         fTransHitHistoRad->Delete();
00082         fTransHitHistoRad=0;
00083     }
00084 
00085     if(fScatterAxis !=0){
00086         fScatterAxis->Delete();
00087         fScatterAxis=0;
00088     }
00089 
00090     if(fLongFit !=0){
00091         fLongFit->Delete();
00092         fLongFit=0;
00093     }
00094 
00095     if(fTransFit !=0){
00096         fTransFit->Delete();
00097         fTransFit=0;
00098     }
00099 
00100 }
00101 
00102 void AngClusterFitAna::Set3DHit(DeqFloat_t &x
00103                                 , DeqFloat_t &y
00104                                 , DeqFloat_t &z
00105                                 , DeqFloat_t &e){
00106     
00107     if(x.size()!=0 && y.size()!=0 && z.size()!=0 && e.size()!=0){
00108 
00109         fX=x;
00110         fY=y;
00111         fZ=z;
00112         fE=e;  
00113     
00114     }
00115 
00116 }
00117 
00118 void AngClusterFitAna::SetAngCluster(Int_t &primShow
00119                                      , DeqDeqInt_t &clusterMap
00120                                      , TVector3 &primDir)
00121 {
00122 
00123     if(clusterMap.size()!=0){
00124 
00125         fClusterMap=clusterMap;
00126         fPrimShow  =primShow;
00127         fPrimDir   =primDir;
00128 
00129     } else fPrimShow=ANtpDefVal::kInt; //Will trigger stopping of 
00130                                        //further computations.
00131     
00132 }
00133 
00134 
00135 static Double_t shwLongFunc(Double_t *x, Double_t *par)
00136 {
00137 
00138     // Longitudinal profile to be fitted to em showers 
00139     //xx=Distance from the vertex in units of X0 
00140     //par[0]=a, par[1]=b, par[2]=E0
00141    
00142     Float_t xx=x[0];
00143     
00144     Double_t lnf = TMath::Log(kRadLengthToPlane*par[2]*par[1])+(par[0]-1)
00145         *TMath::Log(xx*par[1])-par[1]*xx-TMath::LnGamma(par[0]);
00146     
00147     Double_t f = TMath::Exp(lnf);
00148 
00149     return f;
00150 }
00151 
00152 static Double_t shwTransFunc(Double_t *x, Double_t *par)
00153 {
00154     //Transverse fit profile
00155     //x[0]=Radial distance from shower axis in units of Moliere radius
00156     //par[0]=Lambda1 attenuation length for high energy core
00157     //par[1]=Lambda2 attenuation length for shower halo
00158     //par[2]=Relative weight between lambda1 and lambda2
00159     //par[3]=E0 Energy normalization
00160 
00161     Float_t xx=x[0];
00162     Double_t f=par[3]*(TMath::Exp(-TMath::Sqrt(xx/par[0]))
00163                        +par[2]*TMath::Exp(-xx/par[1]));    
00164     return f;
00165 }
00166 
00167 void AngClusterFitAna::Analyze(int evtn
00168                             ,RecRecordImp<RecCandHeader> *srobj)
00169 {
00170 
00171     if (ANtpDefVal::IsDefault(fPrimShow)){
00172         
00173         MSG("AngClusterFitAna",Msg::kInfo)<< "No primary cluster " 
00174                                           << " found in event."
00175                                           << endl << endl;
00176         fAngClusterFit.Reset();
00177         return;
00178     }
00179 
00180     if(fZ.size() < 3 || fZ.size() > 1000) return;  
00181 
00182     if(fPrimDir.X()==0.0 && fPrimDir.Y()==0.0 && fPrimDir.Z()==0) return;
00183     
00184     MSG("AngClusterFitAna",Msg::kDebug)<<"In AngClusterFitAna::Analyze"<<endl;
00185     MSG("AngClusterFitAna",Msg::kDebug)<<"On Snarl "
00186                                       <<srobj->GetHeader().GetSnarl()
00187                                       <<" event "<<evtn<<endl;
00188     
00189     if(srobj==0){
00190         return;
00191     }
00192     if(((dynamic_cast<NtpStRecord *>(srobj))==0)&&
00193        ((dynamic_cast<NtpSRRecord *>(srobj))==0)){
00194         return;
00195     }
00196 
00197     Float_t cosZAngle=fPrimDir.CosTheta(); //Angle between 
00198                                            //Primary Shower dir 
00199                                            //and the Z axis.
00200     
00201     const Int_t kLHistBin=40;
00202     const Int_t kTHistBin=40;
00203     
00204     const Float_t kTHistXLim=static_cast<Float_t>(kTHistBin)*kStripWidth;
00205 
00206     char lhh[100];
00207     char thh[100];
00208     char lhhr[100];
00209     char thhr[100];
00210 
00211     sprintf(lhh ,"LongHitHisto_%d_%d"    ,srobj->GetHeader().GetSnarl(),evtn);
00212     sprintf(lhhr,"LongHitHistoRad_%d_%d" ,srobj->GetHeader().GetSnarl(),evtn);
00213     sprintf(thh ,"TransHitHisto_%d_%d"   ,srobj->GetHeader().GetSnarl(),evtn);
00214     sprintf(thhr,"TransHitHistoRad_%d_%d",srobj->GetHeader().GetSnarl(),evtn);
00215    
00216 
00217     fLongHitHisto  = new TH1F(lhh,"longitudinal energy hit profile"
00218                               ,kLHistBin,0.0,kLHistBin);
00219     
00220     fLongHitHistoRad  = new TH1F(lhhr,"longitudinal energy hit profile(X0)"
00221                                  ,kLHistBin
00222                                  ,0.0
00223                                  ,kLHistBin*kRadLengthToPlane/cosZAngle);
00224    
00225     fTransHitHisto = new TH1F(thh,"trasverse energy hit profile",
00226                               kTHistBin,-kTHistXLim,kTHistXLim);
00227 
00228     fTransHitHistoRad = new TH1F(thhr,"trasverse energy hit profile(Rm)",
00229                                  kTHistBin,0,kTHistBin);
00230 
00231 
00232     TVector3 unitDir; //Unit vector parallel to PrimShow dir
00233     unitDir=fPrimDir.Unit();
00234 
00235     TRotation r,a;
00236 
00237     r.SetToIdentity();
00238     a.SetToIdentity();
00239     
00240     if (TMath::Abs(fPrimDir.Z()) > 0.0){  //Do not rotate Z axis if shower dir
00241                                           //goes along it.
00242 
00243         r.SetZAxis(unitDir); //Rotates the Z axis to be parallel to the 
00244                              //primDir vector
00245        
00246         a=r.Inverse(); //Inverting the rotation means when one
00247                        //applies the transformation to the hit 
00248                        //vector, one obtains the primDir coordinates
00249     }
00250     
00251     TVector3 hitTemp;
00252 
00253     DeqDeqInt_t cluTemp; 
00254     cluTemp.push_back(fClusterMap.at(fPrimShow)); //Get PrimShow cluster
00255 
00256     Float_t rPrime; //Hit Radius around the shower direction. 
00257 
00258     Float_t totalHitEnergy=0.0;
00259 
00260     //Loop over primary shower hits and fill Long and Trans histograms
00261     for (IterDeqDeqInt_t cluIter = cluTemp.begin()
00262              ;cluIter != cluTemp.end()
00263              ; ++cluIter ){
00264         for (IterDeqInt_t hitIter =  cluIter->begin()
00265                  ;hitIter < cluIter->end()
00266                  ; ++hitIter){
00267             
00268             hitTemp.SetX(fX.at(*hitIter));
00269             hitTemp.SetY(fY.at(*hitIter));
00270             hitTemp.SetZ(fZ.at(*hitIter));
00271      
00272             totalHitEnergy+=(fE.at(*hitIter))*kMeuToGeV;
00273 
00274             Float_t nPlane=hitTemp.Z()/(kDetectorPitch);
00275             
00276             Float_t nPlaneDir=nPlane/(cosZAngle);
00277 
00278             Float_t radDist=nPlaneDir*kRadLengthToPlane; //Longitudinal hit 
00279                                                          //distance from 
00280                                                          //vertex in X0 units.
00281 
00282             //Transform hit coordinate system to Zprime
00283             hitTemp.Transform(a);
00284 
00285             rPrime=TMath::Sqrt(hitTemp.X()*hitTemp.X()
00286                                +hitTemp.Y()*hitTemp.Y());
00287 
00288             Float_t rPrimeRad=rPrime/(kMolRadDet); //Transverse hit distance 
00289                                                    //from Zprime in Moliere 
00290                                                    //radius units
00291 
00292             fLongHitHisto    ->Fill(nPlane,fE.at(*hitIter)*kMeuToGeV);
00293             
00294             fLongHitHistoRad ->Fill(radDist,fE.at(*hitIter)*kMeuToGeV);
00295 
00296             fTransHitHisto   ->Fill(hitTemp.X(),fE.at(*hitIter)*kMeuToGeV);
00297 
00298             fTransHitHistoRad->Fill(rPrimeRad,fE.at(*hitIter)*kMeuToGeV);
00299 
00300         }//for (hitIter)
00301     }//for(cluIter)
00302 
00303     //Do the Longitudinal and Transverse fits
00304     FitShower(totalHitEnergy);
00305 
00306     //fill transverse variables
00307     TransVarCalc(fTransHitHisto);
00308 
00309     MSG("AngClusterFitAna",Msg::kDebug)<<"Leaving AngClusterFitAna::Analyze"
00310                                        <<endl;
00311 }
00312 
00313 void AngClusterFitAna::FitShower(Float_t totalHitEnergy)
00314 {
00315         
00316     //Configure fit functions
00317     Int_t lMin=0;
00318     Int_t lMax=100;
00319     Int_t tMin=0;
00320     Int_t tMax=100;
00321     Int_t nParLong=3;
00322     Int_t nParTrans=4;
00323     char lfn[100];
00324     char tfn[100];
00325     sprintf(lfn,"lfit_%d_%d",lMin,lMax);
00326     sprintf(tfn,"tfit_%d_%d",tMin,tMax);
00327     
00328     fLongFit  = new TF1(lfn,shwLongFunc,lMin+0.001,lMax,nParLong);
00329     fTransFit = new TF1(tfn,shwTransFunc,tMin+0.001,tMax,nParTrans);
00330    
00331     fLongFit->SetParNames("a","b","E0");
00332     fLongFit->SetParLimits(0,lMin+0.001,lMax);
00333     fLongFit->SetParLimits(1,0.001,20000);
00334     fLongFit->SetParLimits(2,0.001,20000);
00335 
00336     fTransFit->SetParNames("l1,l2,c12,E0");
00337     fTransFit->SetParLimits(0,tMin+0.001,tMax);
00338     fTransFit->SetParLimits(1,tMin+0.001,tMax);
00339     fTransFit->SetParLimits(2,0.001,50);
00340     fTransFit->SetParLimits(3,0.001,20000);
00341 
00342     //Do the Longitudinal fit
00343     fLongFit->SetParameters(3.,0.5,totalHitEnergy);
00344     fLongHitHistoRad->Fit(fLongFit,"RLQ0+");
00345 
00346     string fitStatus = (string)(gMinuit->fCstatu);
00347     MSG("AngClusterFitAna",Msg::kDebug)<<"LongFit status:  "<<fitStatus<<endl;
00348         
00349     if(fitStatus=="CONVERGED " 
00350        && fLongFit->GetParameter(0)<29.9
00351        && fLongFit->GetNDF()>0
00352        && totalHitEnergy > 0.0){ 
00353         
00354         fAngClusterFit.fACluFitParA     =fLongFit->GetParameter(0);
00355         fAngClusterFit.fACluFitParB     =fLongFit->GetParameter(1);
00356         fAngClusterFit.fACluFitParLongE0=fLongFit->GetParameter(2);
00357         fAngClusterFit.fACluFitShwMax   =GetMaximum(fLongFit);
00358         fAngClusterFit.fACluFitLongChiSq=fLongFit->GetChisquare();
00359         fAngClusterFit.fACluFitLongConv =1;
00360         fAngClusterFit.fACluFitLongNDF
00361              =fAngClusterFit.fACluFitLongChiSq/fLongFit->GetNDF();
00362         fAngClusterFit.fACluFitE0EnergyRatio
00363             =fAngClusterFit.fACluFitParLongE0/totalHitEnergy;
00364     }//if(fitStatus)
00365 
00366     //Do the transverse fit
00367     fTransFit->SetParameters(1.0,1.0,1.0,totalHitEnergy);
00368     fTransHitHistoRad->Fit(fTransFit,"RLQ0+");
00369 
00370     fitStatus = (string)(gMinuit->fCstatu);
00371     MSG("AngClusterFitAna",Msg::kDebug)<<"TransFit status:  "<<fitStatus<<endl;
00372         
00373 
00374     if(fitStatus=="CONVERGED " 
00375        && fTransFit->GetNDF()>0 ){ 
00376         fAngClusterFit.fACluFitParL1     =fTransFit->GetParameter(0);
00377         fAngClusterFit.fACluFitParL2     =fTransFit->GetParameter(1);
00378         fAngClusterFit.fACluFitParC12    =fTransFit->GetParameter(2);
00379         fAngClusterFit.fACluFitParTransE0=fTransFit->GetParameter(3);
00380         fAngClusterFit.fACluFitTransChiSq=fTransFit->GetChisquare();
00381         fAngClusterFit.fACluFitTransConv =1;
00382         fAngClusterFit.fACluFitTransNDF
00383              =fAngClusterFit.fACluFitTransChiSq/fTransFit->GetNDF();
00384     }
00385 
00386     MSG("AngClusterFitAna",Msg::kDebug)<< "total energy " 
00387                                        << totalHitEnergy 
00388                                        << " from long histo " 
00389                                        << fLongHitHistoRad->Integral() 
00390                                        << " from trans histo " 
00391                                        << fTransHitHistoRad->Integral() 
00392                                        << endl;
00393 }
00394 
00395 void AngClusterFitAna::TransVarCalc(TH1F *transHisto)
00396 {
00397     Int_t kTransHistoBin = static_cast <Int_t> (transHisto->GetNbinsX());
00398     
00399     Int_t binMax   =transHisto->GetMaximumBin();
00400     Int_t binVert  =Int_t(((kTransHistoBin-1)/2)+1);
00401 
00402     Float_t peakBin=transHisto->Integral(binMax,binMax);
00403     Float_t vertBin=transHisto->Integral(binVert,binVert);
00404     Float_t total  =transHisto->Integral(1,kTransHistoBin);
00405   
00406     Float_t asymPeak= ANtpDefVal::kFloat;
00407     Float_t asymVert= ANtpDefVal::kFloat;
00408 
00409     MSG("AngClusterFitAna",Msg::kDebug)<<"total:  "
00410                                        <<total<<" binMax: "<<binMax
00411                                        <<" binVert "<<binVert
00412                                        <<"  peakBin: "<<peakBin
00413                                        <<" kTransHistBin "<<kTransHistoBin
00414                                        <<endl;
00415 
00416     MSG("AngClusterFitAna",Msg::kDebug)<<"HistoI1:  "
00417                                        <<transHisto->Integral(binMax+1
00418                                                               ,kTransHistoBin)
00419                                        <<endl;
00420     
00421     MSG("AngClusterFitAna",Msg::kDebug)<<"HistoI1:  "
00422                                        <<transHisto->Integral(1, binMax-1)
00423                                        <<endl;
00424     
00425     MSG("AngClusterFitAna",Msg::kDebug)<<"HistoI1:  "
00426                                        <<transHisto->Integral(binVert+1
00427                                                               ,kTransHistoBin)
00428                                        <<endl;
00429     
00430     MSG("AngClusterFitAna",Msg::kDebug)<<"HistoI1:  "
00431                                        <<transHisto->Integral(1, binVert-1)
00432                                        <<endl;
00433    
00434 
00435 
00436     
00437     if(TMath::Abs(total-peakBin) > 1e-5){
00438         MSG("AngClusterFitAna",Msg::kDebug)<<"Evaluating asymPeak:  "
00439                                            <<TMath::Abs(transHisto
00440                                                         ->Integral(binMax+1
00441                                                         ,kTransHistoBin)
00442                                                         -transHisto
00443                                                         ->Integral(1,binMax-1))
00444                                            <<" divided by "<<(total-peakBin)
00445                                            <<endl;
00446                 
00447         asymPeak=(TMath::Abs(transHisto->Integral(binMax+1,kTransHistoBin)
00448                              -transHisto->Integral(1,binMax-1)))
00449                              /(total-peakBin);
00450     }
00451 
00452     
00453     if(TMath::Abs(total-vertBin) > 1e-5){
00454        MSG("AngClusterFitAna",Msg::kDebug)<<"Evaluating asymVert:  "<<
00455                TMath::Abs(transHisto->Integral(binVert+1,kTransHistoBin)
00456                                             -transHisto->Integral(1,binVert-1))
00457                   <<" divided by "<<(total-peakBin)<<endl;
00458                     
00459         asymVert=(TMath::Abs(transHisto->Integral(binVert+1,kTransHistoBin)
00460                              -transHisto->Integral(1,binVert-1)))
00461                              /(total-vertBin);
00462     }
00463 
00464 
00465     MSG("AngClusterFitAna",Msg::kDebug)<<"asymPeak:  "<<asymPeak
00466                                        <<" asymVert: "<<asymVert<<endl;
00467     
00468     Float_t ratio,molRadPeak=0,molRadVert=0;
00469     
00470     if(total){
00471         for(Int_t i=0; i<=binVert-1;i++){
00472             ratio=transHisto->Integral(binMax-i>0?binMax-i
00473                                        :1,binMax+i<kTransHistoBin
00474                                        ?binMax+i:kTransHistoBin)/total;
00475             if(ratio>0.90){molRadPeak=i+1; break;}
00476         }
00477         for(Int_t i=0; i<=binVert-1;i++){
00478             ratio=transHisto->Integral(binVert-i,binVert+i)/total;
00479             if(ratio>0.90){molRadVert=i+1; break;}
00480         }
00481         
00482     }
00483     
00484     Float_t mean=transHisto->GetMean();
00485     Float_t rms =transHisto->GetRMS();
00486     Float_t skew=0,kurt=0;
00487     Int_t nCount=0;
00488                                                    
00489     MSG("AngClusterFitAna",Msg::kDebug)<<"mean:  "<<mean<<" rms: "<<rms<<endl;
00490         
00491     
00492     for(Int_t i=1;i<=kTransHistoBin;i++){
00493         
00494         if(transHisto->GetBinContent(i)){
00495             skew=skew+(TMath::Power((transHisto->GetBinCenter(i)-mean),3)
00496                        *transHisto->GetBinContent(i));
00497             kurt=kurt+(TMath::Power((transHisto->GetBinCenter(i)-mean),4)
00498                        *transHisto->GetBinContent(i));
00499             nCount++;
00500         }
00501         
00502     }
00503     
00504     if(rms>1e-8 && nCount>1){
00505         skew=skew /(static_cast<Float_t>(nCount-1)*TMath::Power((rms),3.));
00506         kurt=(kurt/(static_cast<Float_t>(nCount-1)*TMath::Power((rms),4.)))-3.;
00507         if(TMath::Abs(skew) > 1.e8 ) skew = ANtpDefVal::kFloat;
00508         if(TMath::Abs(kurt) > 1.e8 ) kurt = ANtpDefVal::kFloat;
00509     }
00510     else{skew= ANtpDefVal::kFloat; kurt= ANtpDefVal::kFloat;}
00511     
00512     //Fill variables
00513     fAngClusterFit.fACluFitAsymPeak=asymPeak;
00514     fAngClusterFit.fACluFitAsymVert=asymVert;
00515     fAngClusterFit.fACluFitMolRadPeak=molRadPeak;
00516     fAngClusterFit.fACluFitMolRadVert=molRadVert;
00517     fAngClusterFit.fACluFitMean=mean;
00518     fAngClusterFit.fACluFitRMS=rms;
00519     fAngClusterFit.fACluFitSkew=skew;
00520     fAngClusterFit.fACluFitKurt=kurt;
00521    
00522     MSG("AngClusterFitAna",Msg::kDebug)
00523         <<"Leaving AngClusterFitAna::TransCalc"
00524         <<endl;
00525 }
00526 
00527 void AngClusterFitAna::Draw(TPad *pad)
00528 {
00529 
00530     if (ANtpDefVal::IsDefault(fPrimShow)){
00531         
00532         MSG("AngClusterFitAna",Msg::kInfo)<< "No primary cluster " 
00533                                           << " found in event."
00534                                           << endl << endl;
00535         fAngClusterFit.Reset();
00536         return;
00537     }
00538 
00539     if(fZ.size() < 3 || fZ.size() > 1000) return;  
00540     if(fPrimDir.X()==0.0 && fPrimDir.Y()==0.0 && fPrimDir.Z()==0) return;
00541 
00542     DeqTPoly scatterCluster;
00543     TPolyMarker3D *scatterTemp;
00544 
00545     DeqDeqInt_t scatTemp; 
00546 
00547     Int_t nClu=0;
00548     Int_t nHit=0;
00549 
00550     //Loop over all cluster hits and fill polymarkers
00551     for (IterDeqDeqInt_t cluIter = fClusterMap.begin()
00552              ;cluIter != fClusterMap.end()
00553              ; ++cluIter ){
00554         scatTemp.push_back(fClusterMap.at(nClu));
00555         scatterTemp = new TPolyMarker3D(scatTemp.size(),1);
00556         
00557         for (IterDeqInt_t hitIter =  cluIter->begin()
00558                  ;hitIter < cluIter->end()
00559                  ; ++hitIter){
00560             
00561             scatterTemp->SetPoint(nHit,fZ.at(*hitIter)
00562                                   ,fX.at(*hitIter)
00563                                   ,fY.at(*hitIter));
00564             nHit++;
00565         }//for (hitIter)
00566         scatterCluster.push_back(scatterTemp);
00567         nHit=0;
00568         nClu++;
00569     }//for(cluIter)
00570 
00571     // Find min and max values;
00572     Float_t minX=static_cast <Float_t> (*min_element(fX.begin(),fX.end()));
00573     Float_t maxX=static_cast <Float_t> (*max_element(fX.begin(),fX.end()));
00574     Float_t minY=static_cast <Float_t> (*min_element(fY.begin(),fY.end()));
00575     Float_t maxY=static_cast <Float_t> (*max_element(fY.begin(),fY.end()));
00576     Float_t minZ=static_cast <Float_t> (*min_element(fZ.begin(),fZ.end()));
00577     Float_t maxZ=static_cast <Float_t> (*max_element(fZ.begin(),fZ.end()));
00578     
00579     if(fScatterAxis) delete fScatterAxis;
00580     fScatterAxis = new TH3F("cluHisto","Scatter plot of all hits"
00581                            , 10,minZ-0.1,maxZ+0.2
00582                            , 10,minX-0.2,maxX+0.2
00583                            , 10,minY-0.2,maxY+0.2);
00584     
00585     fScatterAxis->SetXTitle("Z(m)");
00586     fScatterAxis->SetYTitle("X(m)");
00587     fScatterAxis->SetZTitle("Y(m)");
00588     
00589     TPolyMarker3D *evtVertex;
00590     evtVertex = new TPolyMarker3D(1,8);
00591     evtVertex->SetMarkerColor(1);
00592     evtVertex->SetPoint(0,0.,0.,0.);
00593     
00594     TPolyLine3D *lineCluster;
00595     lineCluster= new TPolyLine3D(2);
00596     
00597     lineCluster->SetPoint(0,0.,0.,0.);
00598     lineCluster->SetPoint(1,fPrimDir.Z(),fPrimDir.X(),fPrimDir.Y());
00599     
00600     pad->cd(1);
00601     
00602     fScatterAxis->Draw();
00603 
00604     Int_t polyColor=1;
00605 
00606     for (IterDeqTPoly polyIter = scatterCluster.begin()
00607              ;polyIter < scatterCluster.end()
00608              ; ++polyIter){
00609         (*polyIter)->SetMarkerStyle(7);
00610         (*polyIter)->SetMarkerColor(polyColor);
00611         (*polyIter)->Draw("sames");
00612         polyColor++;   
00613     }//for(polyIter)
00614 
00615     evtVertex->Draw("sames");
00616     lineCluster->SetLineWidth(2);
00617     lineCluster->Draw("sames");
00618 
00619     pad->cd(3);
00620     fLongHitHistoRad->Draw();
00621     fLongFit->Draw("sames");
00622     pad->cd(2);
00623     
00624     fTransHitHistoRad->Draw();
00625     fTransFit->Draw("sames");
00626     pad->cd(4);
00627     fTransHitHisto->Draw();
00628     pad->cd();
00629     pad->Modified();
00630     
00631     //Clearing ROOT objects
00632     evtVertex->Delete();
00633     scatterTemp->Delete();
00634     lineCluster->Delete();
00635     
00636 }
00637 
00638 void AngClusterFitAna::Reset()
00639 {
00640 
00641 
00642 }
00643 
00644 Double_t AngClusterFitAna::GetMaximum(TF1* efit, Double_t xmin, Double_t xmax)
00645 {
00646    Double_t fXmin = 0.001;
00647    Double_t fXmax = 100;
00648    Double_t fNpx = 40;
00649                                                                                 
00650    Double_t x,y;
00651    if (xmin >= xmax) {xmin = fXmin; xmax = fXmax;}
00652    Double_t dx = (xmax-xmin)/fNpx;
00653    Double_t xxmax = xmin;
00654    Double_t yymax = efit->Eval(xmin+dx);
00655    for (Int_t i=0;i<fNpx;i++) {
00656       x = xmin + (i+0.5)*dx;
00657       y = efit->Eval(x);
00658       if (y > yymax) {xxmax = x; yymax = y;}
00659    }
00660    if (dx < 1.e-9*(fXmax-fXmin)) return TMath::Min(xmax,xxmax);
00661    else return GetMaximum(efit, TMath::Max(xmin,xxmax-dx), TMath::Min(xmax,xxmax+dx));
00662 }
00663