X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=ITS%2FAliITSMultReconstructor.cxx;h=88cb54ae516942ebaf9af7d51aa540e6dc09e43d;hb=cd28560b43434095dbb93239ca70d1ce93e690a9;hp=e23bcbad36cf15c160b9676707c830bb2674dfcf;hpb=7fdf95b00c70f48ef1ac4178771ed45ac5e6af4d;p=u%2Fmrichter%2FAliRoot.git diff --git a/ITS/AliITSMultReconstructor.cxx b/ITS/AliITSMultReconstructor.cxx index e23bcbad36c..88cb54ae516 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" @@ -86,6 +87,7 @@ #include "AliV0.h" #include "AliKFParticle.h" #include "AliKFVertex.h" +#include "AliRefArray.h" //____________________________________________________________________ ClassImp(AliITSMultReconstructor) @@ -93,25 +95,22 @@ 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), @@ -144,22 +143,51 @@ 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()); @@ -189,6 +217,8 @@ fhphiClustersLay1(0){ SetRemoveClustersFromOverlaps(); SetPhiOverlapCut(); SetZetaOverlapCut(); + SetPhiRotationAngle(); + // SetCutPxDrSPDin(); SetCutPxDrSPDout(); @@ -211,16 +241,10 @@ fhphiClustersLay1(0){ 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); @@ -239,32 +263,29 @@ 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), @@ -297,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"); @@ -328,46 +363,42 @@ AliITSMultReconstructor::~AliITSMultReconstructor(){ delete fhphiTracklets; delete fhetaClustersLay1; delete fhphiClustersLay1; - delete[] fUsedClusLay1; - delete[] fUsedClusLay2; + // // delete arrays - for(Int_t i=0; iGetX(),vtx->GetY(),vtx->GetZ()}; - FindTracklets(vtxf); + if(vtx){ + float vtxf[3] = {vtx->GetX(),vtx->GetY(),vtx->GetZ()}; + FindTracklets(vtxf); + } + else { + FindTracklets(0); + } // CreateMultiplicityObject(); } @@ -399,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 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.AddAtAndExpand(cluster,nclLayer++); - 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 @@ -851,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; @@ -868,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 @@ -917,9 +809,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; @@ -932,11 +824,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; @@ -947,7 +839,7 @@ 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.) fOverlapFlagClustersLay[0][iiC1]=kTRUE; // if (distClSameMod<=1.) { // if (distClSameModMin==0. || distClSameModTMath::Pi()) dePhi=2.*TMath::Pi()-dePhi; @@ -981,7 +873,7 @@ 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.) 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) { @@ -1015,7 +1122,7 @@ 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); + FlagTrackClusters(itr); FlagIfSecondary(track,vtx); } FlagV0s(vtx); @@ -1023,22 +1130,22 @@ void AliITSMultReconstructor::ProcessESDTracks() } //____________________________________________________________________ -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; } // } @@ -1189,3 +1296,15 @@ Bool_t AliITSMultReconstructor::CanBeElectron(const AliESDtrack* trc) const 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; + } +}