First version of Raw Data reconstruction. Added appropriate Reconstruct method to...
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALRecPoint.cxx
CommitLineData
ab48128d 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15/* $Id$ */
16//_________________________________________________________________________
70a93198 17// Reconstructed Points for the EMCAL
18// A RecPoint is a cluster of digits
d64c959b 19//*-- Author: Yves Schutz (SUBATECH)
70a93198 20//*-- Author: Dmitri Peressounko (RRC KI & SUBATECH)
21//*-- Author: Heather Gray (LBL) merged AliEMCALRecPoint and AliEMCALTowerRecPoint 02/04
ab48128d 22
23// --- ROOT system ---
1d59832c 24class Riostream;
e52475ed 25#include <TPad.h>
1d59832c 26class TGraph;
27class TPaveText;
e52475ed 28#include <TClonesArray.h>
29#include <TMath.h>
ab48128d 30
31// --- Standard library ---
ab48128d 32
33// --- AliRoot header files ---
1d59832c 34//#include "AliGenerator.h"
35class AliGenerator;
4635df1f 36#include "AliRunLoader.h"
37#include "AliRun.h"
1d59832c 38class AliEMCAL;
4635df1f 39#include "AliEMCALLoader.h"
ab48128d 40#include "AliEMCALGeometry.h"
4635df1f 41#include "AliEMCALHit.h"
ab48128d 42#include "AliEMCALDigit.h"
43#include "AliEMCALRecPoint.h"
ab48128d 44
45ClassImp(AliEMCALRecPoint)
46
ab48128d 47//____________________________________________________________________________
48AliEMCALRecPoint::AliEMCALRecPoint()
18a21c7c 49 : AliRecPoint(),
50 fGeomPtr(0),
51 fClusterType(-1),
52 fCoreEnergy(0),
53 fDispersion(0),
54 fEnergyList(0),
55 fTimeList(0),
56 fAbsIdList(0),
57 fTime(0.),
58 fCoreRadius(10), //HG check this
af5bdd85 59 fDETracksList(0),
18a21c7c 60 fMulParent(0),
61 fMaxParent(0),
62 fParentsList(0),
af5bdd85 63 fDEParentsList(0),
18a21c7c 64 fSuperModuleNumber(0)
ab48128d 65{
66 // ctor
692088ae 67 fMaxTrack = 0 ;
18a21c7c 68 fMulDigit = 0 ;
70a93198 69 fAmp = 0. ;
4635df1f 70
71 AliRunLoader *rl = AliRunLoader::GetRunLoader();
af5bdd85 72 if (rl->GetAliRun() && rl->GetAliRun()->GetDetector("EMCAL"))
73 fGeomPtr = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
74 else
75 fGeomPtr = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaulGeometryName());
14ce0a6e 76 fGeomPtr->GetTransformationForSM(); // Global <-> Local
ab48128d 77}
78
79//____________________________________________________________________________
18a21c7c 80AliEMCALRecPoint::AliEMCALRecPoint(const char * opt)
81 : AliRecPoint(opt),
82 fGeomPtr(0),
83 fClusterType(-1),
84 fCoreEnergy(0),
85 fDispersion(0),
86 fEnergyList(0),
87 fTimeList(0),
88 fAbsIdList(0),
89 fTime(-1.),
90 fCoreRadius(10), //HG check this
af5bdd85 91 fDETracksList(0),
18a21c7c 92 fMulParent(0),
93 fMaxParent(1000),
94 fParentsList(0),
af5bdd85 95 fDEParentsList(0),
18a21c7c 96 fSuperModuleNumber(0)
ab48128d 97{
98 // ctor
ff1e7e2f 99 fMaxTrack = 1000 ;
87cdc3be 100 fMulDigit = 0 ;
70a93198 101 fAmp = 0. ;
af5bdd85 102 fDETracksList = new Float_t[fMaxTrack];
87cdc3be 103 fParentsList = new Int_t[fMaxParent];
af5bdd85 104 fDEParentsList = new Float_t[fMaxParent];
105 for (Int_t i = 0; i < fMaxTrack; i++)
106 fDETracksList[i] = 0;
107 for (Int_t i = 0; i < fMaxParent; i++)
108 fDEParentsList[i] = 0;
18a21c7c 109
4635df1f 110 AliRunLoader *rl = AliRunLoader::GetRunLoader();
af5bdd85 111 if (rl->GetAliRun() && rl->GetAliRun()->GetDetector("EMCAL"))
112 fGeomPtr = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
113 else
114 fGeomPtr = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaulGeometryName());
14ce0a6e 115 fGeomPtr->GetTransformationForSM(); // Global <-> Local
70a93198 116}
18a21c7c 117
70a93198 118//____________________________________________________________________________
18a21c7c 119AliEMCALRecPoint::AliEMCALRecPoint(const AliEMCALRecPoint & rp)
120 : AliRecPoint(rp),
121 fGeomPtr(rp.fGeomPtr),
122 fClusterType(rp.fClusterType),
123 fCoreEnergy(rp.fCoreEnergy),
124 fDispersion(rp.fDispersion),
125 fEnergyList(0),
126 fTimeList(0),
127 fAbsIdList(0),
128 fTime(rp.fTime),
129 fCoreRadius(rp.fCoreRadius),
af5bdd85 130 fDETracksList(0),
18a21c7c 131 fMulParent(rp.fMulParent),
132 fMaxParent(rp.fMaxParent),
133 fParentsList(0),
af5bdd85 134 fDEParentsList(0),
18a21c7c 135 fSuperModuleNumber(rp.fSuperModuleNumber)
0a4cb131 136{
137 //copy ctor
0a4cb131 138 fLambda[0] = rp.fLambda[0];
139 fLambda[1] = rp.fLambda[1];
18a21c7c 140
af5bdd85 141 fEnergyList = new Float_t[rp.fMaxDigit];
142 fTimeList = new Float_t[rp.fMaxDigit];
143 fAbsIdList = new Int_t[rp.fMaxDigit];
0a4cb131 144 for(Int_t i = 0; i < rp.fMulDigit; i++) {
145 fEnergyList[i] = rp.fEnergyList[i];
146 fTimeList[i] = rp.fTimeList[i];
147 fAbsIdList[i] = rp.fAbsIdList[i];
148 }
af5bdd85 149 fDETracksList = new Float_t[rp.fMaxTrack];
150 for(Int_t i = 0; i < rp.fMulTrack; i++) fDETracksList[i] = rp.fDETracksList[i];
151 fParentsList = new Int_t[rp.fMaxParent];
0a4cb131 152 for(Int_t i = 0; i < rp.fMulParent; i++) fParentsList[i] = rp.fParentsList[i];
af5bdd85 153 fDEParentsList = new Float_t[rp.fMaxParent];
154 for(Int_t i = 0; i < rp.fMulParent; i++) fDEParentsList[i] = rp.fDEParentsList[i];
0a4cb131 155
156}
157//____________________________________________________________________________
70a93198 158AliEMCALRecPoint::~AliEMCALRecPoint()
159{
160 // dtor
161 if ( fEnergyList )
162 delete[] fEnergyList ;
85c60a8e 163 if ( fTimeList )
164 delete[] fTimeList ;
e52475ed 165 if ( fAbsIdList )
166 delete[] fAbsIdList ;
af5bdd85 167 if ( fDETracksList)
168 delete[] fDETracksList;
87cdc3be 169 if ( fParentsList)
170 delete[] fParentsList;
af5bdd85 171 if ( fDEParentsList)
172 delete[] fDEParentsList;
70a93198 173}
174
175//____________________________________________________________________________
176void AliEMCALRecPoint::AddDigit(AliEMCALDigit & digit, Float_t Energy)
177{
178 // Adds a digit to the RecPoint
179 // and accumulates the total amplitude and the multiplicity
180
181 if(fEnergyList == 0)
182 fEnergyList = new Float_t[fMaxDigit];
85c60a8e 183 if(fTimeList == 0)
184 fTimeList = new Float_t[fMaxDigit];
e52475ed 185 if(fAbsIdList == 0) {
186 fAbsIdList = new Int_t[fMaxDigit];
14ce0a6e 187 fSuperModuleNumber = fGeomPtr->GetSuperModuleNumber(digit.GetId());
e52475ed 188 }
70a93198 189
190 if ( fMulDigit >= fMaxDigit ) { // increase the size of the lists
191 fMaxDigit*=2 ;
a64a06d6 192 Int_t * tempo = new Int_t[fMaxDigit];
193 Float_t * tempoE = new Float_t[fMaxDigit];
85c60a8e 194 Float_t * tempoT = new Float_t[fMaxDigit];
e52475ed 195 Int_t * tempoId = new Int_t[fMaxDigit];
70a93198 196
197 Int_t index ;
198 for ( index = 0 ; index < fMulDigit ; index++ ){
e52475ed 199 tempo[index] = fDigitsList[index] ;
200 tempoE[index] = fEnergyList[index] ;
85c60a8e 201 tempoT[index] = fTimeList[index] ;
e52475ed 202 tempoId[index] = fAbsIdList[index] ;
70a93198 203 }
204
205 delete [] fDigitsList ;
a64a06d6 206 fDigitsList = new Int_t[fMaxDigit];
70a93198 207
208 delete [] fEnergyList ;
a64a06d6 209 fEnergyList = new Float_t[fMaxDigit];
70a93198 210
85c60a8e 211 delete [] fTimeList ;
212 fTimeList = new Float_t[fMaxDigit];
213
e52475ed 214 delete [] fAbsIdList ;
215 fAbsIdList = new Int_t[fMaxDigit];
216
70a93198 217 for ( index = 0 ; index < fMulDigit ; index++ ){
218 fDigitsList[index] = tempo[index] ;
219 fEnergyList[index] = tempoE[index] ;
85c60a8e 220 fTimeList[index] = tempoT[index] ;
e52475ed 221 fAbsIdList[index] = tempoId[index] ;
70a93198 222 }
223
224 delete [] tempo ;
225 delete [] tempoE ;
85c60a8e 226 delete [] tempoT ;
e52475ed 227 delete [] tempoId ;
70a93198 228 } // if
229
230 fDigitsList[fMulDigit] = digit.GetIndexInList() ;
231 fEnergyList[fMulDigit] = Energy ;
85c60a8e 232 fTimeList[fMulDigit] = digit.GetTime() ;
e52475ed 233 fAbsIdList[fMulDigit] = digit.GetId();
70a93198 234 fMulDigit++ ;
235 fAmp += Energy ;
236
237}
238//____________________________________________________________________________
239Bool_t AliEMCALRecPoint::AreNeighbours(AliEMCALDigit * digit1, AliEMCALDigit * digit2 ) const
240{
241 // Tells if (true) or not (false) two digits are neighbours
242 // A neighbour is defined as being two digits which share a corner
243
e52475ed 244 static Bool_t areNeighbours = kFALSE ;
2bb3725c 245 static Int_t nSupMod=0, nModule=0, nIphi=0, nIeta=0;
246 static int nSupMod1=0, nModule1=0, nIphi1=0, nIeta1=0;
e52475ed 247 static Int_t relid1[2] , relid2[2] ; // ieta, iphi
248 static Int_t rowdiff=0, coldiff=0;
70a93198 249
e52475ed 250 areNeighbours = kFALSE ;
251
2bb3725c 252 fGeomPtr->GetCellIndex(digit1->GetId(), nSupMod,nModule,nIphi,nIeta);
253 fGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, relid1[0],relid1[1]);
e52475ed 254
2bb3725c 255 fGeomPtr->GetCellIndex(digit2->GetId(), nSupMod1,nModule1,nIphi1,nIeta1);
256 fGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod1,nModule1,nIphi1,nIeta1, relid2[0],relid2[1]);
70a93198 257
e52475ed 258 rowdiff = TMath::Abs( relid1[0] - relid2[0] ) ;
259 coldiff = TMath::Abs( relid1[1] - relid2[1] ) ;
70a93198 260
261 if (( coldiff <= 1 ) && ( rowdiff <= 1 ) && (coldiff + rowdiff > 0))
e52475ed 262 areNeighbours = kTRUE ;
ab48128d 263
70a93198 264 return areNeighbours;
265}
266
267//____________________________________________________________________________
268Int_t AliEMCALRecPoint::Compare(const TObject * obj) const
269{
270 // Compares two RecPoints according to their position in the EMCAL modules
271
272 Float_t delta = 1 ; //Width of "Sorting row". If you change this
273 //value (what is senseless) change as well delta in
274 //AliEMCALTrackSegmentMakerv* and other RecPoints...
275 Int_t rv ;
276
277 AliEMCALRecPoint * clu = (AliEMCALRecPoint *)obj ;
278
279 TVector3 locpos1;
280 GetLocalPosition(locpos1);
281 TVector3 locpos2;
282 clu->GetLocalPosition(locpos2);
283
9848d950 284 Int_t rowdif = (Int_t)(TMath::Ceil(locpos1.X()/delta)-TMath::Ceil(locpos2.X()/delta)) ;
70a93198 285 if (rowdif> 0)
286 rv = 1 ;
287 else if(rowdif < 0)
288 rv = -1 ;
289 else if(locpos1.Y()>locpos2.Y())
290 rv = -1 ;
291 else
292 rv = 1 ;
293
294 return rv ;
ab48128d 295}
296
297//____________________________________________________________________________
298Int_t AliEMCALRecPoint::DistancetoPrimitive(Int_t px, Int_t py)
299{
300 // Compute distance from point px,py to a AliEMCALRecPoint considered as a Tmarker
301 // Compute the closest distance of approach from point px,py to this marker.
302 // The distance is computed in pixels units.
70a93198 303 // HG Still need to update -> Not sure what this should achieve
ab48128d 304
305 TVector3 pos(0.,0.,0.) ;
70a93198 306 GetLocalPosition(pos) ;
ab48128d 307 Float_t x = pos.X() ;
70a93198 308 Float_t y = pos.Y() ;
ab48128d 309 const Int_t kMaxDiff = 10;
310 Int_t pxm = gPad->XtoAbsPixel(x);
311 Int_t pym = gPad->YtoAbsPixel(y);
312 Int_t dist = (px-pxm)*(px-pxm) + (py-pym)*(py-pym);
313
314 if (dist > kMaxDiff) return 9999;
315 return dist;
316}
317
318//___________________________________________________________________________
319 void AliEMCALRecPoint::Draw(Option_t *option)
320 {
321 // Draw this AliEMCALRecPoint with its current attributes
322
323 AppendPad(option);
324 }
325
326//______________________________________________________________________________
70a93198 327void AliEMCALRecPoint::ExecuteEvent(Int_t /*event*/, Int_t, Int_t)
ab48128d 328{
329 // Execute action corresponding to one event
330 // This member function is called when a AliEMCALRecPoint is clicked with the locator
331 //
332 // If Left button is clicked on AliEMCALRecPoint, the digits are switched on
333 // and switched off when the mouse button is released.
334
335 // static Int_t pxold, pyold;
336
70a93198 337 /* static TGraph * digitgraph = 0 ;
ab48128d 338 static TPaveText* clustertext = 0 ;
339
340 if (!gPad->IsEditable()) return;
341
342 switch (event) {
343
344
345 case kButton1Down:{
346 AliEMCALDigit * digit ;
88cb7938 347 AliEMCALGeometry * emcalgeom = (AliEMCALGetter::Instance())->EMCALGeometry() ;
ab48128d 348
349 Int_t iDigit;
70a93198 350 Int_t relid[2] ;
ab48128d 351
352 const Int_t kMulDigit=AliEMCALRecPoint::GetDigitsMultiplicity() ;
353 Float_t * xi = new Float_t [kMulDigit] ;
354 Float_t * zi = new Float_t [kMulDigit] ;
355
356 for(iDigit = 0; iDigit < kMulDigit; iDigit++) {
357 Fatal("AliEMCALRecPoint::ExecuteEvent", " -> Something wrong with the code");
358 digit = 0 ; //dynamic_cast<AliEMCALDigit *>((fDigitsList)[iDigit]);
359 emcalgeom->AbsToRelNumbering(digit->GetId(), relid) ;
360 emcalgeom->PosInAlice(relid, xi[iDigit], zi[iDigit]) ;
361 }
362
363 if (!digitgraph) {
364 digitgraph = new TGraph(fMulDigit,xi,zi);
365 digitgraph-> SetMarkerStyle(5) ;
366 digitgraph-> SetMarkerSize(1.) ;
367 digitgraph-> SetMarkerColor(1) ;
368 digitgraph-> Draw("P") ;
369 }
370 if (!clustertext) {
371
372 TVector3 pos(0.,0.,0.) ;
373 GetLocalPosition(pos) ;
374 clustertext = new TPaveText(pos.X()-10,pos.Z()+10,pos.X()+50,pos.Z()+35,"") ;
375 Text_t line1[40] ;
376 Text_t line2[40] ;
377 sprintf(line1,"Energy=%1.2f GeV",GetEnergy()) ;
378 sprintf(line2,"%d Digits",GetDigitsMultiplicity()) ;
379 clustertext ->AddText(line1) ;
380 clustertext ->AddText(line2) ;
381 clustertext ->Draw("");
382 }
383 gPad->Update() ;
9e5d2067 384 Print("") ;
ab48128d 385 delete[] xi ;
386 delete[] zi ;
387 }
388
389break;
390
391 case kButton1Up:
392 if (digitgraph) {
393 delete digitgraph ;
394 digitgraph = 0 ;
395 }
396 if (clustertext) {
397 delete clustertext ;
398 clustertext = 0 ;
399 }
400
401 break;
402
70a93198 403 }*/
404}
405//____________________________________________________________________________
406void AliEMCALRecPoint::EvalAll(Float_t logWeight,TClonesArray * digits)
407{
408 // Evaluates all shower parameters
409
410 EvalLocalPosition(logWeight, digits) ;
500aeccc 411 // printf("eval position done\n");
70a93198 412 EvalElipsAxis(logWeight, digits) ;
500aeccc 413 // printf("eval axis done\n");
70a93198 414 EvalDispersion(logWeight, digits) ;
500aeccc 415 // printf("eval dispersion done\n");
4635df1f 416 //EvalCoreEnergy(logWeight, digits);
1963b290 417 // printf("eval energy done\n");
70a93198 418 EvalTime(digits) ;
500aeccc 419 // printf("eval time done\n");
70a93198 420
87cdc3be 421 EvalPrimaries(digits) ;
500aeccc 422 // printf("eval pri done\n");
87cdc3be 423 EvalParents(digits);
500aeccc 424 // printf("eval parent done\n");
70a93198 425}
426
427//____________________________________________________________________________
428void AliEMCALRecPoint::EvalDispersion(Float_t logWeight, TClonesArray * digits)
429{
430 // Calculates the dispersion of the shower at the origin of the RecPoint
1d46d1f6 431 // in cell units - Nov 16,2006
70a93198 432
1d46d1f6 433 Double_t d = 0., wtot = 0., w = 0.;
af5bdd85 434 Int_t iDigit=0, nstat=0;
70a93198 435 AliEMCALDigit * digit ;
70a93198 436
1d46d1f6 437 // Calculates the dispersion in cell units
438 Double_t etai, phii, etaMean=0.0, phiMean=0.0;
2bb3725c 439 int nSupMod=0, nModule=0, nIphi=0, nIeta=0;
1d46d1f6 440 int iphi=0, ieta=0;
441 // Calculate mean values
70a93198 442 for(iDigit=0; iDigit < fMulDigit; iDigit++) {
443 digit = (AliEMCALDigit *) digits->At(fDigitsList[iDigit]) ;
1963b290 444
1d46d1f6 445 if (fAmp>0 && fEnergyList[iDigit]>0) {
2bb3725c 446 fGeomPtr->GetCellIndex(digit->GetId(), nSupMod,nModule,nIphi,nIeta);
447 fGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi,ieta);
1d46d1f6 448 etai=(Double_t)ieta;
449 phii=(Double_t)iphi;
450 w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
451
452 if(w>0.0) {
453 phiMean += phii*w;
454 etaMean += etai*w;
455 wtot += w;
456 }
457 }
458 }
459 if (wtot>0) {
460 phiMean /= wtot ;
461 etaMean /= wtot ;
462 } else AliError(Form("Wrong weight %f\n", wtot));
70a93198 463
1d46d1f6 464 // Calculate dispersion
465 for(iDigit=0; iDigit < fMulDigit; iDigit++) {
466 digit = (AliEMCALDigit *) digits->At(fDigitsList[iDigit]) ;
467
468 if (fAmp>0 && fEnergyList[iDigit]>0) {
2bb3725c 469 fGeomPtr->GetCellIndex(digit->GetId(), nSupMod,nModule,nIphi,nIeta);
470 fGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi,ieta);
1d46d1f6 471 etai=(Double_t)ieta;
472 phii=(Double_t)iphi;
473 w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
474
475 if(w>0.0) {
476 nstat++;
477 d += w*((etai-etaMean)*(etai-etaMean)+(phii-phiMean)*(phii-phiMean));
e52475ed 478 }
479 }
ab48128d 480 }
70a93198 481
e52475ed 482 if ( wtot > 0 && nstat>1) d /= wtot ;
483 else d = 0. ;
70a93198 484
485 fDispersion = TMath::Sqrt(d) ;
ab48128d 486}
70a93198 487
ab48128d 488//____________________________________________________________________________
70a93198 489void AliEMCALRecPoint::EvalLocalPosition(Float_t logWeight, TClonesArray * digits)
88cb7938 490{
70a93198 491 // Calculates the center of gravity in the local EMCAL-module coordinates
e52475ed 492 // Info("Print", " logWeight %f : cluster energy %f ", logWeight, fAmp); // for testing
70a93198 493
e52475ed 494 AliEMCALDigit * digit;
e52475ed 495 Int_t i=0, nstat=0;
496 Double_t clXYZ[3]={0.,0.,0.}, clRmsXYZ[3]={0.,0.,0.}, xyzi[3], wtot=0., w=0.;
70a93198 497
e52475ed 498 for(Int_t iDigit=0; iDigit<fMulDigit; iDigit++) {
70a93198 499 digit = dynamic_cast<AliEMCALDigit *>(digits->At(fDigitsList[iDigit])) ;
500
14ce0a6e 501 fGeomPtr->RelPosCellInSModule(digit->GetId(), xyzi[0], xyzi[1], xyzi[2]);
e52475ed 502 // printf(" Id %i : Local x,y,z %f %f %f \n", digit->GetId(), xyzi[0], xyzi[1], xyzi[2]);
503
af5bdd85 504 if(logWeight > 0.0) w = TMath::Max( 0., logWeight + TMath::Log( fEnergyList[iDigit] / fAmp ));
505 else w = fEnergyList[iDigit]; // just energy
ab48128d 506
e52475ed 507 if(w>0.0) {
508 wtot += w ;
509 nstat++;
510 for(i=0; i<3; i++ ) {
511 clXYZ[i] += (w*xyzi[i]);
512 clRmsXYZ[i] += (w*xyzi[i]*xyzi[i]);
513 }
514 }
515 }
516 // cout << " wtot " << wtot << endl;
70a93198 517 if ( wtot > 0 ) {
e52475ed 518 // xRMS = TMath::Sqrt(x2m - xMean*xMean);
519 for(i=0; i<3; i++ ) {
520 clXYZ[i] /= wtot;
521 if(nstat>1) {
522 clRmsXYZ[i] /= (wtot*wtot);
523 clRmsXYZ[i] = clRmsXYZ[i] - clXYZ[i]*clXYZ[i];
524 if(clRmsXYZ[i] > 0.0) {
525 clRmsXYZ[i] = TMath::Sqrt(clRmsXYZ[i]);
526 } else clRmsXYZ[i] = 0;
527 } else clRmsXYZ[i] = 0;
528 }
70a93198 529 } else {
e52475ed 530 for(i=0; i<3; i++ ) {
531 clXYZ[i] = clRmsXYZ[i] = -1.;
532 }
70a93198 533 }
e52475ed 534 // clRmsXYZ[i] ??
535 fLocPos.SetX(clXYZ[0]);
536 fLocPos.SetY(clXYZ[1]);
537 fLocPos.SetZ(clXYZ[2]);
70a93198 538
1963b290 539// if (gDebug==2)
500aeccc 540// printf("EvalLocalPosition: eta,phi,r = %f,%f,%f", fLocPos.X(), fLocPos.Y(), fLocPos.Z()) ;
e52475ed 541 fLocPosM = 0 ; // covariance matrix
ab48128d 542}
543
e52475ed 544//void AliEMCALRecPoint::EvalLocalPositionSimple()
545//{ // Weight is proportional of cell energy
546//}
547
70a93198 548//______________________________________________________________________________
549void AliEMCALRecPoint::EvalCoreEnergy(Float_t logWeight, TClonesArray * digits)
550{
551 // This function calculates energy in the core,
4635df1f 552 // i.e. within a radius rad = fCoreEnergy around the center. Beyond this radius
70a93198 553 // in accordance with shower profile the energy deposition
554 // should be less than 2%
1d46d1f6 555 // Unfinished - Nov 15,2006
556 // Distance is calculate in (phi,eta) units
70a93198 557
558 AliEMCALDigit * digit ;
5dee926e 559
70a93198 560 Int_t iDigit;
561
e52475ed 562 if (!fLocPos.Mag()) {
70a93198 563 EvalLocalPosition(logWeight, digits);
564 }
565
1d46d1f6 566 Double_t phiPoint = fLocPos.Phi(), etaPoint = fLocPos.Eta();
567 Double_t eta, phi, distance;
70a93198 568 for(iDigit=0; iDigit < fMulDigit; iDigit++) {
569 digit = (AliEMCALDigit *) ( digits->At(fDigitsList[iDigit]) ) ;
4635df1f 570
1d46d1f6 571 eta = phi = 0.0;
572 fGeomPtr->EtaPhiFromIndex(digit->GetId(),eta, phi) ;
573 phi = phi * TMath::DegToRad();
70a93198 574
1d46d1f6 575 distance = TMath::Sqrt((eta-etaPoint)*(eta-etaPoint)+(phi-phiPoint)*(phi-phiPoint));
70a93198 576 if(distance < fCoreRadius)
577 fCoreEnergy += fEnergyList[iDigit] ;
578 }
579
580}
ab48128d 581//____________________________________________________________________________
70a93198 582void AliEMCALRecPoint::EvalElipsAxis(Float_t logWeight,TClonesArray * digits)
ab48128d 583{
70a93198 584 // Calculates the axis of the shower ellipsoid in eta and phi
1d46d1f6 585 // in cell units
ab48128d 586
1d46d1f6 587 static TString gn(fGeomPtr->GetName());
588
589 Double_t wtot = 0.;
70a93198 590 Double_t x = 0.;
591 Double_t z = 0.;
592 Double_t dxx = 0.;
593 Double_t dzz = 0.;
594 Double_t dxz = 0.;
ab48128d 595
1d46d1f6 596 AliEMCALDigit * digit = 0;
70a93198 597
1d46d1f6 598 Double_t etai , phii, w;
2bb3725c 599 int nSupMod=0, nModule=0, nIphi=0, nIeta=0;
1d46d1f6 600 int iphi=0, ieta=0;
601 for(Int_t iDigit=0; iDigit<fMulDigit; iDigit++) {
70a93198 602 digit = (AliEMCALDigit *) digits->At(fDigitsList[iDigit]) ;
1d46d1f6 603 etai = phii = 0.;
604 if(gn.Contains("SHISH")) {
605 // Nov 15,2006 - use cell numbers as coordinates
606 // Copied for shish-kebab geometry, ieta,iphi is cast as double as eta,phi
607 // We can use the eta,phi(or coordinates) of cell
2bb3725c 608 nSupMod = nModule = nIphi = nIeta = iphi = ieta = 0;
1d46d1f6 609
2bb3725c 610 fGeomPtr->GetCellIndex(digit->GetId(), nSupMod,nModule,nIphi,nIeta);
611 fGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi,ieta);
1d46d1f6 612 etai=(Double_t)ieta;
613 phii=(Double_t)iphi;
614 } else { //
14ce0a6e 615 fGeomPtr->EtaPhiFromIndex(digit->GetId(), etai, phii);
1d46d1f6 616 phii = phii * TMath::DegToRad();
1963b290 617 }
618
1d46d1f6 619 w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
620 // fAmp summed amplitude of digits, i.e. energy of recpoint
621 // Gives smaller value of lambda than log weight
622 // w = fEnergyList[iDigit] / fAmp; // Nov 16, 2006 - try just energy
ff1e7e2f 623
70a93198 624 dxx += w * etai * etai ;
625 x += w * etai ;
626 dzz += w * phii * phii ;
627 z += w * phii ;
1963b290 628
ff1e7e2f 629 dxz += w * etai * phii ;
1963b290 630
70a93198 631 wtot += w ;
632 }
ff1e7e2f 633
70a93198 634 if ( wtot > 0 ) {
635 dxx /= wtot ;
636 x /= wtot ;
637 dxx -= x * x ;
638 dzz /= wtot ;
639 z /= wtot ;
640 dzz -= z * z ;
641 dxz /= wtot ;
642 dxz -= x * z ;
ab48128d 643
70a93198 644 fLambda[0] = 0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ;
645 if(fLambda[0] > 0)
646 fLambda[0] = TMath::Sqrt(fLambda[0]) ;
647 else
648 fLambda[0] = 0;
649
650 fLambda[1] = 0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ;
ff1e7e2f 651
70a93198 652 if(fLambda[1] > 0) //To avoid exception if numerical errors lead to negative lambda.
653 fLambda[1] = TMath::Sqrt(fLambda[1]) ;
654 else
655 fLambda[1]= 0. ;
656 } else {
657 fLambda[0]= 0. ;
658 fLambda[1]= 0. ;
ab48128d 659 }
ff1e7e2f 660
1963b290 661 // printf("Evalaxis: lambdas = %f,%f", fLambda[0],fLambda[1]) ;
ff1e7e2f 662
ab48128d 663}
664
665//______________________________________________________________________________
666void AliEMCALRecPoint::EvalPrimaries(TClonesArray * digits)
667{
af5bdd85 668 // Constructs the list of primary particles (tracks) which
669 // have contributed to this RecPoint and calculate deposited energy
670 // for each track
ab48128d 671
672 AliEMCALDigit * digit ;
af5bdd85 673 Int_t * primArray = new Int_t[fMaxTrack] ;
674 Float_t * dEPrimArray = new Float_t[fMaxTrack] ;
ab48128d 675
676 Int_t index ;
677 for ( index = 0 ; index < GetDigitsMultiplicity() ; index++ ) { // all digits
678 digit = dynamic_cast<AliEMCALDigit *>(digits->At( fDigitsList[index] )) ;
679 Int_t nprimaries = digit->GetNprimary() ;
5c0368b8 680 if ( nprimaries == 0 ) continue ;
ab48128d 681 Int_t jndex ;
682 for ( jndex = 0 ; jndex < nprimaries ; jndex++ ) { // all primaries in digit
683 if ( fMulTrack > fMaxTrack ) {
f792c312 684 fMulTrack = fMaxTrack ;
9859bfc0 685 Error("GetNprimaries", "increase fMaxTrack ") ;
ab48128d 686 break ;
687 }
af5bdd85 688 Int_t newPrimary = digit->GetPrimary(jndex+1);
689 Float_t dEPrimary = digit->GetDEPrimary(jndex+1);
ab48128d 690 Int_t kndex ;
691 Bool_t already = kFALSE ;
692 for ( kndex = 0 ; kndex < fMulTrack ; kndex++ ) { //check if not already stored
af5bdd85 693 if ( newPrimary == primArray[kndex] ){
ab48128d 694 already = kTRUE ;
af5bdd85 695 dEPrimArray[kndex] += dEPrimary;
ab48128d 696 break ;
697 }
698 } // end of check
5c0368b8 699 if ( !already && (fMulTrack < fMaxTrack)) { // store it
af5bdd85 700 primArray[fMulTrack] = newPrimary ;
701 dEPrimArray[fMulTrack] = dEPrimary ;
ab48128d 702 fMulTrack++ ;
703 } // store it
704 } // all primaries in digit
ab48128d 705 } // all digits
706
af5bdd85 707 Int_t *sortIdx = new Int_t[fMulTrack];
708 TMath::Sort(fMulTrack,dEPrimArray,sortIdx);
709 for(index = 0; index < fMulTrack; index++) {
710 fTracksList[index] = primArray[sortIdx[index]] ;
711 fDETracksList[index] = dEPrimArray[sortIdx[index]] ;
712 }
713 delete [] sortIdx;
714 delete [] primArray ;
715 delete [] dEPrimArray ;
ab48128d 716
717}
7ee5c5be 718
87cdc3be 719//______________________________________________________________________________
720void AliEMCALRecPoint::EvalParents(TClonesArray * digits)
721{
722 // Constructs the list of parent particles (tracks) which have contributed to this RecPoint
723
724 AliEMCALDigit * digit ;
af5bdd85 725 Int_t * parentArray = new Int_t[fMaxTrack] ;
726 Float_t * dEParentArray = new Float_t[fMaxTrack] ;
87cdc3be 727
728 Int_t index ;
729 for ( index = 0 ; index < GetDigitsMultiplicity() ; index++ ) { // all digits
730 digit = dynamic_cast<AliEMCALDigit *>(digits->At( fDigitsList[index] )) ;
731 Int_t nparents = digit->GetNiparent() ;
5c0368b8 732 if ( nparents == 0 ) continue ;
87cdc3be 733
734 Int_t jndex ;
735 for ( jndex = 0 ; jndex < nparents ; jndex++ ) { // all primaries in digit
736 if ( fMulParent > fMaxParent ) {
737 fMulTrack = - 1 ;
738 Error("GetNiparent", "increase fMaxParent") ;
739 break ;
740 }
af5bdd85 741 Int_t newParent = digit->GetIparent(jndex+1) ;
742 Float_t newdEParent = digit->GetDEParent(jndex+1) ;
87cdc3be 743 Int_t kndex ;
744 Bool_t already = kFALSE ;
f1d429fd 745 for ( kndex = 0 ; kndex < fMulParent ; kndex++ ) { //check if not already stored
af5bdd85 746 if ( newParent == parentArray[kndex] ){
747 dEParentArray[kndex] += newdEParent;
87cdc3be 748 already = kTRUE ;
749 break ;
750 }
751 } // end of check
5c0368b8 752 if ( !already && (fMulTrack < fMaxTrack)) { // store it
af5bdd85 753 parentArray[fMulParent] = newParent ;
754 dEParentArray[fMulParent] = newdEParent ;
87cdc3be 755 fMulParent++ ;
756 } // store it
757 } // all parents in digit
87cdc3be 758 } // all digits
759
27e2a47c 760 if (fMulParent>0) {
af5bdd85 761 Int_t *sortIdx = new Int_t[fMulParent];
762 TMath::Sort(fMulParent,dEParentArray,sortIdx);
763 for(index = 0; index < fMulParent; index++) {
764 fParentsList[index] = parentArray[sortIdx[index]] ;
765 fDEParentsList[index] = dEParentArray[sortIdx[index]] ;
766 }
767 delete [] sortIdx;
27e2a47c 768 }
87cdc3be 769
af5bdd85 770 delete [] parentArray;
771 delete [] dEParentArray;
87cdc3be 772}
773
ab48128d 774//____________________________________________________________________________
70a93198 775void AliEMCALRecPoint::GetLocalPosition(TVector3 & lpos) const
776{
777 // returns the position of the cluster in the local reference system of ALICE
70a93198 778
779 lpos.SetX(fLocPos.X()) ;
780 lpos.SetY(fLocPos.Y()) ;
781 lpos.SetZ(fLocPos.Z()) ;
782}
783
784//____________________________________________________________________________
ab48128d 785void AliEMCALRecPoint::GetGlobalPosition(TVector3 & gpos) const
786{
787 // returns the position of the cluster in the global reference system of ALICE
70a93198 788 // These are now the Cartesian X, Y and Z
e52475ed 789 // cout<<" geom "<<geom<<endl;
14ce0a6e 790 fGeomPtr->GetGlobal(fLocPos, gpos, fSuperModuleNumber);
70a93198 791}
792
793//____________________________________________________________________________
794Float_t AliEMCALRecPoint::GetMaximalEnergy(void) const
795{
796 // Finds the maximum energy in the cluster
ab48128d 797
70a93198 798 Float_t menergy = 0. ;
799
800 Int_t iDigit;
801
802 for(iDigit=0; iDigit<fMulDigit; iDigit++) {
803
804 if(fEnergyList[iDigit] > menergy)
805 menergy = fEnergyList[iDigit] ;
806 }
807 return menergy ;
ab48128d 808}
809
aad8e277 810//____________________________________________________________________________
70a93198 811Int_t AliEMCALRecPoint::GetMultiplicityAtLevel(Float_t H) const
aad8e277 812{
70a93198 813 // Calculates the multiplicity of digits with energy larger than H*energy
814
815 Int_t multipl = 0 ;
816 Int_t iDigit ;
817 for(iDigit=0; iDigit<fMulDigit; iDigit++) {
818
819 if(fEnergyList[iDigit] > H * fAmp)
820 multipl++ ;
821 }
822 return multipl ;
823}
824
825//____________________________________________________________________________
826Int_t AliEMCALRecPoint::GetNumberOfLocalMax(AliEMCALDigit ** maxAt, Float_t * maxAtEnergy,
827 Float_t locMaxCut,TClonesArray * digits) const
828{
829 // Calculates the number of local maxima in the cluster using fLocalMaxCut as the minimum
830 // energy difference between two local maxima
831
832 AliEMCALDigit * digit ;
833 AliEMCALDigit * digitN ;
834
835 Int_t iDigitN ;
836 Int_t iDigit ;
837
838 for(iDigit = 0; iDigit < fMulDigit; iDigit++)
839 maxAt[iDigit] = (AliEMCALDigit*) digits->At(fDigitsList[iDigit]) ;
840
841 for(iDigit = 0 ; iDigit < fMulDigit; iDigit++) {
842 if(maxAt[iDigit]) {
843 digit = maxAt[iDigit] ;
844
845 for(iDigitN = 0; iDigitN < fMulDigit; iDigitN++) {
846 digitN = (AliEMCALDigit *) digits->At(fDigitsList[iDigitN]) ;
847
848 if ( AreNeighbours(digit, digitN) ) {
849 if (fEnergyList[iDigit] > fEnergyList[iDigitN] ) {
850 maxAt[iDigitN] = 0 ;
851 // but may be digit too is not local max ?
852 if(fEnergyList[iDigit] < fEnergyList[iDigitN] + locMaxCut)
853 maxAt[iDigit] = 0 ;
854 }
855 else {
856 maxAt[iDigit] = 0 ;
857 // but may be digitN too is not local max ?
858 if(fEnergyList[iDigit] > fEnergyList[iDigitN] - locMaxCut)
859 maxAt[iDigitN] = 0 ;
860 }
861 } // if Areneighbours
862 } // while digitN
863 } // slot not empty
864 } // while digit
865
866 iDigitN = 0 ;
867 for(iDigit = 0; iDigit < fMulDigit; iDigit++) {
868 if(maxAt[iDigit] ){
869 maxAt[iDigitN] = maxAt[iDigit] ;
870 maxAtEnergy[iDigitN] = fEnergyList[iDigit] ;
871 iDigitN++ ;
872 }
873 }
874 return iDigitN ;
875}
4635df1f 876
877//____________________________________________________________________________
878Int_t AliEMCALRecPoint::GetPrimaryIndex() const
879{
880 // Get the primary track index in TreeK which deposits the most energy
af5bdd85 881 // in Digits which forms RecPoint.
4635df1f 882
af5bdd85 883 if (fMulTrack)
884 return fTracksList[0];
885 return -12345;
4635df1f 886}
887
70a93198 888//____________________________________________________________________________
889void AliEMCALRecPoint::EvalTime(TClonesArray * digits){
890 // time is set to the time of the digit with the maximum energy
891
892 Float_t maxE = 0;
893 Int_t maxAt = 0;
894 for(Int_t idig=0; idig < fMulDigit; idig++){
895 if(fEnergyList[idig] > maxE){
896 maxE = fEnergyList[idig] ;
897 maxAt = idig;
898 }
899 }
900 fTime = ((AliEMCALDigit*) digits->At(fDigitsList[maxAt]))->GetTime() ;
aad8e277 901
aad8e277 902}
ab48128d 903
904//______________________________________________________________________________
905void AliEMCALRecPoint::Paint(Option_t *)
906{
907 // Paint this ALiRecPoint as a TMarker with its current attributes
908
909 TVector3 pos(0.,0.,0.) ;
910 GetLocalPosition(pos) ;
911 Coord_t x = pos.X() ;
912 Coord_t y = pos.Z() ;
913 Color_t markercolor = 1 ;
914 Size_t markersize = 1. ;
915 Style_t markerstyle = 5 ;
916
917 if (!gPad->IsBatch()) {
918 gVirtualX->SetMarkerColor(markercolor) ;
919 gVirtualX->SetMarkerSize (markersize) ;
920 gVirtualX->SetMarkerStyle(markerstyle) ;
921 }
922 gPad->SetAttMarkerPS(markercolor,markerstyle,markersize) ;
923 gPad->PaintPolyMarker(1,&x,&y,"") ;
924}
70a93198 925
926//______________________________________________________________________________
927Float_t AliEMCALRecPoint::EtaToTheta(Float_t arg) const
928{
929 //Converts Theta (Radians) to Eta(Radians)
930 return (2.*TMath::ATan(TMath::Exp(-arg)));
931}
932
933//______________________________________________________________________________
934Float_t AliEMCALRecPoint::ThetaToEta(Float_t arg) const
935{
936 //Converts Eta (Radians) to Theta(Radians)
937 return (-1 * TMath::Log(TMath::Tan(0.5 * arg)));
938}
261b1065 939
940//____________________________________________________________________________
941void AliEMCALRecPoint::Print(Option_t *) const
942{
943 // Print the list of digits belonging to the cluster
e52475ed 944 return;
261b1065 945 TString message ;
4800667c 946 message = "AliEMCALRecPoint:\n" ;
261b1065 947 message += " digits # = " ;
948 Info("Print", message.Data()) ;
949
950 Int_t iDigit;
951 for(iDigit=0; iDigit<fMulDigit; iDigit++)
952 printf(" %d ", fDigitsList[iDigit] ) ;
e52475ed 953 printf("\n");
954
261b1065 955 Info("Print", " Energies = ") ;
956 for(iDigit=0; iDigit<fMulDigit; iDigit++)
957 printf(" %f ", fEnergyList[iDigit] ) ;
e52475ed 958 printf("\n");
959
960 Info("Print", "\n Abs Ids = ") ;
961 for(iDigit=0; iDigit<fMulDigit; iDigit++)
962 printf(" %i ", fAbsIdList[iDigit] ) ;
963 printf("\n");
964
965 Info("Print", " Primaries ") ;
261b1065 966 for(iDigit = 0;iDigit < fMulTrack; iDigit++)
967 printf(" %d ", fTracksList[iDigit]) ;
e52475ed 968
969 printf("\n Local x %6.2f y %7.2f z %7.1f \n", fLocPos[0], fLocPos[1], fLocPos[2]);
970
85c60a8e 971 message = " ClusterType = %d" ;
972 message += " Multiplicity = %d" ;
261b1065 973 message += " Cluster Energy = %f" ;
974 message += " Core energy = %f" ;
975 message += " Core radius = %f" ;
976 message += " Number of primaries %d" ;
977 message += " Stored at position %d" ;
85c60a8e 978 Info("Print", message.Data(), fClusterType, fMulDigit, fAmp, fCoreEnergy, fCoreRadius, fMulTrack, GetIndexInList() ) ;
261b1065 979}
1d46d1f6 980
981Double_t AliEMCALRecPoint::GetPointEnergy() const
982{
983 static double e;
984 e=0.0;
985 for(int ic=0; ic<GetMultiplicity(); ic++) e += double(fEnergyList[ic]);
986 return e;
987}