--- /dev/null
+// --- ROOT system ---
+#include "TBenchmark.h"
+#include "TROOT.h"
+
+// --- Standard library ---
+#include <iostream.h>
+#include <iomanip.h>
+
+// --- AliRoot header files ---
+#include "AliPHOSClusterizerv2.h"
+#include "AliPHOSGetter.h"
+#include "TFolder.h"
+#include "AliPHOSEvalRecPoint.h"
+#include "AliPHOSRecCpvManager.h"
+#include "AliPHOSRecEmcManager.h"
+
+ClassImp(AliPHOSClusterizerv2)
+
+AliPHOSClusterizerv2::AliPHOSClusterizerv2() : AliPHOSClusterizerv1()
+{}
+
+AliPHOSClusterizerv2::AliPHOSClusterizerv2(const char* File, const char* name):
+ AliPHOSClusterizerv1(File,name)
+{}
+
+void AliPHOSClusterizerv2::GetNumberOfClustersFound(int* numb) const
+{
+
+ AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
+ numb[0] = gime->EmcRecPoints()->GetEntries();
+ numb[1] = gime->CpvRecPoints()->GetEntries();
+}
+
+void AliPHOSClusterizerv2::Exec(Option_t* option)
+{
+
+ if(strstr(option,"tim"))
+ gBenchmark->Start("PHOSClusterizer");
+
+ if(strstr(option,"print"))
+ Print("") ;
+
+ AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
+
+ TFolder* aliceF = (TFolder*)gROOT->FindObjectAny("YSAlice");
+ TFolder* storage = (TFolder*)aliceF->FindObject("WhiteBoard/RecPoints/PHOS");
+ TFolder* wPoolF = storage->AddFolder("SmP","SmartRecPoints for PHOS");
+
+ TObjArray* wPool = new TObjArray(400);
+ wPool->SetName("SmartPoints");
+ wPoolF->Add(wPool);
+ wPoolF->Add(this);
+
+ Int_t nevents = (Int_t) gAlice->TreeE()->GetEntries() ;
+ Int_t ievent ;
+
+ for(ievent = 0; ievent<nevents; ievent++) {
+
+ gAlice->GetEvent(ievent) ;
+ gAlice->SetEvent(ievent) ;
+
+ gime->Event(ievent,"D") ;
+// if(!ReadDigits(ievent)) //reads digits for event fEvent
+// continue;
+
+ cout<<" MakeClusters invoked..";
+ MakeClusters() ;
+ cout<<" done."<<endl;
+
+
+ //SmartRecPoints will communicate with wPool.
+
+ AliPHOSEvalRecPoint* rp=0;
+
+ // CPV reconstruction
+
+ AliPHOSRecCpvManager* recCpv = new AliPHOSRecCpvManager();
+ wPoolF->Add(recCpv);
+
+ for(Int_t iPoint=0; iPoint<gime->CpvRecPoints()->GetEntriesFast(); iPoint++) {
+ rp = new AliPHOSEvalRecPoint(iPoint,AliPHOSEvalRecPoint::cpv);
+ rp->MakeJob();
+ }
+
+ AliPHOSEvalRecPoint pt;
+ pt.UpdateWorkingPool();
+
+ TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
+ Int_t nOldCpv = cpvRecPoints->GetEntries();
+ cpvRecPoints->Delete();
+ cpvRecPoints->Compress();
+
+ for(Int_t i=0; i<wPool->GetEntries(); i++)
+ cpvRecPoints->Add(wPool->At(i));
+
+ wPool->Clear();
+ wPool->Compress();
+
+ wPoolF->Remove(recCpv);
+ delete recCpv;
+
+ cout<<" "<<gime->CpvRecPoints()->GetEntries()<<endl;
+ cout<<" "<<cpvRecPoints->GetEntries()<<" cpvRecPoints."<<endl<<endl;
+
+
+ // Now Emc reconstruction
+
+ AliPHOSRecEmcManager* recEmc = new AliPHOSRecEmcManager();
+ wPoolF->Add(recEmc);
+
+ for(Int_t iPoint=0; iPoint<gime->EmcRecPoints()->GetEntriesFast(); iPoint++) {
+ rp = new AliPHOSEvalRecPoint(iPoint,AliPHOSEvalRecPoint::emc);
+ rp->MakeJob();
+ }
+
+ pt.UpdateWorkingPool();
+
+ TObjArray * emcRecPoints = gime->EmcRecPoints() ;
+ Int_t nOldEmc = emcRecPoints->GetEntries();
+ emcRecPoints->Delete();
+ emcRecPoints->Compress();
+
+ for(Int_t i=0; i<wPool->GetEntries(); i++)
+ emcRecPoints->Add(wPool->At(i));
+
+ wPool->Clear();
+ wPool->Compress();
+
+ wPoolF->Remove(recEmc);
+ delete recEmc;
+
+ cout<<" "<<nOldCpv<<" OLD cpvRecPoints."<<endl;
+ cout<<" "<<gime->CpvRecPoints()->GetEntries()<<endl;
+ cout<<" "<<cpvRecPoints->GetEntries()<<" cpvRecPoints."<<endl<<endl;
+
+ cout<<" "<<nOldEmc<<" OLD emcRecPoints."<<endl;
+ cout<<" "<<gime->EmcRecPoints()->GetEntries()<<endl;
+ cout<<" "<<emcRecPoints->GetEntries()<<" emcRecPoints."<<endl<<endl;
+
+ WriteRecPoints(ievent);
+
+
+ } // loop over events
+
+ if(strstr(option,"tim")) {
+ gBenchmark->Stop("PHOSClusterizer");
+ cout << "AliPHOSClusterizer:" << endl ;
+ cout << " took " << gBenchmark->GetCpuTime("PHOSClusterizer") << " seconds for Clusterizing " << endl;
+ cout << endl ;
+
+ }
+
+}
+
+Int_t AliPHOSClusterizerv2::AreNeighbours(AliPHOSDigit* d1, AliPHOSDigit* d2) const
+{
+ // Points are neighbours if they have common edge.
+ // Points with common vertex are NOT neighbours.
+ // This treatment of neighbourship is the part of
+ // IHEP algorithm of clusterization.
+
+ // Gives the neighbourness of two digits = 0 are not neighbour but continue searching
+ // = 1 are neighbour
+ // = 2 are not neighbour but do not continue searching
+ // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
+ // which is compared to a digit (d2) not yet in a cluster
+
+ const AliPHOSGeometry * geom = AliPHOSGetter::GetInstance()->PHOSGeometry();
+
+ Int_t rv = 0 ;
+
+ Int_t relid1[4] ;
+ geom->AbsToRelNumbering(d1->GetId(), relid1) ;
+
+ Int_t relid2[4] ;
+ geom->AbsToRelNumbering(d2->GetId(), relid2) ;
+
+ if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module and the same PPSD Module
+ Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
+ Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
+
+ if ( ( (coldiff < 1) && (rowdiff <= 1) ) || ( ( coldiff <= 1 ) && ( rowdiff < 1 ) ) ){
+ rv = 1 ;
+ }
+ else {
+ if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1))
+ rv = 2; // Difference in row numbers is too large to look further
+ }
+
+ }
+ else {
+
+ if( (relid1[0] < relid2[0]) || (relid1[1] < relid2[1]) )
+ rv=2 ;
+
+ }
+
+// //Do NOT clusterize upper PPSD // YVK 30.09.2001
+// if( IsInPpsd(d1) && IsInPpsd(d2) &&
+// relid1[1] > 0 &&
+// relid1[1] < geom->GetNumberOfPadsPhi()*geom->GetNumberOfPadsPhi() ) rv = 2 ;
+
+ return rv ;
+
+}
--- /dev/null
+#ifndef AliPHOSClusterizerv2_H
+#define AliPHOSClusterizerv2_H
+
+// --- AliRoot header files ---
+#include "AliPHOSClusterizerv1.h"
+
+class AliPHOSClusterizerv2 : public AliPHOSClusterizerv1 {
+
+public:
+
+ AliPHOSClusterizerv2();
+ AliPHOSClusterizerv2(const char* headerFile, const char* name = 0);
+ ~AliPHOSClusterizerv2() {}
+
+ Int_t AreNeighbours(AliPHOSDigit* d1, AliPHOSDigit* d2) const ;
+ void GetNumberOfClustersFound(int* numb ) const;
+
+ void Exec(Option_t* option);
+
+ ClassDef(AliPHOSClusterizerv2,1)
+};
+
+#endif
--- /dev/null
+// ---- ROOT system ---
+#include "TDirectory.h"
+#include "TBranch.h"
+#include "TTree.h"
+#include "TROOT.h"
+#include "TFolder.h"
+
+// --- AliRoot header files ---
+#include "AliPHOSEvalRecPoint.h"
+#include "AliRun.h"
+#include "AliPHOSGetter.h"
+#include "AliPHOSRecCpvManager.h"
+#include "AliPHOSRecEmcManager.h"
+#include "AliPHOSDigitizer.h"
+
+// --- Standard library ---
+#include <iostream.h>
+
+
+ClassImp(AliPHOSEvalRecPoint)
+
+ AliPHOSEvalRecPoint::AliPHOSEvalRecPoint(): AliPHOSCpvRecPoint()
+{
+ fParent=-333;
+ fChi2Dof=-111;
+ fIsCpv = kTRUE;
+ fIsEmc = kFALSE;
+}
+
+AliPHOSEvalRecPoint::AliPHOSEvalRecPoint(Bool_t cpv, AliPHOSEvalRecPoint* parent) : AliPHOSCpvRecPoint()
+{
+
+ fParent=-333;
+ fChi2Dof=-111;
+
+// fParent=parent;
+ TObjArray* wPool = (TObjArray*)GetWorkingPool();
+ if(!wPool) {
+ cout<<" Couldn't find working pool. Exit."<<endl;
+ exit(1);
+ }
+
+ fParent = wPool->IndexOf((TObject*)parent);
+ fChi2Dof = parent->Chi2Dof();
+
+ if(cpv) {
+ fIsEmc = kFALSE;
+ fIsCpv = kTRUE;
+ }
+ else {
+ fIsEmc = kTRUE;
+ fIsCpv = kFALSE;
+ }
+
+ // Add itself to working pool
+ AddToWorkingPool((TObject*)this);
+
+}
+
+AliPHOSEvalRecPoint::AliPHOSEvalRecPoint(Int_t i, Bool_t cpv) : AliPHOSCpvRecPoint()
+{
+
+ fChi2Dof=-111;
+ fParent=-333;
+
+ AliPHOSEmcRecPoint* rp=0;
+
+ AliPHOSGetter* fGetter = AliPHOSGetter::GetInstance();
+
+ if(cpv) {
+ rp = (AliPHOSCpvRecPoint*)fGetter->CpvRecPoints()->At(i);
+ fIsEmc = kFALSE;
+ fIsCpv = kTRUE;
+ }
+ else {
+ rp = (AliPHOSEmcRecPoint*)fGetter->EmcRecPoints()->At(i);
+ fIsEmc = kTRUE;
+ fIsCpv = kFALSE;
+ }
+
+ Int_t* Digits = rp->GetDigitsList();
+ Float_t* Energies = rp->GetEnergiesList();
+ Int_t nDigits = rp->GetMultiplicity();
+
+ for(Int_t iDigit=0; iDigit<nDigits; iDigit++) {
+ AliPHOSDigit* digit = (AliPHOSDigit*)fGetter->Digits()->At( Digits[iDigit] );
+ Float_t eDigit = Energies[iDigit];
+ this->AddDigit(*digit,eDigit);
+ }
+
+ TVector3 locpos;
+ rp->GetLocalPosition(locpos);
+ fLocPos = locpos;
+
+ // Add itself to working pool
+ AddToWorkingPool((TObject*)this);
+
+}
+
+AliPHOSClusterizer* AliPHOSEvalRecPoint::GetClusterizer()
+{
+ TFolder* aliceF = (TFolder*)gROOT->FindObjectAny("YSAlice");
+ TFolder* wPoolF = (TFolder*)aliceF->FindObject("WhiteBoard/RecPoints/PHOS/SmP");
+ AliPHOSClusterizer* clu = (AliPHOSClusterizer*)wPoolF->FindObject("PHOS:clu-v1");
+ if(!clu) {
+ cout<<" Couldn't find Clusterizer. Exit."<<endl;
+ exit(1);
+ }
+
+ return clu;
+}
+
+Bool_t AliPHOSEvalRecPoint::TooClose(AliPHOSRecPoint* pt) const
+{
+ TVector3 her_pos;
+ TVector3 my_pos;
+ pt->GetLocalPosition(her_pos);
+ this->GetLocalPosition(my_pos);
+ Float_t dx = her_pos.X() - my_pos.X();
+ Float_t dz = her_pos.Z() - my_pos.Z();
+ Float_t dr = TMath::Sqrt(dx*dx + dz*dz);
+
+ if(dr<GetReconstructionManager()->MergeGammasMinDistanceCut())
+ return kTRUE;
+ else
+ return kFALSE;
+
+}
+
+Bool_t AliPHOSEvalRecPoint::NeedToSplit() const
+{
+ return kFALSE;
+}
+
+void AliPHOSEvalRecPoint::DeleteParent()
+{
+ fParent=-333;
+}
+
+void AliPHOSEvalRecPoint::UpdateWorkingPool()
+{
+
+ for(Int_t i=0; i<this->InWorkingPool(); i++) {
+ AliPHOSEvalRecPoint* parent = (AliPHOSEvalRecPoint*)GetFromWorkingPool(i);
+ TObjArray children;
+ Int_t nChild = parent->HasChild(children);
+ for(Int_t iChild=0; iChild<nChild; iChild++)
+ {
+ ((AliPHOSEvalRecPoint*)children.At(iChild))->DeleteParent();
+ }
+
+ if(nChild) {
+ RemoveFromWorkingPool(parent);
+ delete parent;
+ }
+
+ }
+
+ for(Int_t i=0; i<this->InWorkingPool(); i++) {
+ AliPHOSEvalRecPoint* weak = (AliPHOSEvalRecPoint*)GetFromWorkingPool(i);
+ if (weak->KillWeakPoint()) delete weak;
+ }
+
+ for(Int_t i=0; i<this->InWorkingPool(); i++) {
+ AliPHOSEvalRecPoint* close = (AliPHOSEvalRecPoint*)GetFromWorkingPool(i);
+ close->MergeClosePoint();
+ }
+
+ for(Int_t i=0; i<this->InWorkingPool(); i++)
+ ((AliPHOSEvalRecPoint*)AliPHOSEvalRecPoint::GetFromWorkingPool(i))->SetIndexInList(i);
+
+}
+
+Int_t AliPHOSEvalRecPoint::HasChild(TObjArray& children)
+{
+ for(Int_t iChild=0; iChild<InWorkingPool(); iChild++)
+ {
+ AliPHOSEvalRecPoint* child = (AliPHOSEvalRecPoint*)GetFromWorkingPool(iChild);
+ TObject* parent = (TObject*)child->Parent();
+ TObject* me = (TObject*)this;
+ if(parent==me) children.Add(child);
+ }
+
+ return children.GetEntries();
+}
+
+void AliPHOSEvalRecPoint::Init()
+{
+
+ AliPHOSClusterizer* clusterizer = GetClusterizer();
+ if(!clusterizer) {
+ cout<<" Cannot get clusterizer. Exit."<<endl;
+ exit(1);
+ }
+
+ TClonesArray* digits = AliPHOSGetter::GetInstance()->Digits();
+ Float_t LogWeight=0;
+
+ if(this->IsEmc()) {
+ LogWeight = clusterizer->GetEmcLogWeight();
+ }
+ else {
+ LogWeight = clusterizer->GetCpvLogWeight();
+ }
+
+ EvalLocalPosition(LogWeight,digits); // evaluate initial position
+}
+
+
+void AliPHOSEvalRecPoint::MakeJob()
+{
+ // Reconstruction algoritm implementation.
+
+ AliPHOSRecManager* recMng = GetReconstructionManager();
+
+ Init();
+
+ UnfoldLocalMaxima();
+
+ TObjArray children;
+ Int_t nChild = HasChild(children);
+
+ if(!nChild) {
+ EvaluatePosition();
+ if(Chi2Dof()>recMng->OneGamChisqCut()) SplitMergedShowers();
+ }
+
+ for(Int_t iChild=0; iChild<nChild; iChild++) {
+ AliPHOSEvalRecPoint* child = (AliPHOSEvalRecPoint*)children.At(iChild);
+ child->EvaluatePosition();
+
+ if(child->Chi2Dof()>recMng->OneGamChisqCut())
+ child->SplitMergedShowers();
+ }
+
+}
+
+void AliPHOSEvalRecPoint::InitTwoGam(Float_t* gamma1, Float_t* gamma2)
+{
+ //Compute start values for two gamma fit algorithm.
+ // gamma1[0], gamma1[1], gamma1[2] are Energy,x,z of the reconstructed "gammas".
+
+
+ TVector3 lpos; // start point choosing.
+ GetLocalPosition(lpos);
+ Float_t xx = lpos.Z();
+ Float_t yy = lpos.X();
+ Float_t E = GetEnergy();
+
+ cout<<" (x,z,E)[old] = ("<<yy<<","<<xx<<","<<E<<")"<<endl;
+
+// xx = XY(xx/E);
+// yy = XY(yy/E);
+
+
+ Float_t eDigit ;
+ AliPHOSDigit * digit ;
+ Int_t nDigits = GetMultiplicity();
+ Int_t * Digits = GetDigitsList();
+ Float_t * Energies = GetEnergiesList();
+
+ Float_t ix ;
+ Float_t iy ;
+ Int_t relid[4] ;
+
+ Float_t exx = 0;
+ Float_t eyy = 0;
+ Float_t exy = 0;
+ Float_t sqr;
+ Float_t cos2fi = 1.;
+
+ AliPHOSGetter* fGetter = AliPHOSGetter::GetInstance();
+ const AliPHOSGeometry* fGeom = fGetter->PHOSGeometry();
+
+ for(Int_t iDigit = 0 ; iDigit < nDigits ; iDigit ++)
+ {
+ digit = (AliPHOSDigit*)fGetter->Digits()->At( Digits[iDigit] );
+ eDigit = Energies[iDigit];
+ fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
+ fGeom->RelPosInModule(relid, iy, ix);
+
+ Float_t dx = ix - xx;
+ Float_t dy = iy - yy;
+ exx += eDigit*dx*dx;
+ eyy += eDigit*dy*dy;
+ exy += eDigit*dx*dy;
+
+ }
+
+ sqr = TMath::Sqrt(4.*exy*exy + (exx-eyy)*(exx-eyy));
+ Float_t euu = (exx+eyy+sqr)/2.;
+ if(sqr>1.e-10) cos2fi = (exx - eyy)/sqr;
+ Float_t cosfi = TMath::Sqrt((1.+cos2fi)/2.);
+ Float_t sinfi = TMath::Sqrt((1.-cos2fi)/2.);
+ if(exy<0) sinfi = -sinfi;
+
+ Float_t eu3 = 0;
+ for(Int_t iDigit = 0 ; iDigit < nDigits ; iDigit ++)
+ {
+ digit = (AliPHOSDigit*)fGetter->Digits()->At( Digits[iDigit] );
+ eDigit = Energies[iDigit];
+ fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
+ fGeom->RelPosInModule(relid, iy, ix);
+
+ Float_t dx = ix - xx;
+ Float_t dy = iy - yy;
+ Float_t du = dx*cosfi + dy*sinfi;
+ eu3 += eDigit*du*du*du;
+ }
+
+ Float_t c = E*eu3*eu3/(euu*euu*euu)/2.;
+ Float_t sign = 1.;
+ if(eu3<0) sign = -1.;
+ Float_t r = 1.+ c + sign*TMath::Sqrt(2.*c + c*c);
+ Float_t u = 0;
+ if(TMath::Abs(r-1.)>0.1) u = eu3/euu*(r+1.)/(r-1.);
+ if(TMath::Abs(r-1.)<0.1) u = TMath::Sqrt(sqr/E/r)*(1.+r);
+
+ Float_t e2c = E/(1.+r);
+ Float_t e1c = E-e2c;
+ Float_t u1 = -u/(1.+r);
+ Float_t u2 = u+u1;
+
+ Float_t x1c = xx+u1*cosfi;
+ Float_t y1c = yy+u1*sinfi;
+ Float_t x2c = xx+u2*cosfi;
+ Float_t y2c = yy+u2*sinfi;
+
+// printf("e1c -> %f\n",e1c);
+// printf("x1c -> %f\n",x1c);
+// printf("y1c -> %f\n",y1c);
+// printf("e2c -> %f\n",e2c);
+// printf("x2c -> %f\n",x2c);
+// printf("y2c -> %f\n",y2c);
+
+ gamma1[0] = e1c;
+ gamma1[1] = y1c;
+ gamma1[2] = x1c;
+
+ gamma2[0] = e2c;
+ gamma2[1] = y2c;
+ gamma2[2] = x2c;
+
+}
+
+void AliPHOSEvalRecPoint::TwoGam(Float_t* gamma1, Float_t* gamma2)
+{
+ //Fitting algorithm for the two very closed gammas
+ //that merged into the cluster with one maximum.
+ // Starting points gamma1 and gamma2 must be provided before ( see InitTwoGam)
+ //Set chisquare of the fit.
+
+
+ Float_t st0 = GetReconstructionManager()->TwoGamInitialStep();
+ Float_t chmin = GetReconstructionManager()->TwoGamChisqMin();
+ Float_t emin = GetReconstructionManager()->TwoGamEmin();
+ Float_t stpmin = GetReconstructionManager()->TwoGamStepMin();
+ Int_t Niter = GetReconstructionManager()->TwoGamNumOfIterations(); // Number of iterations.
+
+ Float_t chisq = 100.; //Initial chisquare.
+
+ Int_t nadc = GetMultiplicity();
+ if(nadc<3)
+ fChi2Dof= -111.;
+
+ Int_t dof = nadc - 5;
+ if(dof<1) dof=1;
+ Float_t chstop = chmin*dof;
+
+ Float_t ch = 1.e+20;
+ Float_t st = st0;
+ Float_t grx1 = 0.;
+ Float_t gry1 = 0.;
+ Float_t grx2 = 0.;
+ Float_t gry2 = 0.;
+ Float_t gre = 0.;
+ Float_t gr = 1.;
+ Float_t ee1=0,xx1=0,yy1=0,ee2=0,xx2=0,yy2=0,cosi;
+
+ Float_t e1c = gamma1[0];
+ Float_t y1c = gamma1[1];
+ Float_t x1c = gamma1[2];
+
+ Float_t e2c = gamma2[0];
+ Float_t y2c = gamma2[1];
+ Float_t x2c = gamma2[2];
+
+ Float_t E = GetEnergy();
+
+ Float_t eDigit ;
+ AliPHOSDigit * digit ;
+ Int_t nDigits = GetMultiplicity();
+ Int_t * Digits = GetDigitsList();
+ Float_t * Energies = GetEnergiesList();
+
+ Float_t ix ;
+ Float_t iy ;
+ Int_t relid[4] ;
+
+ AliPHOSGetter* fGetter = AliPHOSGetter::GetInstance();
+ const AliPHOSGeometry* fGeom = fGetter->PHOSGeometry();
+
+ for(Int_t Iter=0; Iter<Niter; Iter++)
+ {
+ Float_t chisqc = 0.;
+ Float_t grx1c = 0.;
+ Float_t gry1c = 0.;
+ Float_t grx2c = 0.;
+ Float_t gry2c = 0.;
+ Float_t grec = 0.;
+
+ for(Int_t iDigit = 0 ; iDigit < nDigits ; iDigit ++)
+ {
+ digit = (AliPHOSDigit*)fGetter->Digits()->At( Digits[iDigit] );
+ eDigit = Energies[iDigit];
+ fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
+ fGeom->RelPosInModule(relid, iy, ix);
+
+ Float_t a1,gx1,gy1;
+ Float_t a2,gx2,gy2;
+
+ Float_t dx1 = x1c - ix;
+ Float_t dy1 = y1c - iy;
+// cout<<" Mult "<<nDigits<<" dx1 "<<dx1<<" dy1 "<<dy1<<endl;
+// AG(e1c,dx1,dy1,a1,gx1,gy1);
+ GetReconstructionManager()->AG(e1c,dx1,dy1,a1,gx1,gy1);
+
+ Float_t dx2 = x2c - ix;
+ Float_t dy2 = y2c - iy;
+// cout<<" "<<" dx2 "<<dx2<<" dy2 "<<dy2<<endl;
+// AG(e2c,dx2,dy2,a2,gx2,gy2);
+ GetReconstructionManager()->AG(e2c,dx2,dy2,a2,gx2,gy2);
+
+ Float_t A = a1+a2;
+// Float_t D = Const*A*(1. - A/E);
+// if(D<0) D=0;
+
+// // D = 9.; // ????
+
+// Float_t da = A - eDigit;
+// chisqc += da*da/D;
+// Float_t dd = da/D;
+// dd = dd*(2.-dd*Const*(1.-2.*A/E));
+
+ Float_t dd;
+ chisqc += GetReconstructionManager()->TwoGamChi2(A,eDigit,E,dd);
+
+ grx1c += gx1*dd;
+ gry1c += gy1*dd;
+ grx2c += gx2*dd;
+ gry2c += gy2*dd;
+ grec += (a1/e1c - a2/e2c)*E*dd;
+
+ }
+
+ Float_t grc = TMath::Sqrt(grx1c*grx1c + gry1c*gry1c + grx2c*grx2c + gry2c*gry2c);
+ if(grc<1.e-10) grc=1.e-10;
+
+ Float_t sc = 1. + chisqc/ch;
+ st = st/sc;
+
+ if(chisqc>ch)
+ goto loop20;
+ cosi = (grx1*grx1c + gry1*gry1c + grx2*grx2c + gry2*gry2c + gre*grec)/gr/grc;
+ st = st*sc/(1.4 - cosi);
+
+ ee1 = e1c;
+ xx1 = x1c;
+ yy1 = y1c;
+ ee2 = e2c;
+ xx2 = x2c;
+ yy2 = y2c;
+
+ ch = chisqc;
+
+ if(ch<chstop)
+ goto loop101;
+
+ grx1 = grx1c;
+ gry1 = gry1c;
+ grx2 = grx2c;
+ gry2 = gry2c;
+ gre = grec;
+ gr = grc;
+
+ loop20: ;
+ Float_t step = st*gr;
+
+
+ cout<<" Iteration "<<Iter<<" dof "<<dof<<" chisq/dof "<<ch/dof<<" chstop/dof "<<chstop/dof<<" step " <<step<<" stpmin "<<stpmin<<endl;
+
+
+ if(step<stpmin)
+ goto loop101;
+
+ Float_t de = st*gre*E;
+ e1c = ee1 - de;
+ e2c = ee2 + de;
+
+ if( (e1c>emin) && (e2c>emin) )
+ goto loop25;
+ st = st/2.;
+ goto loop20;
+
+ loop25: ;
+ x1c = xx1 - st*grx1;
+ y1c = yy1 - st*gry1;
+ x2c = xx2 - st*grx2;
+ y2c = yy2 - st*gry2;
+
+
+ }
+
+ loop101:
+
+// if(ch>chisq*(nadc-2)-delch)
+// return ch/dof;
+
+ chisq = ch/dof;
+ gamma1[0] = ee1;
+ gamma1[1] = yy1;
+ gamma1[2] = xx1;
+
+ gamma2[0] = ee2;
+ gamma2[1] = yy2;
+ gamma2[2] = xx2;
+
+ Float_t x1_new = yy1;
+ Float_t z1_new = xx1;
+ Float_t e1_new = ee1;
+ Float_t x2_new = yy2;
+ Float_t z2_new = xx2;
+ Float_t e2_new = ee2;
+
+ cout<<" (x,z,E)[1 fit] = ("<<x1_new<<","<<z1_new<<","<<e1_new<<")"<<endl;
+ cout<<" (x,z,E)[2 fit] = ("<<x2_new<<","<<z2_new<<","<<e2_new<<")"<<endl;
+
+ fChi2Dof = chisq;
+
+}
+
+void AliPHOSEvalRecPoint::UnfoldTwoMergedPoints(Float_t* gamma1, Float_t* gamma2)
+{
+ //Unfold TWO merged rec. points in the case when cluster has only one maximum,
+ //but it's fitting to the one gamma shower is too bad.
+ // Use TwoGam() to estimate the positions and energies of merged points.
+
+
+ Int_t Nmax = 2;
+ Float_t* gamma;
+
+ Int_t* Digits = GetDigitsList();
+ Int_t nDigits = GetMultiplicity();
+ Float_t* Energies = GetEnergiesList();
+ Float_t* eFit = new Float_t[nDigits];
+
+ AliPHOSGetter* fGetter = AliPHOSGetter::GetInstance();
+ const AliPHOSGeometry* fGeom = fGetter->PHOSGeometry();
+
+ for(Int_t iDigit=0; iDigit<nDigits; iDigit++)
+ {
+ AliPHOSDigit* digit = (AliPHOSDigit*)fGetter->Digits()->At( Digits[iDigit] );
+ Int_t relid[4] ;
+ fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
+ Float_t x,z;
+ fGeom->RelPosInModule(relid, x, z);
+
+ Float_t gain = 0.;
+ for(Int_t iMax=0; iMax<Nmax; iMax++)
+ {
+ if(iMax==0)
+ gamma = gamma1;
+ else
+ gamma = gamma2;
+
+ Float_t eMax = gamma[0];
+ Float_t xMax = gamma[1];
+ Float_t zMax = gamma[2];
+
+ Float_t dx = xMax - x;
+ Float_t dz = zMax - z;
+ Float_t Amp,gx,gy;
+ GetReconstructionManager()->AG(eMax,dz,dx,Amp,gx,gy);
+ gain += Amp;
+ }
+
+ eFit[iDigit] = gain;
+ }
+
+
+ for( Int_t iMax=0; iMax<Nmax; iMax++)
+ {
+ if(iMax==0)
+ gamma = gamma1;
+ else
+ gamma = gamma2;
+
+ Float_t eMax = gamma[0];
+ Float_t xMax = gamma[1];
+ Float_t zMax = gamma[2];
+
+ AliPHOSEvalRecPoint* newRP = new AliPHOSEvalRecPoint(IsCPV(),this);
+ TVector3 newpos(xMax,0,zMax);
+ newRP->SetLocalPosition(newpos);
+
+ for( Int_t iDigit = 0 ; iDigit < nDigits ; iDigit ++)
+ {
+ AliPHOSDigit* digit = (AliPHOSDigit*)fGetter->Digits()->At( Digits[iDigit] );
+ Float_t eDigit = Energies[iDigit];
+ Int_t relid[4] ;
+ fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
+ Float_t ix,iz;
+ fGeom->RelPosInModule(relid, ix, iz);
+
+ Float_t dx = xMax - ix;
+ Float_t dz = zMax - iz;
+ Float_t single_shower_gain,gxMax,gyMax;
+ GetReconstructionManager()->AG(eMax,dz,dx,single_shower_gain,gxMax,gyMax);
+ Float_t total_gain = eFit[iDigit];
+ Float_t ratio = single_shower_gain/total_gain;
+// cout<<" ratio -> "<<ratio<<endl;
+ eDigit = eDigit*ratio;
+ newRP->AddDigit(*digit,eDigit);
+ }
+
+ cout<<"======= Split: daughter rec point "<<iMax<<" ================="<<endl;
+ newRP->Print("");
+
+ }
+
+ delete[] eFit;
+ cout<<endl;
+
+
+}
+
+
+void AliPHOSEvalRecPoint::EvaluatePosition()
+{
+ // One gamma fit algorithm.
+ // Set chisq/dof of the fit.
+ // gamma1[0], gamma1[1], gamma1[2] are Energy,x,z of the reconstructed gamma, respectively.
+
+ Int_t nDigits = GetMultiplicity();
+ if(nDigits<2)
+ return;
+
+ Int_t Niter = GetReconstructionManager()->OneGamNumOfIterations(); // number of iterations
+ Float_t St0 = GetReconstructionManager()->OneGamInitialStep(); // initial step
+// const Float_t Stpmin = 0.005;
+ Float_t Stpmin = GetReconstructionManager()->OneGamStepMin();
+ Float_t chmin = GetReconstructionManager()->OneGamChisqMin();
+
+ TVector3 locpos;
+ AliPHOSDigit* digit;
+ Float_t eDigit;
+ Int_t relid[4] ;
+ Float_t ix, iy;
+
+ GetLocalPosition(locpos);
+ Float_t E = GetEnergy();
+ Float_t xc = locpos.Z();
+ Float_t yc = locpos.X();
+ Float_t dx,dy,gx,gy,grxc,gryc;
+ Float_t st = St0;
+ Float_t chisq = 1.e+20;
+ Float_t gr = 1.;
+ Float_t grx = 0.;
+ Float_t gry = 0.;
+ Int_t dof = GetMultiplicity() - 2;
+ if(dof<1)
+ dof = 1;
+ Float_t chstop = chmin*dof;
+ Float_t cosi,x1=0,y1=0;
+ Float_t chisqc;
+
+ AliPHOSGetter* fGetter = AliPHOSGetter::GetInstance();
+ const AliPHOSGeometry* fGeom = fGetter->PHOSGeometry();
+
+ for(Int_t Iter=0; Iter<Niter; Iter++)
+ {
+
+ chisqc = 0.;
+ grxc = 0.;
+ gryc = 0.;
+ for(Int_t iDigit = 0 ; iDigit < nDigits ; iDigit ++)
+ {
+
+ Float_t* Energies = GetEnergiesList();
+ Int_t* Digits = GetDigitsList();
+ digit = (AliPHOSDigit*)fGetter->Digits()->At( Digits[iDigit] );
+ eDigit = Energies[iDigit];
+ fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
+ fGeom->RelPosInModule(relid, iy, ix);
+
+ dx = xc - ix;
+ dy = yc - iy;
+
+ if(!dx) dx=dy;
+ if(!dy) dy=dx;
+
+ Float_t A;
+ GetReconstructionManager()->AG(E,dx,dy,A,gx,gy);
+ Float_t dd;
+ cout<<" (ix iy xc yc dx dy) "<<ix<<" "<<iy<<" "<<xc<<" "<<yc<<" "<<dx<<" "<<dy<<endl;
+ Float_t chi2dg = GetReconstructionManager()->OneGamChi2(A,eDigit,E,dd);
+
+ // Exclude digit with too large chisquare.
+ if(chi2dg > 10) { continue; }
+
+ chisqc += chi2dg;
+ grxc += gx*dd;
+ gryc += gy*dd;
+
+ }
+
+ Float_t grc = TMath::Sqrt(grxc*grxc + gryc*gryc);
+ if(grc<1.e-10) grc=1.e-10;
+ Float_t sc = 1. + chisqc/chisq;
+ cout<<" chisq: "<<chisq<<endl;
+ st = st/sc;
+ if(chisqc>chisq)
+ goto loop20;
+ cosi = (grx*grxc + gry*gryc)/gr/grc;
+ st = st*sc/(1.4 - cosi);
+ x1 = xc;
+ y1 = yc;
+ grx = grxc;
+ gry = gryc;
+ gr = grc;
+ chisq = chisqc;
+
+ if(chisq<chstop)
+ goto loop101;
+
+ loop20: ;
+ Float_t step = st*gr;
+
+
+ cout<<" Iteration "<<Iter<<" dof "<<dof<<" chisq/dof "<<chisq/dof<<" chstop/dof "<<chstop/dof<<" step " <<step<<" Stpmin "<<Stpmin<<endl;
+
+
+ if(step<Stpmin)
+ goto loop101;
+
+ if(step>1.)
+ st = 1./gr;
+ xc = x1 - st*grx;
+ yc = y1 - st*gry;
+ }
+
+
+ loop101:
+ chisq = chisq/dof;
+// if(chisq>Chsqcut)
+// {
+// // TwoGam();
+// }
+
+ Float_t gamma1[3];
+ gamma1[0] = E;
+ gamma1[1] = y1;
+ gamma1[2] = x1;
+
+ TVector3 newpos(gamma1[1],0,gamma1[2]);
+ //SetLocalPosition(newpos);
+
+ fLocPos=newpos;
+ fAmp=E;
+
+// TVector3 pos;
+// RP->GetLocalPosition(pos);
+// cout<<" (x,z)[old] = ("<<x_old<<","<<z_old<<")"<<endl;
+// cout<<" (x,z)[new] = ("<<pos.X()<<","<<pos.Z()<<")"<<endl;
+
+
+// cout<<" gamma[1]: "<<gamma1[1]<<" gamma1[2]: "<<gamma1[2]<<endl;
+
+ fChi2Dof = chisq;
+
+}
+
+
+Bool_t AliPHOSEvalRecPoint::KillWeakPoint()
+{
+ // Remove this point from further procession
+ // if it's energy is too small.
+
+ Float_t Thr0 = GetReconstructionManager()->KillGamMinEnergy();
+
+ if(GetEnergy()<Thr0) {
+ cout<<"+++++++ Killing this rec point ++++++++++"<<endl;
+ RemoveFromWorkingPool(this);
+ return kTRUE;
+ }
+
+ return kFALSE;
+}
+
+
+void AliPHOSEvalRecPoint::SplitMergedShowers()
+{
+ // Split two very close rec. points.
+
+ if(GetMultiplicity()<2)
+ return;
+
+ Float_t gamma1[3];
+ Float_t gamma2[3];
+
+ InitTwoGam(gamma1,gamma2); // start points for fit
+ TwoGam(gamma1,gamma2); // evaluate the positions and energies
+ UnfoldTwoMergedPoints(gamma1,gamma2); // make the job
+
+}
+
+
+void AliPHOSEvalRecPoint::MergeClosePoint()
+{
+
+ AliPHOSGetter* fGetter = AliPHOSGetter::GetInstance();
+// AliPHOSDigitizer* digitizer = fGetter->Digitizer();
+// Float_t fPedestal = digitizer->GetPedestal(); // YVK 30.09.2001
+// Float_t fSlope = digitizer->GetSlope();
+
+ for (Int_t i=0;i<InWorkingPool(); i++)
+ {
+ TObject* obj = GetFromWorkingPool(i);
+ if(!(TObject*)this->IsEqual(obj))
+ {
+ AliPHOSRecPoint* rp = (AliPHOSRecPoint*)obj;
+ if(GetPHOSMod() == rp->GetPHOSMod())
+ {
+ if(TooClose(rp))
+ {
+ cout<<"+++++++ Merging point 1: ++++++"<<endl;
+ this->Print("");
+ cout<<"+++++++ and point 2: ++++++++++"<<endl;
+ ((AliPHOSEvalRecPoint*)rp)->Print("");
+
+ //merge two rec. points
+ TVector3 lpos1;
+ TVector3 lpos2;
+ this->GetLocalPosition(lpos1);
+ rp->GetLocalPosition(lpos2);
+
+ Int_t* Digits = rp->GetDigitsList();
+ Float_t dE = rp->GetEnergy()/(rp->GetEnergy()+this->GetEnergy());
+ Float_t new_x = lpos1.X()*dE + lpos2.X()*(1.-dE);
+ Float_t new_z = lpos1.Z()*dE + lpos2.Z()*(1.-dE);
+ Float_t new_E = rp->GetEnergy()+this->GetEnergy();
+ Float_t* Energies = ((AliPHOSEmcRecPoint*)rp)->GetEnergiesList();
+
+ for(Int_t iDigit=0; iDigit<rp->GetDigitsMultiplicity(); iDigit++)
+ {
+ AliPHOSDigit* digit = (AliPHOSDigit*)fGetter->Digits()->At(Digits[iDigit]);
+ Float_t eDigit = Energies[iDigit];
+ this->AddDigit(*digit,eDigit);
+ }
+
+ TVector3 newpos(new_x,0,new_z);
+ fLocPos = newpos;
+ fAmp = new_E;
+ RemoveFromWorkingPool(rp);
+ delete rp;
+
+ cout<<"++++++ Resulting point: ++++++++"<<endl;
+ this->Print("");
+
+ break;
+ }
+ }
+ }
+ }
+
+}
+
+Int_t AliPHOSEvalRecPoint::UnfoldLocalMaxima()
+{
+ // Make unfolding in the reconstruction point with several local maxima.
+ // Return the number of local maxima.
+
+ // if multiplicity less then 2 - nothing to unfold
+ if(GetMultiplicity()<2) return 1;
+
+ Int_t maxAt[1000];
+ Float_t maxAtEnergy[1000];
+ Float_t LocMaxCut, LogWeight;
+ Int_t relid[4] ;
+ Float_t xMax;
+ Float_t zMax;
+
+// AliPHOSClusterizer* clusterizer = fGetter->Clusterizer("PHOSclu-v1");
+ AliPHOSClusterizer* clusterizer = GetClusterizer();
+ if(!clusterizer) {
+ cout<<" Cannot get clusterizer. Exit."<<endl;
+ exit(1);
+ }
+
+ if(this->IsEmc()) {
+ LocMaxCut = clusterizer->GetEmcLocalMaxCut();
+ LogWeight = clusterizer->GetEmcLogWeight();
+ }
+ else {
+ LocMaxCut = clusterizer->GetCpvLocalMaxCut();
+ LogWeight = clusterizer->GetCpvLogWeight();
+ }
+
+ AliPHOSGetter* fGetter = AliPHOSGetter::GetInstance();
+ const AliPHOSGeometry* fGeom = fGetter->PHOSGeometry();
+ TClonesArray* digits = fGetter->Digits();
+
+ // if number of local maxima less then 2 - nothing to unfold
+ Int_t Nmax = GetNumberOfLocalMax(maxAt,maxAtEnergy,LocMaxCut,digits);
+ if(Nmax<2) return Nmax;
+
+ Int_t* Digits = GetDigitsList();
+ Int_t nDigits = GetMultiplicity();
+ Float_t* Energies = GetEnergiesList();
+ Float_t* eFit = new Float_t[nDigits];
+
+ for(Int_t iDigit=0; iDigit<nDigits; iDigit++)
+ {
+ eFit[iDigit] = 0.;
+ }
+
+ for(Int_t iDigit=0; iDigit<nDigits; iDigit++)
+ {
+
+ AliPHOSDigit* digit = (AliPHOSDigit*)fGetter->Digits()->At( Digits[iDigit] );
+ fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
+ Float_t x,z;
+ fGeom->RelPosInModule(relid, x, z);
+
+ for(Int_t iMax=0; iMax<Nmax; iMax++)
+ {
+ AliPHOSDigit* digitMax = (AliPHOSDigit*)maxAt[iMax];
+ Float_t eMax = maxAtEnergy[iMax];
+ fGeom->AbsToRelNumbering(digitMax->GetId(), relid) ;
+ fGeom->RelPosInModule(relid, xMax, zMax);
+ Float_t dx = xMax - x;
+ Float_t dz = zMax - z;
+ Float_t Amp,gx,gy;
+ GetReconstructionManager()->AG(eMax,dz,dx,Amp,gx,gy);
+// Amp = Amp + 0.5;
+ eFit[iDigit] += Amp;
+ }
+ }
+
+
+ for(Int_t iMax=0; iMax<Nmax; iMax++)
+ {
+ AliPHOSDigit* digitMax = (AliPHOSDigit*)maxAt[iMax];
+ fGeom->AbsToRelNumbering(digitMax->GetId(), relid) ;
+ fGeom->RelPosInModule(relid, xMax, zMax);
+ Float_t eMax = maxAtEnergy[iMax];
+
+ AliPHOSEvalRecPoint* newRP = new AliPHOSEvalRecPoint(IsCPV(),this);
+ newRP->AddDigit(*digitMax,maxAtEnergy[iMax]);
+
+ //Neighbous ( matrix 3x3 around the local maximum)
+ for(Int_t iDigit=0; iDigit<nDigits; iDigit++)
+ {
+ AliPHOSDigit* digit = (AliPHOSDigit*)fGetter->Digits()->At( Digits[iDigit] );
+ Float_t eDigit = Energies[iDigit];
+ fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
+ Float_t ix,iz;
+ fGeom->RelPosInModule(relid, ix, iz);
+
+ if(AreNeighbours(digitMax,digit))
+ {
+ Float_t dx = xMax - ix;
+ Float_t dz = zMax - iz;
+ Float_t single_shower_gain,gxMax,gyMax;
+ GetReconstructionManager()->AG(eMax,dz,dx,single_shower_gain,gxMax,gyMax);
+ Float_t total_gain = eFit[iDigit];
+ Float_t ratio = single_shower_gain/total_gain;
+ cout<<" ratio -> "<<ratio<<endl;
+ eDigit = eDigit*ratio;
+ newRP->AddDigit(*digit,eDigit);
+ }
+ }
+
+ newRP->EvalLocalPosition(LogWeight,digits);
+ cout<<"======= Unfold: daughter rec point "<<iMax<<" ================="<<endl;
+ newRP->Print("");
+ }
+
+// RemoveFromWorkingPool(this);
+
+ delete[] eFit;
+ cout<<endl;
+
+ return Nmax;
+
+}
+
+void AliPHOSEvalRecPoint::PrintPoint(Option_t* opt)
+{
+ AliPHOSCpvRecPoint::Print(opt);
+
+ TVector3 lpos;
+ GetLocalPosition(lpos);
+
+ cout<<" Chi2/dof = "<<Chi2Dof()<<endl;
+ cout<<" Local (x,z) = ("<<lpos.X()<<","<<lpos.Z()<<") in module "<<GetPHOSMod()<<endl;
+// if(GetReconstructionManager())
+// cout<<" Distance of merger = "<<GetReconstructionManager()->MergeGammasMinDistanceCut()<<endl;
+}
+
+AliPHOSRecManager* AliPHOSEvalRecPoint::GetReconstructionManager() const
+{
+
+ TFolder* aliceF = (TFolder*)gROOT->FindObjectAny("YSAlice");
+ TFolder* wPoolF = (TFolder*)aliceF->FindObject("WhiteBoard/RecPoints/PHOS/SmP");
+ AliPHOSRecManager* recMng = (AliPHOSRecManager*)wPoolF->FindObject("AliPHOSRecManager");
+ if(!recMng) {
+ cout<<" Couldn't find Reconstruction Manager. Exit."<<endl;
+ exit(1);
+ }
+
+ return recMng;
+}
+
+
+AliPHOSRecPoint* AliPHOSEvalRecPoint::Parent()
+{
+ if(fParent<0) return NULL;
+ else
+ return (AliPHOSRecPoint*)GetFromWorkingPool(fParent);
+// return fParent;
+}
+
+Float_t AliPHOSEvalRecPoint::Chi2Dof() const
+{
+ return fChi2Dof;
+}
+
+const TObject* AliPHOSEvalRecPoint::GetWorkingPool()
+{
+
+ TFolder* aliceF = (TFolder*)gROOT->FindObjectAny("YSAlice");
+ TFolder* wPoolF = (TFolder*)aliceF->FindObject("WhiteBoard/RecPoints/PHOS/SmP");
+ TObject* wPool = wPoolF->FindObject("SmartPoints");
+ if(!wPool) {
+ cout<<" Couldn't find Working Pool. Exit."<<endl;
+ exit(1);
+ }
+
+ return wPool;
+}
+
+void AliPHOSEvalRecPoint::AddToWorkingPool(TObject* obj)
+{
+ ((TObjArray*)GetWorkingPool())->Add(obj);
+}
+
+TObject* AliPHOSEvalRecPoint::GetFromWorkingPool(Int_t i)
+{
+// return fWorkingPool->At(i);
+ return ((TObjArray*)GetWorkingPool())->At(i);
+}
+
+Int_t AliPHOSEvalRecPoint::InWorkingPool()
+{
+ return ((TObjArray*)GetWorkingPool())->GetEntries();
+}
+
+void AliPHOSEvalRecPoint::RemoveFromWorkingPool(TObject* obj)
+{
+ ((TObjArray*)GetWorkingPool())->Remove(obj);
+ ((TObjArray*)GetWorkingPool())->Compress();
+}
+
+void AliPHOSEvalRecPoint::PrintWorkingPool()
+{
+ ((TObjArray*)GetWorkingPool())->Print("");
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
--- /dev/null
+#ifndef AliPHOSEvalRecPoint_H
+#define AliPHOSEvalRecPoint_H
+
+// --- ROOT system ---
+#include "TObjArray.h"
+
+
+#include "AliPHOSCpvRecPoint.h"
+#include "AliPHOSClusterizer.h"
+#include "AliPHOSRecPoint.h"
+#include "AliPHOSRecManager.h"
+
+class AliPHOSEvalRecPoint : public AliPHOSCpvRecPoint
+{
+public:
+
+ AliPHOSEvalRecPoint();
+ AliPHOSEvalRecPoint(Bool_t cpv,AliPHOSEvalRecPoint* parent);
+ AliPHOSEvalRecPoint(Int_t cluster, Bool_t cpv);
+ virtual ~AliPHOSEvalRecPoint() {}
+
+ Bool_t TooClose(AliPHOSRecPoint* pt) const ;
+ Bool_t NeedToSplit() const ;
+
+ void MergeClosePoint();
+ void SplitMergedShowers();
+ Int_t UnfoldLocalMaxima();
+ void EvaluatePosition();
+ Bool_t KillWeakPoint();
+
+ Bool_t IsEmc(void) const { return fIsEmc; }
+ Bool_t IsCPV(void) const { return fIsCpv; }
+
+ void SetLocalPosition(TVector3& pos) { fLocPos=pos; }
+ void UpdateWorkingPool();
+ void DeleteParent();
+
+ Int_t HasChild(TObjArray& children);
+
+ void MakeJob();
+
+ AliPHOSClusterizer* GetClusterizer();
+ AliPHOSRecManager* GetReconstructionManager() const;
+
+ void PrintPoint(Option_t* opt);
+
+ AliPHOSRecPoint* Parent();
+ Float_t Chi2Dof() const;
+
+ const TObject* GetWorkingPool();
+
+ void AddToWorkingPool(TObject* obj);
+ TObject* GetFromWorkingPool(Int_t i);
+ Int_t InWorkingPool();
+ void RemoveFromWorkingPool(TObject* obj);
+ void PrintWorkingPool();
+
+ enum RecPointType {emc,cpv};
+
+private:
+
+ void Init();
+ void InitTwoGam(Float_t* gamma1, Float_t* gamma2);
+ void TwoGam(Float_t* gamma1, Float_t* gamma2);
+ void UnfoldTwoMergedPoints(Float_t* gamma1, Float_t* gamma2);
+
+private:
+
+ Bool_t fIsEmc;
+ Bool_t fIsCpv;
+ Int_t fParent;
+ Float_t fChi2Dof;
+
+ ClassDef(AliPHOSEvalRecPoint,1)
+
+};
+
+#endif
+
+
+
+
+
+
+
+
+
+
+
+
+
--- /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. *
+ **************************************************************************/
+
+//_________________________________________________________________________
+// Algorythm class to analyze PHOSv1 events:
+// Construct histograms and displays them.
+// IHEP CPV/PHOS reconstruction algorithm used.
+// Use the macro EditorBar.C for best access to the functionnalities
+//*--
+//*-- Author: B. Polichtchouk (IHEP)
+//////////////////////////////////////////////////////////////////////////////
+
+// --- ROOT system ---
+
+#include "TFile.h"
+#include "TH1.h"
+#include "TPad.h"
+#include "TH2.h"
+#include "TParticle.h"
+#include "TClonesArray.h"
+#include "TTree.h"
+#include "TMath.h"
+#include "TCanvas.h"
+#include "TStyle.h"
+
+// --- Standard library ---
+
+#include <iostream.h>
+#include <stdio.h>
+
+// --- AliRoot header files ---
+
+#include "AliPHOSIhepAnalyze.h"
+#include "AliPHOSDigit.h"
+#include "AliPHOSRecParticle.h"
+#include "AliPHOSGetter.h"
+#include "AliPHOSHit.h"
+#include "AliPHOSImpact.h"
+#include "AliPHOSvImpacts.h"
+#include "AliPHOSCpvRecPoint.h"
+#include "AliRun.h"
+#include "AliPHOSGeometry.h"
+#include "AliPHOSEvalRecPoint.h"
+
+ClassImp(AliPHOSIhepAnalyze)
+
+//____________________________________________________________________________
+
+ AliPHOSIhepAnalyze::AliPHOSIhepAnalyze() {}
+
+//____________________________________________________________________________
+
+AliPHOSIhepAnalyze::AliPHOSIhepAnalyze(Text_t * name) : fFileName(name) {}
+
+//____________________________________________________________________________
+void AliPHOSIhepAnalyze::AnalyzeCPV1(Int_t Nevents)
+{
+ //
+ // Analyzes CPV characteristics: resolutions, cluster multiplicity,
+ // cluster lengths along Z and Phi.
+ // Author: Yuri Kharlov
+ // 9 October 2000
+ // Modified by Boris Polichtchouk, 3.07.2001
+ //
+
+ // Book histograms
+
+ TH1F *hDx = new TH1F("hDx" ,"CPV x-resolution@reconstruction",100,-5. , 5.);
+ TH1F *hDz = new TH1F("hDz" ,"CPV z-resolution@reconstruction",100,-5. , 5.);
+ TH1F *hChi2 = new TH1F("hChi2" ,"Chi2/dof of one-gamma fit",30, 0. , 10.);
+ TH1S *hNrp = new TH1S("hNrp" ,"CPV rec.point multiplicity", 21,-0.5,20.5);
+ TH1S *hNrpX = new TH1S("hNrpX","CPV rec.point Phi-length" , 21,-0.5,20.5);
+ TH1S *hNrpZ = new TH1S("hNrpZ","CPV rec.point Z-length" , 21,-0.5,20.5);
+
+ TH1F *hEg = new TH1F("hEg" ,"Energies of impacts",30,0.,6.);
+ TH1F *hEr = new TH1F("hEr" ,"Amplitudes of rec. points",50,0.,20.);
+
+ TList * fCpvImpacts ;
+ TBranch * branchCPVimpacts;
+
+ AliPHOSGetter * please = AliPHOSGetter::GetInstance(GetFileName().Data(),"PHOS");
+ const AliPHOSGeometry * fGeom = please->PHOSGeometry();
+
+ cout << "Start CPV Analysis-1. Resolutions, cluster multiplicity and lengths"<<endl;
+ for ( Int_t ievent=0; ievent<Nevents; ievent++) {
+
+ Int_t nTotalGen = 0;
+ Int_t nChargedGen = 0;
+
+ Int_t ntracks = gAlice->GetEvent(ievent);
+ cout<<" >>>>>>>Event "<<ievent<<".<<<<<<<"<<endl;
+
+ // Get branch of CPV impacts
+ if (! (branchCPVimpacts =gAlice->TreeH()->GetBranch("PHOSCpvImpacts")) ) {
+ cout<<" Couldn't find branch PHOSCpvImpacts. Exit."<<endl;
+ return;
+ }
+
+ // Create and fill arrays of hits for each CPV module
+
+ Int_t nOfModules = fGeom->GetNModules();
+ TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
+ Int_t iModule = 0;
+ for (iModule=0; iModule < nOfModules; iModule++)
+ hitsPerModule[iModule] = new TClonesArray("AliPHOSImpact",100);
+
+ TClonesArray *impacts;
+ AliPHOSImpact *impact;
+ TLorentzVector p;
+
+ // First go through all primary tracks and fill the arrays
+ // of hits per each CPV module
+
+ for (Int_t itrack=0; itrack<ntracks; itrack++) {
+ branchCPVimpacts ->SetAddress(&fCpvImpacts);
+ branchCPVimpacts ->GetEntry(itrack,0);
+
+ for (Int_t iModule=0; iModule < nOfModules; iModule++) {
+ impacts = (TClonesArray *)fCpvImpacts->At(iModule);
+ // Do loop over impacts in the module
+ for (Int_t iImpact=0; iImpact<impacts->GetEntries(); iImpact++) {
+ impact=(AliPHOSImpact*)impacts->At(iImpact);
+ hEg->Fill(impact->GetMomentum().E());
+ TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
+ if(IsCharged(impact->GetPid())) nChargedGen++;
+ new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSImpact(*impact);
+ }
+ }
+ fCpvImpacts->Clear();
+ }
+ for (iModule=0; iModule < nOfModules; iModule++) {
+ Int_t nsum = hitsPerModule[iModule]->GetEntriesFast();
+ printf("CPV module %d has %d impacts\n",iModule,nsum);
+ nTotalGen += nsum;
+ }
+
+ // Then go through reconstructed points and for each find
+ // the closeset hit
+ // The distance from the rec.point to the closest hit
+ // gives the coordinate resolution of the CPV
+
+ please->Event(ievent);
+ TIter nextRP(please->CpvRecPoints()) ;
+ AliPHOSCpvRecPoint *cpvRecPoint ;
+ Float_t xgen, ygen, zgen;
+ while( ( cpvRecPoint = (AliPHOSCpvRecPoint *)nextRP() ) ) {
+
+ Float_t chi2dof = ((AliPHOSEvalRecPoint*)cpvRecPoint)->Chi2Dof();
+ hChi2->Fill(chi2dof);
+ hEr->Fill(cpvRecPoint->GetEnergy());
+
+ TVector3 locpos;
+ cpvRecPoint->GetLocalPosition(locpos);
+ Int_t phosModule = cpvRecPoint->GetPHOSMod();
+ Int_t rpMult = cpvRecPoint->GetMultiplicity();
+ Int_t rpMultX, rpMultZ;
+ cpvRecPoint->GetClusterLengths(rpMultX,rpMultZ);
+ Float_t xrec = locpos.X();
+ Float_t zrec = locpos.Z();
+ Float_t dxmin = 1.e+10;
+ Float_t dzmin = 1.e+10;
+ Float_t r2min = 1.e+10;
+ Float_t r2;
+
+ Int_t nCPVhits = (hitsPerModule[phosModule-1])->GetEntriesFast();
+ Float_t locImpX=1e10,locImpZ=1e10; // local coords of closest impact
+ Float_t gImpX=1e10, gImpZ=1e10,gImpY=1e10; // global coords of closest impact
+ for (Int_t ihit=0; ihit<nCPVhits; ihit++) {
+ impact = (AliPHOSImpact*)(hitsPerModule[phosModule-1])->UncheckedAt(ihit);
+ xgen = impact->X();
+ zgen = impact->Z();
+ ygen = impact->Y();
+
+ //Transform to the local ref.frame
+ const AliPHOSGeometry* geom = please->PHOSGeometry();
+ Float_t phig = geom->GetPHOSAngle(phosModule);
+ Float_t phi = TMath::Pi()/180*phig;
+ Float_t distanceIPtoCPV = geom->GetIPtoOuterCoverDistance() -
+ (geom->GetFTPosition(1)+
+ geom->GetFTPosition(2)+
+ geom->GetCPVTextoliteThickness()
+ )/2;
+ Float_t xoL,yoL,zoL ;
+// xoL = xgen*TMath::Cos(phig)+ygen*TMath::Sin(phig) ;
+// yoL = -xgen*TMath::Sin(phig)+ygen*TMath::Cos(phig) + distanceIPtoCPV;
+ xoL = xgen*TMath::Cos(phi)-ygen*TMath::Sin(phi) ;
+ yoL = xgen*TMath::Sin(phi)+ygen*TMath::Cos(phi) + distanceIPtoCPV;
+ zoL = zgen;
+
+ r2 = TMath::Power((xoL-xrec),2) + TMath::Power((zoL-zrec),2);
+ if ( r2 < r2min ) {
+ r2min = r2;
+ dxmin = xoL - xrec;
+ dzmin = zoL - zrec;
+ locImpX = xoL;
+ locImpZ = zoL;
+ gImpX = xgen;
+ gImpZ = zgen;
+ gImpY = ygen;
+ }
+ }
+ cout<<" Impact global (X,Z,Y) = "<<gImpX<<" "<<gImpZ<<" "<<gImpY<<endl;
+ cout<<" Impact local (X,Z) = "<<locImpX<<" "<<locImpZ<<endl;
+ cout<<" Reconstructed (X,Z) = "<<xrec<<" "<<zrec<<endl;
+ cout<<" dxmin "<<dxmin<<" dzmin "<<dzmin<<endl<<endl;
+ hDx ->Fill(dxmin);
+ hDz ->Fill(dzmin);
+// hDr ->Fill(TMath::Sqrt(r2min));
+ hNrp ->Fill(rpMult);
+ hNrpX->Fill(rpMultX);
+ hNrpZ->Fill(rpMultZ);
+ }
+ delete [] hitsPerModule;
+
+ cout<<"++++ Event "<<ievent<<": total "<<nTotalGen<<" impacts, "
+ <<nChargedGen<<" charged impacts and "<<please->CpvRecPoints()->GetEntries()
+ <<" rec. points."<<endl<<endl;
+ }
+ // Save histograms
+
+ Text_t outputname[80] ;
+ sprintf(outputname,"%s.analyzed",GetFileName().Data());
+ TFile output(outputname,"RECREATE");
+ output.cd();
+
+ // Plot histograms
+
+ TCanvas *cpvCanvas = new TCanvas("Cpv1","CPV analysis-I",20,20,800,600);
+ gStyle->SetOptStat(111111);
+ gStyle->SetOptFit(1);
+ gStyle->SetOptDate(1);
+ cpvCanvas->Divide(3,3);
+
+ cpvCanvas->cd(1);
+ gPad->SetFillColor(10);
+ hNrp->SetFillColor(16);
+ hNrp->Draw();
+
+ cpvCanvas->cd(2);
+ gPad->SetFillColor(10);
+ hNrpX->SetFillColor(16);
+ hNrpX->Draw();
+
+ cpvCanvas->cd(3);
+ gPad->SetFillColor(10);
+ hNrpZ->SetFillColor(16);
+ hNrpZ->Draw();
+
+ cpvCanvas->cd(4);
+ gPad->SetFillColor(10);
+ hDx->SetFillColor(16);
+ hDx->Fit("gaus");
+ hDx->Draw();
+
+ cpvCanvas->cd(5);
+ gPad->SetFillColor(10);
+ hDz->SetFillColor(16);
+ hDz->Fit("gaus");
+ hDz->Draw();
+
+ cpvCanvas->cd(6);
+ gPad->SetFillColor(10);
+ hChi2->SetFillColor(16);
+ hChi2->Draw();
+
+ cpvCanvas->cd(7);
+ gPad->SetFillColor(10);
+ hEg->SetFillColor(16);
+ hEg->Draw();
+
+ cpvCanvas->cd(8);
+ gPad->SetFillColor(10);
+ hEr->SetFillColor(16);
+ hEr->Draw();
+
+ cpvCanvas->Write(0,kOverwrite);
+
+}
+
+
+void AliPHOSIhepAnalyze::AnalyzeEMC1(Int_t Nevents)
+{
+ //
+ // Analyzes Emc characteristics: resolutions, cluster multiplicity,
+ // cluster lengths along Z and Phi.
+ // Author: Boris Polichtchouk, 24.08.2001
+ //
+
+ // Book histograms
+
+ TH1F *hDx = new TH1F("hDx" ,"EMC x-resolution@reconstruction",100,-5. , 5.);
+ TH1F *hDz = new TH1F("hDz" ,"EMC z-resolution@reconstruction",100,-5. , 5.);
+ TH1F *hChi2 = new TH1F("hChi2" ,"Chi2/dof of one-gamma fit",30, 0. , 3.);
+ TH1S *hNrp = new TH1S("hNrp" ,"EMC rec.point multiplicity", 21,-0.5,20.5);
+ TH1S *hNrpX = new TH1S("hNrpX","EMC rec.point Phi-length" , 21,-0.5,20.5);
+ TH1S *hNrpZ = new TH1S("hNrpZ","EMC rec.point Z-length" , 21,-0.5,20.5);
+
+ TList * fEmcImpacts ;
+ TBranch * branchEMCimpacts;
+
+ AliPHOSGetter * please = AliPHOSGetter::GetInstance(GetFileName().Data(),"PHOS");
+ const AliPHOSGeometry * fGeom = please->PHOSGeometry();
+
+ cout << "Start EMC Analysis-1. Resolutions, cluster multiplicity and lengths"<<endl;
+ for ( Int_t ievent=0; ievent<Nevents; ievent++) {
+
+ Int_t nTotalGen = 0;
+
+ Int_t ntracks = gAlice->GetEvent(ievent);
+ cout<<" >>>>>>>Event "<<ievent<<".<<<<<<<"<<endl;
+
+ // Get branch of EMC impacts
+ if (! (branchEMCimpacts =gAlice->TreeH()->GetBranch("PHOSEmcImpacts")) ) {
+ cout<<" Couldn't find branch PHOSEmcImpacts. Exit."<<endl;
+ return;
+ }
+
+ // Create and fill arrays of hits for each EMC module
+
+ Int_t nOfModules = fGeom->GetNModules();
+ TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
+ Int_t iModule = 0;
+ for (iModule=0; iModule < nOfModules; iModule++)
+ hitsPerModule[iModule] = new TClonesArray("AliPHOSImpact",100);
+
+ TClonesArray *impacts;
+ AliPHOSImpact *impact;
+ TLorentzVector p;
+
+ // First go through all primary tracks and fill the arrays
+ // of hits per each EMC module
+
+ for (Int_t itrack=0; itrack<ntracks; itrack++) {
+ branchEMCimpacts ->SetAddress(&fEmcImpacts);
+ branchEMCimpacts ->GetEntry(itrack,0);
+
+ for (Int_t iModule=0; iModule < nOfModules; iModule++) {
+ impacts = (TClonesArray *)fEmcImpacts->At(iModule);
+ // Do loop over impacts in the module
+ for (Int_t iImpact=0; iImpact<impacts->GetEntries(); iImpact++) {
+ impact=(AliPHOSImpact*)impacts->At(iImpact);
+ TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
+ new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSImpact(*impact);
+ }
+ }
+ fEmcImpacts->Clear();
+ }
+ for (iModule=0; iModule < nOfModules; iModule++) {
+ Int_t nsum = hitsPerModule[iModule]->GetEntriesFast();
+ printf("EMC module %d has %d hits\n",iModule,nsum);
+ nTotalGen += nsum;
+ }
+
+ // Then go through reconstructed points and for each find
+ // the closeset hit
+ // The distance from the rec.point to the closest hit
+ // gives the coordinate resolution of the EMC
+
+ please->Event(ievent);
+ TIter nextRP(please->EmcRecPoints()) ;
+ AliPHOSEmcRecPoint *emcRecPoint ;
+ Float_t xgen, ygen, zgen;
+ while( ( emcRecPoint = (AliPHOSEmcRecPoint *)nextRP() ) ) {
+
+ Float_t chi2dof = ((AliPHOSEvalRecPoint*)emcRecPoint)->Chi2Dof();
+ hChi2->Fill(chi2dof);
+
+ TVector3 locpos;
+ emcRecPoint->GetLocalPosition(locpos);
+ Int_t phosModule = emcRecPoint->GetPHOSMod();
+ Int_t rpMult = emcRecPoint->GetMultiplicity();
+ Int_t rpMultX, rpMultZ;
+ ((AliPHOSEvalRecPoint*)emcRecPoint)->GetClusterLengths(rpMultX,rpMultZ);
+ Float_t xrec = locpos.X();
+ Float_t zrec = locpos.Z();
+ Float_t dxmin = 1.e+10;
+ Float_t dzmin = 1.e+10;
+ Float_t r2min = 1.e+10;
+ Float_t r2;
+
+ Int_t nEMChits = (hitsPerModule[phosModule-1])->GetEntriesFast();
+ Float_t locImpX=1e10,locImpZ=1e10; // local coords of closest impact
+ Float_t gImpX=1e10, gImpZ=1e10,gImpY=1e10; // global coords of closest impact
+ for (Int_t ihit=0; ihit<nEMChits; ihit++) {
+ impact = (AliPHOSImpact*)(hitsPerModule[phosModule-1])->UncheckedAt(ihit);
+ xgen = impact->X();
+ zgen = impact->Z();
+ ygen = impact->Y();
+
+
+ //Transform to the local ref.frame
+ const AliPHOSGeometry* geom = please->PHOSGeometry();
+ Float_t phig = geom->GetPHOSAngle(phosModule);
+ Float_t phi = TMath::Pi()/180*phig;
+ Float_t distanceIPtoEMC = geom->GetIPtoCrystalSurface();
+ Float_t xoL,yoL,zoL ;
+// xoL = xgen*TMath::Cos(phig)+ygen*TMath::Sin(phig) ;
+// yoL = -xgen*TMath::Sin(phig)+ygen*TMath::Cos(phig) + distanceIPtoEMC;
+ xoL = xgen*TMath::Cos(phi)-ygen*TMath::Sin(phi) ;
+ yoL = xgen*TMath::Sin(phi)+ygen*TMath::Cos(phi) + distanceIPtoEMC;
+ zoL = zgen;
+
+ r2 = TMath::Power((xoL-xrec),2) + TMath::Power((zoL-zrec),2);
+ if ( r2 < r2min ) {
+ r2min = r2;
+ dxmin = xoL - xrec;
+ dzmin = zoL - zrec;
+ locImpX = xoL;
+ locImpZ = zoL;
+ gImpX = xgen;
+ gImpZ = zgen;
+ gImpY = ygen;
+ }
+ }
+ cout<<" Impact global (X,Z,Y) = "<<gImpX<<" "<<gImpZ<<" "<<gImpY<<endl;
+ cout<<" Impact local (X,Z) = "<<locImpX<<" "<<locImpZ<<endl;
+ cout<<" Reconstructed (X,Z) = "<<xrec<<" "<<zrec<<endl;
+ cout<<" dxmin "<<dxmin<<" dzmin "<<dzmin<<endl<<endl;
+ hDx ->Fill(dxmin);
+ hDz ->Fill(dzmin);
+// hDr ->Fill(TMath::Sqrt(r2min));
+ hNrp ->Fill(rpMult);
+ hNrpX->Fill(rpMultX);
+ hNrpZ->Fill(rpMultZ);
+ }
+ delete [] hitsPerModule;
+
+ cout<<"++++ Event "<<ievent<<": total "<<nTotalGen<<" impacts, "
+ <<please->EmcRecPoints()->GetEntriesFast()<<" Emc rec. points."<<endl<<endl;
+
+ }
+ // Save histograms
+
+ Text_t outputname[80] ;
+ sprintf(outputname,"%s.analyzed",GetFileName().Data());
+ TFile output(outputname,"update");
+ output.cd();
+
+ // Plot histograms
+
+ TCanvas *emcCanvas = new TCanvas("Emc1","EMC analysis-I",20,20,800,400);
+ gStyle->SetOptStat(111111);
+ gStyle->SetOptFit(1);
+ gStyle->SetOptDate(1);
+ emcCanvas->Divide(3,2);
+
+ emcCanvas->cd(1);
+ gPad->SetFillColor(10);
+ hNrp->SetFillColor(16);
+ hNrp->Draw();
+
+ emcCanvas->cd(2);
+ gPad->SetFillColor(10);
+ hNrpX->SetFillColor(16);
+ hNrpX->Draw();
+
+ emcCanvas->cd(3);
+ gPad->SetFillColor(10);
+ hNrpZ->SetFillColor(16);
+ hNrpZ->Draw();
+
+ emcCanvas->cd(4);
+ gPad->SetFillColor(10);
+ hDx->SetFillColor(16);
+ hDx->Fit("gaus");
+ hDx->Draw();
+
+ emcCanvas->cd(5);
+ gPad->SetFillColor(10);
+ hDz->SetFillColor(16);
+ hDz->Fit("gaus");
+ hDz->Draw();
+
+ emcCanvas->cd(6);
+ gPad->SetFillColor(10);
+ hChi2->SetFillColor(16);
+ hChi2->Draw();
+
+ emcCanvas->Write(0,kOverwrite);
+}
+
+//____________________________________________________________________________
+void AliPHOSIhepAnalyze::AnalyzeCPV2(Int_t Nevents)
+{
+ // CPV analysis - part II.
+ // Ratio of combinatoric distances between generated
+ // and reconstructed hits.
+ // Author: Boris Polichtchouk (polishchuk@mx.ihep.su)
+ // 24 March 2001
+
+
+ TH1F* hDrij_cpv_r = new TH1F("Drij_cpv_r","Distance between reconstructed hits in CPV",140,0,50);
+ TH1F* hDrij_cpv_g = new TH1F("Drij_cpv_g","Distance between generated hits in CPV",140,0,50);
+ TH1F* hDrij_cpv_ratio = new TH1F("Drij_cpv_ratio","R_{ij}^{rec}/R_{ij}^{gen} in CPV",140,0,50);
+
+// TH1F* hT0 = new TH1F("hT0","Type of entering particle",20000,-10000,10000);
+
+ hDrij_cpv_r->Sumw2();
+ hDrij_cpv_g->Sumw2();
+ hDrij_cpv_ratio->Sumw2(); //correct treatment of errors
+
+ TList * fCpvImpacts = new TList();
+ TBranch * branchCPVimpacts;
+
+ AliPHOSGetter * please = AliPHOSGetter::GetInstance(GetFileName().Data(),"PHOS");
+ const AliPHOSGeometry * fGeom = please->PHOSGeometry();
+
+ for (Int_t nev=0; nev<Nevents; nev++)
+ {
+ printf("\n=================== Event %10d ===================\n",nev);
+ Int_t ntracks = gAlice->GetEvent(nev);
+ please->Event(nev);
+
+ Int_t nrec_cpv = 0; // Reconstructed points in event
+ Int_t ngen_cpv = 0; // Impacts in event
+
+ // Get branch of CPV impacts
+ if (! (branchCPVimpacts =gAlice->TreeH()->GetBranch("PHOSCpvImpacts")) ) return;
+
+ // Create and fill arrays of hits for each CPV module
+ Int_t nOfModules = fGeom->GetNModules();
+ TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
+ Int_t iModule = 0;
+ for (iModule=0; iModule < nOfModules; iModule++)
+ hitsPerModule[iModule] = new TClonesArray("AliPHOSImpact",100);
+
+ TClonesArray *impacts;
+ AliPHOSImpact *impact;
+
+ for (Int_t itrack=0; itrack<ntracks; itrack++) {
+ branchCPVimpacts ->SetAddress(&fCpvImpacts);
+ cout<<" branchCPVimpacts ->SetAddress(&fCpvImpacts) OK."<<endl;
+ branchCPVimpacts ->GetEntry(itrack,0);
+
+ for (Int_t iModule=0; iModule < nOfModules; iModule++) {
+ impacts = (TClonesArray *)fCpvImpacts->At(iModule);
+ // Do loop over impacts in the module
+ for (Int_t iImpact=0; iImpact<impacts->GetEntries(); iImpact++) {
+ impact=(AliPHOSImpact*)impacts->At(iImpact);
+ TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
+ if(IsCharged(impact->GetPid()))
+ new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSImpact(*impact);
+ }
+ }
+ fCpvImpacts->Clear();
+ }
+
+ for (iModule=0; iModule < nOfModules; iModule++) {
+ Int_t nsum = hitsPerModule[iModule]->GetEntriesFast();
+ printf("CPV module %d has %d hits\n",iModule,nsum);
+
+ AliPHOSImpact* GenHit1;
+ AliPHOSImpact* GenHit2;
+ Int_t irp1,irp2;
+ for(irp1=0; irp1< nsum; irp1++)
+ {
+ GenHit1 = (AliPHOSImpact*)((hitsPerModule[iModule])->At(irp1));
+ for(irp2 = irp1+1; irp2<nsum; irp2++)
+ {
+ GenHit2 = (AliPHOSImpact*)((hitsPerModule[iModule])->At(irp2));
+ Float_t dx = GenHit1->X() - GenHit2->X();
+ Float_t dz = GenHit1->Z() - GenHit2->Z();
+ Float_t dr = TMath::Sqrt(dx*dx + dz*dz);
+ hDrij_cpv_g->Fill(dr);
+// cout<<"(dx dz dr): "<<dx<<" "<<dz<<" "<<endl;
+ }
+ }
+ }
+
+
+ //--------- Combinatoric distance between rec. hits in CPV
+
+ TObjArray* cpvRecPoints = please->CpvRecPoints();
+ nrec_cpv = cpvRecPoints->GetEntriesFast();
+
+ if(nrec_cpv)
+ {
+ AliPHOSCpvRecPoint* RecHit1;
+ AliPHOSCpvRecPoint* RecHit2;
+ TIter next_cpv_rec1(cpvRecPoints);
+ while(TObject* obj1 = next_cpv_rec1() )
+ {
+ TIter next_cpv_rec2(cpvRecPoints);
+ while (TObject* obj2 = next_cpv_rec2())
+ {
+ if(!obj2->IsEqual(obj1))
+ {
+ RecHit1 = (AliPHOSCpvRecPoint*)obj1;
+ RecHit2 = (AliPHOSCpvRecPoint*)obj2;
+ TVector3 locpos1;
+ TVector3 locpos2;
+ RecHit1->GetLocalPosition(locpos1);
+ RecHit2->GetLocalPosition(locpos2);
+ Float_t dx = locpos1.X() - locpos2.X();
+ Float_t dz = locpos1.Z() - locpos2.Z();
+ Float_t dr = TMath::Sqrt(dx*dx + dz*dz);
+ if(RecHit1->GetPHOSMod() == RecHit2->GetPHOSMod())
+ hDrij_cpv_r->Fill(dr);
+ }
+ }
+ }
+ }
+
+ cout<<" Event "<<nev<<". Total of "<<ngen_cpv<<" hits, "<<nrec_cpv<<" rec.points."<<endl;
+
+ delete [] hitsPerModule;
+
+ } // End of loop over events.
+
+
+// hDrij_cpv_g->Draw();
+// hDrij_cpv_r->Draw();
+ hDrij_cpv_ratio->Divide(hDrij_cpv_r,hDrij_cpv_g);
+ hDrij_cpv_ratio->Draw();
+
+// hT0->Draw();
+
+}
+
+
+void AliPHOSIhepAnalyze::CpvSingle(Int_t nevents)
+{
+ // Distributions of coordinates and cluster lengths of rec. points
+ // in the case of single initial particle.
+
+ TH1F* hZr = new TH1F("Zrec","Reconstructed Z distribution",150,-5,5);
+ TH1F* hXr = new TH1F("Xrec","Reconstructed X distribution",150,-14,-2);
+ TH1F *hChi2 = new TH1F("hChi2" ,"Chi2/dof of one-gamma fit",100, 0. , 20.);
+
+ TH1S *hNrp = new TH1S("hNrp" ,"CPV rec.point multiplicity",21,-0.5,20.5);
+ TH1S *hNrpX = new TH1S("hNrpX","CPV rec.point Phi-length" ,21,-0.5,20.5);
+ TH1S *hNrpZ = new TH1S("hNrpZ","CPV rec.point Z-length" ,21,-0.5,20.5);
+
+ AliPHOSGetter* gime = AliPHOSGetter::GetInstance(GetFileName().Data(),"PHOS");
+
+ for(Int_t ievent=0; ievent<nevents; ievent++)
+ {
+ gime->Event(ievent);
+ if(gime->CpvRecPoints()->GetEntriesFast()>1) continue;
+
+ AliPHOSCpvRecPoint* pt = (AliPHOSCpvRecPoint*)(gime->CpvRecPoints())->At(0);
+ if(pt) {
+ TVector3 lpos;
+ pt->GetLocalPosition(lpos);
+ hXr->Fill(lpos.X());
+ hZr->Fill(lpos.Z());
+
+ Int_t rpMult = pt->GetMultiplicity();
+ hNrp->Fill(rpMult);
+ Int_t rpMultX, rpMultZ;
+ pt->GetClusterLengths(rpMultX,rpMultZ);
+ hNrpX->Fill(rpMultX);
+ hNrpZ->Fill(rpMultZ);
+ hChi2->Fill(((AliPHOSEvalRecPoint*)pt)->Chi2Dof());
+ cout<<"+++++ Event "<<ievent<<". (Mult,MultX,MultZ) = "<<rpMult<<" "<<rpMultX<<" "<<rpMultZ<<"+++++"<<endl;
+
+ }
+
+ }
+
+ Text_t outputname[80] ;
+ sprintf(outputname,"%s.analyzed.single",GetFileName().Data());
+ TFile output(outputname,"RECREATE");
+ output.cd();
+
+ // Plot histograms
+ TCanvas *cpvCanvas = new TCanvas("SingleParticle","Single particle events",20,20,800,400);
+ gStyle->SetOptStat(111111);
+ gStyle->SetOptFit(1);
+ gStyle->SetOptDate(1);
+ cpvCanvas->Divide(3,2);
+
+ cpvCanvas->cd(1);
+ gPad->SetFillColor(10);
+ hXr->SetFillColor(16);
+ hXr->Draw();
+
+ cpvCanvas->cd(2);
+ gPad->SetFillColor(10);
+ hZr->SetFillColor(16);
+ hZr->Draw();
+
+ cpvCanvas->cd(3);
+ gPad->SetFillColor(10);
+ hChi2->SetFillColor(16);
+ hChi2->Draw();
+
+ cpvCanvas->cd(4);
+ gPad->SetFillColor(10);
+ hNrp->SetFillColor(16);
+ hNrp->Draw();
+
+ cpvCanvas->cd(5);
+ gPad->SetFillColor(10);
+ hNrpX->SetFillColor(16);
+ hNrpX->Draw();
+
+ cpvCanvas->cd(6);
+ gPad->SetFillColor(10);
+ hNrpZ->SetFillColor(16);
+ hNrpZ->Draw();
+
+ cpvCanvas->Write(0,kOverwrite);
+
+}
+
+void AliPHOSIhepAnalyze::HitsCPV(TClonesArray& hits, Int_t nev)
+{
+ // Cumulative list of charged CPV impacts in event nev.
+
+ TList * fCpvImpacts ;
+ TBranch * branchCPVimpacts;
+
+ AliPHOSGetter * please = AliPHOSGetter::GetInstance(GetFileName().Data(),"PHOS");
+ const AliPHOSGeometry * fGeom = please->PHOSGeometry();
+
+
+ printf("\n=================== Event %10d ===================\n",nev);
+ Int_t ntracks = gAlice->GetEvent(nev);
+ please->Event(nev);
+
+// Int_t nrec_cpv = 0; // Reconstructed points in event // 01.10.2001
+// Int_t ngen_cpv = 0; // Impacts in event
+
+ // Get branch of CPV impacts
+ if (! (branchCPVimpacts =gAlice->TreeH()->GetBranch("PHOSCpvImpacts")) ) return;
+
+ // Create and fill arrays of hits for each CPV module
+ Int_t nOfModules = fGeom->GetNModules();
+ TClonesArray **hitsPerModule = new TClonesArray *[nOfModules];
+ Int_t iModule = 0;
+ for (iModule=0; iModule < nOfModules; iModule++)
+ hitsPerModule[iModule] = new TClonesArray("AliPHOSImpact",100);
+
+ TClonesArray *impacts;
+ AliPHOSImpact *impact;
+
+ for (Int_t itrack=0; itrack<ntracks; itrack++) {
+ branchCPVimpacts ->SetAddress(&fCpvImpacts);
+ cout<<" branchCPVimpacts ->SetAddress(&fCpvImpacts) OK."<<endl;
+ branchCPVimpacts ->GetEntry(itrack,0);
+
+ for (Int_t iModule=0; iModule < nOfModules; iModule++) {
+ impacts = (TClonesArray *)fCpvImpacts->At(iModule);
+ // Do loop over impacts in the module
+ for (Int_t iImpact=0; iImpact<impacts->GetEntries(); iImpact++) {
+ impact=(AliPHOSImpact*)impacts->At(iImpact);
+ TClonesArray &lhits = *(TClonesArray *)hitsPerModule[iModule];
+ new(lhits[hitsPerModule[iModule]->GetEntriesFast()]) AliPHOSImpact(*impact);
+ }
+ }
+ fCpvImpacts->Clear();
+ }
+
+ for (iModule=0; iModule < nOfModules; iModule++) {
+ Int_t nsum = hitsPerModule[iModule]->GetEntriesFast();
+ printf("CPV module %d has %d hits\n",iModule,nsum);
+ }
+
+// TList * fCpvImpacts ;
+// TBranch * branchCPVimpacts;
+// AliPHOSImpact* impact;
+// TClonesArray* impacts;
+
+// AliPHOSGetter * please = AliPHOSGetter::GetInstance(GetFileName().Data(),"PHOS");
+// const AliPHOSGeometry * fGeom = please->PHOSGeometry();
+
+// Int_t ntracks = gAlice->GetEvent(ievent);
+// Int_t nOfModules = fGeom->GetNModules();
+// cout<<" Tracks: "<<ntracks<<" Modules: "<<nOfModules<<endl;
+
+// if (! (branchCPVimpacts =gAlice->TreeH()->GetBranch("PHOSCpvImpacts")) ) return;
+
+// for (Int_t itrack=0; itrack<ntracks; itrack++) {
+// branchCPVimpacts ->SetAddress(&fCpvImpacts);
+// cout<<" branchCPVimpacts ->SetAddress(&fCpvImpacts) OK."<<endl;
+// branchCPVimpacts ->GetEntry(itrack,0);
+// cout<<" branchCPVimpacts ->GetEntry(itrack,0) OK."<<endl;
+
+// for (Int_t iModule=0; iModule < nOfModules; iModule++) {
+// impacts = (TClonesArray *)fCpvImpacts->At(iModule);
+// cout<<" fCpvImpacts->At(iModule) OK."<<endl;
+// // Do loop over impacts in the module
+// for (Int_t iImpact=0; iImpact<impacts->GetEntries(); iImpact++) {
+// impact=(AliPHOSImpact*)impacts->At(iImpact);
+// impact->Print();
+// if(IsCharged(impact->GetPid()))
+// {
+// cout<<" Add charged hit..";
+// new(hits[hits.GetEntriesFast()]) AliPHOSImpact(*impact);
+// cout<<"done."<<endl;
+// }
+// }
+// }
+// fCpvImpacts->Clear();
+// }
+
+// cout<<" PHOS event "<<ievent<<": "<<hits.GetEntries()<<" charged CPV hits."<<endl;
+
+}
+
+
+// void AliPHOSIhepAnalyze::ChargedHitsCPV(TClonesArray* hits, Int_t ievent, Int_t iModule)
+// {
+// // Cumulative list of charged CPV hits in event ievent
+// // in PHOS module iModule.
+
+// HitsCPV(hits,ievent,iModule);
+// TIter next(hits);
+
+// while(AliPHOSCPVHit* cpvhit = (AliPHOSCPVHit*)next())
+// {
+// if(!IsCharged(cpvhit->GetIpart()))
+// {
+// hits->Remove(cpvhit);
+// delete cpvhit;
+// hits->Compress();
+// }
+// }
+
+// cout<<" PHOS module "<<iModule<<": "<<hits->GetEntries()<<" charged CPV hits."<<endl;
+// }
+
+Bool_t AliPHOSIhepAnalyze::IsCharged(Int_t pdg_code)
+{
+ // For HIJING
+ cout<<" pdg_code "<<pdg_code<<endl;
+ if(pdg_code==211 || pdg_code==-211 || pdg_code==321 || pdg_code==-321 || pdg_code==11 || pdg_code==-11 || pdg_code==2212 || pdg_code==-2212) return kTRUE;
+ else
+ return kFALSE;
+}
+//---------------------------------------------------------------------------
+
+
+
+
+
+
--- /dev/null
+#ifndef AliPHOSIhepAnalyze_H
+#define AliPHOSIhepAnalyze_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+
+//_________________________________________________________________________
+// Algorythm class to analyze PHOSv1 events:
+// Construct histograms and displays them.
+// Used the IHEP CPV/PHOS reconstruction algorithm.
+//*--
+//*-- Author : Boris Polichtchouk (IHEP)
+
+// --- ROOT system ---
+#include "TObject.h"
+
+// --- Standard library ---
+
+// --- AliRoot header files ---
+
+class AliPHOSIhepAnalyze : public TObject {
+
+ public:
+
+ AliPHOSIhepAnalyze() ; // ctor
+ AliPHOSIhepAnalyze(Text_t * name) ; // ctor
+
+ void AnalyzeCPV1(Int_t Nevents); // resolutions, mult and cluster lengths for CPV
+ void AnalyzeEMC1(Int_t Nevents); // resolutions, mult and cluster lengths for EMC
+ void AnalyzeCPV2(Int_t Nevents); // delta(gen)/delta(rec) between hits
+ void CpvSingle(Int_t Nevents); // signle particle analysis
+ virtual void HitsCPV(TClonesArray& hits, Int_t event);
+ TString GetFileName() { return fFileName; }
+
+ private:
+
+ Bool_t IsCharged(Int_t pdg_code);
+
+ private:
+
+ TString fFileName;
+
+ClassDef(AliPHOSIhepAnalyze,1) // PHOSv1 event analyzis algorithm
+
+};
+
+#endif // AliPHOSIhepAnalyze_H
+
+
--- /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. *
+ **************************************************************************/
+
+
+//_________________________________________________________________________
+// Class for the management by the CPV reconstruction.
+//
+//*-- Author : Boris Polichtchouk (IHEP, Protvino) 6 Mar 2001
+//
+// --- ROOT system ---
+
+// --- Standard library ---
+
+#include <iostream.h>
+
+// --- AliRoot header files ---
+
+#include "AliPHOSRecCpvManager.h"
+#include "AliPHOS.h"
+#include "AliRun.h"
+#include "AliPHOSGetter.h"
+
+ClassImp(AliPHOSRecCpvManager)
+
+//____________________________________________________________________________
+
+ AliPHOSRecCpvManager::AliPHOSRecCpvManager()
+{
+
+ fOneGamChisqCut = 3.; // If Chi2/dof > fOneGamChisqCut, split point.
+
+ fOneGamInitialStep = 0.00005;
+ fOneGamChisqMin = 1.;
+ fOneGamStepMin = 0.0005;
+ fOneGamNumOfIterations = 50;
+
+ fTwoGamInitialStep = 0.00005;
+ fTwoGamChisqMin = 1.;
+ fTwoGamEmin = 0.1;
+ fTwoGamStepMin = 0.00005;
+ fTwoGamNumOfIterations = 50;
+
+// fThr0 = 0.0285; // Min. energy of rec. point. If E<fThr0, point deleted.
+// fSqdCut = 3.; // Min. distance (in cm) between two rec points.
+
+ fThr0 = 0.;
+ fSqdCut = 0.;
+
+ SetTitle("Cpv Reconstruction Manager");
+}
+
+AliPHOSRecCpvManager::~AliPHOSRecCpvManager(void) {}
+
+Float_t AliPHOSRecCpvManager::Dispersion(Float_t Etot, Float_t Ai, Float_t Ei)
+{
+ //"Dispresion" of energy deposition in the cell.
+ // Etot is the total shower energy, Ai is the
+ // calculated cell response,
+ // Ei is the measured cell response.
+
+ const Float_t Const = 1.5;
+ return Const*Ai*(1.-Ai/Etot);
+}
+
+Float_t AliPHOSRecCpvManager::OneGamChi2(Float_t Ai, Float_t Ei, Float_t Etot, Float_t& Gi)
+{
+ //"Chi2" for one cell.
+ // Etot is the total "shower" energy, Ai is the
+ // calculated cell response,
+ // Ei is the measured cell response.
+
+ const Float_t Const = 1.5;
+
+ Float_t da = Ai - Ei;
+ Float_t D = Const*Ai*(1.-Ai/Etot);
+
+ Float_t dd = da/D;
+ Gi = dd*(2.- dd*Const*(1.-2.*Ai/Etot));
+
+ cout<<" OneGamChi2 (Ai,Ei,Etot,&Gi,chi2) "<<Ai<<" "<<Ei<<" "<<Etot<<" "<<Gi<<" "<<da*da/D<<endl<<endl;
+
+ return da*da/D;
+
+}
+
+Float_t AliPHOSRecCpvManager::TwoGamChi2(Float_t Ai, Float_t Ei, Float_t Etot, Float_t& Gi)
+{
+
+ const Float_t Const = 1.5;
+
+ Float_t da = Ai - Ei;
+ Float_t D = Const*Ai*(1.-Ai/Etot);
+
+ Float_t dd = da/D;
+ Gi = dd*(2.- dd*Const*(1.-2.*Ai/Etot));
+
+ return da*da/D;
+
+}
+
+void AliPHOSRecCpvManager::AG(Float_t Ei, Float_t Xi, Float_t Yi, Float_t& Ai, Float_t& GXi, Float_t& GYi )
+{
+ //Calculates amplitude (Ai) and gradients (GXi, GYi) of CPV pad response.
+ //Integrated response (total "shower energy") is E,
+ //Xi and Yi are the distances along x and y from reference point
+ // to the pad center.
+
+ AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
+ const AliPHOSGeometry* geom = gime->PHOSGeometry();
+ Float_t CelZ = geom->GetPadSizeZ();
+ Float_t CelY = geom->GetPadSizePhi();
+
+// // cout<<"CelZ: "<<CelZ<<" CelY: "<<CelY<<endl;
+
+ Float_t dx = CelZ/2.;
+ Float_t dy = CelY/2.;
+
+// // Float_t x = Xi*CelZ;
+// // Float_t y = Yi*CelZ;
+
+ Float_t x = Xi*CelZ;
+ Float_t y = Yi*CelY;
+
+ Float_t E = Ei;
+
+ Float_t A = Fcml(x+dx,y+dy) - Fcml(x+dx,y-dy) - Fcml(x-dx,y+dy) + Fcml(x-dx,y-dy);
+ Ai = A*E;
+
+
+ Float_t Gx = GradX(x+dx,y+dy) - GradX(x+dx,y-dy) - GradX(x-dx,y+dy) + GradX(x-dx,y-dy);
+ GXi = Gx*E*E;
+
+ Float_t Gy = GradY(x+dx,y+dy) - GradY(x+dx,y-dy) - GradY(x-dx,y+dy) + GradY(x-dx,y-dy);
+ GYi = Gy*E*E;
+
+}
+
+Float_t AliPHOSRecCpvManager::Fcml(Float_t x, Float_t y)
+{
+ //Cumulative function
+
+ const Float_t A = 1.0;
+ const Float_t b = 0.70;
+
+ Float_t Fff = TMath::ATan(x*y/( b*TMath::Sqrt( (b*b) + x*x+y*y)))
+ - TMath::ATan(x*y/(3*b*TMath::Sqrt((3*b)*(3*b) + x*x+y*y)))
+ + TMath::ATan(x*y/(5*b*TMath::Sqrt((5*b)*(5*b) + x*x+y*y)))
+ - TMath::ATan(x*y/(7*b*TMath::Sqrt((7*b)*(7*b) + x*x+y*y)))
+ + TMath::ATan(x*y/(9*b*TMath::Sqrt((9*b)*(9*b) + x*x+y*y)));
+
+ Float_t Fcml = A*Fff/6.2831853071796;
+// cout<<" Fcml: "<<Fcml<<endl;
+ return Fcml;
+
+}
+
+
+Float_t AliPHOSRecCpvManager::GradX(Float_t x, Float_t y)
+{
+
+ const Float_t A = 1.0;
+ const Float_t b = 0.70;
+
+ Float_t skv = b*b + x*x + y*y;
+
+ Float_t Gradient = y*(1.-x*x/skv)* b*TMath::Sqrt(skv)/( b*b*skv+x*x*y*y)
+ - y*(1.-x*x/skv)*3*b*TMath::Sqrt(skv)/((3*b)*(3*b)*skv+x*x*y*y)
+ + y*(1.-x*x/skv)*5*b*TMath::Sqrt(skv)/((5*b)*(5*b)*skv+x*x*y*y)
+ - y*(1.-x*x/skv)*7*b*TMath::Sqrt(skv)/((7*b)*(7*b)*skv+x*x*y*y)
+ + y*(1.-x*x/skv)*9*b*TMath::Sqrt(skv)/((9*b)*(9*b)*skv+x*x*y*y);
+
+ Float_t Grad = A*Gradient/6.2831853071796;
+ return Grad;
+}
+
+
+Float_t AliPHOSRecCpvManager::GradY(Float_t x, Float_t y)
+{
+
+ const Float_t A = 1.0;
+ const Float_t b = 0.70;
+
+ Float_t skv = b*b + x*x + y*y;
+ Float_t Gradient = x*(1.-y*y/skv)* b*TMath::Sqrt(skv)/( b*b*skv+x*x*y*y)
+ - x*(1.-y*y/skv)*3*b*TMath::Sqrt(skv)/((3*b)*(3*b)*skv+x*x*y*y)
+ + x*(1.-y*y/skv)*5*b*TMath::Sqrt(skv)/((5*b)*(5*b)*skv+x*x*y*y)
+ - x*(1.-y*y/skv)*7*b*TMath::Sqrt(skv)/((7*b)*(7*b)*skv+x*x*y*y)
+ + x*(1.-y*y/skv)*9*b*TMath::Sqrt(skv)/((9*b)*(9*b)*skv+x*x*y*y);
+
+ Float_t Grad = A*Gradient/6.2831853071796;
+ return Grad;
+}
+
+
--- /dev/null
+#ifndef AliPHOSRecCpvManager_H
+#define AliPHOSRecCpvManager_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+//_________________________________________________________________________
+// Class for the management of the CPV reconstruction.
+// Author : Boris Polichtchouk (IHEP, Protvino)
+// 6 March 2001
+
+#include "AliPHOSRecManager.h"
+#include "AliPHOSGeometry.h"
+
+class AliPHOSRecCpvManager : public AliPHOSRecManager {
+
+ public:
+
+ AliPHOSRecCpvManager();
+ ~AliPHOSRecCpvManager(void);
+
+
+ void AG(Float_t E, Float_t dx, Float_t dy, Float_t& A, Float_t& grad_x, Float_t& grad_y );
+ Float_t Dispersion(Float_t Etot, Float_t Ai, Float_t Ei);
+
+ Float_t OneGamChi2(Float_t Ai, Float_t Ei, Float_t Etot, Float_t& Gi);
+ Float_t TwoGamChi2(Float_t Ai, Float_t Ei, Float_t Etot, Float_t& Gi);
+
+ Float_t OneGamChisqCut() { return fOneGamChisqCut; }
+ Float_t OneGamInitialStep() { return fOneGamInitialStep; }
+ Float_t OneGamChisqMin() { return fOneGamChisqMin; }
+ Float_t OneGamStepMin() { return fOneGamStepMin; }
+ Int_t OneGamNumOfIterations() { return fOneGamNumOfIterations; }
+
+ Float_t TwoGamInitialStep() { return fTwoGamInitialStep; }
+ Float_t TwoGamChisqMin() { return fTwoGamChisqMin; }
+ Float_t TwoGamEmin() { return fTwoGamEmin; }
+ Float_t TwoGamStepMin() { return fTwoGamStepMin; }
+ Int_t TwoGamNumOfIterations() { return fTwoGamNumOfIterations; }
+
+ Float_t KillGamMinEnergy() { return fThr0; }
+ Float_t MergeGammasMinDistanceCut() { return fSqdCut; }
+
+ void SetTwoPointsMinDistance(Float_t dist) { fSqdCut=dist; }
+ void SetPointMinEnergy(Float_t emin) { fThr0=emin; }
+
+ private:
+
+ Float_t Fcml(Float_t x, Float_t y);
+ Float_t GradX(Float_t x, Float_t y);
+ Float_t GradY(Float_t x, Float_t y);
+
+ private:
+
+ Float_t fOneGamChisqCut;
+
+ Float_t fOneGamInitialStep;
+ Float_t fOneGamChisqMin;
+ Float_t fOneGamStepMin;
+ Int_t fOneGamNumOfIterations;
+
+ Float_t fTwoGamInitialStep;
+ Float_t fTwoGamChisqMin;
+ Float_t fTwoGamEmin;
+ Float_t fTwoGamStepMin;
+ Int_t fTwoGamNumOfIterations;
+
+ Float_t fThr0;
+ Float_t fSqdCut;
+
+ public:
+
+ ClassDef(AliPHOSRecCpvManager,1) // CPV reconstruction management class
+
+} ;
+
+#endif // AliPHOSRecCpvManager_H
+
+
+
--- /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. *
+ **************************************************************************/
+
+
+//_________________________________________________________________________
+// Class for the management by the Emc reconstruction.
+//
+//*-- Author : Boris Polichtchouk (IHEP, Protvino) 6 Mar 2001
+//
+// --- ROOT system ---
+
+// --- Standard library ---
+
+#include <iostream.h>
+
+// --- AliRoot header files ---
+
+#include "AliPHOSRecEmcManager.h"
+#include "AliPHOS.h"
+#include "AliRun.h"
+
+ClassImp(AliPHOSRecEmcManager)
+
+//____________________________________________________________________________
+
+ AliPHOSRecEmcManager::AliPHOSRecEmcManager()
+{
+
+// fOneGamChisqCut = 3.;
+ fOneGamChisqCut = 1.3; // bvp 31.08.2001
+
+ fOneGamInitialStep = 0.00005;
+ fOneGamChisqMin = 1.;
+ fOneGamStepMin = 0.0005;
+ fOneGamNumOfIterations = 50;
+
+ fTwoGamInitialStep = 0.00005;
+ fTwoGamChisqMin = 1.;
+ fTwoGamEmin = 0.1;
+ fTwoGamStepMin = 0.00005;
+ fTwoGamNumOfIterations = 50;
+
+// fThr0 = 0.6;
+ fThr0 = 0.;
+// fSqdCut = 3.;
+// fSqdCut = 0.5; // bvp 31.08.2001
+ fSqdCut = 0.;
+
+ SetTitle("Emc Reconstruction Manager");
+
+}
+
+AliPHOSRecEmcManager::~AliPHOSRecEmcManager(void) {}
+
+Float_t AliPHOSRecEmcManager::Dispersion(Float_t Etot, Float_t Ai, Float_t Ei)
+{
+ //"Dispresion" of energy deposition in the cell.
+ // Etot is the total shower energy, Ai is the
+ // calculated cell response,
+ // Ei is the measured cell response.
+
+ return Ei;
+}
+
+Float_t AliPHOSRecEmcManager::OneGamChi2(Float_t Ai, Float_t Ei, Float_t Etot, Float_t& Gi)
+{
+ // Chi2 used in OneGam (one-gamma fitting).
+ // Gi is d(Chi2)/d(Ai).
+
+ Float_t da = Ai - Ei;
+ Float_t D = Ei; // we assume that sigma(E) = sqrt(E)
+ Gi = 2.*(Ai-Ei)/D;
+
+ return da*da/D;
+
+}
+
+Float_t AliPHOSRecEmcManager::TwoGamChi2(Float_t Ai, Float_t Ei, Float_t Etot, Float_t& Gi)
+{
+
+ Float_t da = Ai - Ei;
+ Float_t D = Ei; // we assume that sigma(E) = sqrt(E)
+ Gi = 2.*(Ai-Ei)/D;
+
+ return da*da/D;
+
+}
+
+void AliPHOSRecEmcManager::AG(Float_t Ei, Float_t Xi, Float_t Yi, Float_t& Ai, Float_t& GXi, Float_t& GYi )
+{
+ //Calculates amplitude (Ai) and gradients (GXi, GYi) of CPV pad response.
+ //Integrated response (total "shower energy") is E,
+ //Xi and Yi are the distances along x and y from reference point
+ // to the pad center.
+ //Shape of the shower is from PHOS TDR.
+
+
+ Float_t r = TMath::Sqrt(Xi*Xi + Yi*Yi);
+ Float_t r4 = r*r*r*r ;
+ Float_t r295 = TMath::Power(r, 2.95) ;
+ Float_t shape = Ei*TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
+ Ai = shape;
+
+
+ //d(shape)/d(Xi)
+ GXi = (-(TMath::Power(TMath::Power(Xi,2) + TMath::Power(Yi,2),2)*
+ ((-0.006077944*Xi*TMath::Power(TMath::Power(Xi,2) + TMath::Power(Yi,2),
+ 0.4750000000000001))/
+ TMath::Power(1 + 0.0652*TMath::Power(TMath::Power(Xi,2) + TMath::Power(Yi,2),1.475),2) -
+ (1.04*Xi*(TMath::Power(Xi,2) + TMath::Power(Yi,2)))/
+ TMath::Power(2.32 + 0.26*TMath::Power(TMath::Power(Xi,2) + TMath::Power(Yi,2),2),2))) -
+ 4*Xi*(TMath::Power(Xi,2) + TMath::Power(Yi,2))*
+ (0.0316/(1 + 0.0652*TMath::Power(TMath::Power(Xi,2) + TMath::Power(Yi,2),1.475)) +
+ 1./(2.32 + 0.26*TMath::Power(TMath::Power(Xi,2) + TMath::Power(Yi,2),2))))/
+ TMath::Power(TMath::E(),TMath::Power(TMath::Power(Xi,2) + TMath::Power(Yi,2),2)*
+ (0.0316/(1 + 0.0652*TMath::Power(TMath::Power(Xi,2) + TMath::Power(Yi,2),1.475)) +
+ 1./(2.32 + 0.26*TMath::Power(TMath::Power(Xi,2) + TMath::Power(Yi,2),2))));
+
+ GXi = GXi*Ei;
+
+ //d(shape)/d(Yi)
+ GYi = (-(TMath::Power(TMath::Power(Xi,2) + TMath::Power(Yi,2),2)*
+ ((-0.006077944*Yi*TMath::Power(TMath::Power(Xi,2) + TMath::Power(Yi,2),
+ 0.4750000000000001))/
+ TMath::Power(1 + 0.0652*TMath::Power(TMath::Power(Xi,2) + TMath::Power(Yi,2),1.475),2) -
+ (1.04*Yi*(TMath::Power(Xi,2) + TMath::Power(Yi,2)))/
+ TMath::Power(2.32 + 0.26*TMath::Power(TMath::Power(Xi,2) + TMath::Power(Yi,2),2),2))) -
+ 4*Yi*(TMath::Power(Xi,2) + TMath::Power(Yi,2))*
+ (0.0316/(1 + 0.0652*TMath::Power(TMath::Power(Xi,2) + TMath::Power(Yi,2),1.475)) +
+ 1./(2.32 + 0.26*TMath::Power(TMath::Power(Xi,2) + TMath::Power(Yi,2),2))))/
+ TMath::Power(TMath::E(),TMath::Power(TMath::Power(Xi,2) + TMath::Power(Yi,2),2)*
+ (0.0316/(1 + 0.0652*TMath::Power(TMath::Power(Xi,2) + TMath::Power(Yi,2),1.475)) +
+ 1./(2.32 + 0.26*TMath::Power(TMath::Power(Xi,2) + TMath::Power(Yi,2),2))));
+
+
+ GYi = GYi*Ei;
+
+}
+
+
+
+
+
+
+
--- /dev/null
+#ifndef AliPHOSRecEmcManager_H
+#define AliPHOSRecEmcManager_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+//_________________________________________________________________________
+// Class for the management by the Emc reconstruction.
+// Author : Boris Polichtchouk (IHEP, Protvino)
+// 6 March 2001
+
+#include "AliPHOSRecManager.h"
+#include "AliPHOSGeometry.h"
+
+class AliPHOSRecEmcManager : public AliPHOSRecManager {
+
+ public:
+
+ AliPHOSRecEmcManager();
+ ~AliPHOSRecEmcManager(void);
+
+
+ void AG(Float_t E, Float_t dx, Float_t dy, Float_t& A, Float_t& grad_x, Float_t& grad_y );
+ Float_t Dispersion(Float_t Etot, Float_t Ai, Float_t Ei);
+
+ Float_t OneGamChi2(Float_t Ai, Float_t Ei, Float_t Etot, Float_t& Gi);
+ Float_t TwoGamChi2(Float_t Ai, Float_t Ei, Float_t Etot, Float_t& Gi);
+
+ Float_t OneGamChisqCut() { return fOneGamChisqCut; }
+ Float_t OneGamInitialStep() { return fOneGamInitialStep; }
+ Float_t OneGamChisqMin() { return fOneGamChisqMin; }
+ Float_t OneGamStepMin() { return fOneGamStepMin; }
+ Int_t OneGamNumOfIterations() { return fOneGamNumOfIterations; }
+
+ Float_t TwoGamInitialStep() { return fTwoGamInitialStep; }
+ Float_t TwoGamChisqMin() { return fTwoGamChisqMin; }
+ Float_t TwoGamEmin() { return fTwoGamEmin; }
+ Float_t TwoGamStepMin() { return fTwoGamStepMin; }
+ Int_t TwoGamNumOfIterations() { return fTwoGamNumOfIterations; }
+
+ Float_t KillGamMinEnergy() { return fThr0; }
+ Float_t MergeGammasMinDistanceCut() { return fSqdCut; }
+
+ void SetTwoPointsMinDistance(Float_t dist) { fSqdCut=dist; }
+ void SetPointMinEnergy(Float_t emin) { fThr0=emin; }
+
+ private:
+
+ Float_t fOneGamChisqCut;
+
+ Float_t fOneGamInitialStep;
+ Float_t fOneGamChisqMin;
+ Float_t fOneGamStepMin;
+ Int_t fOneGamNumOfIterations;
+
+ Float_t fTwoGamInitialStep;
+ Float_t fTwoGamChisqMin;
+ Float_t fTwoGamEmin;
+ Float_t fTwoGamStepMin;
+ Int_t fTwoGamNumOfIterations;
+
+ Float_t fThr0;
+ Float_t fSqdCut;
+
+ public:
+
+ ClassDef(AliPHOSRecEmcManager,1) // Emc reconstruction management class
+
+} ;
+
+#endif // AliPHOSRecEmcManager_H
+
+
+
--- /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. *
+ **************************************************************************/
+
+
+//_________________________________________________________________________
+// Base class for the management by the PHOS reconstruction.
+//
+//*-- Author : Boris Polichtchouk (IHEP, Protvino) 6 Mar 2001
+//
+// --- ROOT system ---
+
+// --- Standard library ---
+
+#include <iostream.h>
+
+// --- AliRoot header files ---
+
+#include "AliPHOSRecManager.h"
+
+ClassImp(AliPHOSRecManager)
+
+ AliPHOSRecManager::AliPHOSRecManager() : TNamed("AliPHOSRecManager","PHOS Reconstruction Manager") {}
+//____________________________________________________________________________
+
+
--- /dev/null
+#ifndef AliPHOSRecManager_H
+#define AliPHOSRecManager_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+//_________________________________________________________________________
+// Base class for the management by the PHOS reconstruction.
+// It contains only virtual member functions
+// which will be implemented for the Emc and CPV reconstruction
+// in the appropriate derived classes.
+// Author : Boris Polichtchouk (IHEP, Protvino)
+// 6 March 2001
+
+#include "TNamed.h"
+
+class AliPHOSRecManager : public TNamed {
+
+ public:
+
+ AliPHOSRecManager();
+ virtual ~AliPHOSRecManager(void) {}
+
+ virtual void AG(Float_t E, Float_t dx, Float_t dy, Float_t& A, Float_t& grad_x, Float_t& grad_y ) = 0;
+ virtual Float_t Dispersion(Float_t Etot, Float_t Ai, Float_t Ei) = 0;
+
+ virtual Float_t OneGamChi2(Float_t Ai, Float_t Ei, Float_t Etot, Float_t& Gi) = 0;
+ virtual Float_t TwoGamChi2(Float_t Ai, Float_t Ei, Float_t Etot, Float_t& Gi) = 0;
+
+ virtual Float_t OneGamChisqCut() = 0 ;
+ virtual Float_t OneGamInitialStep() = 0;
+ virtual Float_t OneGamChisqMin() = 0;
+ virtual Float_t OneGamStepMin() = 0;
+ virtual Int_t OneGamNumOfIterations() = 0;
+
+ virtual Float_t TwoGamInitialStep() = 0;
+ virtual Float_t TwoGamChisqMin() = 0;
+ virtual Float_t TwoGamEmin() = 0;
+ virtual Float_t TwoGamStepMin() = 0;
+ virtual Int_t TwoGamNumOfIterations() = 0;
+
+ virtual Float_t KillGamMinEnergy() = 0;
+ virtual Float_t MergeGammasMinDistanceCut() = 0;
+
+ virtual void SetTwoPointsMinDistance(Float_t dist) = 0;
+ virtual void SetPointMinEnergy(Float_t emin) = 0;
+
+ ClassDef(AliPHOSRecManager,1)
+
+} ;
+
+#endif // AliPHOSRecManager_H
+
+
+
+
+
+
+
+
+
+
+
+
+
AliPHOSGetter.cxx AliPHOSTick.cxx \
AliPHOSQAVirtualCheckable.cxx AliPHOSQAIntCheckable.cxx \
AliPHOSQAFloatCheckable.cxx\
- AliPHOSQAObjectCheckable.cxx AliPHOSQAChecker.cxx AliPHOSQAMeanChecker.cxx AliPHOSQAAlarm.cxx
+ AliPHOSQAObjectCheckable.cxx AliPHOSQAChecker.cxx AliPHOSQAMeanChecker.cxx AliPHOSQAAlarm.cxx \
+ AliPHOSIhepAnalyze.cxx AliPHOSEvalRecPoint.cxx \
+ AliPHOSRecManager.cxx AliPHOSRecCpvManager.cxx AliPHOSRecEmcManager.cxx \
+ AliPHOSClusterizerv2.cxx
# C++ Headers
-include tgt_$(ALICE_TARGET)/Make-depend
-
test:
@echo " ____________________________________________________________ "
@echo " "
@echo " Starting the test of the simulation/reconstruction software. Please don't take the warning messages into account. "
@echo " ____________________________________________________________ "
- @aliroot -b -q "testsim.C(100)" > out
+ @aliroot -b -q testsim.C > out
@aliroot -b -q testsimglobal.C > out
@rm out
@rm testPHOS.root
-test10:
- @echo " ____________________________________________________________ "
- @echo " "
- @echo " Starting the test of the simulation/reconstruction software. Please don't take the warning messages into account. "
- @echo " ____________________________________________________________ "
- @aliroot -b -q "testsim.C(10)" > out
- @aliroot -b -q testsimglobal.C > out
- @rm out
- @rm testPHOS.root
+
+
#pragma link C++ class AliPHOSQAMeanChecker+;
#pragma link C++ class AliPHOSQAAlarm+;
#pragma link C++ class AliPHOSTick+;
+#pragma link C++ class AliPHOSIhepAnalyze+;
+#pragma link C++ class AliPHOSRecManager+;
+#pragma link C++ class AliPHOSRecCpvManager+;
+#pragma link C++ class AliPHOSRecEmcManager+;
+#pragma link C++ class AliPHOSClusterizerv2+;
+#pragma link C++ class AliPHOSEvalRecPoint+;
#endif