DA for calibrating energy by pi0 and MIP and related classes (Hisayuki Torii).
[u/mrichter/AliRoot.git] / PHOS / AliPHOSDATreeCluster.cxx
CommitLineData
12f6dd57 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"
14ClassImp(AliPHOSDATreeCluster)
15//------------------------------------------------------------------------
16AliPHOSDATreeCluster::AliPHOSDATreeCluster(const AliPHOSDATreeCluster& cluster):
17//fEnergy(cluster.fEnergy),fX(cluster.fX),fY(cluster.fY),fZ(cluster.fZ),fNDigits(cluster.fNDigits){
18fEnergy(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//------------------------------------------------------------------------
32AliPHOSDATreeCluster& 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//------------------------------------------------------------------------
55void 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//------------------------------------------------------------------------
67std::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//------------------------------------------------------------------------
79void 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//------------------------------------------------------------------------
92bool 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//------------------------------------------------------------------------
128bool 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//------------------------------------------------------------------------
160bool 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//------------------------------------------------------------------------
172bool 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//------------------------------------------------------------------------
184bool 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//------------------------------------------------------------------------