X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PHOS%2FAliPHOSTrackSegmentMakerv1.cxx;h=50a1c9c1cbfb314617ce867bef3e379fad96d53a;hb=c02178127774b63798ce29ce6c6e25fd2c106e64;hp=00e65a955029638bc1a70ce7c204035f41c81d1c;hpb=6727be7ec8e9058096aae06e1edcc20bb98835a4;p=u%2Fmrichter%2FAliRoot.git diff --git a/PHOS/AliPHOSTrackSegmentMakerv1.cxx b/PHOS/AliPHOSTrackSegmentMakerv1.cxx index 00e65a95502..50a1c9c1cbf 100644 --- a/PHOS/AliPHOSTrackSegmentMakerv1.cxx +++ b/PHOS/AliPHOSTrackSegmentMakerv1.cxx @@ -12,565 +12,616 @@ * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ - +/* $Id$ */ //_________________________________________________________________________ -// Algorithm class to construct track segments connection RecPoints in -// EMCA and Ppsd. Unfolds also the clusters in EMCA. -//*-- Author : D. Peressounko SUBATECH -////////////////////////////////////////////////////////////////////////////// +// 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) +// // --- ROOT system --- - -#include "TObjArray.h" -#include "TClonesArray.h" -#include "TObjectTable.h" - +#include "TROOT.h" +#include "TFile.h" +#include "TTree.h" +#include "TSystem.h" +#include "TBenchmark.h" // --- Standard library --- -#include -#include +#include +#include // --- AliRoot header files --- #include "AliPHOSTrackSegmentMakerv1.h" +#include "AliPHOSClusterizerv1.h" #include "AliPHOSTrackSegment.h" +#include "AliPHOSCpvRecPoint.h" +#include "AliPHOSPpsdRecPoint.h" #include "AliPHOSLink.h" #include "AliPHOSv0.h" #include "AliRun.h" -extern void UnfoldingChiSquare(Int_t &nPar, Double_t *Grad, Double_t & fret, Double_t *x, Int_t iflag) ; - ClassImp( AliPHOSTrackSegmentMakerv1) //____________________________________________________________________________ - AliPHOSTrackSegmentMakerv1::AliPHOSTrackSegmentMakerv1() + AliPHOSTrackSegmentMakerv1::AliPHOSTrackSegmentMakerv1() : AliPHOSTrackSegmentMaker() { // ctor - 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 ; + SetTitle("version 1") ; + SetName("AliPHOSTrackSegmentMaker") ; + fR0 = 10. ; + fEmcFirst = 0 ; + fEmcLast = 0 ; + fCpvFirst = 0 ; + fCpvLast = 0 ; + fPpsdFirst= 0 ; + fPpsdLast = 0 ; + fLinkLowArray = 0 ; + fLinkUpArray = 0 ; + fIsInitialized = kFALSE ; } - //____________________________________________________________________________ - AliPHOSTrackSegmentMakerv1::~AliPHOSTrackSegmentMakerv1() -{ - // dtor - delete fMinuit ; + AliPHOSTrackSegmentMakerv1:: AliPHOSTrackSegmentMakerv1(const char* headerFile, const char* branchTitle): AliPHOSTrackSegmentMaker() +{ + // ctor + SetTitle("version 1") ; + SetName("AliPHOSTrackSegmentMaker") ; + fR0 = 10. ; + fEmcFirst = 0 ; + fEmcLast = 0 ; + fCpvFirst = 0 ; + fCpvLast = 0 ; + fPpsdFirst= 0 ; + fPpsdLast = 0 ; + + fHeaderFileName = headerFile ; + fRecPointsBranchTitle = branchTitle ; + + TFile * file = (TFile*) gROOT->GetFile(fHeaderFileName.Data() ) ; + + if(file == 0){ + file = new TFile(fHeaderFileName.Data(),"update") ; + gAlice = (AliRun *) file->Get("gAlice") ; + } + + AliPHOS * phos = (AliPHOS *) gAlice->GetDetector("PHOS") ; + fGeom = AliPHOSGeometry::GetInstance(phos->GetGeometry()->GetName(),phos->GetGeometry()->GetTitle() ); + + fEmcRecPoints = new TObjArray(200) ; + fCpvRecPoints = new TObjArray(200) ; + fClusterizer = new AliPHOSClusterizerv1() ; + + fTrackSegments = new TClonesArray("AliPHOSTrackSegment",200) ; + + fLinkLowArray = new TClonesArray("AliPHOSLink", 1000); + fLinkUpArray = new TClonesArray("AliPHOSLink", 1000); + + fIsInitialized = kTRUE ; + } //____________________________________________________________________________ -Bool_t AliPHOSTrackSegmentMakerv1::FindFit(AliPHOSEmcRecPoint * emcRP, int * maxAt, Float_t * maxAtEnergy, - Int_t nPar, Float_t * fitparameters) -{ - // Calls TMinuit for fitting cluster with several maxima +void AliPHOSTrackSegmentMakerv1::Init(){ - AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ; - - gMinuit->SetPrintLevel(-1) ; // No Printout - gMinuit->SetFCN(UnfoldingChiSquare) ; // To set the address of the minimization function - gMinuit->SetObjectFit(emcRP) ; // To tranfer pointer to UnfoldingChiSquare - - // filling initial values for fit parameters - AliPHOSDigit * digit ; - - Int_t ierflg = 0; - Int_t index = 0 ; - Int_t nDigits = (Int_t) nPar / 3 ; - - Int_t iDigit ; - - - for(iDigit = 0; iDigit < nDigits; iDigit++){ - digit = (AliPHOSDigit *) maxAt[iDigit]; + if(!fIsInitialized){ + if(fHeaderFileName.IsNull()) + fHeaderFileName = "galice.root" ; + + + TFile * file = (TFile*) gROOT->GetFile(fHeaderFileName.Data() ) ; + + if(file == 0){ + file = new TFile(fHeaderFileName.Data(),"update") ; + gAlice = (AliRun *) file->Get("gAlice") ; + } + + AliPHOS * phos = (AliPHOS *) gAlice->GetDetector("PHOS") ; + fGeom = AliPHOSGeometry::GetInstance(phos->GetGeometry()->GetName(),phos->GetGeometry()->GetTitle() ); - Int_t relid[4] ; - Float_t x ; - Float_t z ; - geom->AbsToRelNumbering(digit->GetId(), relid) ; - geom->RelPosInModule(relid, x, z) ; - Float_t energy = maxAtEnergy[iDigit] ; + fEmcRecPoints = new TObjArray(200) ; + fCpvRecPoints = new TObjArray(200) ; + fClusterizer = new AliPHOSClusterizerv1() ; - 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; - } - } + + fTrackSegments = new TClonesArray("AliPHOSTrackSegment",200) ; - Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly - // depends on it. - Double_t p1 = 1.0 ; - Double_t p2 = 0.0 ; - - gMinuit->mnexcm("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 ; + fLinkLowArray = new TClonesArray("AliPHOSLink", 1000); + fLinkUpArray = new TClonesArray("AliPHOSLink", 1000); + + fIsInitialized = kTRUE ; } - return kTRUE; +} +//____________________________________________________________________________ + AliPHOSTrackSegmentMakerv1::~AliPHOSTrackSegmentMakerv1() +{ + // dtor + if(fLinkLowArray) delete fLinkLowArray ; + if(fLinkUpArray) delete fLinkUpArray ; } //____________________________________________________________________________ -void AliPHOSTrackSegmentMakerv1::FillOneModule(DigitsList * Dl, RecPointsList * emcIn, TObjArray * emcOut, - RecPointsList * ppsdIn, TObjArray * ppsdOutUp, - TObjArray * ppsdOutLow, Int_t & phosmod, Int_t & emcStopedAt, - Int_t & ppsdStopedAt) +void AliPHOSTrackSegmentMakerv1::FillOneModule() { - // Unfold clusters and fill xxxOut arrays with clusters from one PHOS module + // Finds bounds in which clusters from one PHOS module are - AliPHOSEmcRecPoint * emcRecPoint ; - AliPHOSPpsdRecPoint * ppsdRecPoint ; - Int_t index ; + + //First EMC clusters + Int_t totalEmc = fEmcRecPoints->GetEntriesFast() ; + for(fEmcFirst = fEmcLast; (fEmcLast < totalEmc) && + (((AliPHOSRecPoint *) fEmcRecPoints->At(fEmcLast))->GetPHOSMod() == fModule ); + fEmcLast ++) ; + - Int_t nEmcUnfolded = emcIn->GetEntries() ; - for(index = emcStopedAt; index < nEmcUnfolded; index++){ - emcRecPoint = (AliPHOSEmcRecPoint *) (*emcIn)[index] ; + //Now CPV clusters + Int_t totalCpv = fCpvRecPoints->GetEntriesFast() ; - if(emcRecPoint->GetPHOSMod() != phosmod ) - break ; + if(fModule <= fGeom->GetNCPVModules()){ // in CPV geometry - Int_t nMultipl = emcRecPoint->GetMultiplicity() ; - Int_t maxAt[nMultipl] ; - Float_t maxAtEnergy[nMultipl] ; - Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy) ; + for(fCpvFirst = fCpvLast; (fCpvLast < totalCpv) && + (((AliPHOSRecPoint *) fCpvRecPoints->At(fCpvLast))->GetPHOSMod() == fModule ); + fCpvLast ++) ; - 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-- ; - } + fPpsdFirst = fCpvLast ; //To avoid scanning RecPoints between fPpsdFirst and fPpsdLast + fPpsdLast = fCpvLast ; //and to be ready to switch to mixed geometry } - emcStopedAt = index ; - - for(index = ppsdStopedAt; index < ppsdIn->GetEntries(); index++){ - ppsdRecPoint = (AliPHOSPpsdRecPoint *) (*ppsdIn)[index] ; - if(ppsdRecPoint->GetPHOSMod() != phosmod ) - break ; - if(ppsdRecPoint->GetUp() ) - ppsdOutUp->Add(ppsdRecPoint) ; - else - ppsdOutLow->Add(ppsdRecPoint) ; + else{ //in PPSD geometry + fCpvLast = fPpsdLast ; + //Upper layer first + for(fCpvFirst = fCpvLast; (fCpvLast < totalCpv) && + (((AliPHOSPpsdRecPoint *) fCpvRecPoints->At(fCpvLast))->GetPHOSMod() == fModule ) && + (((AliPHOSPpsdRecPoint *) fCpvRecPoints->At(fCpvLast))->GetUp()) ; + fCpvLast ++) ; + + fPpsdLast= fCpvLast ; + for(fPpsdFirst = fPpsdLast; (fPpsdLast < totalCpv) && + (((AliPHOSPpsdRecPoint *) fCpvRecPoints->At(fPpsdLast))->GetPHOSMod() == fModule ) && + (!((AliPHOSPpsdRecPoint *) fCpvRecPoints->At(fPpsdLast))->GetUp()) ; + fPpsdLast ++) ; } - ppsdStopedAt = index ; - - emcOut->Sort() ; - ppsdOutUp->Sort() ; - ppsdOutLow->Sort() ; + } //____________________________________________________________________________ -Float_t AliPHOSTrackSegmentMakerv1::GetDistanceInPHOSPlane(AliPHOSEmcRecPoint * emcclu,AliPHOSPpsdRecPoint * PpsdClu, Bool_t &toofar) +Float_t AliPHOSTrackSegmentMakerv1::GetDistanceInPHOSPlane(AliPHOSEmcRecPoint * emcClu,AliPHOSRecPoint * cpvClu, Bool_t &toofar) { + // Calculates the distance between the EMC RecPoint and the PPSD RecPoint + //clusters are sorted in "rows" and "columns" of width 1 cm + 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() Float_t r = fR0 ; TVector3 vecEmc ; - TVector3 vecPpsd ; - - 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 + TVector3 vecCpv ; + + emcClu->GetLocalPosition(vecEmc) ; + cpvClu->GetLocalPosition(vecCpv) ; + + if(emcClu->GetPHOSMod() == cpvClu->GetPHOSMod()){ + if(vecCpv.X() <= vecEmc.X() + fR0 + 2*delta ){ + + vecCpv = vecCpv - vecEmc ; + r = vecCpv.Mag() ; toofar = kFALSE ; - } // if xPpsd >= xEmc - fDelta + + } // if xPpsd >= xEmc + ... else toofar = kTRUE ; } else toofar = kTRUE ; + + //toofar = kFALSE ; + return r ; } //____________________________________________________________________________ -void AliPHOSTrackSegmentMakerv1::MakeLinks(TObjArray * emcRecPoints, TObjArray * ppsdRecPointsUp, - TObjArray * ppsdRecPointsLow, TClonesArray * linklowArray, - TClonesArray *linkupArray) +void AliPHOSTrackSegmentMakerv1::MakeLinks() { - //Finds distanses (links) between all EMC and PPSD clusters, which are not further from each other than fR0 - - TIter nextEmc(emcRecPoints) ; - Int_t iEmcClu = 0 ; + // Finds distances (links) between all EMC and PPSD clusters, which are not further apart from each other than fR0 - AliPHOSPpsdRecPoint * ppsdlow ; - AliPHOSPpsdRecPoint * ppsdup ; + fLinkUpArray->Clear() ; + fLinkLowArray->Clear() ; + + AliPHOSRecPoint * ppsd ; + AliPHOSRecPoint * 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) ; + Int_t iEmcRP; + for(iEmcRP = fEmcFirst; iEmcRP < fEmcLast; iEmcRP++ ) { + emcclu = (AliPHOSEmcRecPoint *) fEmcRecPoints->At(iEmcRP) ; + + Bool_t toofar ; + Int_t iPpsd ; + for(iPpsd = fPpsdFirst; iPpsd < fPpsdLast;iPpsd++ ) { + ppsd = (AliPHOSRecPoint *) fCpvRecPoints->At(iPpsd) ; + Float_t r = GetDistanceInPHOSPlane(emcclu, ppsd, toofar) ; + if(toofar) break ; - if(r < fR0) - new( (*linklowArray)[iLinkLow++]) AliPHOSLink(r, iEmcClu, iPpsdLow) ; - - iPpsdLow++ ; - + if(r < fR0) + new ((*fLinkLowArray)[iLinkLow++]) AliPHOSLink(r, iEmcRP, iPpsd) ; } - TIter nextPpsdUp(ppsdRecPointsUp ) ; - Int_t iPpsdUp = 0 ; - - while( (ppsdup = (AliPHOSPpsdRecPoint*)nextPpsdUp() ) ) { - Float_t r = GetDistanceInPHOSPlane(emcclu, ppsdup, toofar) ; + Int_t iCpv = 0 ; + for(iCpv = fCpvFirst; iCpv < fCpvLast;iCpv++ ) { + + cpv = (AliPHOSRecPoint *) fCpvRecPoints->At(iCpv) ; + Float_t r = GetDistanceInPHOSPlane(emcclu, cpv, toofar) ; if(toofar) break ; - if(r < fR0) - new( (*linkupArray)[iLinkUp++]) AliPHOSLink(r, iEmcClu, iPpsdUp) ; - - iPpsdUp++ ; - + if(r < fR0) { + new ((*fLinkUpArray)[iLinkUp++]) AliPHOSLink(r, iEmcRP, iCpv) ; + } } - - iEmcClu++ ; - - } // while nextEmC + } - linklowArray->Sort() ; //first links with smallest distances - linkupArray->Sort() ; + fLinkLowArray->Sort() ; //first links with smallest distances + fLinkUpArray->Sort() ; } - + //____________________________________________________________________________ -void AliPHOSTrackSegmentMakerv1::MakePairs(TObjArray * emcRecPoints, TObjArray * ppsdRecPointsUp, - TObjArray * ppsdRecPointsLow, TClonesArray * linklowArray, - TClonesArray * linkupArray, TrackSegmentsList * trsl) +void AliPHOSTrackSegmentMakerv1::MakePairs() { - - // Finds the smallest links and makes pairs of PPSD and EMC clusters with smallest distance - TIter nextLow(linklowArray) ; - TIter nextUp(linkupArray) ; + //Make arrays to mark clusters already chousen + 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 fPpsdFirst) + ppsdExist = new Bool_t[fPpsdLast-fPpsdFirst] ; + for(index = 0; index 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) ; - emcRecPoints->AddAt(nullpointer,linkLow->GetEmc()) ; - ppsdRecPointsLow->AddAt(nullpointer,linkLow->GetPpsd()) ; - - if(ppsdUp) - ppsdRecPointsUp->AddAt(nullpointer,linkUp->GetPpsd()) ; + + if( (emcExist[linkLow->GetEmc()-fEmcFirst]> 0) && ppsdExist[linkLow->GetPpsd()-fPpsdFirst] ){ // RecPoints not removed yet + new ((*fTrackSegments)[fNTrackSegments]) AliPHOSTrackSegment((AliPHOSEmcRecPoint *) fEmcRecPoints->At(linkLow->GetEmc()), + nullpointer, + (AliPHOSPpsdRecPoint *)fCpvRecPoints->At(linkLow->GetPpsd()) ) ; + ((AliPHOSTrackSegment* )fTrackSegments->At(fNTrackSegments))->SetIndexInList(fNTrackSegments); + //replace index of emc to negative and shifted index of TS + emcExist[linkLow->GetEmc()-fEmcFirst] = -2 - fNTrackSegments ; + //mark ppsd as used + ppsdExist[linkLow->GetPpsd()-fPpsdFirst] = kFALSE ; + fNTrackSegments++ ; } } - - 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 ; - } - + while ( (linkUp = (AliPHOSLink *)nextUp() ) ){ + if(emcExist[linkUp->GetEmc()-fEmcFirst] != -1){ //without ppsd Up yet + + if(cpvExist[linkUp->GetPpsd()-fCpvFirst]){ //CPV still exist + + if(emcExist[linkUp->GetEmc()-fEmcFirst] > 0){ //without ppsd Low => create new TS + + new ((* fTrackSegments)[fNTrackSegments]) AliPHOSTrackSegment((AliPHOSEmcRecPoint *) fEmcRecPoints->At(linkUp->GetEmc()) , + (AliPHOSPpsdRecPoint *)fCpvRecPoints->At(linkUp->GetPpsd()), + nullpointer) ; + ((AliPHOSTrackSegment *) fTrackSegments->At(fNTrackSegments))->SetIndexInList(fNTrackSegments); + fNTrackSegments++ ; + } + else{ // append ppsd Up to existing TS + ((AliPHOSTrackSegment *)fTrackSegments->At(-2-emcExist[linkUp->GetEmc()-fEmcFirst]))->SetCpvRecPoint((AliPHOSCpvRecPoint *)fCpvRecPoints->At(linkUp->GetPpsd())); + } + + emcExist[linkUp->GetEmc()-fEmcFirst] = -1 ; //Mark emc that Cpv was found + //mark CPV recpoint as already used + cpvExist[linkUp->GetPpsd()-fCpvFirst] = kFALSE ; + } //if ppsdUp still exist + } + } + + //look through emc recPoints left without CPV/PPSD + if(emcExist){ //if there is emc rec point + Int_t iEmcRP ; + for(iEmcRP = 0; iEmcRP < fEmcLast-fEmcFirst ; iEmcRP++ ){ + if(emcExist[iEmcRP] > 0 ){ + new ((*fTrackSegments)[fNTrackSegments]) AliPHOSTrackSegment((AliPHOSEmcRecPoint *)fEmcRecPoints->At(iEmcRP+fEmcFirst), + nullpointer, + nullpointer ) ; + ((AliPHOSTrackSegment *) fTrackSegments->At(fNTrackSegments))->SetIndexInList(fNTrackSegments); + fNTrackSegments++; + } } - AliPHOSTrackSegment * subtr = new AliPHOSTrackSegment(emc, ppsdUp, ppsdLow ) ; - trsl->Add(subtr) ; - - if(ppsdUp) - ppsdRecPointsUp->AddAt(nullpointer,linkUp->GetPpsd()) ; } - + } //____________________________________________________________________________ -void AliPHOSTrackSegmentMakerv1::MakeTrackSegments(DigitsList * DL, RecPointsList * emcl, - RecPointsList * ppsdl, TrackSegmentsList * trsl) +void AliPHOSTrackSegmentMakerv1::Exec(Option_t * option) { - // main function, does the job + // 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() ; - - while(phosmod <= geom->GetNModules() ){ - - FillOneModule(DL, emcl, emcRecPoints, ppsdl, ppsdRecPointsUp, ppsdRecPointsLow, phosmod, emcStopedAt, ppsdStopedAt) ; + if(! fIsInitialized) Init() ; - MakeLinks(emcRecPoints, ppsdRecPointsUp, ppsdRecPointsLow, linklowArray, linkupArray) ; + if(strstr(option,"tim")) + gBenchmark->Start("PHOSTSMaker"); - MakePairs(emcRecPoints, ppsdRecPointsUp, ppsdRecPointsLow, linklowArray, linkupArray, trsl) ; - - emcRecPoints->Clear() ; - - ppsdRecPointsUp->Clear() ; + Int_t nEvents = (Int_t) gAlice->TreeE()->GetEntries() ; - ppsdRecPointsLow->Clear() ; + for(fEvent = 0;fEvent< nEvents; fEvent++){ + if(!ReadRecPoints()) //reads RecPoints for event fEvent + return; + + for(fModule = 1; fModule <= fGeom->GetNModules() ; fModule++ ){ + + FillOneModule() ; + + MakeLinks() ; + + MakePairs() ; + + } - linkupArray->Clear() ; - - linklowArray->Clear() ; - - phosmod++ ; + WriteTrackSegments() ; + if(strstr(option,"deb")) + PrintTrackSegments(option) ; } - delete emcRecPoints ; - emcRecPoints = 0 ; - - delete ppsdRecPointsUp ; - ppsdRecPointsUp = 0 ; - delete ppsdRecPointsLow ; - ppsdRecPointsLow = 0 ; + if(strstr(option,"tim")){ + gBenchmark->Stop("PHOSTSMaker"); + cout << "AliPHOSTSMaker:" << endl ; + cout << " took " << gBenchmark->GetCpuTime("PHOSTSMaker") << " seconds for making TS " + << gBenchmark->GetCpuTime("PHOSTSMaker")/nEvents << " seconds per event " << endl ; + cout << endl ; + } - delete linkupArray ; - linkupArray = 0 ; - delete linklowArray ; - linklowArray = 0 ; } - //____________________________________________________________________________ -Double_t AliPHOSTrackSegmentMakerv1::ShowerShape(Double_t r) -{ -// If you change this function, change also gradiend 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 ; +void AliPHOSTrackSegmentMakerv1::Print(Option_t * option)const { + if(fIsInitialized){ + cout << "======== AliPHOSTrackSegmentMakerv1 ========" << endl ; + cout << "Making Track segments "<< endl ; + cout << " Headers file: " << fHeaderFileName.Data() << endl ; + cout << " RecPoints branch file name: " <Clear() ; + fCpvRecPoints->Clear() ; + fTrackSegments->Clear() ; + fNTrackSegments = 0 ; + fEmcFirst = 0 ; + fEmcLast = 0 ; + fCpvFirst = 0 ; + fCpvLast = 0 ; + fPpsdFirst= 0 ; + fPpsdLast = 0 ; + + + gAlice->GetEvent(fEvent) ; + + // Get TreeR header from file + char treeName[20]; + sprintf(treeName,"TreeR%d",fEvent); + + if(gAlice->TreeR()==0){ + cout << "Error in AliPHOSTrackSegmentMakerv1 : no "<TreeR()->GetListOfBranches() ; + Int_t ibranch; + Bool_t emcNotFound = kTRUE ; + Bool_t cpvNotFound = kTRUE ; + Bool_t clusterizerNotFound = kTRUE ; - 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[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 ; + for(ibranch = 0;ibranch GetEntries();ibranch++){ + + if(emcNotFound){ + emcBranch=(TBranch *) branches->At(ibranch) ; + if( fRecPointsBranchTitle.CompareTo(emcBranch->GetTitle())==0 ) + if( strcmp(emcBranch->GetName(),emcBranchName) == 0) { + emcNotFound = kFALSE ; + } + } + + if(cpvNotFound){ + cpvBranch=(TBranch *) branches->At(ibranch) ; + if( fRecPointsBranchTitle.CompareTo(cpvBranch->GetTitle())==0 ) + if( strcmp(cpvBranch->GetName(),cpvBranchName) == 0) + cpvNotFound = kFALSE ; + } - 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 * ShowerShape(distance) ; + if(clusterizerNotFound){ + clusterizerBranch = (TBranch *) branches->At(ibranch) ; + if( fRecPointsBranchTitle.CompareTo(clusterizerBranch->GetTitle()) == 0) + if( strcmp(clusterizerBranch->GetName(),cluBranchName) == 0) + clusterizerNotFound = kFALSE ; } + } - 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 * ShowerShape(distance) / efit[iDigit] ; - eDigit = emcEnergies[iDigit] * ratio ; - emcRP->AddDigit( *digit, eDigit ) ; - } + if(clusterizerNotFound || emcNotFound || cpvNotFound){ + cout << "AliPHOSTrackSegmentMakerv1: " << endl ; + cout << " Can't find Branch with RecPoints or Clusterizer " ; + cout << " Do nothing" <SetAddress(&fEmcRecPoints) ; + cpvBranch->SetAddress(&fCpvRecPoints) ; + clusterizerBranch->SetAddress(&fClusterizer) ; + + gAlice->TreeR()->GetEvent(0) ; + + delete emcBranchName; + delete cpvBranchName; + delete cluBranchName; - emcList->Add(emcRP) ; + return kTRUE ; + +} +//____________________________________________________________________________ +void AliPHOSTrackSegmentMakerv1::WriteTrackSegments(){ + char treeName[20]; + sprintf(treeName,"TreeR%d",fEvent); + + + //First, check, if branches already exist + TBranch * tsMakerBranch = 0; + TBranch * tsBranch = 0; + + TObjArray * branches = gAlice->TreeR()->GetListOfBranches() ; + Int_t ibranch; + Bool_t tsMakerNotFound = kTRUE ; + Bool_t tsNotFound = kTRUE ; + + for(ibranch = 0;(ibranch GetEntries())&&(tsMakerNotFound||tsNotFound);ibranch++){ + if(tsMakerNotFound){ + tsMakerBranch=(TBranch *) branches->At(ibranch) ; + if( (strcmp(tsMakerBranch->GetName(),"AliPHOSTrackSegmentMaker") == 0) && + (fTSBranchTitle.CompareTo( tsMakerBranch->GetTitle())==0 )) + tsMakerNotFound = kFALSE ; + } + if(tsNotFound){ + tsBranch=(TBranch *) branches->At(ibranch) ; + if( (strcmp(tsBranch->GetName(),"PHOSTS") == 0) && + (fTSBranchTitle.CompareTo( tsBranch->GetTitle())==0 )) + tsNotFound = kFALSE ; + } } -} -//______________________________________________________________________________ -void UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag) -{ + if(!(tsMakerNotFound && tsNotFound )){ + cout << "AliPHOSTrackSegmentMakerv1 error:"<< endl ; + cout << " Branches PHOSTS and AliPHOSTrackSegementMaker " << endl ; + cout << " with title '"<Getenv("CONFIG_SPLIT_FILE")!=0){ //generating file name + filename = new char[strlen(gAlice->GetBaseFile())+20] ; + sprintf(filename,"%s/PHOS.Reco.root",gAlice->GetBaseFile()) ; + } - AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ; + TDirectory *cwd = gDirectory; + + //First TS + Int_t bufferSize = 32000 ; + tsBranch = gAlice->TreeR()->Branch("PHOSTS",&fTrackSegments,bufferSize); + tsBranch->SetTitle(fTSBranchTitle.Data()); + if (filename) { + tsBranch->SetFile(filename); + TIter next( tsBranch->GetListOfBranches()); + while ((tsBranch=(TBranch*)next())) { + tsBranch->SetFile(filename); + } + cwd->cd(); + } + + //Second -TSMaker + Int_t splitlevel = 0 ; + AliPHOSTrackSegmentMakerv1 * ts = this ; + tsMakerBranch = gAlice->TreeR()->Branch("AliPHOSTrackSegmentMaker","AliPHOSTrackSegmentMakerv1", + &ts,bufferSize,splitlevel); + tsMakerBranch->SetTitle(fTSBranchTitle.Data()); + if (filename) { + tsMakerBranch->SetFile(filename); + TIter next( tsMakerBranch->GetListOfBranches()); + while ((tsMakerBranch=(TBranch*)next())) { + tsMakerBranch->SetFile(filename); + } + cwd->cd(); + } + + gAlice->TreeR()->Fill() ; + gAlice->TreeR()->Write(0,kOverwrite) ; + +} - AliPHOSEmcRecPoint * emcRP = (AliPHOSEmcRecPoint *) gMinuit->GetObjectFit() ; // EmcRecPoint to fit - Int_t * emcDigits = emcRP->GetDigitsList() ; - Float_t * emcEnergies = emcRP->GetEnergiesList() ; - fret = 0. ; - Int_t iparam ; - if(iflag == 2) - for(iparam = 0 ; iparam < nPar ; iparam++) - Grad[iparam] = 0 ; // Will evaluate gradient +//____________________________________________________________________________ +void AliPHOSTrackSegmentMakerv1::PrintTrackSegments(Option_t * option){ - Double_t efit ; + cout << "AliPHOSTrackSegmentMakerv1: " << endl ; + cout << " Found " << fTrackSegments->GetEntriesFast() << " trackSegments " << endl ; - 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(strstr(option,"all")) { // printing found TS + cout << "TrackSegment # " << " EMC RP# " << " CPV RP# " << " PPSD RP#" << endl ; - 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++ ; + Int_t index; + for (index = 0 ; index GetEntriesFast() ; index++) { + AliPHOSTrackSegment * ts = (AliPHOSTrackSegment * )fTrackSegments->At(index) ; + cout<<" "<< setw(4) << ts->GetIndexInList() << " " + <GetEmcIndex()<< " " + <GetCpvIndex()<< " " + <GetPpsdIndex()<< endl ; + } + + cout << "-------------------------------------------------------"<< endl ; } } +//____________________________________________________________________________ +void AliPHOSTrackSegmentMakerv1::SetRecPointsBranch(const char * title){ + //set the title of RecPoints + fRecPointsBranchTitle = title ; + +} +//____________________________________________________________________________ +void AliPHOSTrackSegmentMakerv1::SetTrackSegmentsBranch(const char * title){ + + fTSBranchTitle = title ; +}