X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=ITS%2FAliITSMultReconstructor.cxx;h=899d404b02a0ce8e1f98ca8adce7f6313df7c171;hb=73f560964499de8282b5842b7c38374ef65ebd6f;hp=5bc620179d38b2829e073ab0de46b8392bc68213;hpb=aba7fa71e8c9ab56f61bac1d1b3c4ea70af83a52;p=u%2Fmrichter%2FAliRoot.git diff --git a/ITS/AliITSMultReconstructor.cxx b/ITS/AliITSMultReconstructor.cxx index 5bc620179d3..899d404b02a 100644 --- a/ITS/AliITSMultReconstructor.cxx +++ b/ITS/AliITSMultReconstructor.cxx @@ -13,8 +13,6 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* $Id$ */ - //_________________________________________________________________________ // // Implementation of the ITS-SPD trackleter class @@ -54,24 +52,42 @@ // - update to the new algorithm by Mariella and Jan Fiete // - store also DeltaTheta in the ESD // - less new and delete calls when creating the needed arrays +// +// - RS: to decrease the number of new/deletes the clusters data are stored +// not in float[6] attached to float**, but in 1-D array. +// - RS: Clusters are sorted in Z in roder to have the same numbering as in the ITS reco +// - RS: Clusters used by ESDtrack are flagged, this information is passed to AliMulitiplicity object +// when storing the tracklets and single cluster info +// - MN: first MC label of single clusters stored //_________________________________________________________________________ #include #include #include #include -#include "TArrayI.h" +#include +#include +#include #include "AliITSMultReconstructor.h" #include "AliITSReconstructor.h" -#include "AliITSsegmentationSPD.h" #include "AliITSRecPoint.h" #include "AliITSRecPointContainer.h" #include "AliITSgeom.h" #include "AliITSgeomTGeo.h" +#include "AliITSDetTypeRec.h" +#include "AliESDEvent.h" +#include "AliESDVertex.h" +#include "AliESDtrack.h" +#include "AliMultiplicity.h" #include "AliLog.h" #include "TGeoGlobalMagField.h" #include "AliMagF.h" +#include "AliESDv0.h" +#include "AliV0.h" +#include "AliKFParticle.h" +#include "AliKFVertex.h" +#include "AliRefArray.h" //____________________________________________________________________ ClassImp(AliITSMultReconstructor) @@ -79,25 +95,43 @@ ClassImp(AliITSMultReconstructor) //____________________________________________________________________ AliITSMultReconstructor::AliITSMultReconstructor(): -TObject(), -fClustersLay1(0), -fClustersLay2(0), -fDetectorIndexClustersLay1(0), -fDetectorIndexClustersLay2(0), -fOverlapFlagClustersLay1(0), -fOverlapFlagClustersLay2(0), +fDetTypeRec(0),fESDEvent(0),fTreeRP(0),fTreeRPMix(0), fTracklets(0), fSClusters(0), -fNClustersLay1(0), -fNClustersLay2(0), fNTracklets(0), fNSingleCluster(0), -fPhiWindow(0), -fThetaWindow(0), +fDPhiWindow(0), +fDThetaWindow(0), fPhiShift(0), fRemoveClustersFromOverlaps(0), fPhiOverlapCut(0), fZetaOverlapCut(0), +fPhiRotationAngle(0), +fScaleDTBySin2T(0), +fNStdDev(1.0), +fNStdDevSq(1.0), +// +fCutPxDrSPDin(0.1), +fCutPxDrSPDout(0.15), +fCutPxDz(0.2), +fCutDCArz(0.5), +fCutMinElectronProbTPC(0.5), +fCutMinElectronProbESD(0.1), +fCutMinP(0.05), +fCutMinRGamma(2.), +fCutMinRK0(1.), +fCutMinPointAngle(0.98), +fCutMaxDCADauther(0.5), +fCutMassGamma(0.03), +fCutMassGammaNSigma(5.), +fCutMassK0(0.03), +fCutMassK0NSigma(5.), +fCutChi2cGamma(2.), +fCutChi2cK0(2.), +fCutGammaSFromDecay(-10.), +fCutK0SFromDecay(-10.), +fCutMaxDCA(1.), +// fHistOn(0), fhClustersDPhiAcc(0), fhClustersDThetaAcc(0), @@ -108,15 +142,34 @@ fhDPhiVsDThetaAcc(0), fhetaTracklets(0), fhphiTracklets(0), fhetaClustersLay1(0), -fhphiClustersLay1(0){ - - fNFiredChips[0] = 0; - fNFiredChips[1] = 0; - +fhphiClustersLay1(0), +// + fDPhiShift(0), + fDPhiWindow2(0), + fDThetaWindow2(0), + fPartners(0), + fAssociatedLay1(0), + fMinDists(0), + fBlackList(0), +// + fCreateClustersCopy(0), + fClustersLoaded(0), + fRecoDone(0), + fSPDSeg() +{ + // default c-tor + for (int i=0;i<2;i++) { + fNFiredChips[i] = 0; + fClArr[i] = 0; + for (int j=0;j<2;j++) fUsedClusLay[i][j] = 0; + fDetectorIndexClustersLay[i] = 0; + fOverlapFlagClustersLay[i] = 0; + fNClustersLay[i] = 0; + fClustersLay[i] = 0; + } // Method to reconstruct the charged particles multiplicity with the // SPD (tracklets). - - + SetHistOn(); if(AliITSReconstructor::GetRecoParam()) { @@ -126,6 +179,31 @@ fhphiClustersLay1(0){ SetRemoveClustersFromOverlaps(AliITSReconstructor::GetRecoParam()->GetTrackleterRemoveClustersFromOverlaps()); SetPhiOverlapCut(AliITSReconstructor::GetRecoParam()->GetTrackleterPhiOverlapCut()); SetZetaOverlapCut(AliITSReconstructor::GetRecoParam()->GetTrackleterZetaOverlapCut()); + SetPhiRotationAngle(AliITSReconstructor::GetRecoParam()->GetTrackleterPhiRotationAngle()); + SetNStdDev(AliITSReconstructor::GetRecoParam()->GetTrackleterNStdDevCut()); + SetScaleDThetaBySin2T(AliITSReconstructor::GetRecoParam()->GetTrackleterScaleDThetaBySin2T()); + // + SetCutPxDrSPDin(AliITSReconstructor::GetRecoParam()->GetMultCutPxDrSPDin()); + SetCutPxDrSPDout(AliITSReconstructor::GetRecoParam()->GetMultCutPxDrSPDout()); + SetCutPxDz(AliITSReconstructor::GetRecoParam()->GetMultCutPxDz()); + SetCutDCArz(AliITSReconstructor::GetRecoParam()->GetMultCutDCArz()); + SetCutMinElectronProbTPC(AliITSReconstructor::GetRecoParam()->GetMultCutMinElectronProbTPC()); + SetCutMinElectronProbESD(AliITSReconstructor::GetRecoParam()->GetMultCutMinElectronProbESD()); + SetCutMinP(AliITSReconstructor::GetRecoParam()->GetMultCutMinP()); + SetCutMinRGamma(AliITSReconstructor::GetRecoParam()->GetMultCutMinRGamma()); + SetCutMinRK0(AliITSReconstructor::GetRecoParam()->GetMultCutMinRK0()); + SetCutMinPointAngle(AliITSReconstructor::GetRecoParam()->GetMultCutMinPointAngle()); + SetCutMaxDCADauther(AliITSReconstructor::GetRecoParam()->GetMultCutMaxDCADauther()); + SetCutMassGamma(AliITSReconstructor::GetRecoParam()->GetMultCutMassGamma()); + SetCutMassGammaNSigma(AliITSReconstructor::GetRecoParam()->GetMultCutMassGammaNSigma()); + SetCutMassK0(AliITSReconstructor::GetRecoParam()->GetMultCutMassK0()); + SetCutMassK0NSigma(AliITSReconstructor::GetRecoParam()->GetMultCutMassK0NSigma()); + SetCutChi2cGamma(AliITSReconstructor::GetRecoParam()->GetMultCutChi2cGamma()); + SetCutChi2cK0(AliITSReconstructor::GetRecoParam()->GetMultCutChi2cK0()); + SetCutGammaSFromDecay(AliITSReconstructor::GetRecoParam()->GetMultCutGammaSFromDecay()); + SetCutK0SFromDecay(AliITSReconstructor::GetRecoParam()->GetMultCutK0SFromDecay()); + SetCutMaxDCA(AliITSReconstructor::GetRecoParam()->GetMultCutMaxDCA()); + // } else { SetPhiWindow(); SetThetaWindow(); @@ -133,18 +211,34 @@ fhphiClustersLay1(0){ SetRemoveClustersFromOverlaps(); SetPhiOverlapCut(); SetZetaOverlapCut(); + SetPhiRotationAngle(); + + // + SetCutPxDrSPDin(); + SetCutPxDrSPDout(); + SetCutPxDz(); + SetCutDCArz(); + SetCutMinElectronProbTPC(); + SetCutMinElectronProbESD(); + SetCutMinP(); + SetCutMinRGamma(); + SetCutMinRK0(); + SetCutMinPointAngle(); + SetCutMaxDCADauther(); + SetCutMassGamma(); + SetCutMassGammaNSigma(); + SetCutMassK0(); + SetCutMassK0NSigma(); + SetCutChi2cGamma(); + SetCutChi2cK0(); + SetCutGammaSFromDecay(); + SetCutK0SFromDecay(); + SetCutMaxDCA(); } - - - fClustersLay1 = 0; - fClustersLay2 = 0; - fDetectorIndexClustersLay1 = 0; - fDetectorIndexClustersLay2 = 0; - fOverlapFlagClustersLay1 = 0; - fOverlapFlagClustersLay2 = 0; + // fTracklets = 0; fSClusters = 0; - + // // definition of histograms Bool_t oldStatus = TH1::AddDirectoryStatus(); TH1::AddDirectory(kFALSE); @@ -163,50 +257,85 @@ fhphiClustersLay1(0){ fhphiTracklets = new TH1F("phiTracklets", "phi", 100, 0., 2*TMath::Pi()); fhetaClustersLay1 = new TH1F("etaClustersLay1", "etaCl1", 100,-2.,2.); fhphiClustersLay1 = new TH1F("phiClustersLay1", "phiCl1", 100, 0., 2*TMath::Pi()); - + for (int i=2;i--;) fStoreRefs[i][0] = fStoreRefs[i][1] = kFALSE; TH1::AddDirectory(oldStatus); } //______________________________________________________________________ -AliITSMultReconstructor::AliITSMultReconstructor(const AliITSMultReconstructor &mr) : TObject(mr), -fClustersLay1(mr.fClustersLay1), -fClustersLay2(mr.fClustersLay2), -fDetectorIndexClustersLay1(mr.fDetectorIndexClustersLay1), -fDetectorIndexClustersLay2(mr.fDetectorIndexClustersLay2), -fOverlapFlagClustersLay1(mr.fOverlapFlagClustersLay1), -fOverlapFlagClustersLay2(mr.fOverlapFlagClustersLay2), -fTracklets(mr.fTracklets), -fSClusters(mr.fSClusters), -fNClustersLay1(mr.fNClustersLay1), -fNClustersLay2(mr.fNClustersLay2), -fNTracklets(mr.fNTracklets), -fNSingleCluster(mr.fNSingleCluster), -fPhiWindow(mr.fPhiWindow), -fThetaWindow(mr.fThetaWindow), -fPhiShift(mr.fPhiShift), -fRemoveClustersFromOverlaps(mr.fRemoveClustersFromOverlaps), -fPhiOverlapCut(mr.fPhiOverlapCut), -fZetaOverlapCut(mr.fZetaOverlapCut), -fHistOn(mr.fHistOn), -fhClustersDPhiAcc(mr.fhClustersDPhiAcc), -fhClustersDThetaAcc(mr.fhClustersDThetaAcc), -fhClustersDPhiAll(mr.fhClustersDPhiAll), -fhClustersDThetaAll(mr.fhClustersDThetaAll), -fhDPhiVsDThetaAll(mr.fhDPhiVsDThetaAll), -fhDPhiVsDThetaAcc(mr.fhDPhiVsDThetaAcc), -fhetaTracklets(mr.fhetaTracklets), -fhphiTracklets(mr.fhphiTracklets), -fhetaClustersLay1(mr.fhetaClustersLay1), -fhphiClustersLay1(mr.fhphiClustersLay1) { - // Copy constructor - +AliITSMultReconstructor::AliITSMultReconstructor(const AliITSMultReconstructor &mr) : +AliTrackleter(mr), +fDetTypeRec(0),fESDEvent(0),fTreeRP(0),fTreeRPMix(0), +fTracklets(0), +fSClusters(0), +fNTracklets(0), +fNSingleCluster(0), +fDPhiWindow(0), +fDThetaWindow(0), +fPhiShift(0), +fRemoveClustersFromOverlaps(0), +fPhiOverlapCut(0), +fZetaOverlapCut(0), +fPhiRotationAngle(0), +fScaleDTBySin2T(0), +fNStdDev(1.0), +fNStdDevSq(1.0), +// +fCutPxDrSPDin(0.1), +fCutPxDrSPDout(0.15), +fCutPxDz(0.2), +fCutDCArz(0.5), +fCutMinElectronProbTPC(0.5), +fCutMinElectronProbESD(0.1), +fCutMinP(0.05), +fCutMinRGamma(2.), +fCutMinRK0(1.), +fCutMinPointAngle(0.98), +fCutMaxDCADauther(0.5), +fCutMassGamma(0.03), +fCutMassGammaNSigma(5.), +fCutMassK0(0.03), +fCutMassK0NSigma(5.), +fCutChi2cGamma(2.), +fCutChi2cK0(2.), +fCutGammaSFromDecay(-10.), +fCutK0SFromDecay(-10.), +fCutMaxDCA(1.), +// +fHistOn(0), +fhClustersDPhiAcc(0), +fhClustersDThetaAcc(0), +fhClustersDPhiAll(0), +fhClustersDThetaAll(0), +fhDPhiVsDThetaAll(0), +fhDPhiVsDThetaAcc(0), +fhetaTracklets(0), +fhphiTracklets(0), +fhetaClustersLay1(0), +fhphiClustersLay1(0), +fDPhiShift(0), +fDPhiWindow2(0), +fDThetaWindow2(0), +fPartners(0), +fAssociatedLay1(0), +fMinDists(0), +fBlackList(0), +// +fCreateClustersCopy(0), +fClustersLoaded(0), +fRecoDone(0), +fSPDSeg() + { + // Copy constructor :!!! RS ATTENTION: old c-tor reassigned the pointers instead of creating a new copy -> would crash on delete + AliError("May not use"); } //______________________________________________________________________ AliITSMultReconstructor& AliITSMultReconstructor::operator=(const AliITSMultReconstructor& mr){ // Assignment operator - this->~AliITSMultReconstructor(); - new(this) AliITSMultReconstructor(mr); + if (this != &mr) { + this->~AliITSMultReconstructor(); + new(this) AliITSMultReconstructor(mr); + } return *this; } @@ -225,453 +354,356 @@ AliITSMultReconstructor::~AliITSMultReconstructor(){ delete fhphiTracklets; delete fhetaClustersLay1; delete fhphiClustersLay1; - - // delete arrays - for(Int_t i=0; i>>> RS: this part is equivalent to former AliITSVertexer::FindMultiplicity + // + // see if there is a SPD vertex + Bool_t isVtxOK=kTRUE, isCosmics=kFALSE; + AliESDVertex* vtx = (AliESDVertex*)fESDEvent->GetPrimaryVertexSPD(); + if (!vtx || vtx->GetNContributors()<1) isVtxOK = kFALSE; + if (vtx && strstr(vtx->GetTitle(),"cosmics")) { + isVtxOK = kFALSE; + isCosmics = kTRUE; + } + // + if (!isVtxOK) { + if (!isCosmics) { + AliDebug(1,"Tracklets multiplicity not determined because the primary vertex was not found"); + AliDebug(1,"Just counting the number of cluster-fired chips on the SPD layers"); + } + vtx = 0; + } + if(vtx){ + float vtxf[3] = {vtx->GetX(),vtx->GetY(),vtx->GetZ()}; + FindTracklets(vtxf); + } + else { + FindTracklets(0); + } + // + CreateMultiplicityObject(); } //____________________________________________________________________ void AliITSMultReconstructor::Reconstruct(TTree* clusterTree, Float_t* vtx, Float_t* /* vtxRes*/) { // - // - calls LoadClusterArray that finds the position of the clusters + // RS NOTE - this is old reconstructor invocation, to be used from VertexFinder and in analysis mode + + if (fMult) delete fMult; fMult = 0; + fNClustersLay[0] = 0; + fNClustersLay[1] = 0; + fNTracklets = 0; + fNSingleCluster = 0; + // + if (!clusterTree) { AliError(" Invalid ITS cluster tree !\n"); return; } + // + fESDEvent = 0; + SetTreeRP(clusterTree); + // + FindTracklets(vtx); + // +} + + +//____________________________________________________________________ +void AliITSMultReconstructor::ReconstructMix(TTree* clusterTree, TTree* clusterTreeMix, const Float_t* vtx, Float_t*) +{ + // + // RS NOTE - this is old reconstructor invocation, to be used from VertexFinder and in analysis mode + + if (fMult) delete fMult; fMult = 0; + fNClustersLay[0] = 0; + fNClustersLay[1] = 0; + fNTracklets = 0; + fNSingleCluster = 0; + // + if (!clusterTree) { AliError(" Invalid ITS cluster tree !\n"); return; } + if (!clusterTreeMix) { AliError(" Invalid ITS cluster tree 2nd event !\n"); return; } + // + fESDEvent = 0; + SetTreeRP(clusterTree); + SetTreeRPMix(clusterTreeMix); + // + FindTracklets(vtx); + // +} + + +//____________________________________________________________________ +void AliITSMultReconstructor::FindTracklets(const Float_t *vtx) +{ + // - calls LoadClusterArrays that finds the position of the clusters // (in global coord) + // - convert the cluster coordinates to theta, phi (seen from the - // interaction vertex). + // interaction vertex). Clusters in the inner layer can be now + // rotated for combinatorial studies // - makes an array of tracklets // // After this method has been called, the clusters of the two layers // and the tracklets can be retrieved by calling the Get'er methods. - // reset counters - fNClustersLay1 = 0; - fNClustersLay2 = 0; - fNTracklets = 0; - fNSingleCluster = 0; - // loading the clusters - LoadClusterArrays(clusterTree); + // Find tracklets converging to vertex + // + LoadClusterArrays(fTreeRP,fTreeRPMix); + // flag clusters used by ESD tracks + if (fESDEvent) ProcessESDTracks(); + fRecoDone = kTRUE; - const Double_t pi = TMath::Pi(); - - // dPhi shift is field dependent - // get average magnetic field - Float_t bz = 0; - AliMagF* field = 0; - if (TGeoGlobalMagField::Instance()) - field = dynamic_cast(TGeoGlobalMagField::Instance()->GetField()); - if (!field) - { - AliError("Could not retrieve magnetic field. Assuming no field. Delta Phi shift will be deactivated in AliITSMultReconstructor.") - } - else - bz = TMath::Abs(field->SolenoidField()); - - const Double_t dPhiShift = fPhiShift / 5 * bz; - AliDebug(1, Form("Using phi shift of %f", dPhiShift)); - - const Double_t dPhiWindow2 = fPhiWindow * fPhiWindow; - const Double_t dThetaWindow2 = fThetaWindow * fThetaWindow; - - Int_t* partners = new Int_t[fNClustersLay2]; - Float_t* minDists = new Float_t[fNClustersLay2]; - Int_t* associatedLay1 = new Int_t[fNClustersLay1]; - TArrayI** blacklist = new TArrayI*[fNClustersLay1]; - - for (Int_t i=0; iFill(eta); - fhphiClustersLay1->Fill(fClustersLay1[iC1][1]); - } - } - - // Loop on layer 2 : finding theta, phi and r - for (Int_t iC2=0; iC2 0) { found = 0; + for (Int_t iC1=0; iC1GetSize(); i++) { - if (blacklist[iC1]->At(i) == iC2) { - blacklisted = kTRUE; - break; - } - } - if (blacklisted) continue; - } - - // find the difference in angles - Double_t dTheta = TMath::Abs(fClustersLay2[iC2][0] - fClustersLay1[iC1][0]); - Double_t dPhi = TMath::Abs(fClustersLay2[iC2][1] - fClustersLay1[iC1][1]); - // take into account boundary condition - if (dPhi>pi) dPhi=2.*pi-dPhi; - - if (fHistOn) { - fhClustersDPhiAll->Fill(dPhi); - fhClustersDThetaAll->Fill(dTheta); - fhDPhiVsDThetaAll->Fill(dTheta, dPhi); - } - - dPhi -= dPhiShift; - - // make "elliptical" cut in Phi and Theta! - Float_t d = dPhi*dPhi/dPhiWindow2 + dTheta*dTheta/dThetaWindow2; - - // look for the minimum distance: the minimum is in iC2WithBestDist - if (d<1 && d minDist) { - Int_t oldPartner = partners[iC2WithBestDist]; - partners[iC2WithBestDist] = iC1; - minDists[iC2WithBestDist] = minDist; - - // mark as assigned - associatedLay1[iC1] = 1; - - if (oldPartner != -1) { - // redo partner search for cluster in L0 (oldPartner), putting this one (iC2WithBestDist) on its blacklist - if (blacklist[oldPartner] == 0) { - blacklist[oldPartner] = new TArrayI(1); - } else blacklist[oldPartner]->Set(blacklist[oldPartner]->GetSize()+1); - - blacklist[oldPartner]->AddAt(iC2WithBestDist, blacklist[oldPartner]->GetSize()-1); - - // mark as free - associatedLay1[oldPartner] = 0; - } - } else { - // try again to find a cluster without considering iC2WithBestDist - if (blacklist[iC1] == 0) { - blacklist[iC1] = new TArrayI(1); - } - else - blacklist[iC1]->Set(blacklist[iC1]->GetSize()+1); - - blacklist[iC1]->AddAt(iC2WithBestDist, blacklist[iC1]->GetSize()-1); - } - - } else // cluster has no partner; remove - associatedLay1[iC1] = 2; - } // end of loop over clusters in layer 1 - } - - // Step2: store tracklets; remove used clusters - for (Int_t iC2=0; iC2 TMath::Pi()) - fTracklets[fNTracklets][2] = fTracklets[fNTracklets][2]-2.*TMath::Pi(); - if (fTracklets[fNTracklets][2] < -TMath::Pi()) - fTracklets[fNTracklets][2] = fTracklets[fNTracklets][2]+2.*TMath::Pi(); - - // store the difference between theta1 and theta2 - fTracklets[fNTracklets][3] = fClustersLay1[partners[iC2]][0] - fClustersLay2[iC2][0]; - - if (fHistOn) { - fhClustersDPhiAcc->Fill(fTracklets[fNTracklets][2]); - fhClustersDThetaAcc->Fill(fTracklets[fNTracklets][3]); - fhDPhiVsDThetaAcc->Fill(fTracklets[fNTracklets][3],fTracklets[fNTracklets][2]); - } - - // find label - // if equal label in both clusters found this label is assigned - // if no equal label can be found the first labels of the L1 AND L2 cluster are assigned - Int_t label1 = 0; - Int_t label2 = 0; - while (label2 < 3) { - if ((Int_t) fClustersLay1[partners[iC2]][3+label1] != -2 && (Int_t) fClustersLay1[partners[iC2]][3+label1] == (Int_t) fClustersLay2[iC2][3+label2]) - break; - label1++; - if (label1 == 3) { - label1 = 0; - label2++; +//____________________________________________________________________ +void AliITSMultReconstructor::CreateMultiplicityObject() +{ + // create AliMultiplicity object and store it in the ESD event + // + TBits fastOrFiredMap,firedChipMap; + if (fDetTypeRec) { + fastOrFiredMap = fDetTypeRec->GetFastOrFiredMap(); + firedChipMap = fDetTypeRec->GetFiredChipMap(fTreeRP); + } + // + fMult = new AliMultiplicity(fNTracklets,fNSingleCluster,fNFiredChips[0],fNFiredChips[1],fastOrFiredMap); + fMult->SetMultTrackRefs(kTRUE); + // store some details of reco: + fMult->SetScaleDThetaBySin2T(fScaleDTBySin2T); + fMult->SetDPhiWindow2(fDPhiWindow2); + fMult->SetDThetaWindow2(fDThetaWindow2); + fMult->SetDPhiShift(fDPhiShift); + fMult->SetNStdDev(fNStdDev); + // + fMult->SetFiredChipMap(firedChipMap); + AliITSRecPointContainer* rcont = AliITSRecPointContainer::Instance(); + fMult->SetITSClusters(0,rcont->GetNClustersInLayer(1,fTreeRP)); + for(Int_t kk=2;kk<=6;kk++) fMult->SetITSClusters(kk-1,rcont->GetNClustersInLayerFast(kk)); + // + UInt_t shared[100]; + AliRefArray *refs[2][2] = {{0,0},{0,0}}; + for (int il=2;il--;) + for (int it=2;it--;) // tracklet_clusters->track references to stor + if (fStoreRefs[il][it]) refs[il][it] = new AliRefArray(fNTracklets,0); + // + for (int i=fNTracklets;i--;) { + float* tlInfo = fTracklets[i]; + fMult->SetTrackletData(i,tlInfo); + for (int itp=0;itp<2;itp++) { + for (int ilr=0;ilr<2;ilr++) { + if (!fStoreRefs[ilr][itp]) continue; // nothing to store + int clID = int(tlInfo[ilr ? kClID2:kClID1]); + int nref = fUsedClusLay[ilr][itp]->GetReferences(clID,shared,100); + if (!nref) continue; + else if (nref==1) refs[ilr][itp]->AddReference(i,shared[0]); + else refs[ilr][itp]->AddReferences(i,shared,nref); } } - if (label2 < 3) { - AliDebug(AliLog::kDebug, Form("Found label %d == %d for tracklet candidate %d\n", (Int_t) fClustersLay1[partners[iC2]][3+label1], (Int_t) fClustersLay2[iC2][3+label2], fNTracklets)); - fTracklets[fNTracklets][4] = fClustersLay1[partners[iC2]][3+label1]; - fTracklets[fNTracklets][5] = fClustersLay2[iC2][3+label2]; - } else { - AliDebug(AliLog::kDebug, Form("Did not find label %d %d %d %d %d %d for tracklet candidate %d\n", (Int_t) fClustersLay1[partners[iC2]][3], (Int_t) fClustersLay1[partners[iC2]][4], (Int_t) fClustersLay1[partners[iC2]][5], (Int_t) fClustersLay2[iC2][3], (Int_t) fClustersLay2[iC2][4], (Int_t) fClustersLay2[iC2][5], fNTracklets)); - fTracklets[fNTracklets][4] = fClustersLay1[partners[iC2]][3]; - fTracklets[fNTracklets][5] = fClustersLay2[iC2][3]; - } - - if (fHistOn) { - Float_t eta=fTracklets[fNTracklets][0]; - eta= TMath::Tan(eta/2.); - eta=-TMath::Log(eta); - fhetaTracklets->Fill(eta); - fhphiTracklets->Fill(fTracklets[fNTracklets][1]); - } - - AliDebug(1,Form(" Adding tracklet candidate %d ", fNTracklets)); - AliDebug(1,Form(" Cl. %d of Layer 1 and %d of Layer 2", partners[iC2], iC2)); - fNTracklets++; - - associatedLay1[partners[iC2]] = 1; } - - // Delete the following else if you do not want to save Clusters! - // store the cluster - for (Int_t iC1=0; iC1AttachTracklet2TrackRefs(refs[0][0],refs[0][1],refs[1][0],refs[1][1]); + // + AliRefArray *refsc[2] = {0,0}; + for (int it=2;it--;) if (fStoreRefs[0][it]) refsc[it] = new AliRefArray(fNClustersLay[0]); + for (int i=fNSingleCluster;i--;) { + float* clInfo = fSClusters[i]; + fMult->SetSingleClusterData(i,clInfo); + int clID = int(clInfo[kSCID]); + for (int itp=0;itp<2;itp++) { + if (!fStoreRefs[0][itp]) continue; + int nref = fUsedClusLay[0][itp]->GetReferences(clID,shared,100); + if (!nref) continue; + else if (nref==1) refsc[itp]->AddReference(i,shared[0]); + else refsc[itp]->AddReferences(i,shared,nref); } } + fMult->AttachCluster2TrackRefs(refsc[0],refsc[1]); + fMult->CompactBits(); + // +} - delete[] partners; - delete[] minDists; - - for (Int_t i=0; iFetchClusters(0,itsClusterTree); - if(!rpcont->IsSPDActive()){ - AliWarning("No SPD rec points found, multiplicity not calculated"); - return; - } - Float_t cluGlo[3]={0.,0.,0.}; - - + // + // RS: This method was strongly modified wrt original. In order to have the same numbering + // of clusters as in the ITS reco I had to introduce sorting in Z + // Also note that now the clusters data are stored not in float[6] attached to float**, but in 1-D array + AliDebug(1,Form("Loading clusters and cluster-fired chips for layer %d",il)); + // + fNClustersLay[il] = 0; + fNFiredChips[il] = 0; + for (int i=2;i--;) fStoreRefs[il][i] = kFALSE; + // + AliITSRecPointContainer* rpcont = 0; + static TClonesArray statITSrec("AliITSRecPoint"); + static TObjArray clArr(100); + TBranch* branch = 0; + TClonesArray* itsClusters = 0; + // + if (!fCreateClustersCopy) { + rpcont=AliITSRecPointContainer::Instance(); + itsClusters = rpcont->FetchClusters(0,itsClusterTree); + if(!rpcont->IsSPDActive()){ + AliWarning("No SPD rec points found, multiplicity not calculated"); + return; + } + } + else { + itsClusters = &statITSrec; + branch = itsClusterTree->GetBranch("ITSRecPoints"); + branch->SetAddress(&itsClusters); + if (!fClArr[il]) fClArr[il] = new TClonesArray("AliITSRecPoint",100); + } + // // count clusters // loop over the SPD subdetectors - Int_t nSPDL1 = AliITSgeomTGeo::GetModuleIndex(2,1,1); - for (Int_t iIts=0; iIts < nSPDL1; iIts++) { - itsClusters=rpcont->UncheckedGetClusters(iIts); - fNClustersLay1 += itsClusters->GetEntriesFast(); - } - Int_t nSPDL2=AliITSgeomTGeo::GetModuleIndex(3,1,1); - for (Int_t iIts=nSPDL1; iIts < nSPDL2; iIts++) { - itsClusters=rpcont->UncheckedGetClusters(iIts); - fNClustersLay2 += itsClusters->GetEntriesFast(); - } - - // create arrays - fClustersLay1 = new Float_t*[fNClustersLay1]; - fDetectorIndexClustersLay1 = new Int_t[fNClustersLay1]; - fOverlapFlagClustersLay1 = new Bool_t[fNClustersLay1]; - - fClustersLay2 = new Float_t*[fNClustersLay2]; - fDetectorIndexClustersLay2 = new Int_t[fNClustersLay2]; - fOverlapFlagClustersLay2 = new Bool_t[fNClustersLay2]; - - // no double association allowed - fTracklets = new Float_t*[TMath::Min(fNClustersLay1, fNClustersLay2)]; - fSClusters = new Float_t*[fNClustersLay1]; - - for (Int_t i=0; iUncheckedGetClusters(iIts); - - Int_t nClusters = itsClusters->GetEntriesFast(); - - // number of clusters in each chip of the current module + int nclLayer = 0; + int detMin = TMath::Max(0,AliITSgeomTGeo::GetModuleIndex(il+1,1,1)); + int detMax = AliITSgeomTGeo::GetModuleIndex(il+2,1,1); + for (int idt=detMin;idtUncheckedGetClusters(idt); + else branch->GetEvent(idt); + int nClusters = itsClusters->GetEntriesFast(); + if (!nClusters) continue; Int_t nClustersInChip[5] = {0,0,0,0,0}; - Int_t layer = 0; - - // loop over clusters while(nClusters--) { AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters); - - layer = cluster->GetLayer(); - if (layer>1) continue; - - cluster->GetGlobalXYZ(cluGlo); - Float_t x = cluGlo[0]; - Float_t y = cluGlo[1]; - Float_t z = cluGlo[2]; - - // find the chip for the current cluster - Float_t locz = cluster->GetDetLocalZ(); - Int_t iChip = seg.GetChipFromLocal(0,locz); - nClustersInChip[iChip]++; - - if (layer==0) { - fClustersLay1[fNClustersLay1][0] = x; - fClustersLay1[fNClustersLay1][1] = y; - fClustersLay1[fNClustersLay1][2] = z; - - fDetectorIndexClustersLay1[fNClustersLay1]=iIts; - - for (Int_t i=0; i<3; i++) - fClustersLay1[fNClustersLay1][3+i] = cluster->GetLabel(i); - fNClustersLay1++; - } - if (layer==1) { - fClustersLay2[fNClustersLay2][0] = x; - fClustersLay2[fNClustersLay2][1] = y; - fClustersLay2[fNClustersLay2][2] = z; - - fDetectorIndexClustersLay2[fNClustersLay2]=iIts; - - for (Int_t i=0; i<3; i++) - fClustersLay2[fNClustersLay2][3+i] = cluster->GetLabel(i); - fNClustersLay2++; - } - - }// end of cluster loop - - // get number of fired chips in the current module - - for(Int_t ifChip=0; ifChip<5; ifChip++) { - if(nClustersInChip[ifChip] >= 1) fNFiredChips[layer]++; + if (!cluster) continue; + if (fCreateClustersCopy) cluster = new ((*fClArr[il])[nclLayer]) AliITSRecPoint(*cluster); + clArr.AddAtAndExpand(cluster,nclLayer++); + Int_t chipNo = fSPDSeg.GetChipFromLocal(0,cluster->GetDetLocalZ()); + if(chipNo>=0)nClustersInChip[ chipNo ]++; } - - } // end of its "subdetector" loop - - AliDebug(1,Form("(clusters in layer 1 : %d, layer 2: %d)",fNClustersLay1,fNClustersLay2)); - AliDebug(1,Form("(cluster-fired chips in layer 1 : %d, layer 2: %d)",fNFiredChips[0],fNFiredChips[1])); + for(Int_t ifChip=5;ifChip--;) if (nClustersInChip[ifChip]) fNFiredChips[il]++; + } + // sort the clusters in Z (to have the same numbering as in ITS reco + Float_t *z = new Float_t[nclLayer]; + Int_t *index = new Int_t[nclLayer]; + for (int ic=0;icGetZ(); + TMath::Sort(nclLayer,z,index,kFALSE); + Float_t* clustersLay = new Float_t[nclLayer*kClNPar]; + Int_t* detectorIndexClustersLay = new Int_t[nclLayer]; + Bool_t* overlapFlagClustersLay = new Bool_t[nclLayer]; + // + for (int ic=0;icGetGlobalXYZ( clPar ); + detectorIndexClustersLay[ic] = cluster->GetDetectorIndex(); + overlapFlagClustersLay[ic] = kFALSE; + for (Int_t i=3;i--;) clPar[kClMC0+i] = cluster->GetLabel(i); + } + clArr.Clear(); + delete[] z; + delete[] index; + // + if (fOverlapFlagClustersLay[il]) delete[] fOverlapFlagClustersLay[il]; + fOverlapFlagClustersLay[il] = overlapFlagClustersLay; + // + if (fDetectorIndexClustersLay[il]) delete[] fDetectorIndexClustersLay[il]; + fDetectorIndexClustersLay[il] = detectorIndexClustersLay; + // + for (int it=0;it<2;it++) { + if (fUsedClusLay[il][it]) delete fUsedClusLay[il][it]; + fUsedClusLay[il][it] = new AliRefArray(nclLayer); + } + // + if (fClustersLay[il]) delete[] fClustersLay[il]; + fClustersLay[il] = clustersLay; + fNClustersLay[il] = nclLayer; + // } + //____________________________________________________________________ -void -AliITSMultReconstructor::LoadClusterFiredChips(TTree* itsClusterTree) { - // This method +void AliITSMultReconstructor::LoadClusterFiredChips(TTree* itsClusterTree) { + // This method // - gets the clusters from the cluster tree // - counts the number of (cluster)fired chips @@ -680,9 +712,9 @@ AliITSMultReconstructor::LoadClusterFiredChips(TTree* itsClusterTree) { fNFiredChips[0] = 0; fNFiredChips[1] = 0; - AliITSsegmentationSPD seg; AliITSRecPointContainer* rpcont=AliITSRecPointContainer::Instance(); - TClonesArray* itsClusters=rpcont->FetchClusters(0,itsClusterTree); + TClonesArray* itsClusters=NULL; + rpcont->FetchClusters(0,itsClusterTree); if(!rpcont->IsSPDActive()){ AliWarning("No SPD rec points found, multiplicity not calculated"); return; @@ -697,18 +729,20 @@ AliITSMultReconstructor::LoadClusterFiredChips(TTree* itsClusterTree) { // number of clusters in each chip of the current module Int_t nClustersInChip[5] = {0,0,0,0,0}; Int_t layer = 0; + Int_t ladder=0; + Int_t det=0; + AliITSgeomTGeo::GetModuleId(iIts,layer,ladder,det); + --layer; // layer is from 1 to 6 in AliITSgeomTGeo, but from 0 to 5 here + if(layer<0 || layer >1)continue; // loop over clusters while(nClusters--) { AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters); - - layer = cluster->GetLayer(); - if (layer>1) continue; - + // find the chip for the current cluster Float_t locz = cluster->GetDetLocalZ(); - Int_t iChip = seg.GetChipFromLocal(0,locz); - nClustersInChip[iChip]++; + Int_t iChip = fSPDSeg.GetChipFromLocal(0,locz); + if (iChip>=0) nClustersInChip[iChip]++; }// end of cluster loop @@ -746,9 +780,9 @@ AliITSMultReconstructor::SaveHists() { } //____________________________________________________________________ -void -AliITSMultReconstructor::FlagClustersInOverlapRegions (Int_t iC1, Int_t iC2WithBestDist) { - +void AliITSMultReconstructor::FlagClustersInOverlapRegions (Int_t iC1, Int_t iC2WithBestDist) +{ + // Flags clusters in the overlapping regions Float_t distClSameMod=0.; Float_t distClSameModMin=0.; Int_t iClOverlap =0; @@ -758,23 +792,25 @@ AliITSMultReconstructor::FlagClustersInOverlapRegions (Int_t iC1, Int_t iC2WithB Float_t zproj1=0.; Float_t zproj2=0.; Float_t deZproj=0.; - + Float_t* clPar1 = GetClusterLayer1(iC1); + Float_t* clPar2B = GetClusterLayer2(iC2WithBestDist); // Loop on inner layer clusters - for (Int_t iiC1=0; iiC1TMath::Pi()) dePhi=2.*TMath::Pi()-dePhi; - zproj1=meanRadiusLay1/TMath::Tan(fClustersLay1[iC1][0]); - zproj2=meanRadiusLay1/TMath::Tan(fClustersLay1[iiC1][0]); + zproj1=meanRadiusLay1/TMath::Tan(clPar1[kClTh]); + zproj2=meanRadiusLay1/TMath::Tan(clPar11[kClTh]); deZproj=TMath::Abs(zproj1-zproj2); distClSameMod = TMath::Sqrt(TMath::Power(deZproj/fZetaOverlapCut,2)+TMath::Power(dePhi/fPhiOverlapCut,2)); - if (distClSameMod<=1.) fOverlapFlagClustersLay1[iiC1]=kTRUE; + if (distClSameMod<=1.) fOverlapFlagClustersLay[0][iiC1]=kTRUE; // if (distClSameMod<=1.) { // if (distClSameModMin==0. || distClSameModTMath::Pi()) dePhi=2.*TMath::Pi()-dePhi; - zproj1=meanRadiusLay2/TMath::Tan(fClustersLay2[iC2WithBestDist][0]); - zproj2=meanRadiusLay2/TMath::Tan(fClustersLay2[iiC2][0]); + zproj1=meanRadiusLay2/TMath::Tan(clPar2B[kClTh]); + zproj2=meanRadiusLay2/TMath::Tan(clPar2[kClTh]); deZproj=TMath::Abs(zproj1-zproj2); distClSameMod = TMath::Sqrt(TMath::Power(deZproj/fZetaOverlapCut,2)+TMath::Power(dePhi/fPhiOverlapCut,2)); - if (distClSameMod<=1.) fOverlapFlagClustersLay2[iiC2]=kTRUE; + if (distClSameMod<=1.) fOverlapFlagClustersLay[1][iiC2]=kTRUE; // if (distClSameMod<=1.) { // if (distClSameModMin==0. || distClSameMod(TGeoGlobalMagField::Instance()->GetField()); + if (!field) { + AliError("Could not retrieve magnetic field. Assuming no field. Delta Phi shift will be deactivated in AliITSMultReconstructor."); + } + else bz = TMath::Abs(field->SolenoidField()); + fDPhiShift = fPhiShift / 5 * bz; + AliDebug(1, Form("Using phi shift of %f", fDPhiShift)); + // + if (fPartners) delete[] fPartners; fPartners = new Int_t[fNClustersLay[1]]; + if (fMinDists) delete[] fMinDists; fMinDists = new Float_t[fNClustersLay[1]]; + if (fAssociatedLay1) delete[] fAssociatedLay1; fAssociatedLay1 = new Int_t[fNClustersLay[0]]; + // + if (fBlackList) delete fBlackList; fBlackList = new AliRefArray(fNClustersLay[0]); + // + // Printf("Vertex in find tracklets...%f %f %f",vtx[0],vtx[1],vtx[2]); + for (Int_t i=0; iFill(eta); + fhphiClustersLay1->Fill(clPar[kClPh]); + } + } + } + } + // +} + +//____________________________________________________________________ +Int_t AliITSMultReconstructor::AssociateClusterOfL1(Int_t iC1) +{ + // search association of cluster iC1 of L1 with all clusters of L2 + if (fAssociatedLay1[iC1] != 0) return 0; + Int_t iC2WithBestDist = -1; // reset + Double_t minDist = 2*fNStdDev; // reset + float* clPar1 = GetClusterLayer1(iC1); + for (Int_t iC2=0; iC2IsReferred(iC1,iC2)) continue; + float* clPar2 = GetClusterLayer2(iC2); + // + // find the difference in angles + Double_t dTheta = TMath::Abs(clPar2[kClTh] - clPar1[kClTh]); + Double_t dPhi = TMath::Abs(clPar2[kClPh] - clPar1[kClPh]); + // Printf("detheta %f dephi %f", dTheta,dPhi); + // + if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi; // take into account boundary condition + // + if (fHistOn) { + fhClustersDPhiAll->Fill(dPhi); + fhClustersDThetaAll->Fill(dTheta); + fhDPhiVsDThetaAll->Fill(dTheta, dPhi); + } + Float_t d = CalcDist(dPhi,dTheta,clPar1[kClTh]); // make "elliptical" cut in Phi and Theta! + // look for the minimum distance: the minimum is in iC2WithBestDist + if (d minDist) { + Int_t oldPartner = fPartners[iC2WithBestDist]; + fPartners[iC2WithBestDist] = iC1; + fMinDists[iC2WithBestDist] = minDist; + // + fAssociatedLay1[iC1] = 1; // mark as assigned + // + if (oldPartner != -1) { + // redo partner search for cluster in L0 (oldPartner), putting this one (iC2WithBestDist) on its fBlackList + fBlackList->AddReference(oldPartner,iC2WithBestDist); + fAssociatedLay1[oldPartner] = 0; // mark as free + } + } else { + // try again to find a cluster without considering iC2WithBestDist + fBlackList->AddReference(iC1,iC2WithBestDist); + } + // + } + else fAssociatedLay1[iC1] = 2;// cluster has no partner; remove + // + return 1; +} + +//____________________________________________________________________ +Int_t AliITSMultReconstructor::StoreTrackletForL2Cluster(Int_t iC2) +{ + // build tracklet for cluster iC2 of layer 2 + if (fPartners[iC2] == -1) return 0; + if (fRemoveClustersFromOverlaps) FlagClustersInOverlapRegions (fPartners[iC2],iC2); + // Printf("saving tracklets"); + if (fOverlapFlagClustersLay[0][fPartners[iC2]] || fOverlapFlagClustersLay[1][iC2]) return 0; + float* clPar2 = GetClusterLayer2(iC2); + float* clPar1 = GetClusterLayer1(fPartners[iC2]); + // + Float_t* tracklet = fTracklets[fNTracklets] = new Float_t[kTrNPar]; // RS Add also the cluster id's + // + tracklet[kTrTheta] = clPar1[kClTh]; // use the theta from the clusters in the first layer + tracklet[kTrPhi] = clPar1[kClPh]; // use the phi from the clusters in the first layer + tracklet[kTrDPhi] = clPar1[kClPh] - clPar2[kClPh]; // store the difference between phi1 and phi2 + // + // define dphi in the range [0,pi] with proper sign (track charge correlated) + if (tracklet[kTrDPhi] > TMath::Pi()) tracklet[kTrDPhi] = tracklet[kTrDPhi]-2.*TMath::Pi(); + if (tracklet[kTrDPhi] < -TMath::Pi()) tracklet[kTrDPhi] = tracklet[kTrDPhi]+2.*TMath::Pi(); + // + tracklet[kTrDTheta] = clPar1[kClTh] - clPar2[kClTh]; // store the theta1-theta2 + // + if (fHistOn) { + fhClustersDPhiAcc->Fill(tracklet[kTrDPhi]); + fhClustersDThetaAcc->Fill(tracklet[kTrDTheta]); + fhDPhiVsDThetaAcc->Fill(tracklet[kTrDTheta],tracklet[kTrDPhi]); + } + // + // find label + // if equal label in both clusters found this label is assigned + // if no equal label can be found the first labels of the L1 AND L2 cluster are assigned + Int_t label1=0,label2=0; + while (label2 < 3) { + if ( int(clPar1[kClMC0+label1])!=-2 && int(clPar1[kClMC0+label1])==int(clPar2[kClMC0+label2])) break; + if (++label1 == 3) { label1 = 0; label2++; } + } + if (label2 < 3) { + AliDebug(AliLog::kDebug, Form("Found label %d == %d for tracklet candidate %d\n", + (Int_t) clPar1[kClMC0+label1], (Int_t) clPar1[kClMC0+label2], fNTracklets)); + tracklet[kTrLab1] = tracklet[kTrLab2] = clPar1[kClMC0+label1]; + } else { + AliDebug(AliLog::kDebug, Form("Did not find label %d %d %d %d %d %d for tracklet candidate %d\n", + (Int_t) clPar1[kClMC0], (Int_t) clPar1[kClMC1], (Int_t) clPar1[kClMC2], + (Int_t) clPar2[kClMC0], (Int_t) clPar2[kClMC1], (Int_t) clPar2[kClMC2], fNTracklets)); + tracklet[kTrLab1] = clPar1[kClMC0]; + tracklet[kTrLab2] = clPar2[kClMC0]; + } + // + if (fHistOn) { + Float_t eta = tracklet[kTrTheta]; + eta= TMath::Tan(eta/2.); + eta=-TMath::Log(eta); + fhetaTracklets->Fill(eta); + fhphiTracklets->Fill(tracklet[kTrPhi]); + } + // + tracklet[kClID1] = fPartners[iC2]; + tracklet[kClID2] = iC2; + // + // Printf("Adding tracklet candidate"); + AliDebug(1,Form(" Adding tracklet candidate %d ", fNTracklets)); + AliDebug(1,Form(" Cl. %d of Layer 1 and %d of Layer 2", fPartners[iC2], iC2)); + fNTracklets++; + fAssociatedLay1[fPartners[iC2]] = 1; + // + return 1; +} + +//____________________________________________________________________ +void AliITSMultReconstructor::StoreL1Singles() +{ + // Printf("saving single clusters..."); + for (Int_t iC1=0; iC1GetPrimaryVertexTracks(); + if (!vtx || vtx->GetNContributors()<1) vtx = (AliESDVertex*)fESDEvent->GetPrimaryVertexSPD(); + if (!vtx || vtx->GetNContributors()<1) { + AliDebug(1,"No primary vertex: cannot flag primary tracks"); + return; + } + Int_t ntracks = fESDEvent->GetNumberOfTracks(); + for(Int_t itr=0; itrGetTrack(itr); + if (!track->IsOn(AliESDtrack::kITSin)) continue; // use only tracks propagated in ITS to vtx + FlagTrackClusters(itr); + FlagIfSecondary(track,vtx); + } + FlagV0s(vtx); + // +} + +//____________________________________________________________________ +void AliITSMultReconstructor::FlagTrackClusters(Int_t id) +{ + // RS: flag the SPD clusters of the track if it is useful for the multiplicity estimation + // + const AliESDtrack* track = fESDEvent->GetTrack(id); + Int_t idx[12]; + if ( track->GetITSclusters(idx)<3 ) return; // at least 3 clusters must be used in the fit + Int_t itsType = track->IsOn(AliESDtrack::kITSpureSA) ? 1:0; + + for (int i=6/*AliESDfriendTrack::kMaxITScluster*/;i--;) { // ignore extras: note: i>=6 is for extra clusters + if (idx[i]<0) continue; + int layID= (idx[i] & 0xf0000000) >> 28; + if (layID>1) continue; // SPD only + int clID = (idx[i] & 0x0fffffff); + fUsedClusLay[layID][itsType]->AddReference(clID,id); + fStoreRefs[layID][itsType] = kTRUE; + } + // +} + +//____________________________________________________________________ +void AliITSMultReconstructor::FlagIfSecondary(AliESDtrack* track, const AliVertex* vtx) +{ + // RS: check if the track is primary and set the flag + double cut = (track->HasPointOnITSLayer(0)||track->HasPointOnITSLayer(1)) ? fCutPxDrSPDin:fCutPxDrSPDout; + float xz[2]; + track->GetDZ(vtx->GetX(),vtx->GetY(),vtx->GetZ(), fESDEvent->GetMagneticField(), xz); + if (TMath::Abs(xz[0]*track->P())>cut || TMath::Abs(xz[1]*track->P())>fCutPxDz || + TMath::Abs(xz[0])>fCutDCArz || TMath::Abs(xz[1])>fCutDCArz) + track->SetStatus(AliESDtrack::kMultSec); + else track->ResetStatus(AliESDtrack::kMultSec); +} + +//____________________________________________________________________ +void AliITSMultReconstructor::FlagV0s(const AliESDVertex *vtx) +{ + // flag tracks belonging to v0s + // + const double kK0Mass = 0.4976; + // + AliV0 pvertex; + AliKFVertex vertexKF; + AliKFParticle epKF0,epKF1,pipmKF0,piKF0,piKF1,gammaKF,k0KF; + Double_t mass,massErr,chi2c; + enum {kKFIni=BIT(14)}; + // + double recVtx[3]; + float recVtxF[3]; + vtx->GetXYZ(recVtx); + for (int i=3;i--;) recVtxF[i] = recVtx[i]; + // + int ntracks = fESDEvent->GetNumberOfTracks(); + if (ntracks<2) return; + // + vertexKF.X() = recVtx[0]; + vertexKF.Y() = recVtx[1]; + vertexKF.Z() = recVtx[2]; + vertexKF.Covariance(0,0) = vtx->GetXRes()*vtx->GetXRes(); + vertexKF.Covariance(1,2) = vtx->GetYRes()*vtx->GetYRes(); + vertexKF.Covariance(2,2) = vtx->GetZRes()*vtx->GetZRes(); + // + AliESDtrack *trc0,*trc1; + for (int it0=0;it0GetTrack(it0); + if (trc0->IsOn(AliESDtrack::kMultInV0)) continue; + if (!trc0->IsOn(AliESDtrack::kITSin)) continue; + Bool_t isSAP = trc0->IsPureITSStandalone(); + Int_t q0 = trc0->Charge(); + Bool_t testGamma = CanBeElectron(trc0); + epKF0.ResetBit(kKFIni); + piKF0.ResetBit(kKFIni); + double bestChi2=1e16; + int bestID = -1; + // + for (int it1=it0+1;it1GetTrack(it1); + if (trc1->IsOn(AliESDtrack::kMultInV0)) continue; + if (!trc1->IsOn(AliESDtrack::kITSin)) continue; + if (trc1->IsPureITSStandalone() != isSAP) continue; // pair separately ITS_SA_Pure tracks and TPC/ITS+ITS_SA + if ( (q0+trc1->Charge())!=0 ) continue; // don't pair like signs + // + pvertex.SetParamN(q0<0 ? *trc0:*trc1); + pvertex.SetParamP(q0>0 ? *trc0:*trc1); + pvertex.Update(recVtxF); + if (pvertex.P()fCutMaxDCADauther) continue; + double d = pvertex.GetD(recVtx[0],recVtx[1],recVtx[2]); + if (d>fCutMaxDCA) continue; + double dx=recVtx[0]-pvertex.Xv(), dy=recVtx[1]-pvertex.Yv(); + double rv = TMath::Sqrt(dx*dx+dy*dy); + // + // check gamma conversion hypothesis ----------------------------------------------------------->>> + Bool_t gammaOK = kFALSE; + while (testGamma && CanBeElectron(trc1)) { + if (rv0 ? kPositron:kElectron); + epKF0.SetBit(kKFIni); + } + new(&epKF1) AliKFParticle(*trc1,q0<0 ? kPositron:kElectron); + gammaKF.Initialize(); + gammaKF += epKF0; + gammaKF += epKF1; + gammaKF.SetProductionVertex(vertexKF); + gammaKF.GetMass(mass,massErr); + if (mass>fCutMassGamma || (massErr>0&&(mass>massErr*fCutMassGammaNSigma))) break; + if (gammaKF.GetS()fCutChi2cGamma) break; + gammaOK = kTRUE; + if (chi2c>bestChi2) break; + bestChi2 = chi2c; + bestID = it1; + break; + } + if (gammaOK) continue; + // check gamma conversion hypothesis -----------------------------------------------------------<<< + // check K0 conversion hypothesis ----------------------------------------------------------->>> + while (1) { + if (rv0 ? kPiPlus:kPiMinus); + piKF0.SetBit(kKFIni); + } + new(&piKF1) AliKFParticle(*trc1,q0<0 ? kPiPlus:kPiMinus); + k0KF.Initialize(); + k0KF += piKF0; + k0KF += piKF1; + k0KF.SetProductionVertex(vertexKF); + k0KF.GetMass(mass,massErr); + mass -= kK0Mass; + if (TMath::Abs(mass)>fCutMassK0 || (massErr>0&&(abs(mass)>massErr*fCutMassK0NSigma))) break; + if (k0KF.GetS()fCutChi2cK0) break; + if (chi2c>bestChi2) break; + bestChi2 = chi2c; + bestID = it1; + break; + } + // check K0 conversion hypothesis -----------------------------------------------------------<<< + } + // + if (bestID>=0) { + trc0->SetStatus(AliESDtrack::kMultInV0); + fESDEvent->GetTrack(bestID)->SetStatus(AliESDtrack::kMultInV0); + } + } + // +} + +//____________________________________________________________________ +Bool_t AliITSMultReconstructor::CanBeElectron(const AliESDtrack* trc) const +{ + // check if the track can be electron + Double_t pid[AliPID::kSPECIES]; + if (!trc->IsOn(AliESDtrack::kESDpid)) return kTRUE; + trc->GetESDpid(pid); + return (trc->IsOn(AliESDtrack::kTPCpid)) ? + pid[AliPID::kElectron]>fCutMinElectronProbTPC : + pid[AliPID::kElectron]>fCutMinElectronProbESD; + // +}