Changeset 251 in Idarraga


Ignore:
Timestamp:
Nov 16, 2011, 4:42:27 PM (13 years ago)
Author:
idarraga
Message:
 
Location:
mafalda
Files:
15 edited

Legend:

Unmodified
Added
Removed

  • mafalda/AnalysisCore/MediPixAnalysisCore.h

    r222 r251  
    273273
    274274        fChain->SetBranchAddress("m_frameXC", &m_frameXC, &b_FramesData_m_frameXC);
    275         //fChain->SetBranchAddress("m_lvl1", &m_lvl1, &b_FramesData_m_lvl1);
     275        fChain->SetBranchAddress("m_lvl1", &m_lvl1, &b_FramesData_m_lvl1);
    276276        fChain->SetBranchAddress("m_nEntriesPad", &m_nEntriesPad, &b_FramesData_m_nEntriesPad);
    277277        fChain->SetBranchAddress("m_nHitsInPad", &m_nHitsInPad, &b_FramesData_m_nHitsInPad);

  • mafalda/BlobsFinder/BlobsFinder.cpp

    r222 r251  
    598598        Float_t balanceToMin = 0., balanceToMax = 0.,
    599599                        slopeGuess = 0., intersectionyGuess = 0.,
    600                         chisquare_OverDof = 0.;
     600                        chisquare_OverDof = 0., chisquare = 0., ndf = 0.;
    601601        Int_t minQId = 0, maxQId = 0;
    602602        GetBalance(min_x, max_x, min_y, max_y,
    … …  
    612612        Float_t fitCut = 0.;
    613613
    614         MAFTools::LinearRegression(&fitSlope, &fitCut, &chisquare_OverDof
    615                         ,slopeGuess, intersectionyGuess
    616                         ,oneBlob.GetContentSet());
     614        MAFTools::LinearRegression(&fitSlope, &fitCut, &chisquare_OverDof, &chisquare, &ndf,
     615                        slopeGuess, intersectionyGuess,
     616                        oneBlob.GetContentSet());
    617617
    618618        oneBlob.bP.chisquare_OverDof = chisquare_OverDof;
     619        oneBlob.bP.chisquare = chisquare;
     620        oneBlob.bP.NDF = ndf;
    619621        oneBlob.bP.rotAngle = -1.*TMath::ATan(fitSlope);
    620622        oneBlob.bP.fitSlope = fitSlope;

  • mafalda/BlobsFinder/BlobsFinder.h

    r222 r251  
    8888  Float_t rotAngle;
    8989  Float_t chisquare_OverDof;
     90  Float_t chisquare;
     91  Float_t NDF;
    9092  Float_t fitSlope;
    9193  Float_t fitCut;

  • mafalda/MAFTools/MAFTools.cpp

    r247 r251  
    2929
    3030Bool_t LinearRegression(Float_t * slope, Float_t * cut,
    31                 Float_t * chisquare_OverDof,
     31                Float_t * chisquare_OverDof, Float_t * chisquare, Float_t * ndf,
    3232                Float_t initslope, Float_t initycross,
    3333                set< pair<Int_t, Int_t> > blobContent){
    … …  
    6666
    6767        *chisquare_OverDof = f1->GetChisquare()/(Float_t)f1->GetNDF();
     68        *chisquare = f1->GetChisquare();
     69        *ndf = (Float_t)f1->GetNDF();
    6870        *slope = f1->GetParameter(0);
    6971        *cut = f1->GetParameter(1);

  • mafalda/MAFTools/MAFTools.h

    r247 r251  
    2424 
    2525  Bool_t LinearRegression(Float_t * slope, Float_t * cut,
    26                           Float_t * chisquare_OverDof,
     26                          Float_t * chisquare_OverDof, Float_t * chisquare, Float_t * ndf,
    2727                          Float_t initslope, Float_t initycross,
    2828                          std::set< std::pair<Int_t, Int_t> > blobContent);

  • mafalda/MPXAlgo/MediPixAlgo.cpp

    r240 r251  
    439439//};
    440440
     441// In the future I should write an error handler --> TODO
     442void MediPixAlgo::__FATALMC_ERROR(){
     443        cout << "[FATAL] Request of MC related information on real data ... giving up" << endl;
     444        exit(1);
     445}
     446
    441447Int_t MediPixAlgo::GetPrimaryMCVertex_N(){
     448        if( !IsMCData() ) { __FATALMC_ERROR(); }
    442449        return m_myManager->GetPrimaryMCVertex_N();
    443450}
    444451Double_t MediPixAlgo::GetPrimaryMCVertex_X(int i){
     452        if( !IsMCData() ) { __FATALMC_ERROR(); }
    445453        return m_myManager->GetPrimaryMCVertex_X(i);
    446454}
    447455Double_t MediPixAlgo::GetPrimaryMCVertex_Y(int i){
     456        if( !IsMCData() ) { __FATALMC_ERROR(); }
    448457        return m_myManager->GetPrimaryMCVertex_Y(i);
    449458}
    450459Double_t MediPixAlgo::GetPrimaryMCVertex_Z(int i){
     460        if( !IsMCData() ) { __FATALMC_ERROR(); }
    451461        return m_myManager->GetPrimaryMCVertex_Z(i);
    452462}

  • mafalda/MPXAlgo/MediPixAlgo.h

    r240 r251  
    128128  Double_t GetPrimaryMCVertex_Z(int);
    129129
     130  // error handling --> TODO
     131  void __FATALMC_ERROR();
     132
    130133  /* ****************************************** */
    131134

  • mafalda/OctoCEA/OctoCEA.cpp

    r248 r251  
    1111#include "MAFTools.h"
    1212
     13#include <map>
     14#include <vector>
     15
    1316using namespace MSG;
    1417
    … …  
    1922        m_joinSlopeAngleTolerance = 5; // deg
    2023        m_joinCutPixelTolerance = 20; // pixels
     24        m_minNPixels = 5;
     25        m_minNPixelsClouds = 300;
    2126}
    2227
    … …  
    2732        // This value will be overridden by the configuration since it'll be set up
    2833        //  a few lines below as a configuration value
    29         m_minNPixels = 5;
    3034
    3135        // You will get an ntuple file containing a TTree with the name of this
    … …  
    3640        getMyTree()->Branch("alphaAngle", &m_alphaAngle);
    3741        getMyTree()->Branch("chi2", &m_chi2);
     42        getMyTree()->Branch("chi2Prob", &m_chi2Prob);
     43        getMyTree()->Branch("typeOfRecoFrame",&m_typeOfRecoFrame);
     44
     45        getMyTree()->Branch("assemblyTime", &m_assemblyTime);
     46        getMyTree()->Branch("frameId", &m_frameId);
    3847
    3948        // A configuration value that can be tuned from the Viewer
    4049        RegisterConfigurationValue(&m_minNPixels, "minNPixels");
     50        RegisterConfigurationValue(&m_minNPixelsClouds, "minNPixelsClouds");
    4151        RegisterConfigurationValue(&m_minClusterSize, "minClusterSize");
    4252        RegisterConfigurationValue(&m_joinSlopeAngleTolerance, "joinSlopeAngleTolerance");
    … …  
    4656
    4757void OctoCEA::Execute(){
     58
     59        m_frameId = GetFrameId();
     60        m_typeOfRecoFrame = __RECO_UNDEFINED;
    4861
    4962        // Ask the store gate if the previous algorithm (BlobsFinder --> reponsible for clustering)
    … …  
    6578        vector<Float_t> slopesV;
    6679        vector<Float_t> cutsV;
    67 
     80        vector<Int_t> blobIndexV;
     81
     82        int blobIndex = -1;
    6883        for( ; blobsItr != blobsVector.end() ; blobsItr++)
    6984        {
     85                blobIndex++;
    7086
    7187                // Limit all this to clusters with a minimum size.
    … …  
    7894                if ( (*blobsItr).GetBlobType() >= _MIP ) {
    7995
    80                         vector< pair<Int_t, Int_t> > limitPixels = MAFTools::FindLimitPixels(*blobsItr);
     96                        vector< pair<Float_t, Float_t> > limitPixels = MAFTools::FindLimitPixels(*blobsItr);
    8197
    8298                        float cut = (*blobsItr).bP.fitCut;
    … …  
    86102                        slopesV.push_back(slope);
    87103                        cutsV.push_back(cut);
    88 
    89                         float cx = (*blobsItr).bP.geoCenter_x;
    90                         float hwidthx = (*blobsItr).bP.width_x/10.;
    91                         float alpha = (*blobsItr).bP.rotAngle;
     104                        blobIndexV.push_back(blobIndex);
     105
     106                        //float cx = (*blobsItr).bP.geoCenter_x;
     107                        //float hwidthx = (*blobsItr).bP.width_x/10.;
     108                        //float alpha = (*blobsItr).bP.rotAngle;
    92109
    93110                        //alpha = fabs(alpha);
    … …  
    108125
    109126
    110                         m_clusterSize.push_back(  (*blobsItr).bP.nPixels  );
    111                         m_alphaAngle.push_back(  alpha*180./TMath::Pi()  );
    112                         m_chi2.push_back(  (*blobsItr).bP.chisquare_OverDof  );
    113 
    114 
    115                 }
    116 
    117         }
    118 
     127
     128
     129
     130                }
     131
     132        }
     133
     134
     135        // Map index of the segment with index in the blobs container
     136        map<int, int> matchAssemblyIndexesMap;
    119137
    120138        // Study a possible mip assembly
    … …  
    125143                vector<Float_t>::iterator itrS = slopesV.begin();
    126144                vector<Float_t>::iterator itrC = cutsV.begin();
    127 
    128                 Float_t firstS = *itrS;
    129                 Float_t firstC = *itrC;
    130                 Log << MSG::INFO << " leading --> cut : " << firstC << " | slope : " << firstS << endreq;
    131 
    132                 itrS++; itrC++; // compare starting from the second
     145                vector<int>::iterator itrIndex = blobIndexV.begin();
     146                vector<int>::iterator itrIndex2 = blobIndexV.begin();
     147
     148                Log << MSG::INFO << "List of fragments " << endreq;
     149                int cntr = 0;
    133150                for ( ; itrS != slopesV.end() ; ) {
    134 
    135                         Log << MSG::INFO << " cut : " << *itrC << " | slope : " << *itrS << endreq;
    136 
    137                         if( firstC <= (*itrC) - m_joinCutPixelTolerance &&
    138                                         firstC >= (*itrC) - m_joinCutPixelTolerance ) {
    139                                 Log << MSG::INFO << "   within tolerance " << endreq;
    140                         }
    141 
     151                        Log << MSG::INFO << "[" << cntr << "] " << "cut : " << *itrC << " | slope : " << *itrS;
     152                        Log << MSG::INFO << " | container index : " << *itrIndex2 << endreq;
    142153                        itrS++;
    143154                        itrC++;
    144                 }
    145 
     155                        cntr++;
     156                        itrIndex2++;
     157                }
     158
     159                // Search of a trend.  Rewind iterators.
     160                vector<Float_t>::iterator itrC2 = cutsV.begin();
     161                cntr = 0;
     162                itrIndex2 = blobIndexV.begin();
     163                int cntr2 = 0, matchesCntr = 0;
     164                for (itrC = cutsV.begin() ; itrC != cutsV.end() ; itrC++) {
     165                        cntr2=0;
     166                        itrIndex2 = blobIndexV.begin();
     167                        for (itrC2 = cutsV.begin() ; itrC2 != cutsV.end() ; itrC2++) {
     168                                if(itrC == itrC2) {cntr2++; itrIndex2++; continue;} // don't compare same element
     169                                Log << MSG::INFO << "[" << cntr << "] --> " << "[" << cntr2 << "] ";
     170                                if ( (*itrC2) >= (*itrC) - m_joinCutPixelTolerance &&
     171                                                (*itrC2) <= (*itrC) + m_joinCutPixelTolerance ) {
     172                                        Log << MSG::INFO << " within tolerance ";
     173                                        matchesCntr++;
     174                                        matchAssemblyIndexesMap[cntr] = (*itrIndex);
     175                                        matchAssemblyIndexesMap[cntr2] = (*itrIndex2);
     176                                }
     177                                Log << MSG::INFO << endreq;
     178                                cntr2++;
     179                                itrIndex2++;
     180                        }
     181                        if (matchesCntr > 1) break;
     182                        matchesCntr = 0; // otherwise rewind
     183                        cntr++;
     184                        itrIndex++;
     185                        Log << MSG::INFO << "--------------------" << endreq;
     186                }
     187
     188                if(matchesCntr > 1) {
     189                        Log << MSG::INFO << "Found mip assembly" << endreq;
     190                        map<int, int>::iterator mItr = matchAssemblyIndexesMap.begin();
     191                        for ( ; mItr != matchAssemblyIndexesMap.end() ; mItr++) {
     192                                Log << MSG::INFO << "[" << (*mItr).first << "] -- in container --> " << (*mItr).second << endreq;
     193                        }
     194                }
     195
     196        }
     197
     198        // if an assembly was found
     199        if ( !matchAssemblyIndexesMap.empty() ) {
     200                map<int, int>::iterator mItr = matchAssemblyIndexesMap.begin();
     201                int totalClusterSize = 0;
     202                for ( ; mItr != matchAssemblyIndexesMap.end() ; mItr++) {
     203                        // Mip segment index
     204                        int msI = (*mItr).second;
     205                        blob ms = blobsVector[msI];
     206                        set< pair<Int_t, Int_t> > bc = ms.GetContentSet();
     207                        set< pair<Int_t, Int_t> >::iterator bcI = bc.begin();
     208
     209                        for( ; bcI != bc.end() ; bcI++) {
     210                                //Log << MSG::INFO << GetMatrixElement(*bcI) << endreq;
     211                                m_assemblyTime.push_back( GetMatrixElement(*bcI) ); // saving time
     212                        }
     213                        // characteristics of the assembly
     214                        float alpha = ms.bP.rotAngle;
     215                        totalClusterSize += ms.bP.nPixels;
     216                        m_alphaAngle.push_back(  alpha*180./TMath::Pi()  );
     217                        m_chi2.push_back(  ms.bP.chisquare_OverDof  );
     218                        m_chi2Prob.push_back( TMath::Prob(ms.bP.chisquare, ms.bP.NDF));
     219                }
     220                m_clusterSize.push_back(  totalClusterSize  );
     221
     222                m_typeOfRecoFrame = __RECO_MIP;
     223        }
     224        else // fill it for all in the case of big clusters
     225        {
     226                // look at all blobs here
     227                for(blobsItr = blobsVector.begin() ; blobsItr != blobsVector.end() ; blobsItr++)
     228                {
     229                        set< pair<Int_t, Int_t> > bc = (*blobsItr).GetContentSet();
     230                        if((int)bc.size() >= m_minNPixelsClouds){
     231                                set< pair<Int_t, Int_t> >::iterator bcI = bc.begin();
     232                                for( ; bcI != bc.end() ; bcI++){
     233                                        //Log << MSG::INFO << GetMatrixElement(*bcI) << endreq;
     234                                        m_assemblyTime.push_back( GetMatrixElement(*bcI) ); // saving time
     235                                }
     236
     237                        }
     238                }
     239                m_typeOfRecoFrame = __RECO_CLOUD;
    146240        }
    147241
    … …  
    154248        m_alphaAngle.clear();
    155249        m_chi2.clear();
     250        m_chi2Prob.clear();
     251        m_assemblyTime.clear();
    156252
    157253}
    … …  
    169265 */
    170266OctoClusterContainer::OctoClusterContainer(MediPixAlgo * algo) :
    171 CandidateContainer(algo) {
    172 }
    173 
    174 void OctoClusterContainer::SetBlobType(Int_t id, octoclustertype type){
     267                                                                                CandidateContainer(algo) {
     268}
     269
     270void OctoClusterContainer::SetBlobType(Int_t /*id*/, octoclustertype /*type*/){
    175271        //allBlobs[id].SetBlobType(type);
    176272}

  • mafalda/OctoCEA/OctoCEA.h

    r248 r251  
    1010#define __HEAD_INDX 0
    1111#define __TAIL_INDX 1
     12
     13#define __RECO_UNDEFINED    0
     14#define __RECO_MIP    5
     15#define __RECO_CLOUD 10
    1216
    1317class OctoCEA : public MediPixAlgo {
    … …  
    3135        // Config
    3236        Int_t m_minNPixels;
     37        Int_t m_minNPixelsClouds;
    3338        Int_t m_minClusterSize;
    3439        Float_t m_joinCutPixelTolerance;
    … …  
    3944        vector<Float_t> m_alphaAngle;
    4045        vector<Float_t> m_chi2;
     46        vector<Float_t> m_chi2Prob;
     47
    4148        vector<Int_t> m_clusterSize;
     49        Int_t m_typeOfRecoFrame;
     50
     51        vector<int> m_assemblyTime;
     52        Int_t m_frameId;
    4253
    4354        ClassDef(OctoCEA, 1)

  • mafalda/ShallowAngleFEIX/ShallowAngleFEIX.cpp

    r249 r251  
    1717ShallowAngleFEIX::ShallowAngleFEIX()
    1818:
     19m_minNPixels ( 5 ),
     20m_guardDistanceX ( 1 ),
     21m_guardDistanceY ( 5 ),
    1922m_pixSizeX ( 0.40 ), // FEI3 [mm]
    2023m_pixSizeY ( 0.05 ), // FEI3 [mm]
     24m_dRayBalSearch ( 1.4 ),
    2125m_minXWidth ( 6 ),
    22 m_guardDistanceX ( 1 ),
    23 m_dRayBalSearch ( 1.4 ),
    24 m_guardDistanceY ( 5 ),
    25 m_minNPixels ( 5 )
     26m_clCut ( 0.63 )
    2627{
    2728
    … …  
    5657        RegisterConfigurationValue(&m_pixSizeY, "pixelSizeY");
    5758        RegisterConfigurationValue(&m_dRayBalSearch, "dRayBalSearch");
     59        RegisterConfigurationValue(&m_clCut, "CLCut");
     60
     61        // Load clean map TH2F from file previously processed by WeightTOTMap
     62        TString mapFileName = "ShallowAngleFEIX/MAFOutput_WeightTOTMap.root";
     63        m_cleanMapFile = new TFile(mapFileName, "OPEN");
     64        if(!m_cleanMapFile->IsOpen()){
     65                Log << MSG::ERROR << "Error opening good pixels (C.L.) map --> " << mapFileName << endreq;
     66                Log << MSG::ERROR << "giving up ..." << endreq;
     67                exit(1);
     68        }
     69        m_cleanMapFile->cd();
     70        m_cleanMap = static_cast<TH2F*>( m_cleanMapFile->Get("cutmap") );
     71        m_eqMap = static_cast<TH2F*>( m_cleanMapFile->Get("eqmap") );
     72        m_eqhistofunc = static_cast<TF1*>( m_cleanMapFile->Get( "eqhistofunc" ) );
     73
     74        if( ! m_cleanMap || ! m_eqMap || ! m_eqhistofunc ) {
     75                Log << MSG::ERROR << "Error opening good pixels histogram.  Can't find one of more of the objects with name: "
     76                                << "cutmap, eqmap, eqhistofunc" << endreq;
     77                Log << MSG::ERROR << "giving up ..." << endreq;
     78                exit(1);
     79        }
     80
     81        if ( Log.OutputLevel == MSG::DEBUG ) {
     82                for(int x = 0 ; x < m_cleanMap->GetNbinsX() ; x++){
     83                        for(int y = 0 ; y < m_cleanMap->GetNbinsY() ; y++){
     84                                Log << MSG::DEBUG << TString::Format("%.1f ",m_cleanMap->GetBinContent(x, y));
     85                        }
     86                        Log << MSG::DEBUG << endreq;
     87                }
     88        }
     89        Log << MSG::INFO << "Eq mean = " << m_eqhistofunc->GetParameter(1) << endreq;
     90
     91        Log << MSG::INFO << "map loaded ..." << endreq;
     92
     93        // Performing any other operation on the algorithm's output file from now on
     94        getMyROOTFile()->cd();
    5895
    5996}
    … …  
    76113        vector<blob>::iterator blobsItr = blobsVector.begin(); //allBlobs.begin();
    77114
    78 
    79115        for( ; blobsItr != blobsVector.end() ; blobsItr++)
    80116        {
    … …  
    93129
    94130                if(bT == _MIP) { // only Mips
     131
     132                        // initialize the mip goodness
     133                        m_mipGoodness = 0;
    95134
    96135                        // info about the mip to store
    … …  
    182221                                                << GetMatrixElement(*csItr) << endreq;
    183222                        }
    184                         */
     223                         */
    185224
    186225                        m_clusterHeight = (*blobsItr).GetBlobProperties().width_y;
    … …  
    188227                        int cntr = 0;
    189228
    190                         for(int xi = 0 ; xi <= 30 ; xi++){
     229                        // Fill up with a few zeros to the maximum track size.
     230                        // All vectors should end up by at least one entry with value 0
     231                        for (int xi = 0 ; xi <= 30 ; xi++) {
    191232                                m_totPerSegment.push_back(0.); // add segment initialized at 0
    192233                                m_chargeWeights.push_back(0.); // add segment
    193234                        }
    194235
    195                         for(int x = initx ; x <= endx ; x++){
    196 
    197                                 m_totPerSegment.push_back(0); // add segment initialized at 0
    198                                 m_chargeWeights.push_back(0); // add segment
     236                        for (int x = initx ; x <= endx ; x++) {
     237
     238                                m_totPerSegment.push_back(0); // add TOT entry initialized at 0
     239                                m_chargeWeights.push_back(0); //
    199240
    200241                                for (int y = inity ; y <= endy ; y++) {
    201242
    202                                         Log << MSG::DEBUG << "  --> " << x << " , " << y << " : " << GetMatrixElement(x,y) << endreq;
     243                                        if( m_cleanMap->GetBinContent(x, y) < m_clCut ) { // see if it passes the clCut
     244                                                // in this case don't consider this mip
     245                                                m_mipGoodness = 0;
     246                                        } else {
     247                                                m_mipGoodness = 1;
     248                                        }
     249
     250                                        double correctedTOT = GetWeightedTOT(x, y);
     251                                        Log << MSG::DEBUG << "  --> " << x << " , " << y << " : " << GetMatrixElement(x,y)
     252                                                        << " | corrected : " << correctedTOT << endreq;
     253
     254                                        //m_totPerSegment[cntr] += correctedTOT;
     255                                        //m_chargeWeights[cntr] += correctedTOT; // normalized later
    203256                                        m_totPerSegment[cntr] += GetMatrixElement(x,y);
    204257                                        m_chargeWeights[cntr] += GetMatrixElement(x,y); // normalized later
    … …  
    207260                                cntr++;
    208261                        }
     262
     263                        // analyze totPerSegment and chargeWeights see if some alignement is needed
     264                        // Compare the first entry with the mean of the rest of the track
     265                        int max = m_totPerSegment.size(), nent = 0;
     266                        float firsttot = m_totPerSegment[0], means = 0.;
     267                        // calc mean for the rest
     268                        for(cntr = 1 ; cntr < max ; cntr++) {
     269                                if(m_totPerSegment[cntr] == 0.) break; // end of track
     270                                means += m_totPerSegment[cntr];
     271                                nent++;
     272                        }
     273                        means /= nent;
     274
     275                        /*
     276                        if(firsttot > means*0.7) { // in this case the tracks needs shift
     277                                m_totPerSegment.insert(m_totPerSegment.begin(), 0.);
     278                                m_chargeWeights.insert(m_chargeWeights.begin(), 0.);
     279                        }
     280                        */
    209281
    210282                        // normalize m_chargeWeights to total charge
    … …  
    220292
    221293                                vector<Float_t>::iterator itrW = m_chargeWeights.begin();
    222                                 vector<Int_t>::iterator itrT = m_totPerSegment.begin();
     294                                vector<Float_t>::iterator itrT = m_totPerSegment.begin();
    223295
    224296                                for( ; itrW != m_chargeWeights.end() ; ) {
    … …  
    245317
    246318                        // Fill the output tree of this algorithm
    247                         getMyTree()->Fill();
    248 
     319                        if(m_mipGoodness) {
     320                                getMyTree()->Fill();
     321                        }
     322                        //getMyTree()->Fill();
    249323                        // WARNING ! don't forget to clean up your variables for the next TTree::Fill call
    250324                        m_distanceToCenter.clear();
    … …  
    374448}
    375449
     450bool ShallowAngleFEIX::IsGoodPixel_CheckCL(pair<int, int> p, double cutCL){
     451        return IsGoodPixel_CheckCL(p.first, p.second, cutCL);
     452}
     453
     454bool ShallowAngleFEIX::IsGoodPixel_CheckCL(int x, int y, double){
     455
     456        if(m_cleanMap->GetBinContent(x, y) > 0.) {
     457                return true;
     458        }
     459        return false;
     460}
     461
     462double ShallowAngleFEIX::GetWeightedTOT(int x, int y){
     463
     464        // Use the fit function  m_eqhistofunc
     465        // and the eqmap to find the weight value m_eqMap;
     466        double binContent = m_eqMap->GetBinContent(x,y);
     467        if(binContent == 0) // not present in the equalization map
     468                return 0.;
     469
     470        double totval = (double)GetMatrixElement(x,y);
     471
     472        if (totval == 0.) // in case there is no hit in this location
     473                return 0.;
     474
     475        double mean = m_eqhistofunc->GetParameter(1); // gaus mean
     476        double fraction = binContent/mean;
     477
     478        return totval/fraction; // normalize to mean
     479}
     480
    376481#endif

  • mafalda/ShallowAngleFEIX/ShallowAngleFEIX.h

    r249 r251  
    3232        Int_t GuessNumberOfDeltaRaysFromUnbalance(Int_t &, Int_t &);
    3333        int CheckDiscontinuity(blob, int);
     34        bool IsGoodPixel_CheckCL(int, int, double);
     35        bool IsGoodPixel_CheckCL(pair<int, int>, double);
     36        double GetWeightedTOT(int, int);
    3437
    3538private:
    … …  
    4346        Float_t m_dRayBalSearch;
    4447        Int_t m_minXWidth;
     48        Double_t m_clCut;
     49
     50        // tests
     51        Int_t m_mipGoodness;
    4552
    4653        // for output
    … …  
    5966        AllBlobsContainer * m_reclusteringContainer;
    6067
    61         vector<Int_t> m_totPerSegment; // Total TOT per segment/pixel
     68        vector<Float_t> m_totPerSegment; // Total TOT per segment/pixel
    6269        vector<Float_t> m_chargeWeights; // Charge weight, normalized to the total TOT of the track
     70
     71        // clean map handling. Recover from WeightTOTMap run
     72        TFile * m_cleanMapFile;
     73        TH2F * m_cleanMap;
     74        TH2F * m_eqMap;
     75        TF1 * m_eqhistofunc;
    6376
    6477        ClassDef(ShallowAngleFEIX, 1)

  • mafalda/ShallowAngleFEIX/configuration_algos_ShallowAngleFEIX.txt

    r240 r251  
    1111PRBasicSpecies nInnerPixels 2 int
    1212PRBasicSpecies longGammaMax 3 int
    13 PRBasicSpecies minNPixelsMip 5 int
     13PRBasicSpecies minNPixelsMip 3 int
    1414PRBasicSpecies maxNPixelsCurly 50 int
    15 ShallowAngleFEIX minNPixels 5 int
    16 ShallowAngleFEIX minXWidth 6 int
     15ShallowAngleFEIX minNPixels 3 int
     16ShallowAngleFEIX minXWidth 3 int
    1717ShallowAngleFEIX guardDistanceX 1 int
    1818ShallowAngleFEIX guardDistanceY 5 int

  • mafalda/ShallowAngleFEIX/runAngle.C

    r240 r251  
    1414}
    1515
     16#define __THICKNESS_PLANAR 0.200 // [mm]
     17#define __THICKNESS_3D     0.230 // [mm]
     18
    1619void runAngle(){
    1720
    18   gSystem->Load("/home/idarraga/analysis/root_tools/multipleplots/plots_treeTools.so");
    19   gROOT->ProcessLine(".x ~/styles/AtlasStyle.C");
     21        gSystem->Load("/home/idarraga/analysis/root_tools/multipleplots/plots_treeTools.so");
     22        gROOT->ProcessLine(".x ~/styles/AtlasStyle.C");
     23
     24        GetTree("/home/idarraga/analysis/mafalda/MAFOutput_MPXNtuple_USBPIXI4_22_61548.root",
     25                        "ShallowAngleFEIX", "61548 Planar-LUB4 1000V", 1.0);
     26
     27//      GetTree("/home/idarraga/analysis/mafalda/MAFOutput_MPXNtuple_USBPIXI4_22_61564.root",
     28        //              "ShallowAngleFEIX", "61564 800V", 1.0);
     29
     30//      GetTree("/home/idarraga/analysis/mafalda/MAFOutput_MPXNtuple_USBPIXI4_22_61566.root",
     31        //              "ShallowAngleFEIX", "61566 400V", 1.0);
     32
     33//      GetTree("/home/idarraga/analysis/mafalda/MAFOutput_MPXNtuple_USBPIXI4_22_61569.root",
     34        //              "ShallowAngleFEIX", "61569 600V", 1.0);
     35
     36        GetTree("/home/idarraga/analysis/mafalda/MAFOutput_MPXNtuple_USBPIXI4_22_61573.root",
     37                        "ShallowAngleFEIX", "61573 200V", 1.0);
    2038
    2139
    22   GetTree("/home/idarraga/analysis/mafalda/MAFOutput_MPXNtuple_USBPIXI4_23_61417.root",
    23           "ShallowAngleFEIX", "61417 3D-SCC87 p-5E15", 1.0);
     40        //GetTree("/storage1/idarraga/Testbeam_September_CERN_2011/KEK_20V_200V_4deg/MAFOutput/MAFOutput_MPXNtuple_USBPIXI4_21.root",
     41        //  "ShallowAngleFEIX", "DUT2", 1.0);
    2442
    25     GetTree("/home/idarraga/analysis/mafalda/MAFOutput_MPXNtuple_USBPIXI4_22_61548.root",
    26           "ShallowAngleFEIX", "61548 Planar-LUB4 n-5E15", 1.0);
     43        //GetTree("/storage1/idarraga/Testbeam_September_CERN_2011/KEK_20V_200V_4deg/MAFOutput/MAFOutput_MPXNtuple_USBPIXI4_22.root",
     44        //        "ShallowAngleFEIX", "DUT3", 1.0);
    2745
    28   //GetTree("/storage1/idarraga/Testbeam_September_CERN_2011/KEK_20V_200V_4deg/MAFOutput/MAFOutput_MPXNtuple_USBPIXI4_21.root",
    29 //  "ShallowAngleFEIX", "DUT2", 1.0);
    30 
    31   //GetTree("/storage1/idarraga/Testbeam_September_CERN_2011/KEK_20V_200V_4deg/MAFOutput/MAFOutput_MPXNtuple_USBPIXI4_22.root",
    32   //      "ShallowAngleFEIX", "DUT3", 1.0);
    33 
    34   /*
     46        /*
    3547  GetTree("/storage1/idarraga/Testbeam_July_CERN_2011/KEK_100V_4deg/MAFOutput/MAFOutput_MPXNtuple_USBPIXI4_22_100V.root",
    3648          "ShallowAngleFEIX", "DUT2", 1.0);
    … …  
    3850  GetTree("/storage1/idarraga/Testbeam_July_CERN_2011/KEK_100V_4deg/MAFOutput/MAFOutput_MPXNtuple_USBPIXI4_23_100V.root",
    3951          "ShallowAngleFEIX", "DUT3", 1.0);
    40 */
     52         */
    4153
    42   /////////////////////////////////
    43   // config
    44   OnCanvasOptions opts;
    45   opts.doLegend = true;
    46   opts.constructorStr = "";
    47   opts.drawOpt = "";
    48   opts.markerStyle = 3;
    49   opts.setLogY = false;
    50   opts.drawNormalized = false;
     54        /////////////////////////////////
     55        // config
     56        OnCanvasOptions opts;
     57        opts.doLegend = true;
     58        opts.constructorStr = "";
     59        opts.drawOpt = "";
     60        opts.markerStyle = 3;
     61        opts.setLogY = false;
     62        opts.drawNormalized = true;
    5163
    52   ActionFlags afl;
    53   afl.doIntegral = false;
     64        ActionFlags afl;
     65        afl.doIntegral = false;
    5466
    55   //////////////////////////////////
    56   // changes on style
    57   TStyle * theStyle = (TStyle *)gROOT->FindObject("ATLAS");
    58   float tsize = 0.07;
    59   theStyle->SetTextSize(tsize);
    60   theStyle->SetLabelSize(tsize,"x");
    61   theStyle->SetTitleSize(tsize,"x");
    62   theStyle->SetLabelSize(tsize,"y");
    63   theStyle->SetTitleSize(tsize,"y");
    64   theStyle->SetLabelSize(tsize,"z");
    65   theStyle->SetTitleSize(tsize,"z");
    66   ////////////////////////////////////////////////////
    67  
    68   TCanvas * c11 = new TCanvas("c11", "track multiplicity");
     67        //////////////////////////////////
     68        // changes on style
     69        TStyle * theStyle = (TStyle *)gROOT->FindObject("ATLAS");
     70        float tsize = 0.07;
     71        theStyle->SetTextSize(tsize);
     72        theStyle->SetLabelSize(tsize,"x");
     73        theStyle->SetTitleSize(tsize,"x");
     74        theStyle->SetLabelSize(tsize,"y");
     75        theStyle->SetTitleSize(tsize,"y");
     76        theStyle->SetLabelSize(tsize,"z");
     77        theStyle->SetTitleSize(tsize,"z");
     78        ////////////////////////////////////////////////////
    6979
    70   c11->cd();
    71   opts.doLegend = true;
    72   opts.constructorStr = "(40,0,5)";
    73   TH1 * h_first = (TH1 *) DrawAll("mipLength_X", "clusterHeight < 3 && TMath::Abs(mipThetaAngle) < 0.17", opts, afl);
    74   //h_first->GetXaxis()->SetRangeUser(0, 10);
    75   h_first->GetXaxis()->SetTitle("mip length [mm]");
     80        TCanvas * c11 = new TCanvas("c11", "track multiplicity");
    7681
     82        c11->cd();
     83        opts.doLegend = true;
     84        opts.constructorStr = "(40, 0.5, 3.0)";
     85        TH1 * h_first = (TH1 *) DrawAll("mipLength_X", "clusterHeight < 3 && TMath::Abs(mipThetaAngle) < 0.17", opts, afl);
     86        //h_first->GetXaxis()->SetRangeUser(0, 10);
     87        h_first->GetXaxis()->SetTitle("mip length [mm]");
     88        h_first->GetYaxis()->SetRangeUser(0, 1);
     89
     90        // angle function --> The limits here have to mach
     91        TF1 * f1p = new TF1("anglefPlanar", getangle, 1, 10, 1);
     92        f1p->SetParameter(0, __THICKNESS_PLANAR); // detector thickness
     93        f1p->SetRange(f1p->Eval(3), f1p->Eval(1));
     94
     95        //Double_t maxa = h_first->GetMaximum()*0.95;
     96        Double_t maxa = 0.95;
     97        TGaxis * ax2p = new TGaxis(1, maxa, 3, maxa, "anglefPlanar", 510, "-");
     98        ax2p->SetTitle("shallow angle planar [deg]");
     99        //ax2->SetLabelColor(kRed);
     100        ax2p->Draw("");
     101        ax2p->SetLabelColor(kRed);
     102        ax2p->SetLineColor(kRed);
     103        ax2p->SetTitleColor(kRed);
     104
     105
     106        /*
    77107  // angle function --> The limits here have to mach
    78   TF1 * f1 = new TF1("anglef", getangle, 1.71836, 8.53077, 1);
    79   f1->SetParameter(0, 0.150); // detector thickness
    80   Double_t maxa = h_first->GetMaximum()*0.95;
    81   TGaxis * ax2 = new TGaxis(1, maxa, 5, maxa, "anglef", 510, "-");
    82   ax2->SetTitle("shallow angle [deg]");
     108  TF1 * f13 = new TF1("anglef3D", getangle, 1, 10, 1);
     109  f13->SetParameter(0, __THICKNESS_3D); // detector thickness
     110  f13->SetRange(f13->Eval(3), f13->Eval(1));
     111
     112  //Double_t maxa = h_first->GetMaximum()*0.75;
     113  Double_t maxa = 0.75;
     114  TGaxis * ax23 = new TGaxis(1, maxa, 3, maxa, "anglef3D", 510, "-");
     115  ax23->SetTitle("shallow angle 3D [deg]");
    83116  //ax2->SetLabelColor(kRed);
    84   ax2->Draw("");
     117  ax23->Draw("");
    85118
    86   cout << "1 mil = " << f1->Eval(1.) << " deg" << endl;
    87   cout << "2.75 mil = " << f1->Eval(2.75) << " deg" << endl;
    88   cout << "4 mil = " << f1->Eval(4.) << " deg" << endl;
    89   cout << "5 mil = " << f1->Eval(5.) << " deg" << endl;
     119         */
    90120
    91   c11->Update();
     121        cout << "1.0 mm, planar  = " << f1p->Eval(1) << " deg" << endl;
     122        //cout << "1.0 mm, 3d      = " << f13->Eval(1) << " deg" << endl;
     123
     124        cout << "2.5 mm, planar  = " << f1p->Eval(2.5) << " deg" << endl;
     125        //cout << "2.5 mm, 3d      = " << f13->Eval(2.5) << " deg" << endl;
     126
     127        cout << "5.0 mm, planar  = " << f1p->Eval(5) << " deg" << endl;
     128        //cout << "5.0 mm, 3d      = " << f13->Eval(5) << " deg" << endl;
     129
     130        c11->Update();
    92131
    93132}

  • mafalda/ShallowAngleFEIX/runClusterSize.C

    r241 r251  
    1212
    1313
    14   GetTree("/home/idarraga/analysis/mafalda/MAFOutput_MPXNtuple_USBPIXI4_23_61417.root",
    15           "ShallowAngleFEIX", "61417 3D-SCC87 p-5E15", 1.0);
     14  GetTree("/home/idarraga/analysis/mafalda/MAFOutput_MPXNtuple_USBPIXI4_22_61548.root",
     15          "ShallowAngleFEIX", "LUB4_n-5E15_1000V", 1.0);
    1616
    17   GetTree("/home/idarraga/analysis/mafalda/MAFOutput_MPXNtuple_USBPIXI4_22_61548.root",
    18           "ShallowAngleFEIX", "61548 Planar-LUB4 n-5E15", 1.0);
     17  GetTree("/home/idarraga/analysis/mafalda/MAFOutput_MPXNtuple_USBPIXI4_22_61564.root",
     18          "ShallowAngleFEIX", "LUB4_n-5E15_600V", 1.0);
     19
     20  GetTree("/home/idarraga/analysis/mafalda/MAFOutput_MPXNtuple_USBPIXI4_22_61566.root",
     21          "ShallowAngleFEIX", "LUB4_n-5E15_800V", 1.0);
    1922
    2023  //GetTree("/storage1/idarraga/Testbeam_September_CERN_2011/KEK_20V_200V_4deg/MAFOutput/MAFOutput_MPXNtuple_USBPIXI4_21.root",
    … …  
    5861  ////////////////////////////////////////////////////
    5962 
    60   TCanvas * c11 = new TCanvas("c11", "cluster size");
     63  TCanvas * c11 = new TCanvas("c11", "cluster size X");
    6164
    6265  c11->cd();
    … …  
    6568  TH1 * h_first = (TH1 *) DrawAll("clusterSize_X","clusterHeight < 3 && TMath::Abs(mipThetaAngle) < 0.17", opts, afl);
    6669  //h_first->GetXaxis()->SetRangeUser(0, 10);
    67   h_first->GetXaxis()->SetTitle("Number of sub-clusters in track");
     70  h_first->GetXaxis()->SetTitle("cluster size X [pixels]");
     71
     72
     73  TCanvas * c12 = new TCanvas("c12", "cluster size");
     74  c12->cd();
     75  opts.doLegend = true;
     76  opts.constructorStr = "(30,0,3)"; //mipLength_X==7 && TMath::Abs(mipThetaAngle) < 0.17
     77  TH1 * h_first = (TH1 *) DrawAll("mipLength","clusterHeight < 3 && TMath::Abs(mipThetaAngle) < 0.17", opts, afl);
     78  //h_first->GetXaxis()->SetRangeUser(0, 10);
     79  h_first->GetXaxis()->SetTitle("track length [mm]");
     80
    6881
    6982  c11->Update();
     83  c12->Update();
     84
    7085
    7186}

  • mafalda/ShallowAngleFEIX/runNewAlgo.C

    r240 r251  
    66//#define __THICKNESS 0.320 // [mm]  KEK FE-I3
    77//#define __THICKNESS 0.150 // [mm]    KEK FE-I4 run 50734 to 50744
    8 #define __THICKNESS 0.250 // [mm]    Lub4 run 61548
     8#define __THICKNESS_PLANAR 0.200 // [mm]    Lub4  run 61548
     9#define __THICKNESS_3D     0.230 // [mm]    SCC83 run 61417
    910
    1011#define __CHARGE  10
    … …  
    1718        TChain * T1 = new TChain("ShallowAngleFEIX");
    1819
    19         T1->Add("/home/idarraga/analysis/mafalda/MAFOutput_MPXNtuple_USBPIXI4_23_61417.root");
     20        //T1->Add("/home/idarraga/analysis/mafalda/MAFOutput_MPXNtuple_USBPIXI4_22_61548_weighted.root");
     21        T1->Add("/home/idarraga/analysis/mafalda/MAFOutput_MPXNtuple_USBPIXI4_22_61548.root");
    2022        TString titleT = "";
    2123
    … …  
    3032
    3133        Int_t plotType = __CHARGE;
    32         Int_t ndiv = 7;
    33         Int_t fulldiv = 7;
     34        Int_t ndiv = 5;
     35        Int_t fulldiv = 5;
    3436        Int_t rangeYMax = 14;
    3537
    3638        TCanvas * c10 = new TCanvas();
    3739
    38         TGraphErrors * g1 = GetGraph(T1, "3D", 7, 7, plotType);
     40        TGraphErrors * g1 = GetGraph(T1, "LUB4 - 1000V", 6, 6, plotType, __THICKNESS_PLANAR);
    3941        c10->cd();
    4042        g1->Draw("A*");
    … …  
    5860
    5961        TChain * T2 = new TChain("ShallowAngleFEIX");
    60         T2->Add("/home/idarraga/analysis/mafalda/MAFOutput_MPXNtuple_USBPIXI4_22_61548.root");
     62        //T2->Add("/home/idarraga/analysis/mafalda/MAFOutput_MPXNtuple_USBPIXI4_22_61564_weighted.root");
     63        T2->Add("/home/idarraga/analysis/mafalda/MAFOutput_MPXNtuple_USBPIXI4_22_61564.root");
    6164        //T2->Add("/home/idarraga/analysis/mafalda/MAFOutput_MPXNtuple_AllPix_det_200_shallowAngle_3D.root");
    6265
    63         TGraphErrors * g2 = GetGraph(T2, "Planar", 6, 6, plotType);
     66        TGraphErrors * g2 = GetGraph(T2, "LUB4 600V", 6, 6, plotType, __THICKNESS_PLANAR);
    6467        // fix last point
    6568        c10->cd();
    … …  
    6871        g2->SetMarkerColor(kRed);
    6972
     73        /////////////////////////////////////////////////////////////////////////////////////
     74        TChain * T3 = new TChain("ShallowAngleFEIX");
     75        //T3->Add("/home/idarraga/analysis/mafalda/MAFOutput_MPXNtuple_USBPIXI4_22_61566_weighted.root");
     76        T3->Add("/home/idarraga/analysis/mafalda/MAFOutput_MPXNtuple_USBPIXI4_22_61566.root");
     77        //T2->Add("/home/idarraga/analysis/mafalda/MAFOutput_MPXNtuple_AllPix_det_200_shallowAngle_3D.root");
     78        TGraphErrors * g3 = GetGraph(T3, "LUB4 800V", 6, 6, plotType, __THICKNESS_PLANAR);
     79        // fix last point
     80        c10->cd();
     81        g3->Draw("*");
     82        g3->SetMarkerStyle(22);
     83        g3->SetMarkerColor(kGreen);
     84
     85        /*
     86        /////////////////////////////////////////////////////////////////////////////////////
     87        TChain * T4 = new TChain("ShallowAngleFEIX");
     88        T4->Add("/home/idarraga/analysis/mafalda/MAFOutput_MPXNtuple_USBPIXI4_22_61573.root");
     89        //T2->Add("/home/idarraga/analysis/mafalda/MAFOutput_MPXNtuple_AllPix_det_200_shallowAngle_3D.root");
     90        TGraphErrors * g4 = GetGraph(T4, "LUB4 400V", 4, 4, plotType, __THICKNESS_PLANAR);
     91        // fix last point
     92        c10->cd();
     93        g4->Draw("*");
     94        g4->SetMarkerStyle(22);
     95        g4->SetMarkerColor(kMagenta);
     96*/
    7097
    7198        /////////////////////////////////////////////////////////////////////////////////////
    7299        TLegend * l1 = new TLegend(0.5, 0.7, 0.9, 0.9);
    73         l1->AddEntry(g1, "61417 3D-SCC87 p-5E15", "P");
    74         l1->AddEntry(g2, "61548 Planar-LUB4 n-5E15", "P");
     100        l1->AddEntry(g1, "61548 LUB4 1000V", "P");
     101        l1->AddEntry(g2, "61564 LUB4 600V", "P");
     102        l1->AddEntry(g3, "61566 LUB4 800V", "P");
     103        //l1->AddEntry(g4, "61569 LUB4 400V", "P");
    75104        l1->Draw();
    76105        l1->SetBorderSize(1);
    … …  
    80109}
    81110
    82 TGraphErrors * GetGraph(TChain * T, TString sufx, Int_t __NDIV, Int_t fullDiv, Int_t graphType){
     111TGraphErrors * GetGraph(TChain * T, TString sufx, Int_t __NDIV, Int_t fullDiv, Int_t graphType, Double_t thickness){
    83112
    84113        gROOT->ProcessLine(".x ~/styles/AtlasStyle.C");
    … …  
    106135        TString wordDraw = "";
    107136        TString histo = "";
    108         double partition = double(__THICKNESS)/double(fullDiv);
     137        double partition = double(thickness)/double(fullDiv);
    109138
    110139        int partitionCntr = 0;
    111         //for (int indx = __NDIV-1 ; indx >= 0 ; indx--) {
    112         for (int indx = 0 ; indx < __NDIV ; indx++) {
     140        for (int indx = __NDIV-1 ; indx >= 0 ; indx--) {
     141        //for (int indx = 0 ; indx < __NDIV ; indx++) {
    113142
    114143                wordDraw = "totPerSegment[";
    … …  
    128157                c1->cd(indx+1);
    129158                TString cutS = "clusterHeight < 3 ";
    130                 cutS += " && clusterSize_X == ";
     159                cutS += " && (clusterSize_X == ";
    131160                cutS += __NDIV;
     161                cutS += " || clusterSize_X == ";
     162                cutS += __NDIV-1;
     163                cutS += " )";
    132164                cutS += " && TMath::Abs(mipThetaAngle) < 0.17"; // ~10 deg
    133165                TCut cutC = cutS;
    … …  
    202234        return g;
    203235}
     236
     237//Float_t GetSegment(Float a)
Note: See TracChangeset for help on using the changeset viewer.