+}
+
+//____________________________________________________________________
+void AliITSMultReconstructor::InitAux()
+{
+ // init arrays/parameters for tracklet reconstruction
+
+ // dPhi shift is field dependent, get average magnetic field
+ Float_t bz = 0;
+ AliMagF* field = 0;
+ if (TGeoGlobalMagField::Instance()) field = dynamic_cast<AliMagF*>(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; i<fNClustersLay[1]; i++) {
+ fPartners[i] = -1;
+ fMinDists[i] = 2*fNStdDev;
+ }
+ memset(fAssociatedLay1,0,fNClustersLay[0]*sizeof(Int_t));
+ //
+}
+
+//____________________________________________________________________
+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,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<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;
+}