X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PHOS%2FAliPHOSTrackSegmentMakerv1.cxx;h=ba31b3cc1f3c683d97bf5fdba2511598bf578f93;hb=25312a8ed5d33f78b0423d6e0c7ee3869142528f;hp=d388bd58ffbe1ec40bc5ecba31168b380fa43d01;hpb=88714635380b3cd769507e2f5b0583b7214b1e96;p=u%2Fmrichter%2FAliRoot.git diff --git a/PHOS/AliPHOSTrackSegmentMakerv1.cxx b/PHOS/AliPHOSTrackSegmentMakerv1.cxx index d388bd58ffb..ba31b3cc1f3 100644 --- a/PHOS/AliPHOSTrackSegmentMakerv1.cxx +++ b/PHOS/AliPHOSTrackSegmentMakerv1.cxx @@ -12,606 +12,681 @@ * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ - /* $Id$ */ +/* History of cvs commits: + * + * $Log$ + * Revision 1.93 2007/10/10 09:05:10 schutz + * Changing name QualAss to QA + * + * Revision 1.92 2007/08/28 12:55:08 policheh + * Loaders removed from the reconstruction code (C.Cheshkov) + * + * Revision 1.91 2007/08/07 14:12:03 kharlov + * Quality assurance added (Yves Schutz) + * + * Revision 1.90 2007/07/11 13:43:30 hristov + * New class AliESDEvent, backward compatibility with the old AliESD (Christian) + * + * Revision 1.89 2007/07/03 08:13:04 kharlov + * Bug fix in CPV local coordinates + * + * Revision 1.88 2007/06/27 09:11:07 kharlov + * Bug fix for CPV-EMC distance + * + * Revision 1.87 2007/05/04 14:49:29 policheh + * AliPHOSRecPoint inheritance from AliCluster + * + * Revision 1.86 2007/04/02 15:00:16 cvetan + * No more calls to gAlice in the reconstruction + * + * Revision 1.85 2007/03/28 19:18:15 kharlov + * RecPoints recalculation in TSM removed + * + * Revision 1.84 2007/03/07 07:01:21 hristov + * Fixing copy/paste erro. Additional protections + * + * Revision 1.83 2007/03/06 21:07:37 kharlov + * DP: xz CPV-EMC distance filled to TS + * + * Revision 1.82 2007/03/06 06:54:48 kharlov + * DP:Calculation of cluster properties dep. on vertex added + * + * Revision 1.81 2007/02/05 10:02:40 kharlov + * Module numbering is corrected + * + * Revision 1.80 2006/08/28 10:01:56 kharlov + * Effective C++ warnings fixed (Timur Pocheptsov) + * + * Revision 1.79 2006/04/25 12:41:15 hristov + * Moving non-persistent data to AliESDfriend (Yu.Belikov) + * + * Revision 1.78 2005/11/18 13:04:51 hristov + * Bug fix + * + * Revision 1.77 2005/11/17 23:34:36 hristov + * Corrected logics + * + * Revision 1.76 2005/11/17 22:29:12 hristov + * Faster version, no attempt to match tracks outside the PHOS acceptance + * + * Revision 1.75 2005/11/17 12:35:27 hristov + * Use references instead of objects. Avoid to create objects when they are not really needed + * + * Revision 1.74 2005/07/08 14:01:36 hristov + * Tracking in non-uniform nmagnetic field (Yu.Belikov) + * + * Revision 1.73 2005/05/28 14:19:05 schutz + * Compilation warnings fixed by T.P. + * + */ + //_________________________________________________________________________ // Implementation version 1 of algorithm class to construct PHOS track segments -// Associates EMC and PPSD clusters -// Unfolds the EMC cluster -// -//*-- Author: Dmitri Peressounko (RRC Ki & SUBATECH) +// Track segment for PHOS is list of +// EMC RecPoint + (possibly) CPV RecPoint +// To find TrackSegments we do the following: +// for each EMC RecPoints we look at +// CPV RecPoints in the radious fRcpv. +// If there is such a CPV RecPoint, +// we make "Link" it is just indexes of EMC and CPV RecPoint and distance +// between them in the PHOS plane. +// Then we sort "Links" and starting from the +// least "Link" pointing to the unassined EMC and CPV RecPoints assing them to +// new TrackSegment. +// If there is no CPV RecPoint we make TrackSegment +// consisting from EMC alone. There is no TrackSegments without EMC RecPoint. +//// In principle this class should be called from AliPHOSReconstructor, but +// one can use it as well in standalone mode. +// Use case: +// root [0] AliPHOSTrackSegmentMakerv1 * t = new AliPHOSTrackSegmentMaker("galice.root", "tracksegmentsname", "recpointsname") +// Warning in : object already instantiated +// // reads gAlice from header file "galice.root", uses recpoints stored in the branch names "recpointsname" (default = "Default") +// // and saves recpoints in branch named "tracksegmentsname" (default = "recpointsname") +// root [1] t->ExecuteTask() +// root [3] t->SetTrackSegmentsBranch("max distance 5 cm") +// root [4] t->ExecuteTask("deb all time") +// +//*-- Author: Dmitri Peressounko (RRC Ki & SUBATECH) & Yves Schutz (SUBATECH) // // --- ROOT system --- - -#include "TObjArray.h" -#include "TClonesArray.h" -#include "TObjectTable.h" +#include "TVector3.h" +#include "TTree.h" +#include "TBenchmark.h" // --- Standard library --- - -#include - +#include "Riostream.h" // --- AliRoot header files --- - +#include "AliPHOSGeometry.h" #include "AliPHOSTrackSegmentMakerv1.h" #include "AliPHOSTrackSegment.h" #include "AliPHOSLink.h" -#include "AliPHOSv0.h" -#include "AliRun.h" +#include "AliESDEvent.h" +#include "AliESDtrack.h" +#include "AliPHOSEmcRecPoint.h" +#include "AliPHOSCpvRecPoint.h" +#include "AliLog.h" +#include "AliMagF.h" +#include "AliMagF.h" +#include "AliTracker.h" +#include "AliGeomManager.h" +#include "AliCluster.h" +#include "AliKalmanTrack.h" +#include "AliGlobalQADataMaker.h" +#include "AliVParticle.h" -extern void UnfoldingChiSquare(Int_t &nPar, Double_t *Grad, Double_t & fret, Double_t *x, Int_t iflag) ; ClassImp( AliPHOSTrackSegmentMakerv1) //____________________________________________________________________________ - AliPHOSTrackSegmentMakerv1::AliPHOSTrackSegmentMakerv1() : AliPHOSTrackSegmentMaker() +AliPHOSTrackSegmentMakerv1::AliPHOSTrackSegmentMakerv1() : + AliPHOSTrackSegmentMaker(), + fDefaultInit(kTRUE), + fWrite(kFALSE), + fNTrackSegments(0), + fRcpv(0.f), + fRtpc(0.f), + fVtx(0.f), + fLinkUpArray(0), + fEmcFirst(0), + fEmcLast(0), + fCpvFirst(0), + fCpvLast(0), + fModule(0), + fTrackSegments(NULL) +{ + // default ctor (to be used mainly by Streamer) + InitParameters() ; +} + +//____________________________________________________________________________ +AliPHOSTrackSegmentMakerv1::AliPHOSTrackSegmentMakerv1(AliPHOSGeometry *geom) : + AliPHOSTrackSegmentMaker(geom), + fDefaultInit(kFALSE), + fWrite(kFALSE), + fNTrackSegments(0), + fRcpv(0.f), + fRtpc(0.f), + fVtx(0.f), + fLinkUpArray(0), + fEmcFirst(0), + fEmcLast(0), + fCpvFirst(0), + fCpvLast(0), + fModule(0), + fTrackSegments(NULL) { // ctor + InitParameters() ; + Init() ; + fESD = 0; +} + - fR0 = 4. ; - AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ; - //clusters are sorted in "rows" and "columns" of width geom->GetCrystalSize(0), - fDelta = fR0 + geom->GetCrystalSize(0) ; - fMinuit = new TMinuit(100) ; - fUnfoldFlag = kTRUE ; +AliPHOSTrackSegmentMakerv1::AliPHOSTrackSegmentMakerv1(const AliPHOSTrackSegmentMakerv1 & tsm) : + AliPHOSTrackSegmentMaker(tsm), + fDefaultInit(kFALSE), + fWrite(kFALSE), + fNTrackSegments(0), + fRcpv(0.f), + fRtpc(0.f), + fVtx(0.f), + fLinkUpArray(0), + fEmcFirst(0), + fEmcLast(0), + fCpvFirst(0), + fCpvLast(0), + fModule(0), + fTrackSegments(NULL) +{ + // cpy ctor: no implementation yet + // requested by the Coding Convention + Fatal("cpy ctor", "not implemented") ; } + //____________________________________________________________________________ AliPHOSTrackSegmentMakerv1::~AliPHOSTrackSegmentMakerv1() { // dtor - - delete fMinuit ; + // fDefaultInit = kTRUE if TrackSegmentMaker created by default ctor (to get just the parameters) + if (!fDefaultInit) + delete fLinkUpArray ; + if (fTrackSegments) { + fTrackSegments->Delete(); + delete fTrackSegments; + } } //____________________________________________________________________________ -Bool_t AliPHOSTrackSegmentMakerv1::FindFit(AliPHOSEmcRecPoint * emcRP, int * maxAt, Float_t * maxAtEnergy, - Int_t nPar, Float_t * fitparameters) -{ - // Calls TMinuit to fit the energy distribution of a cluster with several maxima +void AliPHOSTrackSegmentMakerv1::FillOneModule() +{ + // Finds first and last indexes between which + // clusters from one PHOS module are + + //First EMC clusters + Int_t totalEmc = fEMCRecPoints->GetEntriesFast() ; + for(fEmcFirst = fEmcLast; (fEmcLast < totalEmc) && + ((dynamic_cast(fEMCRecPoints->At(fEmcLast)))->GetPHOSMod() == fModule ); + fEmcLast ++) ; + + //Now CPV clusters + Int_t totalCpv = fCPVRecPoints->GetEntriesFast() ; - AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ; + for(fCpvFirst = fCpvLast; (fCpvLast < totalCpv) && + ((dynamic_cast(fCPVRecPoints->At(fCpvLast)))->GetPHOSMod() == fModule ); + fCpvLast ++) ; + +} - gMinuit->SetPrintLevel(-1) ; // No Printout - gMinuit->SetFCN(UnfoldingChiSquare) ; // To set the address of the minimization function - gMinuit->SetObjectFit(emcRP) ; // To tranfer pointer to UnfoldingChiSquare +//____________________________________________________________________________ +void AliPHOSTrackSegmentMakerv1::GetDistanceInPHOSPlane(AliPHOSEmcRecPoint * emcClu, + AliPHOSCpvRecPoint * cpvClu, + Int_t &trackindex, + Float_t &dx, Float_t &dz) const +{ + // Calculates the distance between the EMC RecPoint and the CPV RecPoint + // If no CPV, calculates the distance between the EMC RecPoint and the track + // prolongation to the PHOS module plane. + // Clusters are sorted in "rows" and "columns" of width 1 cm - // filling initial values for fit parameters - AliPHOSDigit * digit ; +// Float_t delta = 1 ; // Width of the rows in sorting of RecPoints (in cm) +// // if you change this value, change it as well in xxxRecPoint::Compare() - Int_t ierflg = 0; - Int_t index = 0 ; - Int_t nDigits = (Int_t) nPar / 3 ; + trackindex = -1; + dx = 999.; + dz = 999.; - Int_t iDigit ; + if(!cpvClu) { + + if(!emcClu) { + return; + } + // *** Start the matching + Int_t nt=fESD->GetNumberOfTracks(); + Int_t iPHOSMod = emcClu->GetPHOSMod() ; + //Calculate actual distance to PHOS module + TVector3 globaPos ; + fGeom->Local2Global(iPHOSMod, 0.,0., globaPos) ; + const Double_t rPHOS = globaPos.Pt() ; //Distance to center of PHOS module + const Double_t kYmax = 72.+10. ; //Size of the module (with some reserve) in phi direction + const Double_t kZmax = 64.+10. ; //Size of the module (with some reserve) in z direction + const Double_t kAlpha0=330./180.*TMath::Pi() ; //First PHOS module angular direction + const Double_t kAlpha= 20./180.*TMath::Pi() ; //PHOS module angular size + Double_t minDistance = 1.e6; + + TVector3 vecEmc ; // Local position of EMC recpoint + emcClu->GetLocalPosition(vecEmc) ; + + Double_t gposTrack[3] ; + Double_t bz = AliTracker::GetBz() ; //B-Field for approximate matching + Double_t b[3]; + for (Int_t i=0; iGetTrack(i); + + // Skip the tracks having "wrong" status (has to be checked/tuned) + ULong_t status = esdTrack->GetStatus(); + if ((status & AliESDtrack::kTPCout) == 0) continue; +// if ((status & AliESDtrack::kTRDout) == 0) continue; +// if ((status & AliESDtrack::kTRDrefit) == 1) continue; + + //Continue extrapolation from TPC outer surface + const AliExternalTrackParam *outerParam=esdTrack->GetOuterParam(); + if (!outerParam) continue; + AliExternalTrackParam t(*outerParam); + + t.GetBxByBz(b) ; + //Direction to the current PHOS module + Double_t phiMod=kAlpha0-kAlpha*iPHOSMod ; + if(!t.Rotate(phiMod)) + continue ; + + Double_t y; // Some tracks do not reach the PHOS + if (!t.GetYAt(rPHOS,bz,y)) continue; // because of the bending + + Double_t z; + if(!t.GetZAt(rPHOS,bz,z)) + continue ; + if (TMath::Abs(z) > kZmax) + continue; // Some tracks miss the PHOS in Z + if(TMath::Abs(y) < kYmax){ + t.PropagateToBxByBz(rPHOS,b); // Propagate to the matching module + //t.CorrectForMaterial(...); // Correct for the TOF material, if needed + t.GetXYZ(gposTrack) ; + TVector3 globalPositionTr(gposTrack) ; + TVector3 localPositionTr ; + fGeom->Global2Local(localPositionTr,globalPositionTr,iPHOSMod) ; + Double_t ddx = vecEmc.X()-localPositionTr.X(); + Double_t ddz = vecEmc.Z()-localPositionTr.Z(); + Double_t d2 = ddx*ddx + ddz*ddz; + if(d2 < minDistance) { + dx = ddx ; + dz = ddz ; + trackindex=i; + minDistance=d2 ; + } + } + } //Scanned all tracks + return ; + } - for(iDigit = 0; iDigit < nDigits; iDigit++){ - digit = (AliPHOSDigit *) maxAt[iDigit]; - Int_t relid[4] ; - Float_t x ; - Float_t z ; - geom->AbsToRelNumbering(digit->GetId(), relid) ; - geom->RelPosInModule(relid, x, z) ; + TVector3 emcGlobal; + fGeom->GetGlobalPHOS((AliPHOSRecPoint*)emcClu,emcGlobal); + + // Radius from IP to current point + Double_t rEMC = TMath::Abs(emcGlobal.Pt()); + + // Extrapolate the global track direction to EMC + // and find the closest track + + Int_t nTracks = fESD->GetNumberOfTracks(); + + AliESDtrack *track; + Double_t xyz[] = {-1,-1,-1}; + Double_t pxyz[3]; + Double_t zEMC,xEMC; + Int_t module; + TVector3 vecP; + TVector3 locClu; + + Float_t minDistance = 1.e6; + Float_t dr; + + for (Int_t iTrack=0; iTrackGetTrack(iTrack); + if (!track->GetXYZAt(rEMC, fESD->GetMagneticField(), xyz)) continue; + + AliDebug(1,Form("Event %d, iTrack: %d, (%.3f,%.3f,%.3f)", + fESD->GetEventNumberInFile(),iTrack,xyz[0],xyz[1],xyz[2])); + + if (track->GetPxPyPzAt(rEMC,fESD->GetMagneticField(),pxyz)) { - Float_t energy = maxAtEnergy[iDigit] ; + vecP.SetXYZ(pxyz[0],pxyz[1],pxyz[2]); + fGeom->ImpactOnEmc(xyz,vecP.Theta(),vecP.Phi(),module,zEMC,xEMC) ; - gMinuit->mnparm(index, "x", x, 0.1, 0, 0, ierflg) ; - index++ ; - if(ierflg != 0){ - cout << "PHOS Unfolding> Unable to set initial value for fit procedure : x = " << x << endl ; - return kFALSE; - } - gMinuit->mnparm(index, "z", z, 0.1, 0, 0, ierflg) ; - index++ ; - if(ierflg != 0){ - cout << "PHOS Unfolding> Unable to set initial value for fit procedure : z = " << z << endl ; - return kFALSE; - } - gMinuit->mnparm(index, "Energy", energy , 0.05*energy, 0., 4.*energy, ierflg) ; - index++ ; - if(ierflg != 0){ - cout << "PHOS Unfolding> Unable to set initial value for fit procedure : energy = " << energy << endl ; - return kFALSE; + if(!module) continue; + AliDebug(1,Form("\t\tTrack hit PHOS! Module: %d, (x,z)=(%.3f,%.3f)",module,xEMC,zEMC)); + + if(emcClu->GetPHOSMod() != module) continue; + + // match track to EMC cluster + emcClu->GetLocalPosition(locClu); + + Float_t delta_x = xEMC - locClu.X(); + Float_t delta_z = zEMC - locClu.Z(); + dr = TMath::Sqrt(delta_x*delta_x + delta_z*delta_z); + AliDebug(1,Form("\tMatch iTrack=%d: (dx,dz)=(%.3f,%.3f)",iTrack,delta_x,delta_z)); + + if(drmnexcm("SET STR", &p2, 0, ierflg) ; // force TgMinuit to reduce function calls - gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient - gMinuit->SetMaxIterations(5); - gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ; // No Warnings - gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize - if(ierflg == 4){ // Minimum not found - cout << "PHOS Unfolding> Fit not converged, cluster abandoned "<< endl ; - return kFALSE ; - } - for(index = 0; index < nPar; index++){ - Double_t err ; - Double_t val ; - gMinuit->GetParameter(index, val, err) ; // Returns value and error of parameter index - fitparameters[index] = val ; - } - return kTRUE; - -} - -//____________________________________________________________________________ -void AliPHOSTrackSegmentMakerv1::FillOneModule(DigitsList * dl, - AliPHOSRecPoint::RecPointsList * emcIn, - TObjArray * emcOut, - AliPHOSRecPoint::RecPointsList * ppsdIn, - TObjArray * ppsdOutUp, - TObjArray * ppsdOutLow, - Int_t & phosmod, - Int_t & emcStopedAt, - Int_t & ppsdStopedAt) -{ - // Unfold clusters and fill xxxOut arrays with clusters from one PHOS module - - AliPHOSEmcRecPoint * emcRecPoint ; - AliPHOSPpsdRecPoint * ppsdRecPoint ; - Int_t index ; + + if(trackindex>=0) + AliDebug(1,Form("\t\tBest match for (xClu,zClu,eClu)=(%.3f,%.3f,%.3f): iTrack=%d, dR=%.3f", + locClu.X(),locClu.Z(),emcClu->GetEnergy(), + trackindex,TMath::Sqrt(dx*dx+dz*dz))); + return; + + Float_t distance2Track = fRtpc ; + + trackindex = -1 ; // closest track within fRCpv + + TVector3 vecEmc ; // Local position of EMC recpoint + TVector3 vecPloc ; // Momentum direction at CPV plain + + //toofar = kTRUE ; + if(emcClu->GetPHOSMod() != cpvClu->GetPHOSMod()){ + dx=999. ; + dz=999. ; + return ; + } + + emcClu->GetLocalPosition(vecEmc) ; - Int_t nEmcUnfolded = emcIn->GetEntries() ; - for(index = emcStopedAt; index < nEmcUnfolded; index++){ - emcRecPoint = (AliPHOSEmcRecPoint *) emcIn->At(index) ; + Double_t xCPV,zCPV ; //EMC-projected coordinates of CPV cluster + TVector3 cpvGlobal; // Global position of the CPV recpoint + fGeom->GetGlobalPHOS((AliPHOSRecPoint*)cpvClu,cpvGlobal); + Double_t vtxCPV[3]={cpvGlobal.X(),cpvGlobal.Y(),cpvGlobal.Z()} ; + Int_t dummyMod ; + + if (fESD == 0x0) { + //if no track information available, assume straight line from IP to emcal + fGeom->ImpactOnEmc(vtxCPV,cpvGlobal.Theta(),cpvGlobal.Phi(),dummyMod,zCPV,xCPV) ; + dx=xCPV - vecEmc.X() ; + dz=zCPV - vecEmc.Z() ; + return ; + } + + //if there is ESD try to correct distance using TPC information on particle direct in CPV + if (fESD != 0x0) { + + Double_t rCPV = cpvGlobal.Pt() ;// Radius from IP to current point + + // Extrapolate the global track direction if any to CPV and find the closest track + Int_t iClosestTrack = -1; + TVector3 inPHOS ; + + for (Int_t iTrack=0; iTrackGetTrack(iTrack); + if (!track->GetXYZAt(rCPV, fESD->GetMagneticField(), xyz)) + continue; //track coord on the cylinder of PHOS radius + if ((TMath::Abs(xyz[0])+TMath::Abs(xyz[1])+TMath::Abs(xyz[2]))<=0) + continue; + //Check if this track hits PHOS + inPHOS.SetXYZ(xyz[0],xyz[1],xyz[2]); + distance2Track = inPHOS.Angle(cpvGlobal) ; + // Find the closest track to the CPV recpoint + if (distance2Track < minDistance) { + minDistance = distance2Track; + iClosestTrack = iTrack; + } + } - if(emcRecPoint->GetPHOSMod() != phosmod ) - break ; - - Int_t nMultipl = emcRecPoint->GetMultiplicity() ; - Int_t * maxAt = new Int_t[nMultipl] ; - Float_t * maxAtEnergy = new Float_t[nMultipl] ; - Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy) ; - - if(nMax <= 1 ) // if cluster is very flat (no pronounced maximum) then nMax = 0 - emcOut->Add(emcRecPoint) ; - else if (fUnfoldFlag) { - UnfoldClusters(dl, emcIn, emcRecPoint, nMax, maxAt, maxAtEnergy, emcOut) ; - emcIn->Remove(emcRecPoint); - emcIn->Compress() ; - nEmcUnfolded-- ; - index-- ; + if (iClosestTrack != -1) { + track = fESD->GetTrack(iClosestTrack); + if (track->GetPxPyPzAt(rCPV, fESD->GetMagneticField(), pxyz)) { // track momentum ibid. + vecP.SetXYZ(pxyz[0],pxyz[1],pxyz[2]); + fGeom->ImpactOnEmc(vtxCPV,vecP.Theta(),vecP.Phi(),dummyMod,zCPV,xCPV) ; + } } - delete[] maxAt ; - delete[] maxAtEnergy ; + if(minDistance < fRtpc ){ + trackindex = iClosestTrack ; + } + } + if(trackindex!=-1){ + // If the closest global track is found, calculate EMC-CPV distance from it + dx=xCPV - vecEmc.X() ; + dz=zCPV - vecEmc.Z() ; } - emcStopedAt = index ; - - for(index = ppsdStopedAt; index < ppsdIn->GetEntries(); index++){ - ppsdRecPoint = (AliPHOSPpsdRecPoint *) ppsdIn->At(index) ; - if(ppsdRecPoint->GetPHOSMod() != phosmod ) - break ; - if(ppsdRecPoint->GetUp() ) - ppsdOutUp->Add(ppsdRecPoint) ; - else - ppsdOutLow->Add(ppsdRecPoint) ; + else{ + // If no global track was found, just take the nearest CPV point + fGeom->ImpactOnEmc(vtxCPV,cpvGlobal.Theta(),cpvGlobal.Phi(),dummyMod,zCPV,xCPV) ; + dx=xCPV - vecEmc.X() ; + dz=zCPV - vecEmc.Z() ; } - ppsdStopedAt = index ; - - emcOut->Sort() ; - ppsdOutUp->Sort() ; - ppsdOutLow->Sort() ; + return ; } //____________________________________________________________________________ -Float_t AliPHOSTrackSegmentMakerv1::GetDistanceInPHOSPlane(AliPHOSEmcRecPoint * emcclu,AliPHOSPpsdRecPoint * PpsdClu, Bool_t &toofar) +void AliPHOSTrackSegmentMakerv1::Init() { - // Calculates the distance between the EMC RecPoint and the PPSD RecPoint - - Float_t r = fR0 ; - - TVector3 vecEmc ; - TVector3 vecPpsd ; + // Make all memory allocations that are not possible in default constructor - emcclu->GetLocalPosition(vecEmc) ; - PpsdClu->GetLocalPosition(vecPpsd) ; - if(emcclu->GetPHOSMod() == PpsdClu->GetPHOSMod()){ - if(vecPpsd.X() >= vecEmc.X() - fDelta ){ - if(vecPpsd.Z() >= vecEmc.Z() - fDelta ){ - AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ; - // Correct to difference in CPV and EMC position due to different distance to center. - // we assume, that particle moves from center - Float_t dCPV = geom->GetIPtoOuterCoverDistance(); - Float_t dEMC = geom->GetIPtoCrystalSurface() ; - dEMC = dEMC / dCPV ; - vecPpsd = dEMC * vecPpsd - vecEmc ; - r = vecPpsd.Mag() ; - } // if zPpsd >= zEmc - fDelta - toofar = kFALSE ; - } // if xPpsd >= xEmc - fDelta - else - toofar = kTRUE ; - } - else - toofar = kTRUE ; - - return r ; + fLinkUpArray = new TClonesArray("AliPHOSLink", 1000); + fTrackSegments = new TClonesArray("AliPHOSTrackSegment",100); + fTrackSegments->SetName("TRACKS"); } //____________________________________________________________________________ -void AliPHOSTrackSegmentMakerv1::MakeLinks(TObjArray * emcRecPoints, TObjArray * ppsdRecPointsUp, - TObjArray * ppsdRecPointsLow, TClonesArray * linklowArray, - TClonesArray *linkupArray) +void AliPHOSTrackSegmentMakerv1::InitParameters() +{ + //Initializes parameters + fRcpv = 10. ; + fRtpc = 4. ; + fEmcFirst = 0 ; + fEmcLast = 0 ; + fCpvFirst = 0 ; + fCpvLast = 0 ; + fLinkUpArray = 0 ; + fWrite = kTRUE ; +} + + +//____________________________________________________________________________ +void AliPHOSTrackSegmentMakerv1::MakeLinks()const { - // Finds distances (links) between all EMC and PPSD clusters, which are not further apart from each other than fR0 - - TIter nextEmc(emcRecPoints) ; - Int_t iEmcClu = 0 ; + // Finds distances (links) between all EMC and CPV clusters, + // which are not further apart from each other than fRcpv + // and sort them in accordance with this distance - AliPHOSPpsdRecPoint * ppsdlow ; - AliPHOSPpsdRecPoint * ppsdup ; + fLinkUpArray->Clear() ; + + AliPHOSCpvRecPoint * cpv ; AliPHOSEmcRecPoint * emcclu ; - - Int_t iLinkLow = 0 ; + Int_t iLinkUp = 0 ; - while( (emcclu = (AliPHOSEmcRecPoint*)nextEmc() ) ) { - Bool_t toofar ; - TIter nextPpsdLow(ppsdRecPointsLow ) ; - Int_t iPpsdLow = 0 ; - - while( (ppsdlow = (AliPHOSPpsdRecPoint*)nextPpsdLow() ) ) { - Float_t r = GetDistanceInPHOSPlane(emcclu, ppsdlow, toofar) ; - - if(toofar) - break ; - if(r < fR0){ - new( (*linklowArray)[iLinkLow++]) AliPHOSLink(r, iEmcClu, iPpsdLow) ; - } - iPpsdLow++ ; - - } - - TIter nextPpsdUp(ppsdRecPointsUp ) ; - Int_t iPpsdUp = 0 ; - - while( (ppsdup = (AliPHOSPpsdRecPoint*)nextPpsdUp() ) ) { - Float_t r = GetDistanceInPHOSPlane(emcclu, ppsdup, toofar) ; - - if(toofar) - break ; - if(r < fR0) { - new( (*linkupArray)[iLinkUp++]) AliPHOSLink(r, iEmcClu, iPpsdUp) ; - } - iPpsdUp++ ; + Int_t iEmcRP; + for(iEmcRP = fEmcFirst; iEmcRP < fEmcLast; iEmcRP++ ) { + emcclu = dynamic_cast(fEMCRecPoints->At(iEmcRP)) ; + + //Bool_t toofar ; + Int_t iCpv = 0 ; + for(iCpv = fCpvFirst; iCpv < fCpvLast;iCpv++ ) { + cpv = dynamic_cast(fCPVRecPoints->At(iCpv)) ; + Int_t track = -1 ; + Float_t dx,dz ; + GetDistanceInPHOSPlane(emcclu, cpv, track,dx,dz) ; + if(TMath::Sqrt(dx*dx+dz*dz) < fRcpv ){ + new ((*fLinkUpArray)[iLinkUp++]) AliPHOSLink(dx, dz, iEmcRP, iCpv, track) ; + } } - - iEmcClu++ ; - - } // while nextEmC + } - linklowArray->Sort() ; //first links with smallest distances - linkupArray->Sort() ; + fLinkUpArray->Sort() ; //first links with smallest distances } - + //____________________________________________________________________________ -void AliPHOSTrackSegmentMakerv1::MakePairs(TObjArray * emcRecPoints, - TObjArray * ppsdRecPointsUp, - TObjArray * ppsdRecPointsLow, - TClonesArray * linklowArray, - TClonesArray * linkupArray, - AliPHOSTrackSegment::TrackSegmentsList * trsl) +void AliPHOSTrackSegmentMakerv1::MakePairs() { - - // Finds the smallest links and makes pairs of PPSD and EMC clusters with smallest distance + // Using the previously made list of "links", we found the smallest link - i.e. + // link with the least distance between EMC and CPV and pointing to still + // unassigned RecParticles. We assign these RecPoints to TrackSegment and + // remove them from the list of "unassigned". + + //Make arrays to mark clusters already chosen + Int_t * emcExist = 0; + if(fEmcLast > fEmcFirst) + emcExist = new Int_t[fEmcLast-fEmcFirst] ; + + Int_t index; + for(index = 0; index fCpvFirst) + cpvExist = new Bool_t[fCpvLast-fCpvFirst] ; + for(index = 0; index (nextUp()) ) ){ - AliPHOSEmcRecPoint * emc ; - AliPHOSPpsdRecPoint * ppsdLow ; - AliPHOSPpsdRecPoint * ppsdUp ; - - AliPHOSRecPoint * nullpointer = 0 ; + if(emcExist[linkUp->GetEmc()-fEmcFirst] != -1){ - while ( (linkLow = (AliPHOSLink *)nextLow() ) ){ - emc = (AliPHOSEmcRecPoint *) emcRecPoints->At(linkLow->GetEmc()) ; - ppsdLow = (AliPHOSPpsdRecPoint *) ppsdRecPointsLow->At(linkLow->GetPpsd()) ; - if( (emc) && (ppsdLow) ){ // RecPoints not removed yet - ppsdUp = 0 ; - - while ( (linkUp = (AliPHOSLink *)nextUp() ) ){ - if(linkLow->GetEmc() == linkUp->GetEmc() ){ - ppsdUp = (AliPHOSPpsdRecPoint *) ppsdRecPointsUp->At(linkUp->GetPpsd()) ; - break ; - } - - } // while nextUp - - nextUp.Reset(); -// AliPHOSTrackSegment * subtr = new AliPHOSTrackSegment(emc, ppsdUp, ppsdLow ) ; -// trsl->Add(subtr) ; - new( (*trsl)[fNTrackSegments] ) AliPHOSTrackSegment(emc, ppsdUp, ppsdLow ) ; - fNTrackSegments++ ; - emcRecPoints->AddAt(nullpointer,linkLow->GetEmc()) ; - ppsdRecPointsLow->AddAt(nullpointer,linkLow->GetPpsd()) ; - - if(ppsdUp) - ppsdRecPointsUp->AddAt(nullpointer,linkUp->GetPpsd()) ; + if(cpvExist[linkUp->GetCpv()-fCpvFirst]){ //CPV still exist + Float_t dx,dz ; + linkUp->GetXZ(dx,dz) ; + new ((* fTrackSegments)[fNTrackSegments]) + AliPHOSTrackSegment(dynamic_cast(fEMCRecPoints->At(linkUp->GetEmc())) , + dynamic_cast(fCPVRecPoints->At(linkUp->GetCpv())) , + linkUp->GetTrack(),dx,dz) ; + (dynamic_cast(fTrackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments); + fNTrackSegments++ ; + emcExist[linkUp->GetEmc()-fEmcFirst] = -1 ; //Mark emc that Cpv was found + //mark CPV recpoint as already used + cpvExist[linkUp->GetCpv()-fCpvFirst] = kFALSE ; + } //if CpvUp still exist } - } - - TIter nextEmc(emcRecPoints) ; - nextEmc.Reset() ; - - while( (emc = (AliPHOSEmcRecPoint*)nextEmc()) ){ //to create pairs if no ppsdlow - ppsdLow = 0 ; - ppsdUp = 0 ; - - while ( (linkUp = (AliPHOSLink *)nextUp() ) ){ - - if(emcRecPoints->IndexOf(emc) == linkUp->GetEmc() ){ - ppsdUp = (AliPHOSPpsdRecPoint *) ppsdRecPointsUp->At(linkUp->GetPpsd()) ; - break ; - } - + } + + //look through emc recPoints left without CPV + if(emcExist){ //if there is emc rec point + Int_t iEmcRP ; + for(iEmcRP = 0; iEmcRP < fEmcLast-fEmcFirst ; iEmcRP++ ){ + if(emcExist[iEmcRP] > 0 ){ + Int_t track = -1 ; + Float_t dx,dz ; + AliPHOSEmcRecPoint *emcclu = dynamic_cast(fEMCRecPoints->At(iEmcRP+fEmcFirst)); + GetDistanceInPHOSPlane(emcclu, 0, track,dx,dz); + if(track<0) + new ((*fTrackSegments)[fNTrackSegments]) AliPHOSTrackSegment(emcclu,nullpointer) ; + else + new ((*fTrackSegments)[fNTrackSegments]) AliPHOSTrackSegment(emcclu,0,track,dx,dz); + (dynamic_cast(fTrackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments); + fNTrackSegments++; + } } -// AliPHOSTrackSegment * subtr = new AliPHOSTrackSegment(emc, ppsdUp, ppsdLow ) ; -// trsl->Add(subtr) ; - new( (*trsl)[fNTrackSegments] ) AliPHOSTrackSegment(emc, ppsdUp, ppsdLow ) ; - fNTrackSegments++ ; - - - if(ppsdUp) - ppsdRecPointsUp->AddAt(nullpointer,linkUp->GetPpsd()) ; } - + delete [] emcExist ; + delete [] cpvExist ; } //____________________________________________________________________________ -void AliPHOSTrackSegmentMakerv1::MakeTrackSegments(DigitsList * dl, - AliPHOSRecPoint::RecPointsList * emcl, - AliPHOSRecPoint::RecPointsList * ppsdl, - AliPHOSTrackSegment::TrackSegmentsList * trsl) +void AliPHOSTrackSegmentMakerv1::Clusters2TrackSegments(Option_t *option) { - // Makes the track segments out of the list of EMC and PPSD Recpoints and stores them in a list - - Int_t phosmod = 1 ; - Int_t emcStopedAt = 0 ; - Int_t ppsdStopedAt = 0 ; - - TObjArray * emcRecPoints = new TObjArray(100) ; // these arrays keep pointers - TObjArray * ppsdRecPointsUp = new TObjArray(100) ; // to RecPoints, which are - TObjArray * ppsdRecPointsLow = new TObjArray(100) ; // kept in TClonesArray's emcl and ppsdl - - - TClonesArray * linklowArray = new TClonesArray("AliPHOSLink", 100); - TClonesArray * linkupArray = new TClonesArray("AliPHOSLink", 100); - - AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ; + // Steering method to perform track segment construction for the current event + // Returns an array with the found track-segments. - while(phosmod <= geom->GetNModules() ){ - - FillOneModule(dl, emcl, emcRecPoints, ppsdl, ppsdRecPointsUp, ppsdRecPointsLow, phosmod, emcStopedAt, ppsdStopedAt) ; - - MakeLinks(emcRecPoints, ppsdRecPointsUp, ppsdRecPointsLow, linklowArray, linkupArray) ; - - MakePairs(emcRecPoints, ppsdRecPointsUp, ppsdRecPointsLow, linklowArray, linkupArray, trsl) ; + if(strstr(option,"tim")) + gBenchmark->Start("PHOSTSMaker"); - emcRecPoints->Clear() ; - - ppsdRecPointsUp->Clear() ; - - ppsdRecPointsLow->Clear() ; - - linkupArray->Clear() ; - - linklowArray->Clear() ; - - phosmod++ ; + if(strstr(option,"print")) { + Print() ; + return ; } - delete emcRecPoints ; - emcRecPoints = 0 ; - - delete ppsdRecPointsUp ; - ppsdRecPointsUp = 0 ; - - delete ppsdRecPointsLow ; - ppsdRecPointsLow = 0 ; - - delete linkupArray ; - linkupArray = 0 ; - - delete linklowArray ; - linklowArray = 0 ; -} - -//____________________________________________________________________________ -Double_t AliPHOSTrackSegmentMakerv1::ShowerShape(Double_t r) -{ - // Shape of the shower (see PHOS TDR) - // If you change this function, change also the gradien evaluation in ChiSquare() - - Double_t r4 = r*r*r*r ; - Double_t r295 = TMath::Power(r, 2.95) ; - Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ; - return shape ; -} + + //Make some initializations + fNTrackSegments = 0 ; + fEmcFirst = 0 ; + fEmcLast = 0 ; + fCpvFirst = 0 ; + fCpvLast = 0 ; -//____________________________________________________________________________ -void AliPHOSTrackSegmentMakerv1::UnfoldClusters(DigitsList * dl, - AliPHOSRecPoint::RecPointsList * emcIn, - AliPHOSEmcRecPoint * iniEmc, - Int_t nMax, - int * maxAt, - Float_t * maxAtEnergy, - TObjArray * emcList) -{ - // Performs the unfolding of a cluster with nMax overlapping showers - // This is time consuming (use the (Un)SetUnfolFlag() ) + fTrackSegments->Clear(); - Int_t nPar = 3 * nMax ; - Float_t * fitparameters = new Float_t[nPar] ; - AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ; + // if(!ReadRecPoints(ievent)) continue; //reads RecPoints for event ievent - Bool_t rv = FindFit(iniEmc, maxAt, maxAtEnergy, nPar, fitparameters) ; - if( !rv ) { - // Fit failed, return and remove cluster - delete[] fitparameters ; - return ; + for(fModule = 1; fModule <= fGeom->GetNModules() ; fModule++ ) { + FillOneModule() ; + MakeLinks() ; + MakePairs() ; } - - Float_t xDigit ; - Float_t zDigit ; - Int_t relid[4] ; - - Int_t nDigits = iniEmc->GetMultiplicity() ; - Float_t xpar ; - Float_t zpar ; - Float_t epar ; - Float_t distance ; - Float_t ratio ; - Float_t * efit = new Float_t[nDigits] ; - Int_t iparam ; - Int_t iDigit ; - - AliPHOSDigit * digit ; - AliPHOSEmcRecPoint * emcRP ; - Int_t * emcDigits = iniEmc->GetDigitsList() ; - Float_t * emcEnergies = iniEmc->GetEnergiesList() ; - - Int_t iRecPoint = emcIn->GetEntries() ; - - for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){ - digit = (AliPHOSDigit *) emcDigits[iDigit]; - geom->AbsToRelNumbering(digit->GetId(), relid) ; - geom->RelPosInModule(relid, xDigit, zDigit) ; - efit[iDigit] = 0; - iparam = 0 ; - while(iparam < nPar ){ - xpar = fitparameters[iparam] ; - zpar = fitparameters[iparam+1] ; - epar = fitparameters[iparam+2] ; - iparam += 3 ; - distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ; - distance = TMath::Sqrt(distance) ; - efit[iDigit] += epar * AliPHOSTrackSegmentMakerv1::ShowerShape(distance) ; - } - } - - iparam = 0 ; - Float_t eDigit ; - - while(iparam < nPar ){ - xpar = fitparameters[iparam] ; - zpar = fitparameters[iparam+1] ; - epar = fitparameters[iparam+2] ; - iparam += 3 ; - new ((*emcIn)[iRecPoint]) AliPHOSEmcRecPoint( iniEmc->GetLogWeightCut(), iniEmc->GetLocMaxCut() ) ; - emcRP = (AliPHOSEmcRecPoint *) emcIn->At(iRecPoint++); - - for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){ - digit = (AliPHOSDigit *) emcDigits[iDigit]; - geom->AbsToRelNumbering(digit->GetId(), relid) ; - geom->RelPosInModule(relid, xDigit, zDigit) ; - distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ; - distance = TMath::Sqrt(distance) ; - ratio = epar * AliPHOSTrackSegmentMakerv1::ShowerShape(distance) / efit[iDigit] ; - eDigit = emcEnergies[iDigit] * ratio ; - emcRP->AddDigit( *digit, eDigit ) ; - } - - emcList->Add(emcRP) ; + if(strstr(option,"deb")) + PrintTrackSegments(option); + if(strstr(option,"tim")){ + gBenchmark->Stop("PHOSTSMaker"); + Info("Exec", "took %f seconds for making TS", + gBenchmark->GetCpuTime("PHOSTSMaker")); } - - delete[] fitparameters ; - delete[] efit ; - } -//______________________________________________________________________________ -void UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag) +//____________________________________________________________________________ +void AliPHOSTrackSegmentMakerv1::Print(const Option_t *)const { - // Calculates th Chi square for the cluster unfolding minimization - // Number of parameters, Gradient, Chi squared, parameters, what to do - - AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ; + // Print TrackSegmentMaker parameters + + TString message("") ; + if( strcmp(GetName(), "") != 0 ) { + message = "\n======== AliPHOSTrackSegmentMakerv1 ========\n" ; + message += "Making Track segments\n" ; + message += "with parameters:\n" ; + message += " Maximal EMC - CPV distance (cm) %f\n" ; + message += "============================================\n" ; + Info("Print", message.Data(),fRcpv) ; + } + else + Info("Print", "AliPHOSTrackSegmentMakerv1 not initialized ") ; +} - AliPHOSEmcRecPoint * emcRP = (AliPHOSEmcRecPoint *) gMinuit->GetObjectFit() ; // EmcRecPoint to fit - Int_t * emcDigits = emcRP->GetDigitsList() ; - Float_t * emcEnergies = emcRP->GetEnergiesList() ; - fret = 0. ; - Int_t iparam ; +//____________________________________________________________________________ +void AliPHOSTrackSegmentMakerv1::PrintTrackSegments(Option_t * option) +{ + // option deb - prints # of found TrackSegments + // option deb all - prints as well indexed of found RecParticles assigned to the TS - if(iflag == 2) - for(iparam = 0 ; iparam < nPar ; iparam++) - Grad[iparam] = 0 ; // Will evaluate gradient + Info("PrintTrackSegments", "Results from TrackSegmentMaker:") ; + printf(" Found %d TrackSegments\n", fTrackSegments->GetEntriesFast() ); - Double_t efit ; - - AliPHOSDigit * digit ; - Int_t iDigit = 0 ; - - while ( (digit = (AliPHOSDigit *)emcDigits[iDigit] )){ - Int_t relid[4] ; - Float_t xDigit ; - Float_t zDigit ; - geom->AbsToRelNumbering(digit->GetId(), relid) ; - geom->RelPosInModule(relid, xDigit, zDigit) ; - - if(iflag == 2){ // calculate gradient - Int_t iParam = 0 ; - efit = 0 ; - while(iParam < nPar ){ - Double_t distance = (xDigit - x[iParam]) * (xDigit - x[iParam]) ; - iParam++ ; - distance += (zDigit - x[iParam]) * (zDigit - x[iParam]) ; - distance = TMath::Sqrt( distance ) ; - iParam++ ; - efit += x[iParam] * AliPHOSTrackSegmentMakerv1::ShowerShape(distance) ; - iParam++ ; - } - Double_t sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E) - iParam = 0 ; - while(iParam < nPar ){ - Double_t xpar = x[iParam] ; - Double_t zpar = x[iParam+1] ; - Double_t epar = x[iParam+2] ; - Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ); - Double_t shape = sum * AliPHOSTrackSegmentMakerv1::ShowerShape(dr) ; - Double_t r4 = dr*dr*dr*dr ; - Double_t r295 = TMath::Power(dr,2.95) ; - Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) + - 0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ; - - Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ; // Derivative over x - iParam++ ; - Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ; // Derivative over z - iParam++ ; - Grad[iParam] += shape ; // Derivative over energy - iParam++ ; - } - } - efit = 0; - iparam = 0 ; - while(iparam < nPar ){ - Double_t xpar = x[iparam] ; - Double_t zpar = x[iparam+1] ; - Double_t epar = x[iparam+2] ; - iparam += 3 ; - Double_t distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ; - distance = TMath::Sqrt(distance) ; - efit += epar * AliPHOSTrackSegmentMakerv1::ShowerShape(distance) ; - } - fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/emcEnergies[iDigit] ; - // Here we assume, that sigma = sqrt(E) - iDigit++ ; + if(strstr(option,"all")) { // printing found TS + printf("TrackSegment # EMC RP# CPV RP#\n") ; + Int_t index; + for (index = 0 ; index GetEntriesFast() ; index++) { + AliPHOSTrackSegment * ts = (AliPHOSTrackSegment * )fTrackSegments->At(index) ; + printf(" %d %d %d \n", ts->GetIndexInList(), ts->GetEmcIndex(), ts->GetCpvIndex() ) ; + } } } +