X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=ITS%2FAliITSMultReconstructor.cxx;h=17ebe5d57b640844557a1a91108e2cd4aa7bacbb;hp=53e4338fd2506951354636eb1e40348053e1b050;hb=fcb289102b9b5323c24bf4ec175fb97fcc0260b8;hpb=1f9831abf71114f584f67d9d98eb3086fdf712a6 diff --git a/ITS/AliITSMultReconstructor.cxx b/ITS/AliITSMultReconstructor.cxx index 53e4338fd25..17ebe5d57b6 100644 --- a/ITS/AliITSMultReconstructor.cxx +++ b/ITS/AliITSMultReconstructor.cxx @@ -58,6 +58,7 @@ // - 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 @@ -66,10 +67,10 @@ #include #include #include +#include #include "AliITSMultReconstructor.h" #include "AliITSReconstructor.h" -#include "AliITSsegmentationSPD.h" #include "AliITSRecPoint.h" #include "AliITSRecPointContainer.h" #include "AliITSgeom.h" @@ -82,6 +83,11 @@ #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) @@ -89,25 +95,44 @@ ClassImp(AliITSMultReconstructor) //____________________________________________________________________ AliITSMultReconstructor::AliITSMultReconstructor(): -fDetTypeRec(0),fESDEvent(0),fTreeRP(0),fUsedClusLay1(0),fUsedClusLay2(0), -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), +fNSingleClusterSPD2(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), @@ -118,22 +143,73 @@ 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), + fBuildRefs(kTRUE), + fStoreSPD2SingleCl(kFALSE), + 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; + fClusterCopyIndex[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()) { + if (AliITSReconstructor::GetRecoParam()) { SetPhiWindow(AliITSReconstructor::GetRecoParam()->GetTrackleterPhiWindow()); SetThetaWindow(AliITSReconstructor::GetRecoParam()->GetTrackleterThetaWindow()); SetPhiShift(AliITSReconstructor::GetRecoParam()->GetTrackleterPhiShift()); SetRemoveClustersFromOverlaps(AliITSReconstructor::GetRecoParam()->GetTrackleterRemoveClustersFromOverlaps()); SetPhiOverlapCut(AliITSReconstructor::GetRecoParam()->GetTrackleterPhiOverlapCut()); SetZetaOverlapCut(AliITSReconstructor::GetRecoParam()->GetTrackleterZetaOverlapCut()); + SetPhiRotationAngle(AliITSReconstructor::GetRecoParam()->GetTrackleterPhiRotationAngle()); + SetNStdDev(AliITSReconstructor::GetRecoParam()->GetTrackleterNStdDevCut()); + SetScaleDThetaBySin2T(AliITSReconstructor::GetRecoParam()->GetTrackleterScaleDThetaBySin2T()); + SetBuildRefs(AliITSReconstructor::GetRecoParam()->GetTrackleterBuildCl2TrkRefs()); + SetStoreSPD2SingleCl(AliITSReconstructor::GetRecoParam()->GetTrackleterStoreSPD2SingleCl()); + // + 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(); @@ -141,17 +217,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); @@ -170,32 +263,51 @@ 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) : AliTrackleter(mr), -fDetTypeRec(0),fESDEvent(0),fTreeRP(0),fUsedClusLay1(0),fUsedClusLay2(0), -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), +fNSingleClusterSPD2(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), @@ -206,7 +318,21 @@ fhDPhiVsDThetaAcc(0), fhetaTracklets(0), fhphiTracklets(0), fhetaClustersLay1(0), -fhphiClustersLay1(0) +fhphiClustersLay1(0), +fDPhiShift(0), +fDPhiWindow2(0), +fDThetaWindow2(0), +fPartners(0), +fAssociatedLay1(0), +fMinDists(0), +fBlackList(0), +// +fCreateClustersCopy(0), +fClustersLoaded(0), +fRecoDone(0), +fBuildRefs(kTRUE), +fStoreSPD2SingleCl(kFALSE), +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"); @@ -237,46 +363,42 @@ AliITSMultReconstructor::~AliITSMultReconstructor(){ delete fhphiTracklets; delete fhetaClustersLay1; delete fhphiClustersLay1; - delete[] fUsedClusLay1; - delete[] fUsedClusLay2; + // // delete arrays - for(Int_t i=0; iGetPrimaryVertexSPD(); - if (!vtx && vtx->GetNContributors()<0) isVtxOK = kFALSE; + if (!vtx || vtx->GetNContributors()<1) isVtxOK = kFALSE; if (vtx && strstr(vtx->GetTitle(),"cosmics")) { isVtxOK = kFALSE; isCosmics = kTRUE; @@ -299,8 +421,13 @@ void AliITSMultReconstructor::Reconstruct(AliESDEvent* esd, TTree* treeRP) } vtx = 0; } - float vtxf[3] = {vtx->GetX(),vtx->GetY(),vtx->GetZ()}; - FindTracklets(vtxf); + if(vtx){ + float vtxf[3] = {static_cast(vtx->GetX()),static_cast(vtx->GetY()),static_cast(vtx->GetZ())}; + FindTracklets(vtxf); + } + else { + FindTracklets(0); + } // CreateMultiplicityObject(); } @@ -308,314 +435,94 @@ void AliITSMultReconstructor::Reconstruct(AliESDEvent* esd, TTree* treeRP) //____________________________________________________________________ void AliITSMultReconstructor::Reconstruct(TTree* clusterTree, Float_t* vtx, Float_t* /* vtxRes*/) { // - // RS NOTE - this is old reconstructor invocation, to be used from VertexFinder + // 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; + fNSingleClusterSPD2 = 0; // - // - calls LoadClusterArray that finds the position of the clusters - // (in global coord) - // - convert the cluster coordinates to theta, phi (seen from the - // interaction vertex). - // - 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. + 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; - fNClustersLay1 = 0; - fNClustersLay2 = 0; + fNClustersLay[0] = 0; + fNClustersLay[1] = 0; fNTracklets = 0; fNSingleCluster = 0; + fNSingleClusterSPD2 = 0; // if (!clusterTree) { AliError(" Invalid ITS cluster tree !\n"); return; } + if (!clusterTreeMix) { AliError(" Invalid ITS cluster tree 2nd event !\n"); return; } // fESDEvent = 0; - fTreeRP = clusterTree; + 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). 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. + + // Find tracklets converging to vertex // - LoadClusterArrays(fTreeRP); + LoadClusterArrays(fTreeRP,fTreeRPMix); // flag clusters used by ESD tracks - ProcessESDTracks(); + if (fESDEvent) ProcessESDTracks(); + fRecoDone = kTRUE; if (!vtx) return; - const Double_t pi = TMath::Pi(); + InitAux(); - // 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(clPar[kClPh]); - } - } - - // Loop on layer 2 : finding theta, phi and r - for (Int_t iC2=0; iC2 0) { found = 0; - - // Step1: find all tracklets allowing double assocation - // Loop on layer 1 - 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(clPar2[kClTh] - clPar1[kClTh]); - Double_t dPhi = TMath::Abs(clPar2[kClPh] - clPar1[kClPh]); - // 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()) tracklet[kTrDPhi] = tracklet[kTrDPhi]-2.*TMath::Pi(); - if (tracklet[kTrDPhi] < -TMath::Pi()) tracklet[kTrDPhi] = tracklet[kTrDPhi]+2.*TMath::Pi(); - - // store the difference between theta1 and theta2 - tracklet[kTrDTheta] = clPar1[kClTh] - clPar2[kClTh]; - - 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; - Int_t label2 = 0; - while (label2 < 3) { - if ((Int_t) clPar1[kClMC0+label1] != -2 && (Int_t) clPar1[kClMC0+label1] == (Int_t) clPar2[kClMC0+label2]) - break; - label1++; - 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] = clPar1[kClMC0+label1]; - tracklet[kTrLab2] = clPar2[kClMC0+label2]; - } 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] = partners[iC2]; - tracklet[kClID2] = iC2; - // - 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; + for (Int_t iC1=0; iC1SetMultTrackRefs( fBuildRefs ); + fMult->SetSPD2SinglesStored(fStoreSPD2SingleCl); + fMult->SetNumberOfSingleClustersSPD2(fNSingleClusterSPD2); + // 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}}; + if (fBuildRefs) { + 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, fUsedClusLay1[int(tlInfo[kClID1])]|fUsedClusLay2[int(tlInfo[kClID2])]); + fMult->SetTrackletData(i,tlInfo); + // + if (!fBuildRefs) continue; // do we need references? + 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 (fBuildRefs) fMult->AttachTracklet2TrackRefs(refs[0][0],refs[0][1],refs[1][0],refs[1][1]); + // + AliRefArray *refsc[2] = {0,0}; + if (fBuildRefs) 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,fUsedClusLay1[int(clInfo[kSCID])]); + fMult->SetSingleClusterData(i,clInfo); + // + if (!fBuildRefs) continue; // do we need references? + int ilr = i>=(fNSingleCluster-fNSingleClusterSPD2) ? 1:0; + int clID = int(clInfo[kSCID]); + for (int itp=0;itp<2;itp++) { + if (!fStoreRefs[ilr][itp]) continue; + int nref = fUsedClusLay[ilr][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); + } } + // + if (fBuildRefs) fMult->AttachCluster2TrackRefs(refsc[0],refsc[1]); fMult->CompactBits(); // } //____________________________________________________________________ -void AliITSMultReconstructor::LoadClusterArrays(TTree* itsClusterTree) +void AliITSMultReconstructor::LoadClusterArrays(TTree* tree, TTree* treeMix) +{ + // load cluster info and prepare tracklets arrays + // + if (AreClustersLoaded()) {AliInfo("Clusters are already loaded"); return;} + LoadClusterArrays(tree,0); + LoadClusterArrays(treeMix ? treeMix:tree,1); + int nmaxT = TMath::Min(fNClustersLay[0], fNClustersLay[1]); + if (fTracklets) delete[] fTracklets; + fTracklets = new Float_t*[nmaxT]; + memset(fTracklets,0,nmaxT*sizeof(Float_t*)); + // + if (fSClusters) delete[] fSClusters; + int nSlots = GetStoreSPD2SingleCl() ? fNClustersLay[0]+fNClustersLay[1] : fNClustersLay[0]; + fSClusters = new Float_t*[nSlots]; + memset(fSClusters,0,nSlots*sizeof(Float_t*)); + // + AliDebug(1,Form("(clusters in layer 1 : %d, layer 2: %d)",fNClustersLay[0],fNClustersLay[1])); + AliDebug(1,Form("(cluster-fired chips in layer 1 : %d, layer 2: %d)",fNFiredChips[0],fNFiredChips[1])); + SetClustersLoaded(); +} + +//____________________________________________________________________ +void AliITSMultReconstructor::LoadClusterArrays(TTree* itsClusterTree, int il) { // This method - // - gets the clusters from the cluster tree + // - gets the clusters from the cluster tree for layer il // - convert them into global coordinates // - store them in the internal arrays // - count the number of cluster-fired chips // - // RS: This method was strongly modified wrt original by Jan Fiete. In order to have the same numbering + // 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,"Loading clusters and cluster-fired chips ..."); - - fNClustersLay1 = 0; - fNClustersLay2 = 0; - fNFiredChips[0] = 0; - fNFiredChips[1] = 0; - - AliITSsegmentationSPD seg; - - AliITSRecPointContainer* rpcont=AliITSRecPointContainer::Instance(); - TClonesArray* itsClusters=rpcont->FetchClusters(0,itsClusterTree); - if(!rpcont->IsSPDActive()){ - AliWarning("No SPD rec points found, multiplicity not calculated"); - return; - } + 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); + delete[] fClusterCopyIndex[il]; + } // // count clusters // loop over the SPD subdetectors - TObjArray clArr(100); - for (int il=0;il<2;il++) { - int nclLayer = 0; - int detMin = AliITSgeomTGeo::GetModuleIndex(il+1,1,1); - int detMax = AliITSgeomTGeo::GetModuleIndex(il+2,1,1); - for (int idt=detMin;idtUncheckedGetClusters(idt); - int nClusters = itsClusters->GetEntriesFast(); - if (!nClusters) continue; - Int_t nClustersInChip[5] = {0,0,0,0,0}; - while(nClusters--) { - AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters); - if (!cluster) continue; - clArr[nclLayer++] = cluster; - nClustersInChip[ seg.GetChipFromLocal(0,cluster->GetDetLocalZ()) ]++; - } - 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]; - Char_t* usedClusLay = new Char_t[nclLayer]; - // - for (int ic=0;icGetGlobalXYZ( clPar ); - detectorIndexClustersLay[ic] = cluster->GetDetectorIndex(); - overlapFlagClustersLay[ic] = kFALSE; - usedClusLay[ic] = 0; - for (Int_t i=3;i--;) clPar[kClMC0+i] = cluster->GetLabel(i); - } - clArr.Clear(); - delete[] z; - delete[] index; - // - if (il==0) { - fClustersLay1 = clustersLay; - fOverlapFlagClustersLay1 = overlapFlagClustersLay; - fDetectorIndexClustersLay1 = detectorIndexClustersLay; - fUsedClusLay1 = usedClusLay; - fNClustersLay1 = nclLayer; + 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}; + while(nClusters--) { + AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters); + 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 ]++; } - else { - fClustersLay2 = clustersLay; - fOverlapFlagClustersLay2 = overlapFlagClustersLay; - fDetectorIndexClustersLay2 = detectorIndexClustersLay; - fUsedClusLay2 = usedClusLay; - fNClustersLay2 = nclLayer; + 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]; + if (fCreateClustersCopy) fClusterCopyIndex[il] = new Int_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); + if (fCreateClustersCopy) fClusterCopyIndex[il][ic] = index[ic]; + } + clArr.Clear(); + delete[] z; + delete[] index; + // + if (fOverlapFlagClustersLay[il]) delete[] fOverlapFlagClustersLay[il]; + fOverlapFlagClustersLay[il] = overlapFlagClustersLay; + // + if (fDetectorIndexClustersLay[il]) delete[] fDetectorIndexClustersLay[il]; + fDetectorIndexClustersLay[il] = detectorIndexClustersLay; + // + if (fBuildRefs) { + for (int it=0;it<2;it++) { + if (fUsedClusLay[il][it]) delete fUsedClusLay[il][it]; + fUsedClusLay[il][it] = new AliRefArray(nclLayer); } } // - // no double association allowed - int nmaxT = TMath::Min(fNClustersLay1, fNClustersLay2); - fTracklets = new Float_t*[nmaxT]; - fSClusters = new Float_t*[fNClustersLay1]; - for (Int_t i=nmaxT;i--;) fTracklets[i] = 0; + if (fClustersLay[il]) delete[] fClustersLay[il]; + fClustersLay[il] = clustersLay; + fNClustersLay[il] = nclLayer; // - 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])); } + //____________________________________________________________________ -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 @@ -760,9 +741,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; @@ -777,18 +758,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 @@ -826,12 +809,10 @@ 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; Float_t meanRadiusLay1 = 3.99335; // average radius inner layer Float_t meanRadiusLay2 = 7.37935; // average radius outer layer; @@ -841,11 +822,11 @@ AliITSMultReconstructor::FlagClustersInOverlapRegions (Int_t iC1, Int_t iC2WithB 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; @@ -856,32 +837,21 @@ AliITSMultReconstructor::FlagClustersInOverlapRegions (Int_t iC1, Int_t iC2WithB 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.) { -// if (distClSameModMin==0. || distClSameModTMath::Pi()) dePhi=2.*TMath::Pi()-dePhi; @@ -890,21 +860,227 @@ AliITSMultReconstructor::FlagClustersInOverlapRegions (Int_t iC1, Int_t iC2WithB 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.) { -// 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,clPar2[kClTh] - clPar1[kClTh],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 = (AliESDVertex*)fESDEvent->GetPrimaryVertexSPD(); - if (!vtx) { + if (!vtx || vtx->GetNContributors()<1) vtx = (AliESDVertex*)fESDEvent->GetPrimaryVertexSPD(); + if (!vtx || vtx->GetNContributors()<1) { AliDebug(1,"No primary vertex: cannot flag primary tracks"); return; } @@ -923,46 +1100,189 @@ void AliITSMultReconstructor::ProcessESDTracks() for(Int_t itr=0; itrGetTrack(itr); if (!track->IsOn(AliESDtrack::kITSin)) continue; // use only tracks propagated in ITS to vtx - FlagTrackClusters(track); - FlagIfPrimary(track,vtx); + FlagTrackClusters(itr); + FlagIfSecondary(track,vtx); } + FlagV0s(vtx); // } //____________________________________________________________________ -void AliITSMultReconstructor::FlagTrackClusters(const AliESDtrack* track) +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 - // - char mark = track->IsOn(AliESDtrack::kITSpureSA) ? kITSSAPBit : kITSTPCBit; - char *uClus[2] = {fUsedClusLay1,fUsedClusLay2}; - for (int i=AliESDfriendTrack::kMaxITScluster;i--;) { - // note: i>=6 is for extra clusters + 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); - uClus[layID][clID] |= mark; + fUsedClusLay[layID][itsType]->AddReference(clID,id); + fStoreRefs[layID][itsType] = kTRUE; } // } //____________________________________________________________________ -void AliITSMultReconstructor::FlagIfPrimary(AliESDtrack* track, const AliVertex* vtx) +void AliITSMultReconstructor::FlagIfSecondary(AliESDtrack* track, const AliVertex* vtx) { // RS: check if the track is primary and set the flag - const double kPDCASPD1 = 0.1; - const double kPDCASPD0 = 0.3; - // - double cut = (track->HasPointOnITSLayer(0)||track->HasPointOnITSLayer(1)) ? kPDCASPD1 : kPDCASPD0; - // in principle, the track must already have been propagated to vertex - /* - Double_t dzRec[2]={0,0}, covdzRec[3]; - track->PropagateToDCA(vtx, fESDEvent->GetMagneticField(), 3.0, dzRec, covdzRec); - */ - double dist = track->GetD(vtx->GetX(),vtx->GetY(),fESDEvent->GetMagneticField()); - if (TMath::Abs(dist*track->P())SetStatus(AliESDtrack::kMultPrimary); + 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; + // +} + +//____________________________________________________________________ +AliITSRecPoint* AliITSMultReconstructor::GetRecPoint(Int_t lr, Int_t n) const +{ + // return a cluster of lr corresponding to orderer cluster index n + if (fClArr[lr] && fClusterCopyIndex[lr] && nAt(fClusterCopyIndex[lr][n]); + else { + AliError("To access the clusters SetCreateClustersCopy should have been called"); + return 0; + } }