From 7284b2b2c8a961f4b4260499e9217194cc6b2f54 Mon Sep 17 00:00:00 2001 From: masera Date: Fri, 3 Jul 2009 20:58:10 +0000 Subject: [PATCH] New algorithm in Multiplicity reconstructor + several mods in the TrackleterSPDEff class, in the ITS RecoParam and in the reconstructor (D. Elia) --- ITS/AliITSMultReconstructor.cxx | 415 +++++++++++++++++--------------- ITS/AliITSMultReconstructor.h | 17 +- ITS/AliITSRecoParam.cxx | 9 +- ITS/AliITSRecoParam.h | 17 +- ITS/AliITSReconstructor.cxx | 4 +- ITS/AliITSTrackleterSPDEff.cxx | 32 +-- ITS/AliITSTrackleterSPDEff.h | 10 +- 7 files changed, 265 insertions(+), 239 deletions(-) diff --git a/ITS/AliITSMultReconstructor.cxx b/ITS/AliITSMultReconstructor.cxx index 26dffa83ba9..c6f0769ed60 100644 --- a/ITS/AliITSMultReconstructor.cxx +++ b/ITS/AliITSMultReconstructor.cxx @@ -15,48 +15,50 @@ /* $Id$ */ -//____________________________________________________________________ -// -// AliITSMultReconstructor - find clusters in the pixels (theta and -// phi) and tracklets. +//_________________________________________________________________________ // -// These can be used to extract charged particles multiplcicity from the ITS. +// Implementation of the ITS-SPD trackleter class // -// A tracklet consist of two ITS clusters, one in the first pixel -// layer and one in the second. The clusters are associates if the -// differencies in Phi (azimuth) and Zeta (longitudinal) are inside -// a fiducial volume. In case of multiple candidates it is selected the -// candidate with minimum distance in Phi. -// The parameter AssociationChoice allows to control if two clusters -// in layer 2 can be associated to the same cluster in layer 1 or not. -// (TRUE means double associations exluded; default = TRUE) +// Retrieves clusters in the pixels (theta and phi) and tracklets. +// These can be used to extract charged particle multiplcicity from the ITS. // -// Two methods return the number of traklets and the number of clusters -// in the first SPD layer (GetNTracklets GetNSingleClusters) +// A tracklet consist of two ITS clusters, one in the first pixel layer and +// one in the second. The clusters are associates if the differencies in +// Phi (azimuth) and Theta (polar) are within fiducial values. +// In case of multiple candidates it is selected the candidate with minimum +// distance. // -// ----------------------------------------------------------------- -// -// NOTE: The cuts on phi and zeta depend on the interacting system (p-p -// or Pb-Pb). Please, check the file AliITSMultReconstructor.h and be -// sure that SetPhiWindow and SetZetaWindow are defined accordingly. +// Two methods return the number of traklets and the number of unassociated +// clusters (i.e. not used in any tracklet) in the first SPD layer +// (GetNTracklets and GetNSingleClusters) +// +// The cuts on phi and theta depend on the interacting system (p-p or Pb-Pb) +// and can be set via AliITSRecoParam class +// (SetPhiWindow and SetThetaWindow) // -// Author : Tiziano Virgili -// -// Recent updates (D. Elia, INFN Bari): -// - multiple association forbidden (fOnlyOneTrackletPerC2 = kTRUE) +// Origin: Tiziano Virgili +// +// Current support and development: +// Domenico Elia, Maria Nicassio (INFN Bari) +// Domenico.Elia@ba.infn.it, Maria.Nicassio@ba.infn.it +// +// Most recent updates: +// - multiple association forbidden (fOnlyOneTrackletPerC2 = kTRUE) // - phi definition changed to ALICE convention (0,2*TMath::pi()) // - cluster coordinates taken with GetGlobalXYZ() // - fGeometry removed // - number of fired chips on the two layers // - option to avoid duplicates in the overlaps (fRemoveClustersFromOverlaps) // - options and fiducial cuts via AliITSRecoParam -// -//____________________________________________________________________ +// - move from DeltaZeta to DeltaTheta cut (M. Nicassio) +// - update to the new algo (M. Nicassio, J.F.Grosse-Oetringhaus) +//_________________________________________________________________________ #include #include #include #include +#include "TArrayI.h" #include "AliITSMultReconstructor.h" #include "AliITSReconstructor.h" @@ -64,6 +66,8 @@ #include "AliITSRecPoint.h" #include "AliITSgeom.h" #include "AliLog.h" +//#include "TGeoGlobalMagField.h" +//#include "AliMagF.h" //____________________________________________________________________ ClassImp(AliITSMultReconstructor) @@ -80,28 +84,22 @@ fOverlapFlagClustersLay1(0), fOverlapFlagClustersLay2(0), fTracklets(0), fSClusters(0), -fAssociationFlag(0), fNClustersLay1(0), fNClustersLay2(0), fNTracklets(0), fNSingleCluster(0), -fOnlyOneTrackletPerC2(0), fPhiWindow(0), -fZetaWindow(0), +fThetaWindow(0), fRemoveClustersFromOverlaps(0), fPhiOverlapCut(0), fZetaOverlapCut(0), fHistOn(0), fhClustersDPhiAcc(0), fhClustersDThetaAcc(0), -fhClustersDZetaAcc(0), fhClustersDPhiAll(0), fhClustersDThetaAll(0), -fhClustersDZetaAll(0), fhDPhiVsDThetaAll(0), fhDPhiVsDThetaAcc(0), -fhDPhiVsDZetaAll(0), -fhDPhiVsDZetaAcc(0), fhetaTracklets(0), fhphiTracklets(0), fhetaClustersLay1(0), @@ -117,16 +115,14 @@ fhphiClustersLay1(0){ SetHistOn(); if(AliITSReconstructor::GetRecoParam()) { - SetOnlyOneTrackletPerC2(AliITSReconstructor::GetRecoParam()->GetTrackleterOnlyOneTrackletPerC2()); SetPhiWindow(AliITSReconstructor::GetRecoParam()->GetTrackleterPhiWindow()); - SetZetaWindow(AliITSReconstructor::GetRecoParam()->GetTrackleterZetaWindow()); + SetThetaWindow(AliITSReconstructor::GetRecoParam()->GetTrackleterThetaWindow()); SetRemoveClustersFromOverlaps(AliITSReconstructor::GetRecoParam()->GetTrackleterRemoveClustersFromOverlaps()); SetPhiOverlapCut(AliITSReconstructor::GetRecoParam()->GetTrackleterPhiOverlapCut()); SetZetaOverlapCut(AliITSReconstructor::GetRecoParam()->GetTrackleterZetaOverlapCut()); } else { - SetOnlyOneTrackletPerC2(); SetPhiWindow(); - SetZetaWindow(); + SetThetaWindow(); SetRemoveClustersFromOverlaps(); SetPhiOverlapCut(); SetZetaOverlapCut(); @@ -141,41 +137,31 @@ fhphiClustersLay1(0){ fOverlapFlagClustersLay2 = new Bool_t[300000]; fTracklets = new Float_t*[300000]; fSClusters = new Float_t*[300000]; - fAssociationFlag = new Bool_t[300000]; for(Int_t i=0; i<300000; i++) { fClustersLay1[i] = new Float_t[6]; fClustersLay2[i] = new Float_t[6]; - fTracklets[i] = new Float_t[5]; + fTracklets[i] = new Float_t[6]; fSClusters[i] = new Float_t[2]; fOverlapFlagClustersLay1[i] = kFALSE; fOverlapFlagClustersLay2[i] = kFALSE; - fAssociationFlag[i] = kFALSE; } // definition of histograms - fhClustersDPhiAcc = new TH1F("dphiacc", "dphi", 100,0.,0.1); + fhClustersDPhiAcc = new TH1F("dphiacc", "dphi", 100,-0.1,0.1); fhClustersDPhiAcc->SetDirectory(0); fhClustersDThetaAcc = new TH1F("dthetaacc","dtheta",100,-0.1,0.1); fhClustersDThetaAcc->SetDirectory(0); - fhClustersDZetaAcc = new TH1F("dzetaacc","dzeta",100,-1.,1.); - fhClustersDZetaAcc->SetDirectory(0); - fhDPhiVsDZetaAcc = new TH2F("dphiVsDzetaacc","",100,-1.,1.,100,0.,0.1); - fhDPhiVsDZetaAcc->SetDirectory(0); - fhDPhiVsDThetaAcc = new TH2F("dphiVsDthetaAcc","",100,-0.1,0.1,100,0.,0.1); + fhDPhiVsDThetaAcc = new TH2F("dphiVsDthetaAcc","",100,-0.1,0.1,100,-0.1,0.1); fhDPhiVsDThetaAcc->SetDirectory(0); fhClustersDPhiAll = new TH1F("dphiall", "dphi", 100,0.0,0.5); fhClustersDPhiAll->SetDirectory(0); - fhClustersDThetaAll = new TH1F("dthetaall","dtheta",100,-0.5,0.5); + fhClustersDThetaAll = new TH1F("dthetaall","dtheta",100,0.0,0.5); fhClustersDThetaAll->SetDirectory(0); - fhClustersDZetaAll = new TH1F("dzetaall","dzeta",100,-5.,5.); - fhClustersDZetaAll->SetDirectory(0); - fhDPhiVsDZetaAll = new TH2F("dphiVsDzetaall","",100,-5.,5.,100,0.,0.5); - fhDPhiVsDZetaAll->SetDirectory(0); - fhDPhiVsDThetaAll = new TH2F("dphiVsDthetaAll","",100,-0.5,0.5,100,0.,0.5); + fhDPhiVsDThetaAll = new TH2F("dphiVsDthetaAll","",100,0.,0.5,100,0.,0.5); fhDPhiVsDThetaAll->SetDirectory(0); fhetaTracklets = new TH1F("etaTracklets", "eta", 100,-2.,2.); @@ -198,28 +184,22 @@ fOverlapFlagClustersLay1(mr.fOverlapFlagClustersLay1), fOverlapFlagClustersLay2(mr.fOverlapFlagClustersLay2), fTracklets(mr.fTracklets), fSClusters(mr.fSClusters), -fAssociationFlag(mr.fAssociationFlag), fNClustersLay1(mr.fNClustersLay1), fNClustersLay2(mr.fNClustersLay2), fNTracklets(mr.fNTracklets), fNSingleCluster(mr.fNSingleCluster), -fOnlyOneTrackletPerC2(mr.fOnlyOneTrackletPerC2), fPhiWindow(mr.fPhiWindow), -fZetaWindow(mr.fZetaWindow), +fThetaWindow(mr.fThetaWindow), fRemoveClustersFromOverlaps(mr.fRemoveClustersFromOverlaps), fPhiOverlapCut(mr.fPhiOverlapCut), fZetaOverlapCut(mr.fZetaOverlapCut), fHistOn(mr.fHistOn), fhClustersDPhiAcc(mr.fhClustersDPhiAcc), fhClustersDThetaAcc(mr.fhClustersDThetaAcc), -fhClustersDZetaAcc(mr.fhClustersDZetaAcc), fhClustersDPhiAll(mr.fhClustersDPhiAll), fhClustersDThetaAll(mr.fhClustersDThetaAll), -fhClustersDZetaAll(mr.fhClustersDZetaAll), fhDPhiVsDThetaAll(mr.fhDPhiVsDThetaAll), fhDPhiVsDThetaAcc(mr.fhDPhiVsDThetaAcc), -fhDPhiVsDZetaAll(mr.fhDPhiVsDZetaAll), -fhDPhiVsDZetaAcc(mr.fhDPhiVsDZetaAcc), fhetaTracklets(mr.fhetaTracklets), fhphiTracklets(mr.fhphiTracklets), fhetaClustersLay1(mr.fhetaClustersLay1), @@ -243,14 +223,10 @@ AliITSMultReconstructor::~AliITSMultReconstructor(){ // delete histograms delete fhClustersDPhiAcc; delete fhClustersDThetaAcc; - delete fhClustersDZetaAcc; delete fhClustersDPhiAll; delete fhClustersDThetaAll; - delete fhClustersDZetaAll; delete fhDPhiVsDThetaAll; delete fhDPhiVsDThetaAcc; - delete fhDPhiVsDZetaAll; - delete fhDPhiVsDZetaAcc; delete fhetaTracklets; delete fhphiTracklets; delete fhetaClustersLay1; @@ -271,18 +247,15 @@ AliITSMultReconstructor::~AliITSMultReconstructor(){ delete [] fOverlapFlagClustersLay2; delete [] fTracklets; delete [] fSClusters; - - delete [] fAssociationFlag; } //____________________________________________________________________ -void -AliITSMultReconstructor::Reconstruct(TTree* clusterTree, Float_t* vtx, Float_t* /* vtxRes*/) { +void AliITSMultReconstructor::Reconstruct(TTree* clusterTree, Float_t* vtx, Float_t* /* vtxRes*/) { // // - calls LoadClusterArray that finds the position of the clusters // (in global coord) // - convert the cluster coordinates to theta, phi (seen from the - // interaction vertex). The third coordinate is used for .... + // interaction vertex). // - makes an array of tracklets // // After this method has been called, the clusters of the two layers @@ -292,10 +265,26 @@ AliITSMultReconstructor::Reconstruct(TTree* clusterTree, Float_t* vtx, Float_t* fNClustersLay1 = 0; fNClustersLay2 = 0; fNTracklets = 0; - fNSingleCluster = 0; + fNSingleCluster = 0; + // loading the clusters LoadClusterArrays(clusterTree); + const Double_t pi = TMath::Pi(); + 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; i 0) { + //Printf("Iteration..."); + 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 - Float_t dTheta = fClustersLay2[iC2][0] - fClustersLay1[iC1][0]; - Float_t dPhi = TMath::Abs(fClustersLay2[iC2][1] - fClustersLay1[iC1][1]); + 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>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi; + if (dPhi>pi) dPhi=2.*pi-dPhi; - // find the difference in z (between linear projection from layer 1 - // and the actual point: Dzeta= z1/r1*r2 -z2) - Float_t r2 = fClustersLay2[iC2][2]/TMath::Cos(fClustersLay2[iC2][0]); - Float_t dZeta = fClustersLay1[iC1][2]*r2 - fClustersLay2[iC2][2]; + // Account for B-field shifting dphi ? ... + // Double_t pos[3]={0.,0.,0.}; + // Double_t B[3]={0.,0.,0.}; + // TGeoGlobalMagField::Instance()->Field(pos,B); + // if (B[2]=!0) dPhi-=0.005; // field dependent if (fHistOn) { - fhClustersDPhiAll->Fill(dPhi); + fhClustersDPhiAll->Fill(dPhi); fhClustersDThetaAll->Fill(dTheta); - fhClustersDZetaAll->Fill(dZeta); fhDPhiVsDThetaAll->Fill(dTheta, dPhi); - fhDPhiVsDZetaAll->Fill(dZeta, dPhi); } - // make "elliptical" cut in Phi and Zeta! - Float_t d = TMath::Sqrt(TMath::Power(dPhi/fPhiWindow,2) + TMath::Power(dZeta/fZetaWindow,2)); - - if (d>1) continue; - - //look for the minimum distance: the minimum is in iC2WithBestDist - if (TMath::Sqrt(dZeta*dZeta+(r2*dPhi*r2*dPhi)) < distmin ) { - distmin=TMath::Sqrt(dZeta*dZeta + (r2*dPhi*r2*dPhi)); - dPhimin = dPhi; - dThetamin = dTheta; - dZetamin = dZeta; + // make "elliptical" cut in Phi and Theta! + Float_t d = TMath::Power(dPhi/fPhiWindow,2) + TMath::Power(dTheta/fThetaWindow,2); + + // look for the minimum distance: the minimum is in iC2WithBestDist + if (d<1 && dFill(dPhimin); - fhClustersDThetaAcc->Fill(dThetamin); - fhClustersDZetaAcc->Fill(dZetamin); - fhDPhiVsDThetaAcc->Fill(dThetamin, dPhimin); - fhDPhiVsDZetaAcc->Fill(dZetamin, dPhimin); - } - - if (fOnlyOneTrackletPerC2) fAssociationFlag[iC2WithBestDist] = kTRUE; // flag the association - - // store the tracklet - - // use the theta from the clusters in the first layer - fTracklets[fNTracklets][0] = fClustersLay1[iC1][0]; - // use the phi from the clusters in the first layer - fTracklets[fNTracklets][1] = fClustersLay1[iC1][1]; - // store the difference between phi1 and phi2 - fTracklets[fNTracklets][2] = fClustersLay1[iC1][1] - fClustersLay2[iC2WithBestDist][1]; - - // define dphi in the range [0,pi] with proper sign (track charge correlated) - if (fTracklets[fNTracklets][2] > 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(); - - // 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[iC1][3+label1] != -2 && (Int_t) fClustersLay1[iC1][3+label1] == (Int_t) fClustersLay2[iC2WithBestDist][3+label2]) - break; - label1++; - if (label1 == 3) - { - label1 = 0; - label2++; - } - } + if (minDist<1) { // This means that a cluster in layer 2 was found that matches with iC1 + + if (minDists[iC2WithBestDist] > 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 (best1) 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[iC1] = 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); + } - if (label2 < 3) - { - AliDebug(AliLog::kDebug, Form("Found label %d == %d for tracklet candidate %d\n", (Int_t) fClustersLay1[iC1][3+label1], (Int_t) fClustersLay2[iC2WithBestDist][3+label2], fNTracklets)); - fTracklets[fNTracklets][3] = fClustersLay1[iC1][3+label1]; - fTracklets[fNTracklets][4] = fClustersLay2[iC2WithBestDist][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[iC1][3], (Int_t) fClustersLay1[iC1][4], (Int_t) fClustersLay1[iC1][5], (Int_t) fClustersLay2[iC2WithBestDist][3], (Int_t) fClustersLay2[iC2WithBestDist][4], (Int_t) fClustersLay2[iC2WithBestDist][5], fNTracklets)); - fTracklets[fNTracklets][3] = fClustersLay1[iC1][3]; - fTracklets[fNTracklets][4] = fClustersLay2[iC2WithBestDist][3]; - } + } 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; iC2Fill(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", iC1, - iC2WithBestDist)); - fNTracklets++; + if (partners[iC2] == -1) continue; - if (fRemoveClustersFromOverlaps) FlagClustersInOverlapRegions (iC1,iC2WithBestDist); + //Printf("Tracklet7: %d %d %f %d %d %f %f %f", partners[iC2], iC2, minDists[iC2], (Int_t) fClustersLay1[partners[iC2]][3], (Int_t) fClustersLay2[iC2][3], fClustersLay1[partners[iC2]][0], fClustersLay2[iC2][0], fClustersLay1[partners[iC2]][1]); + + + if (fOverlapFlagClustersLay1[partners[iC2]] || fOverlapFlagClustersLay2[iC2]) continue; + if (fRemoveClustersFromOverlaps) FlagClustersInOverlapRegions (partners[iC2],iC2); + + // use the theta from the clusters in the first layer + fTracklets[fNTracklets][0] = fClustersLay1[partners[iC2]][0]; + // use the phi from the clusters in the first layer + fTracklets[fNTracklets][1] = fClustersLay1[partners[iC2]][1]; + // store the difference between phi1 and phi2 + fTracklets[fNTracklets][2] = fClustersLay1[partners[iC2]][1] - fClustersLay2[iC2][1]; + + // define dphi in the range [0,pi] with proper sign (track charge correlated) + if (fTracklets[fNTracklets][2] > 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]); } - // Delete the following else if you do not want to save Clusters! + // 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++; + } + } + 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]); + } - else { // This means that the cluster has not been associated + 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++; - // store the cluster - + associatedLay1[partners[iC2]] = 1; + } + + // Delete the following else if you do not want to save Clusters! + // store the cluster + for (Int_t iC1=0; iC1Write(); fhClustersDThetaAll->Write(); - fhClustersDZetaAll->Write(); fhDPhiVsDThetaAll->Write(); - fhDPhiVsDZetaAll->Write(); fhClustersDPhiAcc->Write(); fhClustersDThetaAcc->Write(); - fhClustersDZetaAcc->Write(); fhDPhiVsDThetaAcc->Write(); - fhDPhiVsDZetaAcc->Write(); fhetaTracklets->Write(); fhphiTracklets->Write(); diff --git a/ITS/AliITSMultReconstructor.h b/ITS/AliITSMultReconstructor.h index 943e241197b..a00e1934ddd 100644 --- a/ITS/AliITSMultReconstructor.h +++ b/ITS/AliITSMultReconstructor.h @@ -14,11 +14,9 @@ // // A tracklet consist of two ITS clusters, one in the first pixel // layer and one in the second. The clusters are associates if the -// differencies in Phi (azimuth) and Zeta (longitudinal) are inside +// differencies in Phi (azimuth) and Theta (polar) are inside // a fiducial volume. In case of multiple candidates it is selected the // candidate with minimum distance in Phi. -// The boolean fOnlyOneTrackletPerC2 allows to control if two clusters -// in layer 2 can be associated to the same cluster in layer 1 or not. // ///////////////////////////////////////////////////////////////////////// @@ -41,9 +39,8 @@ public: void FlagClustersInOverlapRegions(Int_t ic1,Int_t ic2); // Following members are set via AliITSRecoParam - void SetOnlyOneTrackletPerC2(Bool_t b = kTRUE) {fOnlyOneTrackletPerC2 = b;} void SetPhiWindow(Float_t w=0.08) {fPhiWindow=w;} - void SetZetaWindow(Float_t w=1.) {fZetaWindow=w;} + void SetThetaWindow(Float_t w=0.025) {fThetaWindow=w;} void SetRemoveClustersFromOverlaps(Bool_t b = kFALSE) {fRemoveClustersFromOverlaps = b;} void SetPhiOverlapCut(Float_t w=0.005) {fPhiOverlapCut=w;} void SetZetaOverlapCut(Float_t w=0.05) {fZetaOverlapCut=w;} @@ -77,7 +74,6 @@ protected: Float_t** fTracklets; // tracklets Float_t** fSClusters; // single clusters (unassociated) - Bool_t* fAssociationFlag; // flag for the associations Int_t fNClustersLay1; // Number of clusters (Layer1) Int_t fNClustersLay2; // Number of clusters (Layer2) @@ -86,9 +82,8 @@ protected: Short_t fNFiredChips[2]; // Number of fired chips in the two SPD layers // Following members are set via AliITSRecoParam - Bool_t fOnlyOneTrackletPerC2; // Allow only one tracklet per cluster in the outer layer Float_t fPhiWindow; // Search window in phi - Float_t fZetaWindow; // Search window in eta + Float_t fThetaWindow; // Search window in eta Bool_t fRemoveClustersFromOverlaps; // Option to skip clusters in the overlaps Float_t fPhiOverlapCut; // Fiducial window in phi for overlap cut Float_t fZetaOverlapCut; // Fiducial window in eta for overlap cut @@ -97,15 +92,11 @@ protected: TH1F* fhClustersDPhiAcc; // Phi2 - Phi1 for tracklets TH1F* fhClustersDThetaAcc; // Theta2 - Theta1 for tracklets - TH1F* fhClustersDZetaAcc; // z2 - z1projected for tracklets TH1F* fhClustersDPhiAll; // Phi2 - Phi1 all the combinations TH1F* fhClustersDThetaAll; // Theta2 - Theta1 all the combinations - TH1F* fhClustersDZetaAll; // z2 - z1projected all the combinations TH2F* fhDPhiVsDThetaAll; // 2D plot for all the combinations TH2F* fhDPhiVsDThetaAcc; // same plot for tracklets - TH2F* fhDPhiVsDZetaAll; // 2d plot for all the combination - TH2F* fhDPhiVsDZetaAcc; // same plot for tracklets TH1F* fhetaTracklets; // Pseudorapidity distr. for tracklets TH1F* fhphiTracklets; // Azimuthal (Phi) distr. for tracklets @@ -115,7 +106,7 @@ protected: void LoadClusterArrays(TTree* tree); - ClassDef(AliITSMultReconstructor,4) + ClassDef(AliITSMultReconstructor,5) }; #endif diff --git a/ITS/AliITSRecoParam.cxx b/ITS/AliITSRecoParam.cxx index 8daa29dd991..26e85d08662 100644 --- a/ITS/AliITSRecoParam.cxx +++ b/ITS/AliITSRecoParam.cxx @@ -119,7 +119,9 @@ fUseTrackletsPlaneEff(kFALSE), fMCTrackletsPlaneEff(kFALSE), fBkgTrackletsPlaneEff(kFALSE), fTrackleterPhiWindowL1(0.10), +fTrackleterPhiWindowL2(0.07), fTrackleterZetaWindowL1(0.6), +fTrackleterZetaWindowL2(0.4), fUpdateOncePerEventPlaneEff(kTRUE), fMinContVtxPlaneEff(3), fIPlanePlaneEff(0), @@ -157,9 +159,8 @@ fUseSDDCorrectionMaps(kTRUE), fUseSDDClusterSizeSelection(kFALSE), fMinClusterChargeSDD(0.), fUseChargeMatchingInClusterFinderSSD(kTRUE), -fTrackleterOnlyOneTrackletPerC2(kTRUE), fTrackleterPhiWindow(0.08), -fTrackleterZetaWindow(1.00), +fTrackleterThetaWindow(0.025), fTrackleterRemoveClustersFromOverlaps(kFALSE), fTrackleterPhiOverlapCut(0.005), fTrackleterZetaOverlapCut(0.05), @@ -606,8 +607,8 @@ AliITSRecoParam *AliITSRecoParam::GetPlaneEffParam(Int_t i) param->SetIPlanePlaneEff(i); param->SetComputePlaneEff(kTRUE,kFALSE); param->SetUseTrackletsPlaneEff(kTRUE); - param->SetTrackleterPhiWindow(0.07); - param->SetTrackleterZetaWindow(0.4); + param->SetTrackleterPhiWindowL2(0.07); + param->SetTrackleterZetaWindowL2(0.4); param->SetTrackleterPhiWindowL1(0.10); param->SetTrackleterZetaWindowL1(0.6); param->SetUpdateOncePerEventPlaneEff(kTRUE); diff --git a/ITS/AliITSRecoParam.h b/ITS/AliITSRecoParam.h index a0e4556291b..03ad9d1d87a 100644 --- a/ITS/AliITSRecoParam.h +++ b/ITS/AliITSRecoParam.h @@ -219,8 +219,12 @@ class AliITSRecoParam : public AliDetectorRecoParam Bool_t GetBkgTrackletsPlaneEff() const {return fBkgTrackletsPlaneEff;} void SetTrackleterPhiWindowL1(Float_t w=0.10) {fTrackleterPhiWindowL1=w; return;} Float_t GetTrackleterPhiWindowL1() const {return fTrackleterPhiWindowL1;} + void SetTrackleterPhiWindowL2(Float_t w=0.07) {fTrackleterPhiWindowL2=w; return;} + Float_t GetTrackleterPhiWindowL2() const {return fTrackleterPhiWindowL2;} void SetTrackleterZetaWindowL1(Float_t w=0.6) {fTrackleterZetaWindowL1=w; return;} Float_t GetTrackleterZetaWindowL1() const {return fTrackleterZetaWindowL1;} + void SetTrackleterZetaWindowL2(Float_t w=0.40) {fTrackleterZetaWindowL2=w; return;} + Float_t GetTrackleterZetaWindowL2() const {return fTrackleterZetaWindowL2;} void SetUpdateOncePerEventPlaneEff(Bool_t use=kTRUE) {fUpdateOncePerEventPlaneEff=use; return;} Bool_t GetUpdateOncePerEventPlaneEff() const {return fUpdateOncePerEventPlaneEff;} void SetMinContVtxPlaneEff(Int_t n=3) {fMinContVtxPlaneEff=n; return;} @@ -323,12 +327,10 @@ class AliITSRecoParam : public AliDetectorRecoParam Bool_t GetUseCosmicRunShiftsSSD() const { return fUseCosmicRunShiftsSSD; } // SPD Tracklets (D. Elia) - void SetTrackleterOnlyOneTrackletPerC2(Bool_t use= kTRUE) {fTrackleterOnlyOneTrackletPerC2=use; return; } - Bool_t GetTrackleterOnlyOneTrackletPerC2() const { return fTrackleterOnlyOneTrackletPerC2; } void SetTrackleterPhiWindow(Float_t w=0.08) {fTrackleterPhiWindow=w;} - void SetTrackleterZetaWindow(Float_t w=1.) {fTrackleterZetaWindow=w;} + void SetTrackleterThetaWindow(Float_t w=0.025) {fTrackleterThetaWindow=w;} Float_t GetTrackleterPhiWindow() const {return fTrackleterPhiWindow;} - Float_t GetTrackleterZetaWindow() const {return fTrackleterZetaWindow;} + Float_t GetTrackleterThetaWindow() const {return fTrackleterThetaWindow;} void SetTrackleterRemoveClustersFromOverlaps(Bool_t use=kTRUE) { fTrackleterRemoveClustersFromOverlaps=use; return; } Bool_t GetTrackleterRemoveClustersFromOverlaps() const { return fTrackleterRemoveClustersFromOverlaps; } void SetTrackleterPhiOverlapCut(Float_t w=0.005) {fTrackleterPhiOverlapCut=w;} @@ -509,7 +511,9 @@ class AliITSRecoParam : public AliDetectorRecoParam Bool_t fMCTrackletsPlaneEff; // flag to enable the use of MC info for corrections (SPD PlaneEff using tracklets) Bool_t fBkgTrackletsPlaneEff; // flag to evaluate background instead of normal use (SPD PlaneEff using tracklets) Float_t fTrackleterPhiWindowL1; // Search window in phi for inner layer (1) (SPD PlaneEff using tracklets) + Float_t fTrackleterPhiWindowL2; // Search window in phi for outer layer (2) (SPD PlaneEff using tracklets) Float_t fTrackleterZetaWindowL1; // Search window in zeta for inner layer (1) (SPD PlaneEff using tracklets) + Float_t fTrackleterZetaWindowL2; // Search window in zeta for outer layer (2) (SPD PlaneEff using tracklets) Bool_t fUpdateOncePerEventPlaneEff; // option to update chip efficiency once/event (to avoid doubles) Int_t fMinContVtxPlaneEff; // min number of contributors to ESD vtx for SPD PlaneEff using tracklets Int_t fIPlanePlaneEff; // index of the plane (in the range [-1,5]) to study the efficiency (-1 ->Tracklets) @@ -557,9 +561,8 @@ class AliITSRecoParam : public AliDetectorRecoParam Bool_t fUseChargeMatchingInClusterFinderSSD; // SSD // SPD Tracklets (D. Elia) - Bool_t fTrackleterOnlyOneTrackletPerC2; // Allow only one tracklet per cluster in the outer layer Float_t fTrackleterPhiWindow; // Search window in phi - Float_t fTrackleterZetaWindow; // Search window in eta + Float_t fTrackleterThetaWindow; // Search window in eta Bool_t fTrackleterRemoveClustersFromOverlaps; // Option to skip clusters in the overlaps Float_t fTrackleterPhiOverlapCut; // Fiducial window in phi for overlap cut Float_t fTrackleterZetaOverlapCut; // Fiducial window in eta for overlap cut @@ -592,7 +595,7 @@ class AliITSRecoParam : public AliDetectorRecoParam Bool_t fAlignFilterFillQANtuples; // fill QA ntuples - ClassDef(AliITSRecoParam,24) // ITS reco parameters + ClassDef(AliITSRecoParam,25) // ITS reco parameters }; #endif diff --git a/ITS/AliITSReconstructor.cxx b/ITS/AliITSReconstructor.cxx index f18da628fc7..102b276ae6a 100644 --- a/ITS/AliITSReconstructor.cxx +++ b/ITS/AliITSReconstructor.cxx @@ -133,8 +133,8 @@ AliTracker* AliITSReconstructor::CreateTrackleter() const if(GetRecoParam()->GetBkgTrackletsPlaneEff()) spdtrackleter->SetReflectClusterAroundZAxisForLayer(1,kTRUE); if(GetRecoParam()->GetMCTrackletsPlaneEff()) spdtrackleter->SetMC(); spdtrackleter->SetHistOn(); - spdtrackleter->SetPhiWindow(GetRecoParam()->GetTrackleterPhiWindow()); - spdtrackleter->SetZetaWindow(GetRecoParam()->GetTrackleterZetaWindow()); + spdtrackleter->SetPhiWindowL2(GetRecoParam()->GetTrackleterPhiWindowL2()); + spdtrackleter->SetZetaWindowL2(GetRecoParam()->GetTrackleterZetaWindowL2()); spdtrackleter->SetPhiWindowL1(GetRecoParam()->GetTrackleterPhiWindowL1()); spdtrackleter->SetZetaWindowL1(GetRecoParam()->GetTrackleterZetaWindowL1()); if(GetRecoParam()->GetUpdateOncePerEventPlaneEff()) spdtrackleter->SetUpdateOncePerEventPlaneEff(); diff --git a/ITS/AliITSTrackleterSPDEff.cxx b/ITS/AliITSTrackleterSPDEff.cxx index b0d1cba64f9..873b2d54e70 100644 --- a/ITS/AliITSTrackleterSPDEff.cxx +++ b/ITS/AliITSTrackleterSPDEff.cxx @@ -63,8 +63,8 @@ fNClustersLay1(0), fNClustersLay2(0), fNTracklets(0), fOnlyOneTrackletPerC2(0), -fPhiWindow(0), -fZetaWindow(0), +fPhiWindowL2(0), +fZetaWindowL2(0), fPhiOverlapCut(0), fZetaOverlapCut(0), fHistOn(0), @@ -131,8 +131,8 @@ fhClustersInModuleLay2(0) { // default constructor // from AliITSMultReconstructor - SetPhiWindow(); - SetZetaWindow(); + SetPhiWindowL2(); + SetZetaWindowL2(); SetOnlyOneTrackletPerC2(); fClustersLay1 = new Float_t*[300000]; fClustersLay2 = new Float_t*[300000]; @@ -175,8 +175,8 @@ fNClustersLay1(mr.fNClustersLay1), fNClustersLay2(mr.fNClustersLay2), fNTracklets(mr.fNTracklets), fOnlyOneTrackletPerC2(mr.fOnlyOneTrackletPerC2), -fPhiWindow(mr.fPhiWindow), -fZetaWindow(mr.fZetaWindow), +fPhiWindowL2(mr.fPhiWindowL2), +fZetaWindowL2(mr.fZetaWindowL2), fPhiOverlapCut(mr.fPhiOverlapCut), fZetaOverlapCut(mr.fZetaOverlapCut), fHistOn(mr.fHistOn), @@ -484,8 +484,8 @@ AliITSTrackleterSPDEff::Reconstruct(AliStack *pStack, TTree *tRef) { } // make "elliptical" cut in Phi and Zeta! - Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindow/fPhiWindow + - dZeta*dZeta/fZetaWindow/fZetaWindow); + Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindowL2/fPhiWindowL2 + + dZeta*dZeta/fZetaWindowL2/fZetaWindowL2); if (d>1) continue; @@ -1344,7 +1344,7 @@ void AliITSTrackleterSPDEff::PrintAscii(ostream *os)const{ // none. // Return: // none. - *os << fPhiWindowL1 <<" "<< fZetaWindowL1 << " " << fPhiWindow <<" "<< fZetaWindow + *os << fPhiWindowL1 <<" "<< fZetaWindowL1 << " " << fPhiWindowL2 <<" "<< fZetaWindowL2 << " " << fOnlyOneTrackletPerC1 << " " << fOnlyOneTrackletPerC2 << " " << fUpdateOncePerEventPlaneEff << " " << fReflectClusterAroundZAxisForLayer0 << " " << fReflectClusterAroundZAxisForLayer1; @@ -1378,7 +1378,7 @@ void AliITSTrackleterSPDEff::ReadAscii(istream *is){ // none. Bool_t tmp= fMC; - *is >> fPhiWindowL1 >> fZetaWindowL1 >> fPhiWindow >> fZetaWindow + *is >> fPhiWindowL1 >> fZetaWindowL1 >> fPhiWindowL2 >> fZetaWindowL2 >> fOnlyOneTrackletPerC1 >> fOnlyOneTrackletPerC2 >> fUpdateOncePerEventPlaneEff >> fReflectClusterAroundZAxisForLayer0 >> fReflectClusterAroundZAxisForLayer1; @@ -1457,8 +1457,8 @@ void AliITSTrackleterSPDEff::SavePredictionMC(TString filename) const { TH1F* cuts = new TH1F("cuts", "list of cuts", 10, 0, 10); // TH1I containing cuts cuts->SetBinContent(1,fPhiWindowL1); cuts->SetBinContent(2,fZetaWindowL1); - cuts->SetBinContent(3,fPhiWindow); - cuts->SetBinContent(4,fZetaWindow); + cuts->SetBinContent(3,fPhiWindowL2); + cuts->SetBinContent(4,fZetaWindowL2); cuts->SetBinContent(5,fOnlyOneTrackletPerC1); cuts->SetBinContent(6,fOnlyOneTrackletPerC2); cuts->SetBinContent(7,fUpdateOncePerEventPlaneEff); @@ -1548,8 +1548,8 @@ void AliITSTrackleterSPDEff::ReadPredictionMC(TString filename) { TH1F *cuts = (TH1F*)mcfile->Get("cuts"); fPhiWindowL1=(Float_t)cuts->GetBinContent(1); fZetaWindowL1=(Float_t)cuts->GetBinContent(2); - fPhiWindow=(Float_t)cuts->GetBinContent(3); - fZetaWindow=(Float_t)cuts->GetBinContent(4); + fPhiWindowL2=(Float_t)cuts->GetBinContent(3); + fZetaWindowL2=(Float_t)cuts->GetBinContent(4); fOnlyOneTrackletPerC1=(Bool_t)cuts->GetBinContent(5); fOnlyOneTrackletPerC2=(Bool_t)cuts->GetBinContent(6); fUpdateOncePerEventPlaneEff=(Bool_t)cuts->GetBinContent(7); @@ -1887,8 +1887,8 @@ for(Int_t iref=0;iref