AMS-02 Offline software: ECAL Reconstruction

1. EcalHitR class

The code of this class is contained in the files root.h and root.C. The class is composed by a few functions and 12 public attributes.


      Failed to include URL http://ams.cern.ch/AMS/TkTemp/GBATCH/html/classEcalHitR.html Can't connect to ams.cern.ch:443

unsigned int Status
//Statusword
32 bits to define the classification of this hit; For the meaning of each bit see the files AMS/include/amsdbc.h and AMS/CC/amsdbc.C.

Int Idsoft
//4digits number SPPC = SuperLayer/PM/subCell 1:9/1:36/1:4
The very strange definition of this variable follows: in decimal base, first digit (from right to left) defines the sub cell; second and third digits define the Photomultiplier; the fourth digit define the Super Layer.

Int Proj
//projection (0-x,1-y)

Int Plane
//ECAL plane number (0,...17)

Int Cell
//ECAL Cell number (0,...71)

float Edep
// ECAL measured energy (MeV)
In the file ecaldbc.h the array with the max capacity of ADC is defined:

const integer ECADCMX[3]={3500,4095,4095};

the first value refers to low gain overflow, the second to the high gain and the third to the dynode.
In case of real data and run number less than 1265450278 the first value becomes 2000 due to an error in DSP.
If ADC value, register in a channel (radc H/L/Dyn gain), is gt (ECADCMC[i] + pedestal value), the channel is saturated.
Then we have four scenarios to compute the correct edep value in the channel:

1)if radcH >0 and overflowHG==0 then we use this value (fadc);
2)if radcL > 3*(sigma of pedestal) and overflowLG==0 then fadc = (radcL – h2l_offset)*h2l_ratio;
3)if overflowHG==1 and radcL < 3*(sigma of pedestal) then fadc = radcH;
4)if overflowLG==1 then hen fadc = (radcL – h2l_offset)*h2l_ratio;

For the high to low offset and high to low ratio values see files described below.
After that the gain correction is applied multiplying by 1/gain() (see below) and then adc value is converted in MeV multiplying for adc to MeV converter factor:

edep = (fadc/gain())*adc2mev();

The value of this factor is defined in the following files contained in AMSDataDir/v4.00 folder (if PGTRACK no def!!!):
EcalAnorRD.1167606001 (1.15)
EcalAnorRD.1265206478 (0.4615) for data starting from TB-2010 period.

The corresponding file for MC simulation is: EcalAnorMC.1 (1.15);

In the case of PGTRACK on, in the folder AMSDataDir/v5.00, there are only two files:
EcalAnorRD.1167606001 for the real data;
EcalAnorMC.1 for the MC data;
copy of the previous files.

The final value for edep (in Mev) is obtained adding to previous value the PM saturation correction (see below/now this contribution is not activated) and the attenuation correction (see below).


float EdCorr
//ECAL Pmsaturation1-correction (MeV) added to Edep

If the dynode signal (radc Dyn) is not saturated (see above) and radc Dyn > 3*(sigma of pedestal) we can try to use this value to correct the saturated signal in the cell. First the radc Dyn is converted in Mev:

dyn edep = radc Dyn * adc2mev()*an2dyr()*pmrgain()

where an2dyr() is the Anode/Dynode ratio (see below).
Then an estimated of energy deposit in each PM's saturated cell can be calculated:

est_edep = (dyn edep – sum of energy deposits in not saturated cell)/number of saturated cell

Moreover a saturation correction coefficient can be derived from the total charge (qmeas) deposited in the PM:

qmeas (pC) = Sum on the PM cells (edep*adc2q)

where edep is the signal in each PM cell and adc2q (adc to pC conversion factor) is 0.045;
the function pmsatf1(...,qmeas) (see files AMS/include/ecaldbc.h and AMS/CC/ecaldbc.C) returns the coefficient (saturf).
So we defined the saturation correction as follows:
EdCorr = edep*( saturf – 1)

Warning: now this correction is not activated!!!

float AttCorr
// Attenuation Correction applied (w/r to center of ecal) (MeV)
To compute the attenuation correction to be added to energy deposit in each sub cell, two attenuation factors are defined:

direct attenuation factor – attfdir = ((1-fastf)*exp(-pmd/lslow) + fastf*exp(-pmd/lfast))
where:
fastf is the percentage of short component;
pmd is the distance between PM and point of energy deposition;
lslow is the attenuation length for the long component;
lfast is the attenuation length for the short component;

reflected attenuation factor – attfrfl = ((1-fastf)*exp(-(2*hflen-pmd)/lslow)+fastf*exp(-(2*hflen-pmd)/lfast))*((1-fastf)*exp(-2*hflen/lslow)+fastf*exp(-2*hflen/lfast))*fendrf
where:
hflen is the half length of a fiber;
fendrf is the end cut reflection factor.

So we define attenuation factor:
attf = 0.5*(attfdir + attfrfl);

the correction for the energy deposit is:
AttCorr = edep*(attf0/attf -1)
where
attf0 is the attenuatione w/r to centere of the fiber
Please note:
the code for this calculation is contained in ecalrec.h and ecalrec.C files;
the value of the factor in attfdir and attfrfl formulas are defined in the following file contained in AMSDataDir/v4.00 folder (if PGTRACK no def!!!):
EcalFiatRD.1167606001
EcalFiatRD.1265206478 for data starting from TB-2010 period

the legenda to read these files is reported below:

//--------------------------------------------------------------------------
PM fibers attenuation length (fast comp)
PM(1) PM(2) ............................... PM(36)    <-- SuperLayer=1
..................................................         
PM(1) PM(2) ............................... PM(36)    <-- SuperLayer=9
 //----------------------------------------------------------------------------   
PM fibers attenuation length (slow comp)
PM(1) PM(2) ............................... PM(36)    <-- SuperLayer=1
..................................................         
PM(1) PM(2) ............................... PM(36)    <-- SuperLayer=9
//----------------------------------------------------------------------------   
PM fibers fast comp fraction
PM(1) PM(2) ............................... PM(36)    <-- SuperLayer=1
..................................................         
PM(1) PM(2) ............................... PM(36)    <-- SuperLayer=9

The fiber end-cut reflection factor is defined in the job.C file:
fendrf = 0.3;

There is also the corresponding file for the MC data (in the same folder):
EcalFiatMC.1
In the case of PGTRACK on, in the folder AMSDataDir/v5.00, there are two files:
EcalFiatRD.1167606001 for the real data;
EcalFiatMC.1 for the MC data;
copy of the previous files!!!

float Coo[3]
\\ECAL Coo (cm)

The ECAL geometry is defined in AMS/include/ecaldbc.h and AMS/CC/ecaldbc.C and in the following files contained in AMSDataDir/v4.00 folder (if PGTRACK no def!!!):

EcalAlignAss1.dat
EcalAlignMC.dat
EcalAlignPMag.dat
EcalAlignPreAss.dat
EcalAlignSpace.dat

float ADC[3]
\\ECAL (ADC-Ped) for Hi/Low/Dynode channels %BR float Ped[3]
ECAL\\Pedestals

The pedestal table is contained in the following files contained in AMSDataDir/v4.00 folder (if PGTRACK no def!!!):

eclp_df_mc.dat for MC data
eclp_df_rl.dat for real data

the legenda to read these files is reported below:

Ped      PM(1)SC(1) ..... PM(36)SC(1)
Sig       PM(1)SC(1) ..... PM(36)SC(1)
Status  PM(1)SC(1) ..... PM(36)SC(1)
...........................................................          ====> Super Layer 1 High Gain
Ped      PM(1)SC(4) ..... PM(36)SC(4)
Sig       PM(1)SC(4) ..... PM(36)SC(4)
Status  PM(1)SC(4) ..... PM(36)SC(4)
............................................................
............................................................
Ped      PM(1)SC(1) ..... PM(36)SC(1)
Sig       PM(1)SC(1) ..... PM(36)SC(1)
Status  PM(1)SC(1) ..... PM(36)SC(1)
...........................................................          ====> Super Layer 9 High Gain
Ped      PM(1)SC(4) ..... PM(36)SC(4)
Sig       PM(1)SC(4) ..... PM(36)SC(4)
Status  PM(1)SC(4) ..... PM(36)SC(4)
............................................................
...........................................................
...........................................................
Ped      PM(1)SC(1) ..... PM(36)SC(1)
Sig       PM(1)SC(1) ..... PM(36)SC(1)
Status  PM(1)SC(1) ..... PM(36)SC(1)
...........................................................          ====> Super Layer 1 Low Gain
Ped      PM(1)SC(4) ..... PM(36)SC(4)
Sig       PM(1)SC(4) ..... PM(36)SC(4)
Status  PM(1)SC(4) ..... PM(36)SC(4)
............................................................
............................................................
Ped      PM(1)SC(1) ..... PM(36)SC(1)
Sig       PM(1)SC(1) ..... PM(36)SC(1)
Status  PM(1)SC(1) ..... PM(36)SC(1)
...........................................................          ====> Super Layer 9 Low Gain
Ped      PM(1)SC(4) ..... PM(36)SC(4)
Sig       PM(1)SC(4) ..... PM(36)SC(4)
Status  PM(1)SC(4) ..... PM(36)SC(4)
...........................................................
..........................................................
.........................................................
Ped      PM(1)SC(1) ..... PM(36)SC(1)
Sig       PM(1)SC(1) ..... PM(36)SC(1)
Status  PM(1)SC(1) ..... PM(36)SC(1)
...........................................................          ====> Super Layer 1 Dyn
Ped      PM(1)SC(4) ..... PM(36)SC(4)
Sig       PM(1)SC(4) ..... PM(36)SC(4)
Status  PM(1)SC(4) ..... PM(36)SC(4)
............................................................
............................................................
Ped      PM(1)SC(1) ..... PM(36)SC(1)
Sig       PM(1)SC(1) ..... PM(36)SC(1)
Status  PM(1)SC(1) ..... PM(36)SC(1)
...........................................................          ====> Super Layer 9 Dyn
Ped      PM(1)SC(4) ..... PM(36)SC(4)
Sig       PM(1)SC(4) ..... PM(36)SC(4)
Status  PM(1)SC(4) ..... PM(36)SC(4)

The corresponding files are contained also in AMSDataDir/v5.00 folder (PGTRACK version).

float Gain
\\ 1/gain (!)
roduct of two factor: pmscgain[i] x pmrgain (where “i” range from 0 to 4);

pmrgain : PM relative (to ref. PM) gain (if A=(sum of 4 PM sub cells) pmrgain = A/Aref);
pmscgain[i] : relative (to averaged) gain of 4 PM sub cells (high Gain chain (average_of_four=1 !!!);

the information about PMs gain are contained in the following files of AMS software project contained in AMSDataDir/v4.00 folder (if PGTRACK no def!!!):

EcalRlgaRD.1167606001 for the data runs starting from 1211903814 run (Date of the first event : Tue May 27 17:56:54 2008);
EcalRlgaRD.1265206478 for the data runs starting from 1265520268 run (Date of the first event : Sun Feb 7 06:24:28 2010), e.g. TB 2010 runs;

the legenda to read these files is reported below:

DBstatus ( good/bad -> 0/1) for 4 sub cells each as HL (H -> Hchain, L -> Lchain) :
PM(1)SC(1) PM(2)SC(1) ... PM(36)SC(1) \
.................................................................           <-- SuperLayer=1
 PM(1)SC(4) PM(2)SC(4) ... PM(36)SC(4)   /
 PM(1)Dyn   PM(2)Dyn    ...  PM(36)Dyn   
/
 .................................................................   

PM(1)SC(1) PM(2)SC(1) ... PM(36)SC(1)  \      
.................................................................          <-- SuperLayer=9
PM(1)SC(4) PM(2)SC(4) ... PM(36)SC(4)  /    
PM(1)Dyn   PM(2)Dyn    ...  PM(36)Dyn   / 
//----------------------------------------------------------------------------   
PM relative gain(relative to some ref.PM) (e.g. ratio of Anode 4 x sub cell (4pixels)
 signal of Given and Ref. PM's) :   
PM(1) PM(2) ... PM(36)    <-- SuperLayer=1
........................................         
PM(1) PM(2) ... PM(36)    <-- SuperLayer=9
//----------------------------------------------------------------------------   
PM sub cell (pixel) relative(to averaged over 4 sub cells) gain :
PM(1)SC(1) PM(2)SC(1) ... PM(36)SC(1)  \  
.................................................................       <-- SuperLayer=1
PM(1)SC(4) PM(2)SC(4) ... PM(36)SC(4)  /
.................................................................    
PM(1)SC(1) PM(2)SC(1) ... PM(36)SC(1)  \      
.................................................................       <-- SuperLayer=9
PM(1)SC(4) PM(2)SC(4) ... PM(36)SC(4)  /    
//----------------------------------------------------------------------------   
PM sub cell (pixel) High/Low-channel gain ratio for each of 4 sub cells) :
PM(1)SC(1) PM(2)SC(1) ... PM(36)SC(1)  \  
.................................................................       <-- SuperLayer=1
PM(1)SC(4) PM(2)SC(4) ... PM(36)SC(4)  /
.................................................................    
PM(1)SC(1) PM(2)SC(1) ... PM(36)SC(1)  \      
.................................................................        <-- SuperLayer=9
PM(1)SC(4) PM(2)SC(4) ... PM(36)SC(4)  /    
//----------------------------------------------------------------------------   
PM sub cell (pixel) Low-channel offset for each of 4 sub cells) :
PM(1)SC(1) PM(2)SC(1) ... PM(36)SC(1)  \  
................................................................       <-- SuperLayer=1
PM(1)SC(4) PM(2)SC(4) ... PM(36)SC(4)  /
.................................................................    
PM(1)SC(1) PM(2)SC(1) ... PM(36)SC(1)  \      
.................................................................        <-- SuperLayer=9
PM(1)SC(4) PM(2)SC(4) ... PM(36)SC(4)  /    
//----------------------------------------------------------------------------   
PM Anode (sum of 4sc)/Dynode ratio :
PM(1) PM(2) ... PM(36)    <-- SuperLayer=1
.................................................         
.PM(1) PM(2) ... PM(36)    <-- SuperLayer=9
//----------------------------------------------------------------------------  

Please NOTE: the offset values are defined only for data runs starting from 1265520268.
there are also the corresponding files for the MC data (in the same folder):
EcalRlgaMC.1 for date starting from 2007/01/01
EcalRlgaMC.2 for date starting from 2008/01/01
In the case of PGTRACK on, in the folder AMSDataDir/v5.00, there are only two files:
EcalRlgaRD.1167606001 for the real data;
EcalRlgaMC.1 for the MC data;
copy of the previous files (without definition of offset)!!!

2. EcalClusterR class

The code of this class is contained in the files root.h, root.C, ecalrec.h and ecalrec.C. The class is composed by 7 functions,10 pubblic attributes, 1 protected attribute and 1 private attribute. There are also two friends classes.


      Failed to include URL http://ams.cern.ch/AMS/TkTemp/GBATCH/html/classEcalClusterR.html Can't connect to ams.cern.ch:443

2.1. Cluster Building Procedure

The building procedure of a hit cluster is contained in the file ecalrec.C. In particular the friend class Ecal1DCluster comprises the method Ecal1DCluster::build(int) used to identify hit clusters in each calorimeter plane.

  1. Search for max energy deposit - The building procedure starts with a loop on planes; for each plane the cell with the max energy deposit (emax) is identified between the cells previously classified like NO-BAD; a pointer (*imax) to Ecal hit corresponding to this deposit is declared. In addition also a pointer to each Ecal hit, corrisponding to a no-bad cell, is defined and the energy deposits of no bad cells are memorized in an array (adc[cell num]).
    Remark: the energy deposit means energy reconstructed (in MeV) in a cell.
  2. Check for Bad Channels - If emax>0 in a plane, the bad channels check is performed (in that plane). First of all mincell and maxcell are fixed by the following algorithm:
    int mincell=0
    int maxcell=72
    
    for (int i=mincell;i<maxcell;i++)
    {
         if(adc[i]<0 || (adc[i]==0 && isBad()))mincell++;
         else break;
     }
    
    for (int i=maxcell-1;i>=0;i--)
    {
         if(adc[i]<0 || (adc[i]==0 && isBad()))maxcell--;
         else break;
    }
    
    Then in the interval [mincell,maxcell] the bad channels are identified and a recovery procedure is attempted. The two channels on the left and on the right of a bad channel are checked:
    for (int i=mincell;i<maxcell;i++)
    {
       if((adc[i]<0 || (adc[i]==0 && isBad())))
           {
               bool bl=i-1>=0;
               if(bl){
                   bl =bl && isBad() && adc[i-1]<=0;
              }
              bool br= i+1<Maxcell;
              if(br){
                    br =br && isBad() && adc[i+1]<=0;
              }
    }
    
    if both are good, the mean edep value of these two channels is assigned to the bad channel and this channel is classified with the flags LEAK + RECOVERED; else if only one (left or right) is good, the half edep value of that good channel is assigned to bad channel and also in this case it is classified with the flags LEAK + RECOVERED; a new pointer to recovered ecal hit are declared:
    if(!bl && !br){
         adc[i]=(adc[i-1]+adc[i+1])/2;
            if(adc[i]){
                status[i] |= LEAK;
                int st=statusa[i] | RECOVERED; 
                AMSEcalHit *ptrh[i]=new AMSEcalHit(st, ...);
            }
    }
    
    else if(!bl){
       adc[i]=(adc[i-1])/2;
       if(adc[i]){
           statusa[i] | = LEAK;
           int st=statusa[i]  | RECOVERED;
           AMSEcalHit * ptrh[i]=new AMSEcalHit(st,...);
       }
    }
    
    else if(!br){
       adc[i]=(adc[i+1])/2;
       if(adc[i]){
             statusa[i] |= LEAK;
             int st=statusa[i] | RECOVERED;
             AMSEcalHit * ptrh[i]=new AMSEcalHit(st,...);
        }
    }
    
    else{
       statusa[i] |= CATLEAK;
       adc[i]=0;
    }
    
    otherwise the bad channels is classified as CATASTROPHIC LEAK and excluded from the cluster building procedure.
  3. Star Cluster Search - First of all, a energy threshold (Thr) for a cluster hit is fixed (for each plane):
    Thr=max( min(emax,Thr1DSeed)),emax*Thr1DRSeed);
    
    where Thr1DSeed=10 and Thr1DRSeed=0.18.
    Also a reference value is set:
     ref= - FLT_MAX;
    Then a loop on plane cell in the interval [mincell,maxcell] starts; if energy deposit in a cell is greater than reference value, the reference value is updated:
    if(adc[i]>Thr)
           ref=adc[i];
    
    else if the signals in two consecutive cells (with indexes i and i+1) is less than reference, the center of a new cluster is found and it is set at cell number i-1; a procedure to identify the cluster edges is performed:
    for (int i=mincell;i<maxcell+1;i++){
        if(adc[i]<ref && adc[i+1]<ref){
    
    // two cases are defined for energy deposit in a cell and in the two consecutives : WIDE and NEAR
    
              if( adc[i]< adc[i+1] && adc[i+1]> Thr && adc[i+2]>adc[i+1]){
                   status | =WIDE;
              } 
              else if(adc[i+1]<adc[i+2] && adc[i+2]>Thr)
                          status |= NEAR;
              }
    
              int center=i-1; // cluster center fixed
              int left=mincell;
              int right=maxcell-1;
    
    //search for left edge of the cluster
              for(int k=center-1;k>=mincell;k--){
                        if(adc[k]<=0){
                             left=k+1;
                              break;
                        }
               }
    
    //search for right edge of the cluster
              for(int k=center+1;k<maxcell;k++){
                        if((adc[k]>adc[k-1] && adc[k]>Thr && adc[k+1]>adc[k]) || (adc[k]==0) ){
                            right=k-1;
                            break;
                        }  
               }
    
  4. Getting Coordinates - The energy center of gravity (c.o.g) is calculated in a core 5 cells wide around the cluster center:
      
       ileft=center-CoreSize ; //=2
       if( ileft<left)ileft=left;
       iright=center+CoreSize;
       if(iright>right)iright=right;
       ilr=min(iright-center,center-ileft);
       ileft=center-ilr;
       iright=center+ilr; 
       
       w=0;
       cx=0;
       cz=0;
       for(int k=ileft;k<=iright;k++){
           if(k==center+1 && (status & WIDE)){
              e=adc[k]/2;  // in the WIDE case the weight of right end is half energy deposit
           }
           else
              e=adc[k];
            
           cx += CellCooX/Y*e;
           cz += CellCooZ*e;
           w+=e;
       }
              
       if(w){
          cx/=w;
          cz/=w;
       }
      
       AMSPoint coo(0,0,cz);  // the cluster coordinates are memorized  
       coo[proj]=cx;
    
  5. Getting Energies - The sum of the total energy (ec) in the cluster is then calculated and also the energies within intervals around cluster center wide respectively 5, 7 and 9 cells (ec5,ec7,ec9) are also calculated. If a cluster has an edge (left or right) on the edge of the plane (mincell or maxcell, see above) and the distance from the center is within a “leakage distance”, an estimate of possible leakage (leak) is calculated. Moreover also the sum (dead) of all hits classified LEAK in the cluster is performed. Finally the c.o.g. of all cluster and the rms are calculated; each single cluster hit is classified USED.
    integer ir=right-center;
    integer il=center-left;
     
    if(left==mincell && il<Cl1DLeakSize && il<ir)   //Leak Size = 9
       {
       
         for(int k=center+il+1;k<=center+min(Cl1DLeakSize,ir);k++)
                 {
                    leak+=adc[k];   // Estimate of left leak
                 }  
        }
    
    else if(right==maxcell-1 && ir<il && ir<Cl1DLeakSize)
            {
            
            for(int k=max(left,center-Cl1DLeakSize);k<=center-ir-1;k++)
                 {
                   leak+=adc[k]; // Estimate of right leak
                 }
             }
    
    for(int k=left;k<=right;k++)
        {
          if(ptrh[k])  (ptrh[k])->setstatus(USED);   //hit classified USED
          
         if(statusa[k] & LEAK && adc[k]>0)
            {
                dead+=adc[k];    // Estimate of total signals in dead channel
            }
    
           if(k==center+1 && (status & WIDE))
             {
                 e=adc[k]/2;  //in the WIDE case the weight is always half deposit 
              }
           else e=adc[k];
    
          if(abs(k-center)<2)
             {
                 ec3+=e;  //energy in 5 cell
              }
           if(abs(k-center)<3)
             {
               ec5+=e; //energy in 7 cell
              }
          if(abs(k-center)<5)
             {
               ec9+=e; //energy in 9 cell
             }
          
         ec += e; //total energy
         x   += CellCoo(ipl,k,0)*e; 
         x2 += CellCoo(ipl,k,0)*CellCoo(ipl,k,0)*e;  
         }
         if(ec)
          {
            x/=ec; //center of gravity
            x2/=ec; 
            x2=sqrt(fabs(x2-x*x)); //RMS
           }
    
The cluster building procedure then ends; the new cluster is archived. A new search starts.

ALERT! Remark: the definitions of constants used in this procedure (Cluster Core Size,...) are in job.C file; at the end of the cluster searching in a plane, the hits in the interval [mincell,maxcell], not used to build a cluster, are included in a JUNK cluster.

3. Ecal2DClusterR class

The code of this class is contained in the files root.h,root.C, ecalrec.h and ecalrec.C. The class is composed by 6 functions, 6 pubblic attributes and 1 protected attribute. There are also two friends classes.


      Failed to include URL http://ams.cern.ch/AMS/TkTemp/GBATCH/html/classEcal2DClusterR.html Can't connect to ams.cern.ch:443

3.2. 2D Cluster Building Procedure

The building procedure of a 2D cluster is contained in the file ecalrec.C. In particular the friend class AMSEcal2DCluster comprises the method build(int) used to identify 2D clusters in each calorimeter projection.

  1. Search for max energy 1D cluster in each view - The building procedure starts with a loop on ecal projections; for each view the 1D cluster with the max energy is identified between the clusters previously classified like GOOD; a pointer ( *pshmax) to this cluster is declared; in addition also a pointer to first 1D (GOOD) cluster for each plane in the projection is memorized in an array (p1d[plane]).
    .....
    for (int proj=0;proj<2;proj++){
        Ecal1DCluster *pshmax=0;
        
       do{
        pshmax=0;
        Ecal1DCluster *ptr=(Ecal1DCluster*)AMSEvent::gethead()->
                                   getheadC("Ecal1DCluster",proj,1);
         bool newplane=true;   
         Ecal1DCluster *p1d[Maxrow];
    .....
         while(ptr){
          if(ptr->Good()){
           if(!pshmax || pshmax->getEnergy()<ptr->getEnergy()) pshmax=ptr;
            if(newplane){
             p1d[ptr->getplane()]=ptr;
             newplane=false;
            }    
          }
          if(ptr->testlast())newplane=true;
         ptr=ptr->next();
         }
    ......
    
  2. Search for clusters matching with pshmax - After that, two procedures to search, for each view, clusters matching in energy with pshmax start; the first procedure scans the planes after the max:
     Ecal1DCluster *p1c[Maxrow];
      ......
      p1c[pshmax->getplane()]=pshmax;
      Ecal1DCluster *plast=pshmax;
      for(int ipl=pshmax->getplane()+1;ipl<ECALDBc::GetLayersNo();ipl++){
           Ecal1DCluster *pcan=0;
           for(Ecal1DCluster *p=p1d[ipl];p;p=p->next()){
             if(p->Good() && p->Distance(plast)<ECREFFKEY.Thr2DMax){
              pcan=p->EnergyMatch(plast,pcan);
             }
             if(p->testlast()) break;
           }
           p1c[ipl]=pcan;
           if(pcan) plast=pcan;   
        }
    
    starting from the first 1D cluster ( p1c[plane]) , all clusters in each plane are checked; if the distance from the cluster identified in the previous plane (at the loop beginning this cluster is phsmax) is less than Thr2DMax and the Cluster Energy is similar to it, the cluster is memorized and the scan of next plane stars;otherwise the second procedure scans the planes before the max:
       plast=pshmax;
         for(int ipl=pshmax->getplane()-1;ipl>=0;ipl--){
           Ecal1DCluster *pcan=0;
           for(Ecal1DCluster *p=p1d[ipl];p;p=p->next()){
             if(p->Good() && p->Distance(plast)<ECREFFKEY.Thr2DMax){
              pcan=p->EnergyMatch(plast,pcan);
             }
             if(p->testlast())break;
          }
           p1c[ipl]=pcan;
           if(pcan)plast=pcan;   
         }.
    
    The match of a new cluster *p depends on the Distance(...) between *p and *plast and on the EnergyMatch(...) between them.
    *Distance(Ecal1DCluster p1Dcl)* is a method of Ecal1DCluster class; it returns the tangent of the angle between the line joining the two clusters ( *p and *p1Dcl) and ECAL Z-axis (see ecalrec.C file):
     number Ecal1DCluster::Distance(Ecal1DCluster *p1Dcl){
      if(_proj==p->getproj()){
       number dz=_Coo[2]-p->getcoo()[2];
       if(dz){
        return fabs(_Coo[_proj]-p->getcoo()[_proj])/fabs(dz);
       }
      }
      return FLT_MAX;
     }
    
    ECREFFKEY.Thr2DMax is the maximum value permitted for the Distance and the default is 1.2 (see job.C file); it can be modified in input card.
    Also EnergyMatch(Ecal1DCluster *pmatch, Ecal1DCluster *pbest) is a method of Ecal1DCluster class; it returns *p itself if the difference between Energies (absolute value) of *p cluster and *pmatch cluster is less than the difference between *pmatch and *pbest (see ecalrec.h):
    Ecal1DCluster* EnergyMatch(Ecal1DCluster *pmatch, Ecal1DCluster *pbest)
               {
                                  return (pbest && fabs(_EnergyC-pmatch->_EnergyC)>
                                   fabs(pmatch->_EnergyC-pbest->_EnergyC))?pbest:this;
                }.
    
    At the end of the two procedures, an array of 1D clusters ( p1c[iplane]) is built (in a view).
  3. Linear fit and cluster 2D declaration - After the building procedures, a linear fit is performed; if we succeed and the chi2 value is less than ECREFFKEY.Chi22DMax, a new 2D Cluster *pcl is declared; the 1D clusters contained in this clusters are classified as USED:
        number chi2,t0,tantz;
         integer tot;
         bool reset=false;
         bool suc=StrLineFit(p1c,ECALDBc::GetLayersNo(),proj,reset,NULL,tot,chi2,t0,tantz);
          if(suc && chi2<ECREFFKEY.Chi22DMax){
          for(int ipl=0;ipl<ECALDBc::GetLayersNo();ipl++){
           if(p1c[ipl])p1c[ipl]->setstatus(AMSDBc::USED);
          }
          AMSEcal2DCluster *pcl=new AMSEcal2DCluster(proj,tot,p1c,t0,tantz,chi2);
          pcl->_Fit();
          AMSEvent::gethead()->addnext(AMSID("Ecal2DCluster",0),pcl);
          }
         else pshmax->setstatus(AMSDBc::DELETED);
         }
       
       }while(pshmax);
    
    StrLineFit(...) is method of AMSEcal2DCluster class; it tries to perform a linear fit using a chi2 minimization algorithm.

    ECREFFKEY.Chi22DMax is the maximum value permitted for the Chi2 and the default is 1000 (see job.C file); it can be modified in input card.
    *_Fit(...)* is a method of AMSEcal2DCluster class; it performs the calculation of EnergyC, SideLeak and DeadLeak by summing the corresponding values from the 1D clusters composing the cluster; it computes also Energy, Energy3C-5C-9C by summing all hits in the 2D cluster.
    When a 2D cluster building procedure ends, a new one starts from the first point (search for the max energy cluster) within the set of the 1D clusters not yet used.
  4. Orphaned Clusters - At the end of search for 2D clusters, some orphaned 1D clusters are not used.A loop on these clusters is performed and if the distance between a 2D clusters and a 1D cluster is less than TransShowerSize2D (divided by cos of 2D cluster polar angle), the 1D cluster is added to the 2D cluster:
    while(ptr){
        if(!ptr->checkstatus(AMSDBc::BAD) && !ptr->checkstatus(AMSDBc::USED) &&
        !ptr->checkstatus(AMSDBc::JUNK))
       { 
                    AMSEcal2DCluster *p2d=(AMSEcal2DCluster*)AMSEvent::gethead()->
                                                                getheadC("Ecal2DCluster",0,1);
                    number dist=FLT_MAX;
                    AMSEcal2DCluster *p2dc=0;
                    while(p2d)
                             {
                                 if(p2d->getproj()==proj)
                                        {
                                            number intercep=p2d->getcoo()+ptr->getcoo()[2]*p2d->gettan();
                                            if(fabs(intercep-ptr->getcoo()[proj])<dist)
                                                {
                                                   dist=fabs(intercep-ptr->getcoo()[proj]);
                                                    p2dc=p2d;
                                             }
                                       }
                              p2d=p2d->next();
                          }
                         if(p2dc)
                           {
                                 if(dist<ECREFFKEY.TransShowerSize2D/fabs(cos(atan(p2dc->gettan()))))
                                     p2dc->_AddOrphane(ptr);
             
                           }
                 }
             ptr=ptr->next();
       }
    
    AddOrphane(...) is a method of AMSEcal2DCluster class (see ecalrec.C); it is similar to _Fit() method (see above), but it computes also _OrpLeak value (contribution to the 2D cluster energy from orphaned clusters);
    *ECREFFKEY.TransShowerSize2D* is the transversal shower size and the default is 10 (see job.C file); it can be modified in input card.

4. EcalShowerR class

The code of this class is contained in the files root.h,root.C, ecalrec.h and ecalrec.C.


      Failed to include URL http://ams.cern.ch/AMS/TkTemp/GBATCH/html/classEcalShowerR.html Can't connect to ams.cern.ch:443

  1. Shower Building.
  2. Add Orphans.

Shower Building

Energy Reconsruction: EnergyE . In order to obtain the particle energy, the shower deposited energy (EnergyD), already corrected for the light attenaution lenghts and the lateral leakage as above explained, must be furtherly corrected for the "Anode efficiency" and the "Rear Leakage".

  1. Anode Efficiency. The energy reconstructed by ECAL shows to be dependent on the particle impact point: this is mainly due to geometrical inefficiency at the anode edges. The correction is done using “S1S3 ratios” : S1 is the sum of most energetic cell in each layer, S3 is the sum of S1 plus the cells adjacent to the most energetic cell in each layer. Two S1S3 ratios are computed: one for X side (S1S3X) and one for Y side (S1S3Y).
    The corrected energy, Ecorr, is obtained from the S1S3 ratio via the following formula (constant parameters are tuned on TB):

    Ecorr = Ecorrx + Ecorry;
    Ecorrx = Edepx / [1 + (EFFmax -EFFmin)/(S1S3max - S1S3min)*(S1S3X -S1S3max)];
    where
    - Edepx is the deposited energy in the side X;
    - EFFmax is the best efficiency ( = 1 );
    - EFFmin is the minimum efficiency ( = 0.922);
    - S1S3max are the S1S3X value where the function reaches the maximum efficiency ( = 0.753);
    - S1S3min are the S1S3X value where the function reaches the minimum efficiency ( = 0.529).
    Similar way for Y side with:
    - Edepy is the deposited energy in the side Y;
    - EFFmax is the best efficiency ( = 1 );
    - EFFmin is the minimum efficiency ( = 0.933);
    - S1S3max is the S1S3Y value where the function reaches the maximum efficiency ( = 0.743);
    - S1S3min is theS1S3Y value where the function reaches the minimum efficiency ( = 0.517).

  2. Rear Leakage. A "rear leakage" correction is appllied to take into account the missing energy due to the longitudinal leakage. We can apply a correction factor to the deposited energy corrected for anode efficiency (Ecorr) to obtain the particle energy. This correction factor is tuned from TB and it is a quadratic function of the fraction of energy deposited in the last 2 layers (Elast/Ecorr) :

    Eparticle = [1 + alpha + beta * Elast/Ecorr + gamma * (Elast/Ecorr)^2] * Ecorr.
    - Eparticle is the particle energy;
    - Elast is the energy in the last two layers.
    The term alpha is not costant and depends, weakly, on the particle energy: the behaviour was studied on TB and fitted with a straight line; the variation is of the order of 5% between 6 and 250 GeV.
    In the reconstruction software an iterative process is used to compute alpha: at first step alpha is assumed constant and Eparticle is computed. Using this Eparticle we obatin a new value of alpha:

    EnergyE = [1 + alpha(Eparticle) + beta * Elast/Ecorr + gamma * (Elast/Ecorr)^2] * Ecorr.

Add Orphans

OrpLeak variable The orphan clusters evaluated during the 2D clustering process with the method AddOrphane() , described in previous section, are added to the shower before the leakage correction.

The OrpLeak variable is then multiplied by a magic number (0.9984) and divided by EnergyC, while in OrpLeakPI this factor is not applied and the orphanes energy is divided by EnergyE.

Edit | Attach | Watch | Print version | History: r13 < r12 < r11 < r10 < r9 | Backlinks | Raw View | WYSIWYG | More topic actions
Topic revision: r13 - 2014-11-20 - TWikiAdminUser
 
    • Cern Search Icon Cern Search
    • TWiki Search Icon TWiki Search
    • Google Search Icon Google Search

    AMS All webs login

This site is powered by the TWiki collaboration platform Powered by PerlCopyright & 2008-2019 by the contributing authors. All material on this collaboration platform is the property of the contributing authors.
Ideas, requests, problems regarding TWiki? Send feedback