PHOS/AliPHOSDATreeCluster.cxx

   1 // --
   2 // --
   3 // Implementation for TTree output in PHOS DA
   4 // for calibrating energy by pi0 and MIP.
   5 // --
   6 // -- Author: Hisayuki Torii (Hiroshima Univ.)
   7 // --
   8
   9 #include <math.h>
  10 #include <Rtypes.h>
  11 #include <iostream>
  12 #include "AliPHOSDATreeDigit.h"
  13 #include "AliPHOSDATreeCluster.h"
  14 ClassImp(AliPHOSDATreeCluster)
  15 //------------------------------------------------------------------------
  16 AliPHOSDATreeCluster::AliPHOSDATreeCluster(const AliPHOSDATreeCluster& cluster):
  17 //fEnergy(cluster.fEnergy),fX(cluster.fX),fY(cluster.fY),fZ(cluster.fZ),fNDigits(cluster.fNDigits){
  18 fEnergy(cluster.fEnergy),fRow(cluster.fRow),fCol(cluster.fCol),fNDigits(cluster.fNDigits),fDigits(0){
  19   // Copy Constructor
  20
  21   if( fNDigits > 0 ){
  22     fDigits = new AliPHOSDATreeDigit[fNDigits];
  23     int ndigits = fNDigits;
  24     while( ndigits-- ){
  25       fDigits[ndigits] = cluster.fDigits[ndigits];
  26     }
  27   } else {
  28     fDigits = 0;
  29   }
  30 }
  31 //------------------------------------------------------------------------
  32 AliPHOSDATreeCluster& AliPHOSDATreeCluster::operator=(const AliPHOSDATreeCluster& cluster){
  33   // Copy Operator
  34
  35   if( fNDigits> 0 ) delete[] fDigits;
  36   fEnergy = cluster.fEnergy;
  37   fNDigits = cluster.fNDigits;
  38   fRow = cluster.fRow;
  39   fCol = cluster.fCol;
  40   //fX = cluster.fX;
  41   //fY = cluster.fY;
  42   //fZ = cluster.fZ;
  43   if( fNDigits > 0 ){
  44     fDigits = new AliPHOSDATreeDigit[fNDigits];
  45     int ndigits = fNDigits;
  46     while( ndigits-- ){
  47       fDigits[ndigits] = cluster.fDigits[ndigits];
  48     }
  49   } else {
  50     fDigits = 0;
  51   }
  52   return *this;
  53 }
  54 //------------------------------------------------------------------------
  55 void AliPHOSDATreeCluster::Print(char* opt){
  56   // Print out
  57   std::cout<<" AliPHOSDATreeCluster:: Energy="<<fEnergy<<" NDigits="<<fNDigits
  58            <<" (row,col)=("<<fRow<<","<<fCol<<")"<<std::endl;
  59     //<<" (x,y,z)=("<<fX<<","<<fY<<","<<fZ<<")"<<std::endl;
  60   int ndigits = fNDigits;
  61   while( ndigits-- ){
  62     std::cout<<"   -->["<<ndigits<<"] : ";
  63     fDigits[ndigits].Print(opt);
  64   }
  65 }
  66 //------------------------------------------------------------------------
  67 std::ostream& operator<<(std::ostream& out, const AliPHOSDATreeCluster& cluster){
  68   // Print out
  69   out<<" AliPHOSDATreeCluster:: Energy="<<cluster.fEnergy<<" NDigits="<<cluster.fNDigits
  70      <<" (row,col)=("<<cluster.fRow<<","<<cluster.fCol<<")"<<std::endl;
  71   int ndigits = cluster.fNDigits;
  72   while( ndigits-- ){
  73     out<<"   -->["<<ndigits<<"] : "<<cluster.fDigits[ndigits];
  74     if( ndigits!=0 ) out<<std::endl;
  75   }
  76   return out;
  77 }
  78 //------------------------------------------------------------------------
  79 void AliPHOSDATreeCluster::Reset(){
  80   // Reset information
  81
  82   fEnergy = 0;
  83   fRow = -100;
  84   fCol = -100;
  85   //fX = 0;
  86   //fY = 0;
  87   //fZ = 0;
  88   if( fNDigits> 0 ) delete[] fDigits;
  89   fNDigits = 0;
  90 }
  91 //------------------------------------------------------------------------
  92 bool AliPHOSDATreeCluster::Append(AliPHOSDATreeDigit& digit){
  93   // Add digit information and sum all energy
  94   //
  95   if(! digit.IsValid() ){
  96     std::cout<<" AliPHOSDATreeCluster::Append():: Error!! Digit is not valid.."<<std::endl;
  97     return false;
  98   }
  99   AliPHOSDATreeDigit* newfDigits = new AliPHOSDATreeDigit[fNDigits+1];
 100   bool bsearching = true;
 101   int ndigit = fNDigits;
 102   while( ndigit-- ){
 103     if( fDigits[ndigit].GetAbsId() == digit.GetAbsId() ){
 104       std::cout<<" AliPHOSDATreeCluster::Append():: Error!! The channel already exist."<<std::endl;
 105       std::cout<<" Add "<<digit<<std::endl;
 106       std::cout<<" into *this"<<*this<<std::endl;
 107       delete[] newfDigits;
 108       return false;
 109     }
 110     if( fDigits[ndigit].GetEnergy() < digit.GetEnergy() ){
 111       newfDigits[ndigit+1] = fDigits[ndigit];
 112     } else {
 113       if( bsearching ) {
 114         bsearching = false;
 115         newfDigits[ndigit+1] = digit;
 116       }
 117       newfDigits[ndigit] = fDigits[ndigit];
 118     }
 119   }
 120   if( bsearching ) newfDigits[0] = digit;
 121   if( fNDigits>0 ) delete[] fDigits;
 122   fNDigits++;
 123   fDigits = newfDigits;
 124   fEnergy += digit.GetEnergy();
 125   return true;
 126 }
 127 //------------------------------------------------------------------------
 128 bool AliPHOSDATreeCluster::Append(AliPHOSDATreeCluster& cluster){
 129   // Add another cluster information and sum all energy
 130   //
 131   AliPHOSDATreeDigit* newfDigits = new AliPHOSDATreeDigit[fNDigits+cluster.fNDigits];
 132   int ndigits1 = fNDigits;
 133   int ndigits2 = cluster.fNDigits;
 134   int ndigitsall = ndigits1 + ndigits2;
 135   while( ndigitsall-- ){
 136     //std::cout<<" ------ ndigits1:"<<ndigits1<<" ndigits2:"<<ndigits2<<std::endl;
 137     if( ndigits1 && ndigits2 ){
 138       if( fDigits[ndigits1-1].GetEnergy() < cluster.fDigits[ndigits2-1].GetEnergy() ){
 139         newfDigits[ndigitsall] = fDigits[--ndigits1];
 140       } else {
 141         newfDigits[ndigitsall]= cluster.fDigits[--ndigits2];
 142       }
 143     } else if ( ndigits1 && ndigits2==0 ){
 144       newfDigits[ndigitsall] = fDigits[--ndigits1];
 145     } else if ( ndigits2 && ndigits1==0 ){
 146       newfDigits[ndigitsall]= cluster.fDigits[--ndigits2];
 147     } else {
 148       std::cout<<" AliPHOSDATreeCluster::Append() Something wrong.. "<<std::endl;
 149       return false;
 150     }
 151   }
 152   if(fNDigits>0) delete[] fDigits;
 153   fDigits = newfDigits;
 154   fNDigits += cluster.fNDigits;
 155   fEnergy += cluster.GetEnergy();
 156   return true;
 157
 158 }
 159 //------------------------------------------------------------------------
 160 bool AliPHOSDATreeCluster::IsNeighbor(const AliPHOSDATreeDigit& digit) const{
 161   // Check wether the given digit is neighboring to this cluster.
 162   // Return true if yes.
 163
 164   bool status = false;
 165   int ndigits = fNDigits;
 166   while( ndigits-- && !status ){
 167     status = digit.IsNeighbor(fDigits[ndigits]);
 168   }
 169   return status;
 170 }
 171 //------------------------------------------------------------------------
 172 bool AliPHOSDATreeCluster::IsNeighbor(const AliPHOSDATreeCluster& cluster) const{
 173   // Check wether the given cluster is neighboring to this cluster.
 174   // Return true if yes.
 175
 176   bool status = false;
 177   int ndigits = fNDigits;
 178   while( ndigits-- && !status ){
 179     status = cluster.IsNeighbor(fDigits[ndigits]);
 180   }
 181   return status;
 182 }
 183 //------------------------------------------------------------------------
 184 bool AliPHOSDATreeCluster::CalculateProperty(){
 185   // Calculate the hit position
 186   // (calculation of dispersion is not valid)
 187
 188   fCol = 0;
 189   fRow = 0;
 190   float totweight = 0;
 191   float weight;
 192   int ndigits = fNDigits;
 193   while( ndigits-- ){
 194     weight = log(fDigits[ndigits].GetEnergy()/fEnergy) + 4.5; //4.5 is for PHOS
 195     //std::cout<<" AliPHOSDATreeCluster::CalculateProperty() DEBUG: ndigits="<<ndigits<<" weight="<<weight<<std::endl;
 196     if( weight > 0 ){
 197       totweight += weight;
 198       fRow += fDigits[ndigits].GetRow() * weight;
 199       fCol += fDigits[ndigits].GetCol() * weight;
 200     }
 201   }
 202   //std::cout<<" AliPHOSDATreeCluster::CalculateProperty() DEBUG: totweight="<<totweight<<std::endl;
 203   if( totweight > 0 ){
 204     fRow /= totweight;
 205     fCol /= totweight;
 206   } else {
 207     fRow = 0;
 208     fCol = 0;
 209   }
 210   /*
 211   float disp = 0;
 212   if( totweight > ){
 213     ndigits = fNDigits;
 214     while( ndigits-- ){
 215       weight = log(fDigits[ndigits].GetEnergy()/fEnergy) + 4.5; //4.5 is for PHOS
 216       disp += weight * ( (fDigits[ndigits].GetRow()-fRow)*(fDigits[ndigits].GetRow()-fRow) +
 217                          (fDigits[ndigits].GetCol()-fCol)*(fDigits[ndigits].GetCol()-fCol) );
 218     }
 219     disp /= totweight;
 220   }
 221   */
 222   return true;
 223 }
 224 //------------------------------------------------------------------------