+//____________________________________________________________________
+void AliITSMultReconstructor::ClusterPos2Angles(const Float_t *vtx)
+{
+ // convert cluster coordinates to angles wrt vertex
+ for (int ilr=0;ilr<2;ilr++) {
+ for (Int_t iC=0; iC<fNClustersLay[ilr]; iC++) {
+ float* clPar = GetClusterOfLayer(ilr,iC);
+ CalcThetaPhi(clPar[kClTh]-vtx[0],clPar[kClPh]-vtx[1],clPar[kClZ]-vtx[2],clPar[kClTh],clPar[kClPh]);
+ if (ilr==0) {
+ clPar[kClPh] = clPar[kClPh] + fPhiRotationAngle; // rotation of inner layer for comb studies
+ if (fHistOn) {
+ Float_t eta = clPar[kClTh];
+ eta= TMath::Tan(eta/2.);
+ eta=-TMath::Log(eta);
+ fhetaClustersLay1->Fill(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; iC2<fNClustersLay[1]; iC2++) {
+ //
+ if (fBlackList->IsReferred(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<fNStdDev && d<minDist) { minDist=d; iC2WithBestDist = iC2; }
+ }
+ //
+ if (minDist<fNStdDev) { // This means that a cluster in layer 2 was found that matches with iC1
+ //
+ if (fMinDists[iC2WithBestDist] > 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; iC1<fNClustersLay[0]; iC1++) {
+ float* clPar1 = GetClusterLayer1(iC1);
+ if (fAssociatedLay1[iC1]==2||fAssociatedLay1[iC1]==0) {
+ fSClusters[fNSingleCluster] = new Float_t[kClNPar];
+ fSClusters[fNSingleCluster][kSCTh] = clPar1[kClTh];
+ fSClusters[fNSingleCluster][kSCPh] = clPar1[kClPh];
+ fSClusters[fNSingleCluster][kSCLab] = clPar1[kClMC0];
+ fSClusters[fNSingleCluster][kSCID] = iC1;
+ AliDebug(1,Form(" Adding a single cluster %d (cluster %d of layer 1)",
+ fNSingleCluster, iC1));
+ fNSingleCluster++;
+ }
+ }
+ //
+}
+
+//____________________________________________________________________
+void AliITSMultReconstructor::ProcessESDTracks()
+{
+ // Flag the clusters used by ESD tracks
+ // Flag primary tracks to be used for multiplicity counting
+ //
+ if (!fESDEvent) return;
+ AliESDVertex* vtx = (AliESDVertex*)fESDEvent->GetPrimaryVertexTracks();
+ if (!vtx || vtx->GetNContributors()<1) vtx = (AliESDVertex*)fESDEvent->GetPrimaryVertexSPD();
+ if (!vtx || vtx->GetNContributors()<1) {
+ AliDebug(1,"No primary vertex: cannot flag primary tracks");
+ return;
+ }
+ Int_t ntracks = fESDEvent->GetNumberOfTracks();
+ for(Int_t itr=0; itr<ntracks; itr++) {
+ AliESDtrack* track = fESDEvent->GetTrack(itr);
+ if (!track->IsOn(AliESDtrack::kITSin)) continue; // use only tracks propagated in ITS to vtx
+ FlagTrackClusters(itr);
+ FlagIfSecondary(track,vtx);
+ }
+ FlagV0s(vtx);
+ //
+}
+
+//____________________________________________________________________
+void AliITSMultReconstructor::FlagTrackClusters(Int_t id)
+{
+ // RS: flag the SPD clusters of the track if it is useful for the multiplicity estimation
+ //
+ const AliESDtrack* track = fESDEvent->GetTrack(id);
+ Int_t idx[12];
+ if ( track->GetITSclusters(idx)<3 ) return; // at least 3 clusters must be used in the fit
+ Int_t itsType = track->IsOn(AliESDtrack::kITSpureSA) ? 1:0;
+
+ for (int i=6/*AliESDfriendTrack::kMaxITScluster*/;i--;) { // ignore extras: note: i>=6 is for extra clusters
+ if (idx[i]<0) continue;
+ int layID= (idx[i] & 0xf0000000) >> 28;
+ if (layID>1) continue; // SPD only
+ int clID = (idx[i] & 0x0fffffff);
+ fUsedClusLay[layID][itsType]->AddReference(clID,id);
+ fStoreRefs[layID][itsType] = kTRUE;
+ }
+ //
+}
+
+//____________________________________________________________________
+void AliITSMultReconstructor::FlagIfSecondary(AliESDtrack* track, const AliVertex* vtx)
+{
+ // RS: check if the track is primary and set the flag
+ double cut = (track->HasPointOnITSLayer(0)||track->HasPointOnITSLayer(1)) ? fCutPxDrSPDin:fCutPxDrSPDout;
+ float xz[2];
+ track->GetDZ(vtx->GetX(),vtx->GetY(),vtx->GetZ(), fESDEvent->GetMagneticField(), xz);
+ if (TMath::Abs(xz[0]*track->P())>cut || TMath::Abs(xz[1]*track->P())>fCutPxDz ||
+ TMath::Abs(xz[0])>fCutDCArz || TMath::Abs(xz[1])>fCutDCArz)
+ track->SetStatus(AliESDtrack::kMultSec);
+ else track->ResetStatus(AliESDtrack::kMultSec);
+}
+
+//____________________________________________________________________
+void AliITSMultReconstructor::FlagV0s(const AliESDVertex *vtx)
+{
+ // flag tracks belonging to v0s
+ //
+ const double kK0Mass = 0.4976;
+ //
+ AliV0 pvertex;
+ AliKFVertex vertexKF;
+ AliKFParticle epKF0,epKF1,pipmKF0,piKF0,piKF1,gammaKF,k0KF;
+ Double_t mass,massErr,chi2c;
+ enum {kKFIni=BIT(14)};
+ //
+ double recVtx[3];
+ float recVtxF[3];
+ vtx->GetXYZ(recVtx);
+ for (int i=3;i--;) recVtxF[i] = recVtx[i];
+ //
+ int ntracks = fESDEvent->GetNumberOfTracks();
+ if (ntracks<2) return;
+ //
+ vertexKF.X() = recVtx[0];
+ vertexKF.Y() = recVtx[1];
+ vertexKF.Z() = recVtx[2];
+ vertexKF.Covariance(0,0) = vtx->GetXRes()*vtx->GetXRes();
+ vertexKF.Covariance(1,2) = vtx->GetYRes()*vtx->GetYRes();
+ vertexKF.Covariance(2,2) = vtx->GetZRes()*vtx->GetZRes();
+ //
+ AliESDtrack *trc0,*trc1;
+ for (int it0=0;it0<ntracks;it0++) {
+ trc0 = fESDEvent->GetTrack(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;it1<ntracks;it1++) {
+ trc1 = fESDEvent->GetTrack(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()<fCutMinP) continue;
+ if (pvertex.GetV0CosineOfPointingAngle()<fCutMinPointAngle) continue;
+ if (pvertex.GetDcaV0Daughters()>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 (rv<fCutMinRGamma) break;
+ if (!epKF0.TestBit(kKFIni)) {
+ new(&epKF0) AliKFParticle(*trc0,q0>0 ? 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()<fCutGammaSFromDecay) break;
+ gammaKF.SetMassConstraint(0.,0.001);
+ chi2c = (gammaKF.GetNDF()!=0) ? gammaKF.GetChi2()/gammaKF.GetNDF() : 1000;
+ if (chi2c>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 (rv<fCutMinRK0) break;
+ if (!piKF0.TestBit(kKFIni)) {
+ new(&piKF0) AliKFParticle(*trc0,q0>0 ? 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()<fCutK0SFromDecay) break;
+ k0KF.SetMassConstraint(kK0Mass,0.001);
+ chi2c = (k0KF.GetNDF()!=0) ? k0KF.GetChi2()/k0KF.GetNDF() : 1000;
+ if (chi2c>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;
+ //
+}