--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+/* $Id$ */
+//_________________________________________________________________________
+// Base Class for EMCAL Reconstructed Points
+// Why should I put meaningless comments
+// just to satisfy
+// the code checker
+//*-- Author: Gines Martinez (SUBATECH)
+
+// --- ROOT system ---
+#include "TPad.h"
+#include "TClonesArray.h"
+
+// --- Standard library ---
+#include <iostream.h>
+#include <stdio.h>
+
+// --- AliRoot header files ---
+
+#include "AliEMCALGeometry.h"
+#include "AliEMCALDigit.h"
+#include "AliEMCALRecPoint.h"
+#include "AliEMCALGetter.h"
+
+ClassImp(AliEMCALRecPoint)
+
+
+//____________________________________________________________________________
+AliEMCALRecPoint::AliEMCALRecPoint()
+ : AliRecPoint()
+{
+ // ctor
+
+ fMaxTrack = 200 ;
+ fTheta = fPhi = 0. ;
+ fEMCALArm = 1;
+
+}
+
+//____________________________________________________________________________
+AliEMCALRecPoint::AliEMCALRecPoint(const char * opt) : AliRecPoint(opt)
+{
+ // ctor
+
+ fMaxTrack = 200 ;
+ fTheta = fPhi = 0. ;
+ fEMCALArm = 1;
+
+}
+
+//____________________________________________________________________________
+Int_t AliEMCALRecPoint::DistancetoPrimitive(Int_t px, Int_t py)
+{
+ // Compute distance from point px,py to a AliEMCALRecPoint considered as a Tmarker
+ // Compute the closest distance of approach from point px,py to this marker.
+ // The distance is computed in pixels units.
+
+ TVector3 pos(0.,0.,0.) ;
+ GetLocalPosition( pos) ;
+ Float_t x = pos.X() ;
+ Float_t y = pos.Z() ;
+ const Int_t kMaxDiff = 10;
+ Int_t pxm = gPad->XtoAbsPixel(x);
+ Int_t pym = gPad->YtoAbsPixel(y);
+ Int_t dist = (px-pxm)*(px-pxm) + (py-pym)*(py-pym);
+
+ if (dist > kMaxDiff) return 9999;
+ return dist;
+}
+
+//___________________________________________________________________________
+ void AliEMCALRecPoint::Draw(Option_t *option)
+ {
+ // Draw this AliEMCALRecPoint with its current attributes
+
+ AppendPad(option);
+ }
+
+//______________________________________________________________________________
+void AliEMCALRecPoint::ExecuteEvent(Int_t event, Int_t px, Int_t py)
+{
+ // Execute action corresponding to one event
+ // This member function is called when a AliEMCALRecPoint is clicked with the locator
+ //
+ // If Left button is clicked on AliEMCALRecPoint, the digits are switched on
+ // and switched off when the mouse button is released.
+
+ // static Int_t pxold, pyold;
+
+ static TGraph * digitgraph = 0 ;
+ static TPaveText* clustertext = 0 ;
+
+ if (!gPad->IsEditable()) return;
+
+ switch (event) {
+
+
+ case kButton1Down:{
+ AliEMCALDigit * digit ;
+ AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
+ AliEMCALGeometry * emcalgeom = const_cast<AliEMCALGeometry*>(gime->EMCALGeometry());
+
+ Int_t iDigit;
+ Int_t relid[4] ;
+
+ const Int_t kMulDigit=AliEMCALRecPoint::GetDigitsMultiplicity() ;
+ Float_t * xi = new Float_t [kMulDigit] ;
+ Float_t * zi = new Float_t [kMulDigit] ;
+
+ for(iDigit = 0; iDigit < kMulDigit; iDigit++) {
+ Fatal("AliEMCALRecPoint::ExecuteEvent", " -> Something wrong with the code");
+ digit = 0 ; //dynamic_cast<AliEMCALDigit *>((fDigitsList)[iDigit]);
+ emcalgeom->AbsToRelNumbering(digit->GetId(), relid) ;
+ emcalgeom->PosInAlice(relid, xi[iDigit], zi[iDigit]) ;
+ }
+
+ if (!digitgraph) {
+ digitgraph = new TGraph(fMulDigit,xi,zi);
+ digitgraph-> SetMarkerStyle(5) ;
+ digitgraph-> SetMarkerSize(1.) ;
+ digitgraph-> SetMarkerColor(1) ;
+ digitgraph-> Draw("P") ;
+ }
+ if (!clustertext) {
+
+ TVector3 pos(0.,0.,0.) ;
+ GetLocalPosition(pos) ;
+ clustertext = new TPaveText(pos.X()-10,pos.Z()+10,pos.X()+50,pos.Z()+35,"") ;
+ Text_t line1[40] ;
+ Text_t line2[40] ;
+ sprintf(line1,"Energy=%1.2f GeV",GetEnergy()) ;
+ sprintf(line2,"%d Digits",GetDigitsMultiplicity()) ;
+ clustertext ->AddText(line1) ;
+ clustertext ->AddText(line2) ;
+ clustertext ->Draw("");
+ }
+ gPad->Update() ;
+ Print() ;
+ delete[] xi ;
+ delete[] zi ;
+ }
+
+break;
+
+ case kButton1Up:
+ if (digitgraph) {
+ delete digitgraph ;
+ digitgraph = 0 ;
+ }
+ if (clustertext) {
+ delete clustertext ;
+ clustertext = 0 ;
+ }
+
+ break;
+
+ }
+}
+//____________________________________________________________________________
+void AliEMCALRecPoint::EvalAll(Float_t logWeight,TClonesArray * digits) {
+ //evaluates (if necessary) all RecPoint data members
+
+ EvalPrimaries(digits) ;
+}
+
+//____________________________________________________________________________
+void AliEMCALRecPoint::EvalEMCALArm(AliEMCALDigit * digit)
+{
+ // Returns the EMCAL module in which the RecPoint is found
+
+
+ if( fEMCALArm == 0){
+ Int_t relid[4] ;
+
+ AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
+ AliEMCALGeometry * emcalgeom = const_cast<AliEMCALGeometry*>(gime->EMCALGeometry());
+
+ emcalgeom->AbsToRelNumbering(digit->GetId(), relid) ;
+ fEMCALArm = relid[0];
+ }
+}
+
+//______________________________________________________________________________
+void AliEMCALRecPoint::EvalPrimaries(TClonesArray * digits)
+{
+ // Constructs the list of primary particles (tracks) which have contributed to this RecPoint
+
+ AliEMCALDigit * digit ;
+ Int_t * tempo = new Int_t[fMaxTrack] ;
+
+ Int_t index ;
+ for ( index = 0 ; index < GetDigitsMultiplicity() ; index++ ) { // all digits
+ digit = dynamic_cast<AliEMCALDigit *>(digits->At( fDigitsList[index] )) ;
+ Int_t nprimaries = digit->GetNprimary() ;
+ Int_t * newprimaryarray = new Int_t[nprimaries] ;
+ Int_t ii ;
+ for ( ii = 0 ; ii < nprimaries ; ii++)
+ newprimaryarray[ii] = digit->GetPrimary(ii+1) ;
+
+ Int_t jndex ;
+ for ( jndex = 0 ; jndex < nprimaries ; jndex++ ) { // all primaries in digit
+ if ( fMulTrack > fMaxTrack ) {
+ fMulTrack = - 1 ;
+ cout << "AliEMCALRecPoint::GetNprimaries ERROR > increase fMaxTrack " << endl ;
+ break ;
+ }
+ Int_t newprimary = newprimaryarray[jndex] ;
+ Int_t kndex ;
+ Bool_t already = kFALSE ;
+ for ( kndex = 0 ; kndex < fMulTrack ; kndex++ ) { //check if not already stored
+ if ( newprimary == tempo[kndex] ){
+ already = kTRUE ;
+ break ;
+ }
+ } // end of check
+ if ( !already) { // store it
+ tempo[fMulTrack] = newprimary ;
+ fMulTrack++ ;
+ } // store it
+ } // all primaries in digit
+ delete newprimaryarray ;
+ } // all digits
+
+
+ fTracksList = new Int_t[fMulTrack] ;
+ for(index = 0; index < fMulTrack; index++)
+ fTracksList[index] = tempo[index] ;
+
+ delete tempo ;
+
+}
+//____________________________________________________________________________
+void AliEMCALRecPoint::GetGlobalPosition(TVector3 & gpos) const
+{
+ // returns the position of the cluster in the global reference system of ALICE
+ // and the uncertainty on this position
+
+ AliEMCALGeometry * emcalgeom = AliEMCALGetter::GetInstance()->EMCALGeometry();
+ gpos.SetX(fPhi) ;
+ gpos.SetY(emcalgeom->GetIPDistance() + emcalgeom->GetAirGap()) ;
+ gpos.SetZ(fTheta) ;
+
+}
+
+
+//______________________________________________________________________________
+void AliEMCALRecPoint::Paint(Option_t *)
+{
+ // Paint this ALiRecPoint as a TMarker with its current attributes
+
+ TVector3 pos(0.,0.,0.) ;
+ GetLocalPosition(pos) ;
+ Coord_t x = pos.X() ;
+ Coord_t y = pos.Z() ;
+ Color_t markercolor = 1 ;
+ Size_t markersize = 1. ;
+ Style_t markerstyle = 5 ;
+
+ if (!gPad->IsBatch()) {
+ gVirtualX->SetMarkerColor(markercolor) ;
+ gVirtualX->SetMarkerSize (markersize) ;
+ gVirtualX->SetMarkerStyle(markerstyle) ;
+ }
+ gPad->SetAttMarkerPS(markercolor,markerstyle,markersize) ;
+ gPad->PaintPolyMarker(1,&x,&y,"") ;
+}
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+/* $Id$ */
+
+//_________________________________________________________________________
+// RecPoint implementation for EMCAL-EMC
+// An TowerRecPoint is a cluster of digits
+//*--
+//*-- Author: Dmitri Peressounko (RRC KI & SUBATECH)
+
+
+// --- ROOT system ---
+#include "TPad.h"
+#include "TH2.h"
+#include "TMath.h"
+#include "TCanvas.h"
+
+// --- Standard library ---
+
+#include <iostream.h>
+
+// --- AliRoot header files ---
+
+ #include "AliGenerator.h"
+#include "AliEMCALGeometry.h"
+#include "AliEMCALTowerRecPoint.h"
+#include "AliRun.h"
+#include "AliEMCALGetter.h"
+
+ClassImp(AliEMCALTowerRecPoint)
+
+//____________________________________________________________________________
+AliEMCALTowerRecPoint::AliEMCALTowerRecPoint() : AliEMCALRecPoint()
+{
+ // ctor
+
+ fMulDigit = 0 ;
+ fAmp = 0. ;
+ fCoreEnergy = 0 ;
+ fEnergyList = 0 ;
+ fTime = -1. ;
+ fLocPos.SetX(1000000.) ; //Local position should be evaluated
+
+}
+
+//____________________________________________________________________________
+AliEMCALTowerRecPoint::AliEMCALTowerRecPoint(const char * opt) : AliEMCALRecPoint(opt)
+{
+ // ctor
+
+ fMulDigit = 0 ;
+ fAmp = 0. ;
+ fCoreEnergy = 0 ;
+ fEnergyList = 0 ;
+ fTime = -1. ;
+ fLocPos.SetX(1000000.) ; //Local position should be evaluated
+
+}
+
+//____________________________________________________________________________
+AliEMCALTowerRecPoint::~AliEMCALTowerRecPoint()
+{
+ // dtor
+
+ if ( fEnergyList )
+ delete[] fEnergyList ;
+}
+
+//____________________________________________________________________________
+void AliEMCALTowerRecPoint::AddDigit(AliEMCALDigit & digit, Float_t Energy)
+{
+ // Adds a digit to the RecPoint
+ // and accumulates the total amplitude and the multiplicity
+
+ if(fEnergyList == 0)
+ fEnergyList = new Float_t[fMaxDigit];
+
+ if ( fMulDigit >= fMaxDigit ) { // increase the size of the lists
+ fMaxDigit*=2 ;
+ Int_t * tempo = new ( Int_t[fMaxDigit] ) ;
+ Float_t * tempoE = new ( Float_t[fMaxDigit] ) ;
+
+ Int_t index ;
+ for ( index = 0 ; index < fMulDigit ; index++ ){
+ tempo[index] = fDigitsList[index] ;
+ tempoE[index] = fEnergyList[index] ;
+ }
+
+ delete [] fDigitsList ;
+ fDigitsList = new ( Int_t[fMaxDigit] ) ;
+
+ delete [] fEnergyList ;
+ fEnergyList = new ( Float_t[fMaxDigit] ) ;
+
+ for ( index = 0 ; index < fMulDigit ; index++ ){
+ fDigitsList[index] = tempo[index] ;
+ fEnergyList[index] = tempoE[index] ;
+ }
+
+ delete [] tempo ;
+ delete [] tempoE ;
+ } // if
+
+ fDigitsList[fMulDigit] = digit.GetIndexInList() ;
+ fEnergyList[fMulDigit] = Energy ;
+ fMulDigit++ ;
+ fAmp += Energy ;
+
+ // EvalEMCALMod(&digit) ;
+}
+
+//____________________________________________________________________________
+Bool_t AliEMCALTowerRecPoint::AreNeighbours(AliEMCALDigit * digit1, AliEMCALDigit * digit2 ) const
+{
+ // Tells if (true) or not (false) two digits are neighbors
+
+ Bool_t aren = kFALSE ;
+
+ AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
+ AliEMCALGeometry * phosgeom = (AliEMCALGeometry*)gime->EMCALGeometry();
+
+ Int_t relid1[4] ;
+ phosgeom->AbsToRelNumbering(digit1->GetId(), relid1) ;
+
+ Int_t relid2[4] ;
+ phosgeom->AbsToRelNumbering(digit2->GetId(), relid2) ;
+
+ Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
+ Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
+
+ if (( coldiff <= 1 ) && ( rowdiff <= 1 ) && (coldiff + rowdiff > 0))
+ aren = kTRUE ;
+
+ return aren ;
+}
+
+//____________________________________________________________________________
+Int_t AliEMCALTowerRecPoint::Compare(const TObject * obj) const
+{
+ // Compares two RecPoints according to their position in the EMCAL modules
+
+ Float_t delta = 1 ; //Width of "Sorting row". If you changibg this
+ //value (what is senseless) change as vell delta in
+ //AliEMCALTrackSegmentMakerv* and other RecPoints...
+ Int_t rv ;
+
+ AliEMCALTowerRecPoint * clu = (AliEMCALTowerRecPoint *)obj ;
+
+
+ Int_t phosmod1 = GetEMCALArm() ;
+ Int_t phosmod2 = clu->GetEMCALArm() ;
+
+ TVector3 locpos1;
+ GetLocalPosition(locpos1) ;
+ TVector3 locpos2;
+ clu->GetLocalPosition(locpos2) ;
+
+ if(phosmod1 == phosmod2 ) {
+ Int_t rowdif = (Int_t)TMath::Ceil(locpos1.X()/delta)-(Int_t)TMath::Ceil(locpos2.X()/delta) ;
+ if (rowdif> 0)
+ rv = 1 ;
+ else if(rowdif < 0)
+ rv = -1 ;
+ else if(locpos1.Z()>locpos2.Z())
+ rv = -1 ;
+ else
+ rv = 1 ;
+ }
+
+ else {
+ if(phosmod1 < phosmod2 )
+ rv = -1 ;
+ else
+ rv = 1 ;
+ }
+
+ return rv ;
+}
+//______________________________________________________________________________
+void AliEMCALTowerRecPoint::ExecuteEvent(Int_t event, Int_t px, Int_t py) const
+{
+
+ // Execute action corresponding to one event
+ // This member function is called when a AliEMCALRecPoint is clicked with the locator
+ //
+ // If Left button is clicked on AliEMCALRecPoint, the digits are switched on
+ // and switched off when the mouse button is released.
+
+
+ // AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
+// if(!gime) return ;
+// AliEMCALGeometry * emcalgeom = (AliEMCALGeometry*)gime->EMCALGeometry();
+
+// static TGraph * digitgraph = 0 ;
+
+// if (!gPad->IsEditable()) return;
+
+// TH2F * histo = 0 ;
+// TCanvas * histocanvas ;
+
+// const TClonesArray * digits = gime->Digits() ;
+
+// switch (event) {
+
+// case kButton1Down: {
+// AliEMCALDigit * digit ;
+// Int_t iDigit;
+// Int_t relid[4] ;
+
+// const Int_t kMulDigit = AliEMCALTowerRecPoint::GetDigitsMultiplicity() ;
+// Float_t * xi = new Float_t[kMulDigit] ;
+// Float_t * zi = new Float_t[kMulDigit] ;
+
+// // create the histogram for the single cluster
+// // 1. gets histogram boundaries
+// Float_t ximax = -999. ;
+// Float_t zimax = -999. ;
+// Float_t ximin = 999. ;
+// Float_t zimin = 999. ;
+
+// for(iDigit=0; iDigit<kMulDigit; iDigit++) {
+// digit = (AliEMCALDigit *) digits->At(fDigitsList[iDigit]) ;
+// emcalgeom->AbsToRelNumbering(digit->GetId(), relid) ;
+// emcalgeom->RelPosInModule(relid, xi[iDigit], zi[iDigit]);
+// if ( xi[iDigit] > ximax )
+// ximax = xi[iDigit] ;
+// if ( xi[iDigit] < ximin )
+// ximin = xi[iDigit] ;
+// if ( zi[iDigit] > zimax )
+// zimax = zi[iDigit] ;
+// if ( zi[iDigit] < zimin )
+// zimin = zi[iDigit] ;
+// }
+// ximax += emcalgeom->GetCrystalSize(0) / 2. ;
+// zimax += emcalgeom->GetCrystalSize(2) / 2. ;
+// ximin -= emcalgeom->GetCrystalSize(0) / 2. ;
+// zimin -= emcalgeom->GetCrystalSize(2) / 2. ;
+// Int_t xdim = (int)( (ximax - ximin ) / emcalgeom->GetCrystalSize(0) + 0.5 ) ;
+// Int_t zdim = (int)( (zimax - zimin ) / emcalgeom->GetCrystalSize(2) + 0.5 ) ;
+
+// // 2. gets the histogram title
+
+// Text_t title[100] ;
+// sprintf(title,"Energy=%1.2f GeV ; Digits ; %d ", GetEnergy(), GetDigitsMultiplicity()) ;
+
+// if (!histo) {
+// delete histo ;
+// histo = 0 ;
+// }
+// histo = new TH2F("cluster3D", title, xdim, ximin, ximax, zdim, zimin, zimax) ;
+
+// Float_t x, z ;
+// for(iDigit=0; iDigit<kMulDigit; iDigit++) {
+// digit = (AliEMCALDigit *) digits->At(fDigitsList[iDigit]) ;
+// emcalgeom->AbsToRelNumbering(digit->GetId(), relid) ;
+// emcalgeom->RelPosInModule(relid, x, z);
+// histo->Fill(x, z, fEnergyList[iDigit] ) ;
+// }
+
+// if (!digitgraph) {
+// digitgraph = new TGraph(kMulDigit,xi,zi);
+// digitgraph-> SetMarkerStyle(5) ;
+// digitgraph-> SetMarkerSize(1.) ;
+// digitgraph-> SetMarkerColor(1) ;
+// digitgraph-> Paint("P") ;
+// }
+
+// // Print() ;
+// histocanvas = new TCanvas("cluster", "a single cluster", 600, 500) ;
+// histocanvas->Draw() ;
+// histo->Draw("lego1") ;
+
+// delete[] xi ;
+// delete[] zi ;
+
+// break;
+// }
+
+// case kButton1Up:
+// if (digitgraph) {
+// delete digitgraph ;
+// digitgraph = 0 ;
+// }
+// break;
+
+// }
+}
+
+//____________________________________________________________________________
+void AliEMCALTowerRecPoint::EvalDispersion(Float_t logWeight,TClonesArray * digits)
+{
+ // Calculates the dispersion of the shower at the origine of the RecPoint
+
+ Float_t d = 0. ;
+ Float_t wtot = 0. ;
+
+ AliEMCALDigit * digit ;
+
+ AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
+ AliEMCALGeometry * emcalgeom = (AliEMCALGeometry*)gime->EMCALGeometry();
+
+
+ // Calculates the center of gravity in the local EMCAL-module coordinates
+
+ Int_t iDigit;
+ Int_t relid[4] ;
+
+ if (!fTheta || !fPhi ) {
+ for(iDigit=0; iDigit<fMulDigit; iDigit++) {
+ digit = dynamic_cast<AliEMCALDigit *>(digits->At(fDigitsList[iDigit])) ;
+
+ Float_t thetai ;
+ Float_t phii ;
+ emcalgeom->AbsToRelNumbering(digit->GetId(), relid) ;
+ emcalgeom->PosInAlice(relid, thetai, phii);
+ Float_t w = TMath::Max( 0., logWeight + TMath::Log( fEnergyList[iDigit] / fAmp ) ) ;
+ fTheta = fTheta + thetai * w ;
+ fPhi += (phii * w );
+ wtot += w ;
+ }
+
+ fTheta /= wtot ;
+ fPhi /= wtot ;
+ }
+
+ const Float_t kDeg2Rad = TMath::Pi() / static_cast<Double_t>(180) ;
+
+ Float_t cyl_radius = emcalgeom->GetIPDistance()+emcalgeom->GetAirGap() ;
+ Float_t x = cyl_radius * TMath::Cos(fPhi * kDeg2Rad ) ;
+ Float_t y = cyl_radius * TMath::Cos(fPhi * kDeg2Rad ) ;
+ Float_t z = cyl_radius * TMath::Tan(fTheta * kDeg2Rad ) ;
+
+// Calculates the dispersion in coordinates
+ wtot = 0.;
+ for(iDigit=0; iDigit < fMulDigit; iDigit++) {
+ digit = (AliEMCALDigit *) digits->At(fDigitsList[iDigit]) ;
+ Float_t thetai = 0. ;
+ Float_t phii = 0.;
+ emcalgeom->AbsToRelNumbering(digit->GetId(), relid) ;
+ emcalgeom->PosInAlice(relid, thetai, phii);
+ Float_t xi = cyl_radius * TMath::Cos(phii * kDeg2Rad ) ;
+ Float_t yi = cyl_radius * TMath::Sin(phii * kDeg2Rad ) ;
+ Float_t zi = cyl_radius * TMath::Tan(thetai * kDeg2Rad ) ;
+
+ Float_t w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
+ d += w*((xi-x)*(xi-x) + (yi-y)*(yi-y)+ (zi-z)*(zi-z) ) ;
+ wtot+=w ;
+ }
+
+ d /= wtot ;
+
+ fDispersion = TMath::Sqrt(d) ;
+
+}
+//______________________________________________________________________________
+void AliEMCALTowerRecPoint::EvalCoreEnergy(Float_t logWeight, TClonesArray * digits)
+{
+ // This function calculates energy in the core,
+ // i.e. within a radius rad = 3cm around the center. Beyond this radius
+ // in accordance with shower profile the energy deposition
+ // should be less than 2%
+
+ Float_t coreRadius = 10. ;
+
+ AliEMCALDigit * digit ;
+ Int_t relid[4] ;
+ Float_t wtot = 0. ;
+
+ AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
+ AliEMCALGeometry * emcalgeom = (AliEMCALGeometry*)gime->EMCALGeometry();
+
+ Int_t iDigit;
+
+ if (!fTheta || !fPhi ) {
+ for(iDigit=0; iDigit<fMulDigit; iDigit++) {
+ digit = dynamic_cast<AliEMCALDigit *>(digits->At(fDigitsList[iDigit])) ;
+
+ Float_t thetai ;
+ Float_t phii ;
+ emcalgeom->AbsToRelNumbering(digit->GetId(), relid) ;
+ emcalgeom->PosInAlice(relid, thetai, phii);
+ Float_t w = TMath::Max( 0., logWeight + TMath::Log( fEnergyList[iDigit] / fAmp ) ) ;
+ fTheta = fTheta + thetai * w ;
+ fPhi += (phii * w );
+ wtot += w ;
+ }
+
+ fTheta /= wtot ;
+ fPhi /= wtot ;
+ }
+
+ const Float_t kDeg2Rad = TMath::Pi() / static_cast<Double_t>(180) ;
+
+ Float_t cyl_radius = emcalgeom->GetIPDistance()+emcalgeom->GetAirGap() ;
+ Float_t x = cyl_radius * TMath::Cos(fPhi * kDeg2Rad ) ;
+ Float_t y = cyl_radius * TMath::Cos(fPhi * kDeg2Rad ) ;
+ Float_t z = cyl_radius * TMath::Tan(fTheta * kDeg2Rad ) ;
+
+ for(iDigit=0; iDigit < fMulDigit; iDigit++) {
+ digit = (AliEMCALDigit *) ( digits->At(fDigitsList[iDigit]) ) ;
+ Int_t relid[4] ;
+ Float_t thetai = 0. ;
+ Float_t phii = 0. ;
+ emcalgeom->AbsToRelNumbering(digit->GetId(), relid) ;
+ emcalgeom->PosInAlice(relid, thetai, phii);
+
+ Float_t xi = cyl_radius * TMath::Cos(phii * kDeg2Rad ) ;
+ Float_t yi = cyl_radius * TMath::Sin(phii * kDeg2Rad ) ;
+ Float_t zi = cyl_radius * TMath::Tan(thetai * kDeg2Rad ) ;
+
+ Float_t distance = TMath::Sqrt((xi-x)*(xi-x)+(yi-y)*(yi-y)+(zi-z)*(zi-z)) ;
+ if(distance < coreRadius)
+ fCoreEnergy += fEnergyList[iDigit] ;
+ }
+
+}
+
+//____________________________________________________________________________
+void AliEMCALTowerRecPoint::EvalElipsAxis(Float_t logWeight,TClonesArray * digits)
+{
+ // Calculates the axis of the shower ellipsoid
+
+ Double_t wtot = 0. ;
+ Double_t x = 0.;
+ Double_t z = 0.;
+ Double_t dxx = 0.;
+ Double_t dzz = 0.;
+ Double_t dxz = 0.;
+
+ AliEMCALDigit * digit ;
+
+ AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
+ AliEMCALGeometry * emcalgeom = (AliEMCALGeometry*)gime->EMCALGeometry();
+
+ Int_t iDigit;
+ const Float_t kDeg2Rad = TMath::Pi() / static_cast<Double_t>(180) ;
+
+ Float_t cyl_radius = emcalgeom->GetIPDistance()+emcalgeom->GetAirGap() ;
+
+ for(iDigit=0; iDigit<fMulDigit; iDigit++) {
+ digit = (AliEMCALDigit *) digits->At(fDigitsList[iDigit]) ;
+ Int_t relid[4] ;
+ Float_t thetai = 0. ;
+ Float_t phii = 0. ;
+ emcalgeom->AbsToRelNumbering(digit->GetId(), relid) ;
+ emcalgeom->PosInAlice(relid, thetai, phii);
+ Double_t w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
+ Float_t xi = cyl_radius * TMath::Cos(fPhi * kDeg2Rad ) ;
+ Float_t zi = cyl_radius * TMath::Tan(fTheta * kDeg2Rad ) ;
+ dxx += w * xi * xi ;
+ x += w * xi ;
+ dzz += w * zi * zi ;
+ z += w * zi ;
+ dxz += w * xi * zi ;
+ wtot += w ;
+ }
+ dxx /= wtot ;
+ x /= wtot ;
+ dxx -= x * x ;
+ dzz /= wtot ;
+ z /= wtot ;
+ dzz -= z * z ;
+ dxz /= wtot ;
+ dxz -= x * z ;
+
+// //Apply correction due to non-perpendicular incidence
+// Double_t CosX ;
+// Double_t CosZ ;
+// AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
+// AliEMCALGeometry * emcalgeom = (AliEMCALGeometry*)gime->EMCALGeometry();
+ // Double_t DistanceToIP= (Double_t ) emcalgeom->GetIPDistance() ;
+
+// CosX = DistanceToIP/TMath::Sqrt(DistanceToIP*DistanceToIP+x*x) ;
+// CosZ = DistanceToIP/TMath::Sqrt(DistanceToIP*DistanceToIP+z*z) ;
+
+// dxx = dxx/(CosX*CosX) ;
+// dzz = dzz/(CosZ*CosZ) ;
+// dxz = dxz/(CosX*CosZ) ;
+
+
+ fLambda[0] = 0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ;
+ if(fLambda[0] > 0)
+ fLambda[0] = TMath::Sqrt(fLambda[0]) ;
+
+ fLambda[1] = 0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ;
+ if(fLambda[1] > 0) //To avoid exception if numerical errors lead to negative lambda.
+ fLambda[1] = TMath::Sqrt(fLambda[1]) ;
+ else
+ fLambda[1]= 0. ;
+}
+
+//____________________________________________________________________________
+void AliEMCALTowerRecPoint::EvalAll(Float_t logWeight, TClonesArray * digits )
+{
+ // Evaluates all shower parameters
+
+ AliEMCALRecPoint::EvalAll(logWeight,digits) ;
+ EvalGlobalPosition(logWeight, digits) ;
+ EvalElipsAxis(logWeight, digits) ;
+ EvalDispersion(logWeight, digits) ;
+ EvalCoreEnergy(logWeight, digits);
+ EvalTime(digits) ;
+}
+
+//____________________________________________________________________________
+void AliEMCALTowerRecPoint::EvalGlobalPosition(Float_t logWeight, TClonesArray * digits)
+{
+ // Calculates the center of gravity in the local EMCAL-module coordinates
+ Float_t wtot = 0. ;
+
+ Int_t relid[4] ;
+
+ AliEMCALDigit * digit ;
+
+ AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
+ AliEMCALGeometry * emcalgeom = static_cast<AliEMCALGeometry*>(gime->EMCALGeometry());
+ Int_t iDigit;
+
+ for(iDigit=0; iDigit<fMulDigit; iDigit++) {
+ digit = dynamic_cast<AliEMCALDigit *>(digits->At(fDigitsList[iDigit])) ;
+
+ Float_t thetai ;
+ Float_t phii ;
+ emcalgeom->AbsToRelNumbering(digit->GetId(), relid) ;
+ emcalgeom->PosInAlice(relid, thetai, phii);
+ Float_t w = TMath::Max( 0., logWeight + TMath::Log( fEnergyList[iDigit] / fAmp ) ) ;
+ fTheta = fTheta + thetai * w ;
+ fPhi += (phii * w );
+ wtot += w ;
+ }
+
+ fTheta /= wtot ;
+ fPhi /= wtot ;
+
+ fLocPos.SetX(0.) ;
+ fLocPos.SetY(0.) ;
+ fLocPos.SetZ(0.) ;
+
+ fLocPosM = 0 ;
+}
+
+//____________________________________________________________________________
+Float_t AliEMCALTowerRecPoint::GetMaximalEnergy(void) const
+{
+ // Finds the maximum energy in the cluster
+
+ Float_t menergy = 0. ;
+
+ Int_t iDigit;
+
+ for(iDigit=0; iDigit<fMulDigit; iDigit++) {
+
+ if(fEnergyList[iDigit] > menergy)
+ menergy = fEnergyList[iDigit] ;
+ }
+ return menergy ;
+}
+
+//____________________________________________________________________________
+Int_t AliEMCALTowerRecPoint::GetMultiplicityAtLevel(const Float_t H) const
+{
+ // Calculates the multiplicity of digits with energy larger than H*energy
+
+ Int_t multipl = 0 ;
+ Int_t iDigit ;
+ for(iDigit=0; iDigit<fMulDigit; iDigit++) {
+
+ if(fEnergyList[iDigit] > H * fAmp)
+ multipl++ ;
+ }
+ return multipl ;
+}
+
+//____________________________________________________________________________
+Int_t AliEMCALTowerRecPoint::GetNumberOfLocalMax(Int_t * maxAt, Float_t * maxAtEnergy,
+ Float_t locMaxCut,TClonesArray * digits) const
+{
+ // Calculates the number of local maxima in the cluster using fLocalMaxCut as the minimum
+ // energy difference between two local maxima
+
+ AliEMCALDigit * digit ;
+ AliEMCALDigit * digitN ;
+
+
+ Int_t iDigitN ;
+ Int_t iDigit ;
+
+ for(iDigit = 0; iDigit < fMulDigit; iDigit++)
+ maxAt[iDigit] = (Int_t) digits->At(fDigitsList[iDigit]) ;
+
+
+ for(iDigit = 0 ; iDigit < fMulDigit; iDigit++) {
+ if(maxAt[iDigit] != -1) {
+ digit = (AliEMCALDigit *) maxAt[iDigit] ;
+
+ for(iDigitN = 0; iDigitN < fMulDigit; iDigitN++) {
+ digitN = (AliEMCALDigit *) digits->At(fDigitsList[iDigitN]) ;
+
+ if ( AreNeighbours(digit, digitN) ) {
+ if (fEnergyList[iDigit] > fEnergyList[iDigitN] ) {
+ maxAt[iDigitN] = -1 ;
+ // but may be digit too is not local max ?
+ if(fEnergyList[iDigit] < fEnergyList[iDigitN] + locMaxCut)
+ maxAt[iDigit] = -1 ;
+ }
+ else {
+ maxAt[iDigit] = -1 ;
+ // but may be digitN too is not local max ?
+ if(fEnergyList[iDigit] > fEnergyList[iDigitN] - locMaxCut)
+ maxAt[iDigitN] = -1 ;
+ }
+ } // if Areneighbours
+ } // while digitN
+ } // slot not empty
+ } // while digit
+
+ iDigitN = 0 ;
+ for(iDigit = 0; iDigit < fMulDigit; iDigit++) {
+ if(maxAt[iDigit] != -1){
+ maxAt[iDigitN] = maxAt[iDigit] ;
+ maxAtEnergy[iDigitN] = fEnergyList[iDigit] ;
+ iDigitN++ ;
+ }
+ }
+ return iDigitN ;
+}
+//____________________________________________________________________________
+void AliEMCALTowerRecPoint::EvalTime(TClonesArray * digits){
+
+ Float_t maxE = 0;
+ Int_t maxAt = 0;
+ for(Int_t idig=0; idig < fMulDigit; idig++){
+ if(fEnergyList[idig] > maxE){
+ maxE = fEnergyList[idig] ;
+ maxAt = idig;
+ }
+ }
+ fTime = ((AliEMCALDigit*) digits->At(fDigitsList[maxAt]))->GetTime() ;
+
+}
+//____________________________________________________________________________
+void AliEMCALTowerRecPoint::Print(Option_t * option)
+{
+ // Print the list of digits belonging to the cluster
+
+ cout << "AliEMCALTowerRecPoint: " << endl ;
+
+ Int_t iDigit;
+ cout << " digits # = " ;
+ for(iDigit=0; iDigit<fMulDigit; iDigit++)
+ cout << fDigitsList[iDigit] << " " ;
+ cout << endl ;
+
+ cout << " Energies = " ;
+ for(iDigit=0; iDigit<fMulDigit; iDigit++)
+ cout << fEnergyList[iDigit] << " ";
+ cout << endl ;
+
+ cout << " Primaries " ;
+ for(iDigit = 0;iDigit < fMulTrack; iDigit++)
+ cout << fTracksList[iDigit] << " " << endl ;
+
+ cout << " Multiplicity = " << fMulDigit << endl ;
+ cout << " Cluster Energy = " << fAmp << endl ;
+ cout << " Number of primaries " << fMulTrack << endl ;
+ cout << " Stored at position " << GetIndexInList() << endl ;
+
+}
+
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff