loon: AngClusterAna.cxx Source File

00001 
00002 
00003 
00004 
00005 
00006 
00007 
00008 
00009 
00010 
00011 
00012 
00013 #include <iostream>
00014 #include <algorithm>
00015 #include <TCanvas.h>
00016 #include <TFile.h>
00017 #include <TH2F.h>
00018 #include <TH1F.h>
00019 #include "StandardNtuple/NtpStRecord.h"
00020 #include "CandNtupleSR/NtpSRRecord.h"
00021 #include "CandNtupleSR/NtpSREvent.h"
00022 #include "CandNtupleSR/NtpSRTrack.h"
00023 #include "CandNtupleSR/NtpSRStrip.h"
00024 #include "MessageService/MsgService.h"
00025 #include "NueAna/AngClusterAna.h"
00026 #include "AnalysisNtuples/ANtpDefaultValue.h"
00027 
00028 #include "NueAna/NueAnaTools/SntpHelpers.h"
00029 #include "VertexFinder/NtpVtxFinder/NtpVtxFinder.h"
00030 
00031 
00032 
00033 CVSID("$Id: AngClusterAna.cxx,v 1.24 2007/06/06 15:35:30 boehm Exp $");
00034 
00035 using std::cout;
00036 using std::endl;
00037 
00038 ClassImp(AngClusterAna)
00039 
00040 const Float_t kDetectorPitch=0.05949; //Munits::meter
00041 const Float_t kStripWidth=0.041; //Munits::meters
00042 
00043 const Float_t kMeuToGeV=1.0;  //Conversion factor from MEU to GeV
00044 
00045 //......................................................................
00046 
00047 AngClusterAna::AngClusterAna(AngCluster &ac):
00048     fAngCluster(ac)
00049 {
00050 
00051 
00052 }
00053 
00054 AngClusterAna::~AngClusterAna()
00055 {
00056 
00057 }
00058 
00059 void AngClusterAna::Analyze(int evtn
00060                             ,RecRecordImp<RecCandHeader> *srobj)
00061 {
00062 
00063     MSG("AngClusterAna",Msg::kInfo)<<"In AngClusterAna::Analyze"<<endl;
00064     MSG("AngClusterAna",Msg::kDebug)<<"On Snarl "<<srobj->GetHeader().GetSnarl()
00065                                  <<" event "<<evtn<<endl;
00066 
00067     if(srobj==0){
00068         return;
00069     }
00070     if(((dynamic_cast<NtpStRecord *>(srobj))==0)&&
00071        ((dynamic_cast<NtpSRRecord *>(srobj))==0)){
00072         return;
00073     }
00074 
00075 
00076     if(fZ.size() >=3 && fZ.size() <= 1000){
00077 
00078         //Find Clusters
00079         FindCluster();
00080 
00081         //Calculate cluster variables and fill branch.
00082         ClusterVarCalc();
00083     }
00084 
00085     WeightedEnergy(evtn, srobj);
00086 }
00087 
00088 void AngClusterAna::Set3DHit(DeqFloat_t &x
00089                              , DeqFloat_t &y
00090                              , DeqFloat_t &z
00091                              , DeqFloat_t &e){
00092 
00093     if(x.size()!=0 && y.size()!=0 && z.size()!=0 && e.size()!=0){
00094 
00095         fX=x;
00096         fY=y;
00097         fZ=z;
00098         fE=e;  
00099     
00100     }
00101 
00102 }
00103 
00104 void AngClusterAna::GetAngCluster(Int_t &primShow
00105                                   , DeqDeqInt_t &clusterMap
00106                                   , TVector3 &primDir)
00107 {
00108     
00109     if(fClusterMap.size()==0){
00110         MSG("AngClusterAna",Msg::kWarning)<< "Primary cluster not found for " 
00111                                        << "event, stopping any further"
00112                                        << " clustering related processing"
00113                                        << endl;            
00114         return;
00115     }
00116 
00117     if(fPrimDir.X()==0.0 && fPrimDir.Y()==0.0 && fPrimDir.Z()==0) return;
00118 
00119     //Reference relevant quantities.
00120     primShow  =fPrimShow;
00121     clusterMap=fClusterMap;
00122     primDir   =fPrimDir;
00123 }
00124 
00125 
00126 
00127 void AngClusterAna::FindCentroidBlob(TMatrixD & grid
00128                                   ,const Int_t nBinX
00129                                   ,const Int_t nBinY
00130                                   ,Int_t indX
00131                                   ,Int_t indY
00132                                   ,TMatrixD & blobCentro)
00133 {
00134     //Recursively find centroid agglomerations in the centroid histogram
00135     // No approximations are used, all centroid aggregates are found.
00136     
00137     //Make sure we did not go over the edges.    
00138     if (indX < 0 || indY < 0 || indX >= nBinX || indY >= nBinY)
00139         return;
00140     
00141     //Finish recursion if the next grid point is empty.
00142     if(TMath::Nint(grid(indX,indY)) != 1)
00143         return;
00144     
00145     //else
00146     grid(indX,indY) = 0.0;    //Mark as used
00147     blobCentro(indX,indY)=1.0;
00148   
00149     MSG("AngClusterAna",Msg::kDebug)<< "Adding point(" << indX
00150                                     <<"," << indY << ") to blob " 
00151                                     << endl;
00152   
00153     FindCentroidBlob(grid, nBinX, nBinY, indX - 1, indY, blobCentro);
00154     FindCentroidBlob(grid, nBinX, nBinY, indX + 1, indY, blobCentro);
00155     FindCentroidBlob(grid, nBinX, nBinY, indX, indY - 1, blobCentro);
00156     FindCentroidBlob(grid, nBinX, nBinY, indX, indY + 1, blobCentro);
00157     
00158 }
00159 
00160 void AngClusterAna::FindCluster()
00161 {
00162 
00163    
00164     //Array of TVector3 hit objects containing
00165     //normal and XZY and ZYX CS rotations  
00166     std::vector<TVector3> hitVecArrayX(fZ.size());
00167     std::vector<TVector3> hitVecArrayY(fZ.size());
00168     
00169     //Vectors to be used in clustering
00170     VecFloat_t cosX;
00171     VecFloat_t cosY;
00172     cosX.assign(fZ.size(),0.0);
00173     cosY.assign(fZ.size(),0.0);
00174     
00175     fTotalEventHitEnergy=0.0;
00176 
00177     for (UInt_t n=0; n < hitVecArrayX.size(); n++){ //Loop over all 3D Hits
00178         
00179         (hitVecArrayX.at(n)).SetZ(fX.at(n));
00180         (hitVecArrayX.at(n)).SetY(fY.at(n));
00181         (hitVecArrayX.at(n)).SetX(fZ.at(n));
00182 
00183         (hitVecArrayY.at(n)).SetX(fX.at(n));
00184         (hitVecArrayY.at(n)).SetZ(fY.at(n));
00185         (hitVecArrayY.at(n)).SetY(fZ.at(n));
00186 
00187         cosX.at(n)=(hitVecArrayX.at(n)).CosTheta();
00188         cosY.at(n)=(hitVecArrayY.at(n)).CosTheta();
00189 
00190         //Calculate total 3D hit event energy  
00191         fTotalEventHitEnergy+=fE.at(n)*kMeuToGeV;
00192         
00193     }
00194 
00195     //============================================================
00196     //Start of cluster finding algorithm (nearest-neighbor method)
00197     //Do it in CosTheta X-Y space for now
00198     
00199     TMatrixD distHit(fZ.size(),fZ.size());   //Holds all distances 
00200                                              //between hits
00201     TMatrixD usedDist(fZ.size(),fZ.size());  //Keep track of already 
00202                                              //used distances
00203     
00204     VecFloat_t usedHit; //keep track of already used hits
00205     usedHit.assign(fZ.size(),0.0);
00206     
00207     Int_t nMin;  //index of smallest inter-hit distance
00208     Int_t kMin;
00209     Float_t minDist=999999.;
00210 
00211     //Calculate all inter-hit distances ommitting reciprocals
00212     for (UInt_t n=0; n < fZ.size(); n++){
00213         for (UInt_t k=0; k < fZ.size(); k++){
00214             distHit(n,k) =0.0;
00215             usedDist(n,k)=0.0;
00216             if (k > n){
00217 
00218                 distHit(n,k)= TMath::Sqrt((cosX.at(k)-cosX.at(n))
00219                                           *(cosX.at(k)-cosX.at(n))
00220                                           +(cosY.at(k)-cosY.at(n))
00221                                           *(cosY.at(k)-cosY.at(n)));
00222                
00223                 if (distHit(n,k) < minDist){
00224                     minDist=distHit(n,k);
00225                     nMin=n;
00226                     kMin=k;
00227                 }
00228                 
00229             } //if
00230         }//for(k)
00231     }//for(n)
00232     
00233     MSG("AngClusterAna",Msg::kDebug)<< "minDist[" 
00234                                     << nMin << "][" << kMin << "] = " 
00235                                     << distHit(nMin,kMin) << endl;
00236 
00237 
00238     Float_t radius=0.5; //limit radius for nearest neighbor 
00239 
00240     DeqInt_t aggHitSize; //Keeps track of the size of each aggregate 
00241 
00242     DeqDeqInt_t aggHitArray; //Container filled with all the 
00243                              //hit aggregate sets.
00244 
00245     UInt_t aggHitSizeMax=0;
00246 
00247     Int_t ahaIndex=0;
00248     Int_t ahaIndexMax=0;
00249 
00250     for (UInt_t n=0; n < fZ.size(); n++){
00251         for (UInt_t k=(n+1); k < fZ.size(); k++){
00252             if (k > n){
00253                 //Take the (n,k) pair as an aggregate seed 
00254                 if (TMath::Nint(usedDist(n,k)) == 0
00255                     && distHit(n,k) < radius
00256                     && distHit(n,k) > 0){
00257                     
00258                     usedDist(n,k)=1.0;
00259                     
00260                     DeqInt_t aggHit; //aggregate hit array
00261                                        
00262                     aggHit.push_back(n); //Store first hit from seed
00263                     aggHit.push_back(k); //Store second hit
00264                     
00265                     //Look for neighbor hits above the upper hit of 
00266                     //the seed pair
00267                     for (UInt_t kTemp=(n+1); kTemp < fZ.size(); kTemp++){
00268                         if (TMath::Nint(usedDist(n,kTemp)) == 0){
00269                             usedDist(n,kTemp)=1.0;
00270                             if (kTemp != k 
00271                                 && distHit(n,kTemp) > 0.0 
00272                                 && distHit(n,kTemp) < radius){
00273                                 
00274                                 aggHit.push_back(kTemp);
00275                                 usedDist(n,kTemp)=1.0;
00276 
00277                             }//if(kTemp !=k)
00278                                                        
00279                         }else continue; //if (usedDist(n,kTemp)) 
00280                   
00281                     }//for(kTemp)
00282                     
00283                     //Look for neighbor hits below the lower hit of 
00284                     //the seed pair
00285                     for (UInt_t nTemp=0; nTemp < k; nTemp++){
00286                         if (TMath::Nint(usedDist(nTemp,k)) == 0){
00287                             usedDist(nTemp,k)=1.0;
00288                             if (nTemp != n 
00289                                 && distHit(nTemp,k) > 0.0 
00290                                 && distHit(nTemp,k) < radius){
00291                                 
00292                                 aggHit.push_back(nTemp);
00293                                 usedDist(nTemp,k)=1.0;
00294 
00295                             }//if(nTemp !=n)
00296 
00297                         }else continue; //if (usedDist(nTemp,k)) 
00298                   
00299                     }//for(nTemp)
00300 
00301                     aggHitSize.push_back(aggHit.size()); //Record aggreg. size
00302                     aggHitArray.push_back(aggHit);       //Store hit aggregate
00303 
00304                     if (aggHit.size() > aggHitSizeMax){
00305                         aggHitSizeMax=aggHit.size();
00306                         ahaIndexMax=ahaIndex;
00307                     }
00308                     ahaIndex++;
00309                 }//if(distHit < radius)
00310                 
00311             }//if(k > n)
00312 
00313         }//for(k)
00314         
00315     }//for(n)
00316     
00317     
00318     MSG("AngClusterAna",Msg::kInfo)<< "Found "
00319                                    << aggHitArray.size() 
00320                                    << " aggregates."
00321                                    << endl;
00322     
00323     MSG("AngClusterAna",Msg::kInfo)<< "Largest aggregate contains "
00324                                    << aggHitSizeMax 
00325                                    << " Hits and has Index "
00326                                    << ahaIndexMax
00327                                    << endl;
00328  
00329     //Cluster determination:
00330     //Eliminate all overlapping aggregates in favor of the largest one.
00331     //Method: Calculate and histogram aggregate centroids. Run a 
00332     //recursion blob finder function on the histogram and determine 
00333     //high centroid density region bounds. From all the centroids
00334     //consistent with a region boundary, choose the one 
00335     //corresponding to the aggregate with most hits and call 
00336     //the aggregate a cluster.
00337 
00338     const Int_t kBinCentroX=10;
00339     const Int_t kBinCentroY=10;
00340 
00341     const Float_t kMinCentroX =-1.0;
00342     const Float_t kMaxCentroX = 1.0;
00343     const Float_t kMinCentroY =-1.0;
00344     const Float_t kMaxCentroY = 1.0;
00345 
00346     Float_t centSumX=0.0;
00347     Float_t centSumY=0.0;
00348     Float_t centValX=0.0;
00349     Float_t centValY=0.0;
00350     
00351     // Write out centroid histogram for debugging purposes 
00352     Bool_t kDebug=0;
00353     
00354     TH2F *centroHisto;
00355     centroHisto = new TH2F("centroHisto","CosX vs CosY centroid histogram"
00356                            ,kBinCentroX,kMinCentroX,kMaxCentroX
00357                            ,kBinCentroY,kMinCentroY,kMaxCentroY);
00358 
00359     centroidX.clear();
00360     centroidY.clear();
00361 
00362     for (IterDeqDeqInt_t ahaIter = aggHitArray.begin()
00363              ;ahaIter != aggHitArray.end()
00364              ; ++ahaIter ){
00365         centSumX=0.0;
00366         centSumY=0.0;
00367         for (IterDeqInt_t ahIter = ahaIter->begin()
00368                  ;ahIter < ahaIter->end()
00369                  ; ++ahIter){
00370             
00371             centSumX+=cosX.at((*ahIter));
00372             centSumY+=cosY.at((*ahIter)); 
00373         }//for(ahIter)
00374         
00375         centValX=centSumX/(static_cast <Float_t> (ahaIter->size()));
00376         centValY=centSumY/(static_cast <Float_t> (ahaIter->size()));
00377         
00378         centroidX.push_back(centValX);
00379         centroidY.push_back(centValY);
00380         
00381         centroHisto->Fill(centValX,centValY);
00382         
00383     }//for(ahaIter)
00384 
00385         
00386     TMatrixD gridHisto(kBinCentroX,kBinCentroY); //Centroid Grid
00387                                                  //for blob finding 
00388 
00389     const Float_t kMinCentroid=1.0; //Minimum number of centroids 
00390                                     //to be found in a  bin that 
00391                                     // will be part of a blob;
00392     
00393     if (kMinCentroid < 1.0){
00394         MSG("AngClusterAna",Msg::kFatal)<< "Attempting to cluster "
00395                                         << "empty space. Check "
00396                                         << "value of kMinCentroid"
00397                                         << endl;
00398     }
00399     
00400     Float_t centroStepX
00401         =(kMaxCentroX-kMinCentroX)/static_cast <Float_t> (kBinCentroX);
00402     Float_t centroStepY
00403         =(kMaxCentroY-kMinCentroY)/static_cast <Float_t> (kBinCentroY);
00404 
00405     //Map each centroid to bin position.
00406     VecInt_t centroBinX;
00407     centroBinX.assign(aggHitSize.size(),0);
00408     VecInt_t centroBinY;
00409     centroBinY.assign(aggHitSize.size(),0);
00410 
00411     //Populate blob grid.
00412     //Reminder: ROOT bin indexing starts at 1
00413     for (Int_t cBinX=0; cBinX < kBinCentroX; cBinX++){
00414       for (Int_t cBinY=0; cBinY < kBinCentroY; cBinY++){
00415 
00416           gridHisto(cBinX,cBinY)=0.0;
00417           
00418           //Find bin coordinates for each centroid
00419           for(UInt_t ah=0; ah < aggHitArray.size(); ah++){
00420 
00421               if (centroidX.at(ah) > (centroHisto->GetXaxis()
00422                                       ->GetBinLowEdge(cBinX+1)) 
00423                   && centroidX.at(ah) < centroStepX
00424                                         +(centroHisto
00425                                           ->GetXaxis()
00426                                           ->GetBinLowEdge(cBinX+1)) 
00427                   && centroidY.at(ah) > (centroHisto->GetYaxis()
00428                                          ->GetBinLowEdge(cBinY+1)) 
00429                   && centroidY.at(ah) < centroStepY
00430                                         +(centroHisto
00431                                           ->GetYaxis()
00432                                           ->GetBinLowEdge(cBinY+1))){
00433               
00434                   centroBinX.at(ah)=cBinX;
00435                   centroBinY.at(ah)=cBinY;
00436 
00437               }//if(centroid)
00438                   
00439           }//for(ah)
00440 
00441 
00442           //Ignore bins with low centroid density
00443           Int_t cBin=centroHisto->GetBin(cBinX+1,cBinY+1);
00444           if (centroHisto->GetBinContent(cBin) >= kMinCentroid ){
00445           
00446               gridHisto(cBinX,cBinY)=1.0;
00447 
00448           }//if(content)
00449       }//for(cBinY)
00450     }//for(cBinX)
00451 
00452     TMatrixD blobCentro(kBinCentroX,kBinCentroY); //Contains all centroids
00453                                                   //on a given blob
00454 
00455     DeqTMatrixD blobArray; //Contains all the found blobs 
00456 
00457     //Find blobs recursively in the centroid histogram
00458     for (Int_t nX=0; nX < kBinCentroX; nX++){
00459         for (Int_t nY=0; nY < kBinCentroY; nY++){
00460             blobCentro(nX,nY)=0.0;
00461             if(TMath::Nint(gridHisto(nX,nY))==1.0){
00462                 
00463                 FindCentroidBlob(gridHisto
00464                                  ,kBinCentroX
00465                                  ,kBinCentroY
00466                                  ,nX
00467                                  ,nY
00468                                  ,blobCentro);
00469 
00470                 blobArray.push_back(blobCentro);
00471                 if(kDebug) blobCentro.Print();
00472                 blobCentro.Zero();
00473             }//if(gridHisto)
00474         }//for(nX)
00475     }//for(nY)
00476 
00477     
00478     Int_t sizeMax=-1;   //keep track of max agg size within a blob
00479     Int_t cluSelect=-1; //index of max agg size in a blob
00480 
00481     VecInt_t cluAggMap; //Correspondance between final cluster and
00482                         //largest Hit aggregate for that cluster
00483     cluAggMap.assign(aggHitArray.size(),0);
00484 
00485     Int_t cluSizeMax=0;   //keep track of maximum size over all blobs
00486     Int_t cluPrimShow=-1; //index of largest aggregate over all blobs
00487     Int_t aggIndexMax=-1;
00488 
00489     fClusterSize.clear(); //Prepare for the new event.
00490     fClusterMap.clear();
00491 
00492     Int_t blobIndex=0;
00493 
00494     //Find the largest aggregate for each blob and call that a cluster
00495     for (IterDeqTMatrixD baIter = blobArray.begin()
00496              ;baIter != blobArray.end()
00497              ; ++baIter ){
00498 
00499         sizeMax=-1;
00500         cluSelect=-1;
00501 
00502         for (Int_t i = 0; i < kBinCentroX; i++){
00503             for (Int_t j = 0; j < kBinCentroY; j++){
00504                 if(TMath::Nint((*baIter)(i,j))==1){
00505 
00506                     for(UInt_t ah=0; ah < aggHitArray.size(); ah++){
00507 
00508                         if (centroBinX.at(ah)==i && centroBinY.at(ah)==j){
00509 
00510                             if(aggHitSize.at(ah) > sizeMax){
00511                                 sizeMax=aggHitSize.at(ah);
00512                                 cluSelect=ah;
00513                             }//if(aggHitSize)
00514 
00515                         }//if(centroBin)
00516 
00517                     }//for(ah)
00518 
00519                 }//if(*baIter(i,j)==1)
00520             }//for(i)
00521         }//for(i)
00522 
00523         if(cluSelect >=0){
00524             cluAggMap.at(blobIndex)=cluSelect;
00525 
00526             //Fill cluster size member deque;
00527             fClusterSize.push_back(aggHitSize.at(cluSelect));
00528 
00529             //Fill cluster map member container
00530             fClusterMap.push_back(aggHitArray.at(cluSelect));
00531 
00532             //Determination of largest cluster index ("primary shower")
00533             if (aggHitSize.at(cluSelect) > cluSizeMax){
00534                 cluSizeMax=aggHitSize.at(cluSelect);
00535                 cluPrimShow=blobIndex;
00536                 aggIndexMax=cluSelect;
00537             }
00538 
00539             MSG("AngClusterAna",Msg::kDebug)<< "Cluster "
00540                                        << blobIndex
00541                                        << " largest agg is "
00542                                        << cluSelect
00543                                        << " with "
00544                                        << sizeMax
00545                                        << " Hits"
00546                                        << endl;
00547 
00548             blobIndex++;
00549 
00550         } else{
00551             MSG("AngClusterAna",Msg::kWarning)<< "Unable to find "
00552                                            << " largest cluster."
00553                                            << " Quitting."
00554                                            << endl << endl;
00555             break;
00556         }
00557 
00558     }//for(baIter)
00559 
00560 
00561 
00562     if(cluPrimShow >= 0){
00563         fPrimShow=cluPrimShow;
00564         MSG("AngClusterAna",Msg::kInfo)<< "Primary Cluster "
00565                                        << "found in blob "
00566                                        << cluPrimShow << "/"
00567                                        << blobArray.size()
00568                                        << " contains "
00569                                        << aggHitSize.at
00570                                        (cluAggMap.at(cluPrimShow))
00571                                        << " hits and was found in aggregate "
00572                                        << aggIndexMax
00573                                        << endl;
00574     }else {
00575 
00576         MSG("AngClusterAna",Msg::kInfo)<< "No primary cluster "
00577                                        << " found in event."
00578                                        << endl << endl;
00579         if(centroHisto)
00580             delete centroHisto;
00581 
00582         fPrimShow=ANtpDefVal::kInt;
00583 
00584         return;
00585     }
00586 
00587 
00588     //Find all high energy hits close to the vertex (low angular resolution)
00589     DeqInt_t highEnergyHits;
00590 
00591     for (UInt_t n=0; n < fZ.size(); n++){
00592         if(fZ.at(n)/kDetectorPitch > 0. && fZ.at(n)/kDetectorPitch < 4.){
00593             if(TMath::Sqrt(fX.at(n)*fX.at(n)+fY.at(n)*fY.at(n))
00594                /kStripWidth < 4.){
00595 
00596                 if(fE.at(n) > 5.0){
00597                     highEnergyHits.push_back(n);
00598                     //Add energy of close hits to total event energy
00599                     fTotalEventHitEnergy+=fE.at(n)*kMeuToGeV;
00600                     MSG("AngClusterAna",Msg::kDebug)<< "(X,Y,Z,E) HighHits = "
00601                                                    << "("
00602                                                    << fX.at(n) << ","
00603                                                    << fY.at(n) << ","
00604                                                    << fZ.at(n) << ","
00605                                                    << fE.at(n) << ")"
00606                                                    << endl;
00607                 }//for(fE.at(n))
00608             }//for(fX.at(n))
00609         } //for(fZ.at(n))
00610     }// for(n)
00611 
00612     //Update the primary shower cluster with the high energy hits above
00613     DeqDeqInt_t clusterMapTemp;
00614     DeqInt_t    hitTemp;
00615     Int_t       cluIndexTemp=0;
00616 
00617     for (IterDeqDeqInt_t cluIter = fClusterMap.begin()
00618              ;cluIter != fClusterMap.end()
00619              ; ++cluIter ){
00620 
00621         hitTemp=*cluIter;
00622 
00623         if (cluIndexTemp==fPrimShow){
00624 
00625             for (UInt_t n=0; n < highEnergyHits.size(); n++){
00626                 hitTemp.push_back(highEnergyHits.at(n));
00627             }//for(n)
00628 
00629         }//if(fPrimShow)
00630 
00631         clusterMapTemp.push_back(hitTemp);
00632         cluIndexTemp++;
00633 
00634     }//for(cluIter)
00635 
00636     //Fill member variables
00637     fClusterMap.clear();
00638     fClusterMap=clusterMapTemp;
00639 
00640     fClusterSize.at(fPrimShow)+=highEnergyHits.size();
00641 
00642     fNCluster=fClusterMap.size();
00643 
00644     Int_t cluIndex=0;
00645 
00646     fClusterID.clear();
00647 
00648     fClusterID.assign(fZ.size(),-1); //hits outside any cluster get -1.
00649 
00650     //Build the clusterID vector (returns the cluster one hit belongs to.)
00651     for (IterDeqDeqInt_t cluIter = fClusterMap.begin()
00652              ;cluIter != fClusterMap.end()
00653              ; ++cluIter ){
00654         for (IterDeqInt_t hitIter = cluIter->begin()
00655                  ;hitIter < cluIter->end()
00656                  ; ++hitIter){
00657 
00658             fClusterID.at((*hitIter))=cluIndex;
00659 
00660         }//for (hitIter)
00661 
00662         cluIndex++;
00663 
00664     }//for(cluIter)
00665 
00666 
00667 
00668     //Delete centroid histogram
00669     if(centroHisto)
00670         delete centroHisto;
00671 
00672 }
00673 
00674 void AngClusterAna::ClusterVarCalc()
00675 {
00676     if (ANtpDefVal::IsDefault(fPrimShow)){
00677 
00678         MSG("AngClusterAna",Msg::kInfo)<< "No primary cluster "
00679                                        << " found in event."
00680                                        << endl << endl;
00681         fAngCluster.Reset();
00682         return;
00683     }
00684 
00685     //Get primary cluster centroid coordinates
00686     Float_t cosXCent=centroidX.at(fPrimShow);
00687     Float_t cosYCent=centroidY.at(fPrimShow);
00688 
00689     Float_t temp = 1.0-(cosXCent*cosXCent)-(cosYCent*cosYCent);
00690     if(temp < 0) temp = 0;
00691     Float_t cosZCent =TMath::Sqrt(temp);
00692 
00693     //Define histograms
00694     Int_t zBin =21;
00695     Int_t rBin =21;
00696     Float_t zMin =-0.05;
00697     Float_t rMin =-0.05;
00698     Float_t zMax =1.05;
00699     Float_t rMax =1.05;
00700 
00701     TH2F *hRZPrime;
00702     hRZPrime = new TH2F("RZPrime","r vs zprime"
00703                         , zBin, zMin, zMax
00704                         , rBin, rMin, rMax);
00705 
00706     TH2F *hRZ;
00707     hRZ = new TH2F("RZ","r vs z"
00708                    ,zBin, zMin, zMax
00709                    ,rBin, rMin, rMax);
00710 
00711     TH1F *hZPrime;
00712     hZPrime = new TH1F("ZPrime","ELongPrime"
00713                      ,zBin, zMin, zMax);
00714     
00715     TH1F *hRPrime;
00716     hRPrime = new TH1F("RPrime","ETransPrime"
00717                        ,rBin, rMin, rMax);
00718     
00719     TH1F *hZAxis;
00720     hZAxis = new TH1F("ZAxis","ELong"
00721                     ,zBin, zMin, zMax);
00722     
00723     TH1F *hRAxis;
00724     hRAxis = new TH1F("RAxis","ETrans"
00725                       ,rBin, rMin, rMax);
00726 
00727     
00728     //Iterate over all hits in primary cluster and 
00729     //calculate XYZ position centroid
00730     Float_t zPrime;
00731     Float_t distVertexSq;
00732 
00733     Float_t rPrime;
00734     Float_t rAxis;
00735 
00736     Float_t centroX=0;
00737     Float_t centroY=0;
00738     Float_t centroZ=0;
00739 
00740     //Needed to calculate cluster centroid weighted by distance to vertex.
00741     Float_t centroWeightX=0;
00742     Float_t centroWeightY=0;
00743     Float_t centroWeightZ=0;
00744 
00745     Float_t  totalHitEnergy=0.0;
00746 
00747     Float_t epsilon=0.0000001;
00748 
00749     DeqDeqInt_t cluTemp; 
00750     cluTemp.push_back(fClusterMap.at(fPrimShow)); //Select the 
00751                                                   //primary cluster
00752   
00753     for (IterDeqDeqInt_t cluIter = cluTemp.begin()
00754              ;cluIter != cluTemp.end()
00755              ; ++cluIter ){
00756         for (IterDeqInt_t hitIter =  cluIter->begin()
00757                  ;hitIter < cluIter->end()
00758                  ; ++hitIter){
00759             
00760             zPrime=fX.at(*hitIter)*cosXCent
00761                 +fY.at(*hitIter)*cosYCent
00762                 +fZ.at(*hitIter)*cosZCent;
00763             
00764             distVertexSq=fX.at(*hitIter)*fX.at(*hitIter)
00765                 +fY.at(*hitIter)*fY.at(*hitIter)
00766                 +fZ.at(*hitIter)*fZ.at(*hitIter);
00767             
00768             rPrime=TMath::Sqrt(
00769                 TMath::Abs(distVertexSq-(zPrime*zPrime)));
00770             
00771             rAxis=TMath::Sqrt(
00772                 TMath::Abs(distVertexSq-fZ.at(*hitIter)*fZ.at(*hitIter)));
00773 
00774             hRZPrime->Fill(zPrime,rPrime,fE.at(*hitIter));
00775             hRZ     ->Fill(fZ.at(*hitIter),rAxis,fE.at(*hitIter));
00776          
00777             hZPrime ->Fill(zPrime,fE.at(*hitIter));
00778             hRPrime ->Fill(rPrime,fE.at(*hitIter));
00779 
00780             hZAxis  ->Fill(fZ.at(*hitIter),fE.at(*hitIter));
00781             hRAxis  ->Fill(rAxis,fE.at(*hitIter));
00782 
00783             centroX+=fX.at(*hitIter);
00784             centroY+=fY.at(*hitIter);
00785             centroZ+=fZ.at(*hitIter);
00786 
00787             //Total Energy of primary cluster
00788             totalHitEnergy+=(fE.at(*hitIter))*kMeuToGeV;
00789 
00790             MSG("AngClusterAna",Msg::kDebug)<< "(X,Y,Z) = " << "(" 
00791                                             << fX.at(*hitIter) << ","
00792                                             << fY.at(*hitIter) << ","
00793                                             << fZ.at(*hitIter) << ")"
00794                                             << endl;
00795             
00796             MSG("AngClusterAna",Msg::kDebug)<< "centroTemp =" << "("
00797                                             << centroX << ","
00798                                             << centroY << ","
00799                                             << centroZ << ")"
00800                                             << endl;
00801            
00802             centroWeightX+=fX.at(*hitIter)*fX.at(*hitIter);
00803             centroWeightY+=fY.at(*hitIter)*fY.at(*hitIter);
00804             centroWeightZ+=fZ.at(*hitIter)*fZ.at(*hitIter);
00805             
00806 
00807         }//for (hitIter)
00808     }//for(cluIter)
00809 
00810     MSG("AngClusterAna",Msg::kDebug)<< "number of hits = " 
00811                                     << fClusterSize.at(fPrimShow)
00812                                     << endl;
00813         
00814     if (centroX == 0.0) centroX+=epsilon;
00815     if (centroY == 0.0) centroY+=epsilon;
00816     if (centroZ == 0.0) centroZ+=epsilon;
00817 
00818     centroWeightX/=centroX;
00819     centroWeightY/=centroY;
00820     centroWeightZ/=centroZ;
00821 
00822     centroX/=(fClusterSize.at(fPrimShow));
00823     centroY/=(fClusterSize.at(fPrimShow));
00824     centroZ/=(fClusterSize.at(fPrimShow));
00825   
00826     //Fill branch variables
00827     fAngCluster.fACluRmsShwAxis=hRPrime->GetRMS();
00828     fAngCluster.fACluRmsZAxis  =hRAxis->GetRMS();
00829     //     fAngCluster.fACluShwDirX   =centroWeightX;
00830     //     fAngCluster.fACluShwDirY   =centroWeightY;
00831     //     fAngCluster.fACluShwDirZ   =centroWeightZ;
00832 
00833     fAngCluster.fACluPrimEnergy     =totalHitEnergy;
00834     if (fTotalEventHitEnergy > 0.0)
00835         fAngCluster.fACluPrimEnergyRatio=totalHitEnergy/fTotalEventHitEnergy;
00836 
00837     fAngCluster.fACluShwDirX   =centroX;
00838     fAngCluster.fACluShwDirY   =centroY;
00839     fAngCluster.fACluShwDirZ   =centroZ;
00840 
00841     //Fill primary direction vector.
00842     fPrimDir.SetX(centroX);
00843     fPrimDir.SetY(centroY);
00844     fPrimDir.SetZ(centroZ);
00845 
00846     //Print variables
00847     MSG("AngClusterAna",Msg::kInfo)<< "fACluRmsShwAxis = "
00848                                    << fAngCluster.fACluRmsShwAxis
00849                                    << " fACluRmsZAxis = "
00850                                    << fAngCluster.fACluRmsZAxis
00851                                    << " fACluShwDirX = "
00852                                    << fAngCluster.fACluShwDirX
00853                                    << " fACluShwDirY = "
00854                                    << fAngCluster.fACluShwDirY
00855                                    << " fACluShwDirZ = "
00856                                    << fAngCluster.fACluShwDirZ
00857                                    << endl;
00858 
00859     //delete histograms
00860     if(hRZPrime)
00861         delete hRZPrime;
00862 
00863     if(hRZ)
00864         delete hRZ;
00865 
00866      if(hZPrime)
00867          delete hZPrime;
00868 
00869      if(hRPrime)
00870          delete hRPrime;
00871 
00872      if(hZAxis)
00873          delete hZAxis;
00874 
00875      if(hRAxis)
00876          delete hRAxis;
00877 
00878 }
00879 
00880 void AngClusterAna::WeightedEnergy(int evtn, RecRecordImp<RecCandHeader> *srobj){
00881   if(srobj==0){
00882     return;
00883   }
00884   if(((dynamic_cast<NtpStRecord *>(srobj))==0)&&
00885      ((dynamic_cast<NtpSRRecord *>(srobj))==0)){
00886     return;
00887   }
00888 
00889   MSG("ShwfitAna",Msg::kDebug)<<"In ShwfitAna::Analyze"<<endl;
00890   MSG("ShwfitAna",Msg::kDebug)<<"On Snarl "<<srobj->GetHeader().GetSnarl()
00891                               <<" event "<<evtn<<endl;
00892 
00893 //  Reset(srobj->GetHeader().GetSnarl(),evtn);
00894 
00895   NtpSREvent *event = SntpHelpers::GetEvent(evtn,srobj);
00896 
00897   if(!event){
00898       MSG("ShwfitAna",Msg::kError)<<"Couldn't get event "<<evtn
00899                                    <<" from Snarl "<<srobj->GetHeader().GetSnarl()<<endl;
00900       return;
00901    }
00902 
00903   Float_t cent_x = fAngCluster.fACluShwDirX;
00904   Float_t cent_y = fAngCluster.fACluShwDirY;
00905   Float_t cent_z = fAngCluster.fACluShwDirZ;
00906   
00907   Float_t cent_u = TMath::Sqrt(2)*(cent_x-cent_y);
00908   Float_t cent_v = TMath::Sqrt(2)*(cent_x+cent_y);
00909 
00910   Int_t MaxStripU = 0;
00911   Int_t MaxStripV = 0;
00912   Float_t MaxStripU_ph = 0;
00913   Float_t MaxStripV_ph = 0;
00914 
00915   Float_t Weighted_ph0 = 0.0;
00916   Float_t Weighted_ph1 = 0.0;
00917   Float_t Weighted_ph2 = 0.0;
00918   Float_t Weighted_ph3 = 0.0;
00919 
00920   Int_t vtxPlane = event->vtx.plane;
00921   Float_t vtxU = event->vtx.u;
00922   Float_t vtxV = event->vtx.v;
00923 
00924   if(ReleaseType::IsCedar(release)){
00925     NtpStRecord* st = dynamic_cast<NtpStRecord *>(srobj);
00926     NtpVtxFinder vtxf(evtn, st);
00927     if(vtxf.FindVertex() > 0){
00928        vtxPlane = vtxf.VtxPlane();
00929        vtxU = vtxf.VtxU();
00930        vtxV = vtxf.VtxV();
00931     }
00932   }
00933 
00934   for(int i=0;i<event->nstrip;i++){
00935     Int_t index = SntpHelpers::GetStripIndex(i,event);
00936     NtpSRStrip *strip = SntpHelpers::GetStrip(index,srobj);
00937     if(!strip){
00938       continue;
00939     }
00940     if(!evtstp0mip){
00941        MSG("AngCluster",Msg::kError)<<"No mip strip information"<<endl;
00942        continue;
00943     }
00944 
00945 
00946     Float_t Strip_ph_meu = evtstp0mip[index] + evtstp1mip[index];
00947     //(strip->ph0.sigcor+strip->ph1.sigcor)/sigcormeu;
00948 
00949     if(strip->plane>=event->vtx.plane && strip->plane<=event->vtx.plane+(cent_z/0.0354) ){
00950       if(strip->planeview==PlaneView::kU){
00951         if(strip->tpos-event->vtx.u+cent_u*(strip->plane-event->vtx.plane)/(cent_z/0.0354)<.0205){
00952           if( Strip_ph_meu>MaxStripU_ph ){MaxStripU_ph=Strip_ph_meu; MaxStripU = i;}
00953         }
00954       }else if(strip->planeview==PlaneView::kV){
00955         if(strip->tpos-event->vtx.v+cent_v*(strip->plane-event->vtx.plane)/(cent_z/0.0354)<.0205){
00956           if( Strip_ph_meu>MaxStripV_ph ){MaxStripV_ph=Strip_ph_meu; MaxStripV = i;}
00957         }
00958       }
00959     }
00960   }
00961 
00962   Int_t Main_index_u = SntpHelpers::GetStripIndex(MaxStripU,event);
00963   NtpSRStrip *Main_strip_u = SntpHelpers::GetStrip(Main_index_u,srobj);
00964   Int_t Main_index_v = SntpHelpers::GetStripIndex(MaxStripV,event);
00965   NtpSRStrip *Main_strip_v = SntpHelpers::GetStrip(Main_index_v,srobj);
00966 
00967   for(int i=0;i<event->nstrip;i++){
00968     Int_t index = SntpHelpers::GetStripIndex(i,event);
00969     NtpSRStrip *strip = SntpHelpers::GetStrip(index,srobj);
00970     if(!strip){
00971       continue;
00972     }
00973     if(!evtstp0mip){
00974        MSG("MSTCalcAna",Msg::kError)<<"No mip strip information"<<endl;
00975        continue;
00976     }
00977 
00978     float charge = evtstp0mip[index] + evtstp1mip[index];;
00979     //  = (strip->ph0.sigcor+strip->ph1.sigcor)/sigcormeu
00980 
00981 
00982     Float_t Strip_ph_meu0 = charge;
00983     Float_t Strip_ph_meu1 = charge;
00984     Float_t Strip_ph_meu2 = charge;
00985     Float_t Strip_ph_meu3 = charge;
00986 
00987     if(strip->planeview==PlaneView::kU){
00988       if(Main_strip_u->plane-strip->plane>0){
00989         Strip_ph_meu0*=TMath::Sqrt(2*TMath::Abs(Main_strip_u->plane-strip->plane)*0.0354) + 40*TMath::Sqrt(TMath::Abs(Main_strip_u->tpos-strip->tpos));
00990         Strip_ph_meu1*=TMath::Sqrt(2*TMath::Abs(Main_strip_u->plane-strip->plane)*0.0354) + 20*TMath::Abs(Main_strip_u->tpos-strip->tpos);
00991         Strip_ph_meu2*=TMath::Sqrt( (5*TMath::Abs(Main_strip_u->plane-strip->plane)*0.0354) + 40*TMath::Abs(Main_strip_u->tpos-strip->tpos) );
00992         Strip_ph_meu3*=TMath::Sqrt( (2*TMath::Abs(Main_strip_u->plane-strip->plane)*0.0354) + 40*TMath::Abs(Main_strip_u->tpos-strip->tpos) );
00993 
00994       }else{
00995         Strip_ph_meu0*=TMath::Sqrt(TMath::Abs(Main_strip_u->plane-strip->plane)*0.0354) + 40*TMath::Sqrt(TMath::Abs(Main_strip_u->tpos-strip->tpos));
00996         Strip_ph_meu1*=TMath::Sqrt(TMath::Abs(Main_strip_u->plane-strip->plane)*0.0354) + 20*TMath::Abs(Main_strip_u->tpos-strip->tpos);
00997         Strip_ph_meu2*=TMath::Sqrt( (TMath::Abs(Main_strip_u->plane-strip->plane)*0.0354) + 40*TMath::Abs(Main_strip_u->tpos-strip->tpos) );
00998         Strip_ph_meu3*=TMath::Sqrt( (TMath::Abs(Main_strip_u->plane-strip->plane)*0.0354) + 40*TMath::Abs(Main_strip_u->tpos-strip->tpos) );
00999 
01000       }
01001 
01002     }else if(strip->planeview==PlaneView::kV){
01003       if(Main_strip_v->plane-strip->plane>0){
01004         Strip_ph_meu0*=TMath::Sqrt(2*TMath::Abs(Main_strip_v->plane-strip->plane)*0.0354) + 40*TMath::Sqrt(TMath::Abs(Main_strip_v->tpos-strip->tpos));
01005         Strip_ph_meu1*=TMath::Sqrt(2*TMath::Abs(Main_strip_v->plane-strip->plane)*0.0354) + 20*TMath::Abs(Main_strip_v->tpos-strip->tpos);
01006         Strip_ph_meu2*=TMath::Sqrt( (5*TMath::Abs(Main_strip_v->plane-strip->plane)*0.0354) + 40*TMath::Abs(Main_strip_v->tpos-strip->tpos) );
01007         Strip_ph_meu3*=TMath::Sqrt( (2*TMath::Abs(Main_strip_v->plane-strip->plane)*0.0354) + 40*TMath::Abs(Main_strip_v->tpos-strip->tpos) );
01008 
01009       }else{
01010         Strip_ph_meu0*=TMath::Sqrt(TMath::Abs(Main_strip_v->plane-strip->plane)*0.0354) + 40*TMath::Sqrt(TMath::Abs(Main_strip_v->tpos-strip->tpos));
01011         Strip_ph_meu1*=TMath::Sqrt(TMath::Abs(Main_strip_v->plane-strip->plane)*0.0354) + 20*TMath::Abs(Main_strip_v->tpos-strip->tpos);
01012         Strip_ph_meu2*=TMath::Sqrt( (TMath::Abs(Main_strip_v->plane-strip->plane)*0.0354) + 40*TMath::Abs(Main_strip_v->tpos-strip->tpos) );
01013         Strip_ph_meu3*=TMath::Sqrt( (TMath::Abs(Main_strip_v->plane-strip->plane)*0.0354) + 40*TMath::Abs(Main_strip_v->tpos-strip->tpos) );
01014       }
01015 
01016     }
01017 
01018     Weighted_ph0+=Strip_ph_meu0;
01019     Weighted_ph1+=Strip_ph_meu1;
01020     Weighted_ph2+=Strip_ph_meu2;
01021     Weighted_ph3+=Strip_ph_meu3;
01022   }
01023 
01024   if(MaxStripU==0 && MaxStripV==0){Weighted_ph0=-100;Weighted_ph1=-100;Weighted_ph2=-100;Weighted_ph3=-100;}
01025   fAngCluster.weightedPH0 = Weighted_ph0;
01026   fAngCluster.weightedPH1 = Weighted_ph1;
01027   fAngCluster.weightedPH2 = Weighted_ph2;
01028   fAngCluster.weightedPH3 = Weighted_ph3;
01029 
01030 }
01031 
01032 
01033 void AngClusterAna::Reset()
01034 {
01035 
01036 
01037 }