+ //
+ Int_t negativepair[30000];
+ Int_t cnegative[3000];
+ Int_t cused1[3000];
+ Int_t positivepair[30000];
+ Int_t cpositive[3000];
+ Int_t cused2[3000];
+ for (Int_t i=0;i<3000;i++) {cnegative[i]=0; cused1[i]=0;}
+ for (Int_t i=0;i<3000;i++) {cpositive[i]=0; cused2[i]=0;}
+ for (Int_t i=0;i<30000;i++){negativepair[i]=positivepair[i]=0;}
+ static Short_t pairs[1000][1000];
+ memset(pairs,0,sizeof(Short_t)*1000000);
+// Short_t ** pairs = new Short_t*[1000];
+// for (Int_t i=0; i<1000; i++) {
+// pairs[i] = new Short_t[1000];
+// memset(pairs[i],0,sizeof(Short_t)*1000);
+// }
+ //
+ // find available pairs
+ //
+ for (Int_t i=0; i<np; i++) {
+ Float_t yp=pos[i].GetY()*fYpitchSSD;
+ if (pos[i].GetQ()<3) continue;
+ for (Int_t j=0; j<nn; j++) {
+ if (neg[j].GetQ()<3) continue;
+ Float_t yn=neg[j].GetY()*fYpitchSSD;
+ Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
+ Float_t yt=yn + tann*zt;
+ zt-=fHlSSD; yt-=fHwSSD;
+ if (TMath::Abs(yt)<fHwSSD+0.01)
+ if (TMath::Abs(zt)<fHlSSD+0.01*(neg[j].GetNd()+pos[i].GetNd())) {
+ negativepair[i*10+cnegative[i]] =j; //index
+ positivepair[j*10+cpositive[j]] =i;
+ cnegative[i]++; //counters
+ cpositive[j]++;
+ pairs[i][j]=100;
+ }
+ }
+ }
+ //
+ for (Int_t i=0; i<np; i++) {
+ Float_t yp=pos[i].GetY()*fYpitchSSD;
+ if (pos[i].GetQ()<3) continue;
+ for (Int_t j=0; j<nn; j++) {
+ if (neg[j].GetQ()<3) continue;
+ if (cpositive[j]&&cnegative[i]) continue;
+ Float_t yn=neg[j].GetY()*fYpitchSSD;
+ Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
+ Float_t yt=yn + tann*zt;
+ zt-=fHlSSD; yt-=fHwSSD;
+ if (TMath::Abs(yt)<fHwSSD+0.1)
+ if (TMath::Abs(zt)<fHlSSD+0.15) {
+ if (cnegative[i]==0) pos[i].SetNd(100); // not available pair
+ if (cpositive[j]==0) neg[j].SetNd(100); // not available pair
+ negativepair[i*10+cnegative[i]] =j; //index
+ positivepair[j*10+cpositive[j]] =i;
+ cnegative[i]++; //counters
+ cpositive[j]++;
+ pairs[i][j]=100;
+ }
+ }
+ }
+ //
+ Float_t lp[5];
+ Int_t milab[10];
+ Double_t ratio;
+
+ //
+ // sign gold tracks
+ //
+ for (Int_t ip=0;ip<np;ip++){
+ Float_t ybest=1000,zbest=1000,qbest=0;
+ //
+ // select gold clusters
+ if ( (cnegative[ip]==1) && cpositive[negativepair[10*ip]]==1){
+ Float_t yp=pos[ip].GetY()*fYpitchSSD;
+ Int_t j = negativepair[10*ip];
+ ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
+ //
+ Float_t yn=neg[j].GetY()*fYpitchSSD;
+ Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
+ Float_t yt=yn + tann*zt;
+ zt-=fHlSSD; yt-=fHwSSD;
+ ybest=yt; zbest=zt;
+ qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
+ lp[0]=-(-ybest+fYshift[fI]);
+ lp[1]= -zbest+fZshift[fI];
+ lp[2]=0.0025*0.0025; //SigmaY2
+ lp[3]=0.110*0.110; //SigmaZ2
+
+ lp[4]=qbest; //Q
+ for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
+ for (Int_t ilab=0;ilab<3;ilab++){
+ milab[ilab] = pos[ip].GetLabel(ilab);
+ milab[ilab+3] = neg[j].GetLabel(ilab);
+ }
+ //
+ CheckLabels2(milab);
+ milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
+ Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fI]};
+ AliITSclusterV2 * cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info);
+ ncl++;
+ cl2->SetChargeRatio(ratio);
+ cl2->SetType(1);
+ pairs[ip][j]=1;
+ if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
+ cl2->SetType(2);
+ pairs[ip][j]=2;
+ }
+ cused1[ip]++;
+ cused2[j]++;
+ }
+ }
+
+ for (Int_t ip=0;ip<np;ip++){
+ Float_t ybest=1000,zbest=1000,qbest=0;
+ //
+ //
+ // select "silber" cluster
+ if ( cnegative[ip]==1 && cpositive[negativepair[10*ip]]==2){
+ Int_t in = negativepair[10*ip];
+ Int_t ip2 = positivepair[10*in];
+ if (ip2==ip) ip2 = positivepair[10*in+1];
+ Float_t pcharge = pos[ip].GetQ()+pos[ip2].GetQ();
+ if (TMath::Abs(pcharge-neg[in].GetQ())<10){
+ //
+ // add first pair
+ if (pairs[ip][in]==100){ //
+ Float_t yp=pos[ip].GetY()*fYpitchSSD;
+ Float_t yn=neg[in].GetY()*fYpitchSSD;
+ Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
+ Float_t yt=yn + tann*zt;
+ zt-=fHlSSD; yt-=fHwSSD;
+ ybest =yt; zbest=zt;
+ qbest =pos[ip].GetQ();
+ lp[0]=-(-ybest+fYshift[fI]);
+ lp[1]= -zbest+fZshift[fI];
+ lp[2]=0.0025*0.0025; //SigmaY2
+ lp[3]=0.110*0.110; //SigmaZ2
+
+ lp[4]=qbest; //Q
+ for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
+ for (Int_t ilab=0;ilab<3;ilab++){
+ milab[ilab] = pos[ip].GetLabel(ilab);
+ milab[ilab+3] = neg[in].GetLabel(ilab);
+ }
+ //
+ CheckLabels2(milab);
+ ratio = (pos[ip].GetQ()-neg[in].GetQ())/(pos[ip].GetQ()+neg[in].GetQ());
+ milab[3]=(((ip<<10) + in)<<10) + idet; // pos|neg|det
+ Int_t info[3] = {pos[ip].GetNd(),neg[in].GetNd(),fNlayer[fI]};
+ AliITSclusterV2 * cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info);
+ ncl++;
+ cl2->SetChargeRatio(ratio);
+ cl2->SetType(5);
+ pairs[ip][in] = 5;
+ if ((pos[ip].GetNd()+neg[in].GetNd())>6){ //multi cluster
+ cl2->SetType(6);
+ pairs[ip][in] = 6;
+ }
+ }
+ //
+ // add second pair
+
+ // if (!(cused1[ip2] || cused2[in])){ //
+ if (pairs[ip2][in]==100){
+ Float_t yp=pos[ip2].GetY()*fYpitchSSD;
+ Float_t yn=neg[in].GetY()*fYpitchSSD;
+ Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
+ Float_t yt=yn + tann*zt;
+ zt-=fHlSSD; yt-=fHwSSD;
+ ybest =yt; zbest=zt;
+ qbest =pos[ip2].GetQ();
+ lp[0]=-(-ybest+fYshift[fI]);
+ lp[1]= -zbest+fZshift[fI];
+ lp[2]=0.0025*0.0025; //SigmaY2
+ lp[3]=0.110*0.110; //SigmaZ2
+
+ lp[4]=qbest; //Q
+ for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
+ for (Int_t ilab=0;ilab<3;ilab++){
+ milab[ilab] = pos[ip2].GetLabel(ilab);
+ milab[ilab+3] = neg[in].GetLabel(ilab);
+ }
+ //
+ CheckLabels2(milab);
+ ratio = (pos[ip2].GetQ()-neg[in].GetQ())/(pos[ip2].GetQ()+neg[in].GetQ());
+ milab[3]=(((ip2<<10) + in)<<10) + idet; // pos|neg|det
+ Int_t info[3] = {pos[ip2].GetNd(),neg[in].GetNd(),fNlayer[fI]};
+ AliITSclusterV2 *cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info);
+ ncl++;
+ cl2->SetChargeRatio(ratio);
+ cl2->SetType(5);
+ pairs[ip2][in] =5;
+ if ((pos[ip2].GetNd()+neg[in].GetNd())>6){ //multi cluster
+ cl2->SetType(6);
+ pairs[ip2][in] =6;
+ }
+ }
+ cused1[ip]++;
+ cused1[ip2]++;
+ cused2[in]++;
+ }
+ }
+ }
+
+ //
+ for (Int_t jn=0;jn<nn;jn++){
+ if (cused2[jn]) continue;
+ Float_t ybest=1000,zbest=1000,qbest=0;
+ // select "silber" cluster
+ if ( cpositive[jn]==1 && cnegative[positivepair[10*jn]]==2){
+ Int_t ip = positivepair[10*jn];
+ Int_t jn2 = negativepair[10*ip];
+ if (jn2==jn) jn2 = negativepair[10*ip+1];
+ Float_t pcharge = neg[jn].GetQ()+neg[jn2].GetQ();
+ //
+ if (TMath::Abs(pcharge-pos[ip].GetQ())<10){
+ //
+ // add first pair
+ // if (!(cused1[ip]||cused2[jn])){
+ if (pairs[ip][jn]==100){
+ Float_t yn=neg[jn].GetY()*fYpitchSSD;
+ Float_t yp=pos[ip].GetY()*fYpitchSSD;
+ Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
+ Float_t yt=yn + tann*zt;
+ zt-=fHlSSD; yt-=fHwSSD;
+ ybest =yt; zbest=zt;
+ qbest =neg[jn].GetQ();
+ lp[0]=-(-ybest+fYshift[fI]);
+ lp[1]= -zbest+fZshift[fI];
+ lp[2]=0.0025*0.0025; //SigmaY2
+ lp[3]=0.110*0.110; //SigmaZ2
+
+ lp[4]=qbest; //Q
+ for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
+ for (Int_t ilab=0;ilab<3;ilab++){
+ milab[ilab] = pos[ip].GetLabel(ilab);
+ milab[ilab+3] = neg[jn].GetLabel(ilab);
+ }
+ //
+ CheckLabels2(milab);
+ ratio = (pos[ip].GetQ()-neg[jn].GetQ())/(pos[ip].GetQ()+neg[jn].GetQ());
+ milab[3]=(((ip<<10) + jn)<<10) + idet; // pos|neg|det
+ Int_t info[3] = {pos[ip].GetNd(),neg[jn].GetNd(),fNlayer[fI]};
+ AliITSclusterV2 * cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info);
+ ncl++;
+ cl2->SetChargeRatio(ratio);
+ cl2->SetType(7);
+ pairs[ip][jn] =7;
+ if ((pos[ip].GetNd()+neg[jn].GetNd())>6){ //multi cluster
+ cl2->SetType(8);
+ pairs[ip][jn]=8;
+ }
+ }
+ //
+ // add second pair
+ // if (!(cused1[ip]||cused2[jn2])){
+ if (pairs[ip][jn2]==100){
+ Float_t yn=neg[jn2].GetY()*fYpitchSSD;
+ Double_t yp=pos[ip].GetY()*fYpitchSSD;
+ Double_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
+ Double_t yt=yn + tann*zt;
+ zt-=fHlSSD; yt-=fHwSSD;
+ ybest =yt; zbest=zt;
+ qbest =neg[jn2].GetQ();
+ lp[0]=-(-ybest+fYshift[fI]);
+ lp[1]= -zbest+fZshift[fI];
+ lp[2]=0.0025*0.0025; //SigmaY2
+ lp[3]=0.110*0.110; //SigmaZ2
+
+ lp[4]=qbest; //Q
+ for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
+ for (Int_t ilab=0;ilab<3;ilab++){
+ milab[ilab] = pos[ip].GetLabel(ilab);
+ milab[ilab+3] = neg[jn2].GetLabel(ilab);
+ }
+ //
+ CheckLabels2(milab);
+ ratio = (pos[ip].GetQ()-neg[jn2].GetQ())/(pos[ip].GetQ()+neg[jn2].GetQ());
+ milab[3]=(((ip<<10) + jn2)<<10) + idet; // pos|neg|det
+ Int_t info[3] = {pos[ip].GetNd(),neg[jn2].GetNd(),fNlayer[fI]};
+ AliITSclusterV2* cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info);
+ ncl++;
+ cl2->SetChargeRatio(ratio);
+ pairs[ip][jn2]=7;
+ cl2->SetType(7);
+ if ((pos[ip].GetNd()+neg[jn2].GetNd())>6){ //multi cluster
+ cl2->SetType(8);
+ pairs[ip][jn2]=8;
+ }
+ }
+ cused1[ip]++;
+ cused2[jn]++;
+ cused2[jn2]++;
+ }
+ }
+ }
+
+ for (Int_t ip=0;ip<np;ip++){
+ Float_t ybest=1000,zbest=1000,qbest=0;
+ //
+ // 2x2 clusters
+ //
+ if ( (cnegative[ip]<5) && cpositive[negativepair[10*ip]]<5){
+ Float_t minchargediff =4.;
+ Int_t j=-1;
+ for (Int_t di=0;di<cnegative[ip];di++){
+ Int_t jc = negativepair[ip*10+di];
+ Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
+ if (TMath::Abs(chargedif)<minchargediff){
+ j =jc;
+ minchargediff = TMath::Abs(chargedif);
+ }
+ }
+ if (j<0) continue; // not proper cluster
+ Int_t count =0;
+ for (Int_t di=0;di<cnegative[ip];di++){
+ Int_t jc = negativepair[ip*10+di];
+ Float_t chargedif = pos[ip].GetQ()-neg[jc].GetQ();
+ if (TMath::Abs(chargedif)<minchargediff+3.) count++;
+ }
+ if (count>1) continue; // more than one "proper" cluster for positive
+ //
+ count =0;
+ for (Int_t dj=0;dj<cpositive[j];dj++){
+ Int_t ic = positivepair[j*10+dj];
+ Float_t chargedif = pos[ic].GetQ()-neg[j].GetQ();
+ if (TMath::Abs(chargedif)<minchargediff+3.) count++;
+ }
+ if (count>1) continue; // more than one "proper" cluster for negative
+
+ Int_t jp = 0;
+
+ count =0;
+ for (Int_t dj=0;dj<cnegative[jp];dj++){
+ Int_t ic = positivepair[jp*10+dj];
+ Float_t chargedif = pos[ic].GetQ()-neg[jp].GetQ();
+ if (TMath::Abs(chargedif)<minchargediff+4.) count++;
+ }
+ if (count>1) continue;
+ if (pairs[ip][j]<100) continue;
+ //
+ //almost gold clusters
+ Float_t yp=pos[ip].GetY()*fYpitchSSD;
+ Float_t yn=neg[j].GetY()*fYpitchSSD;
+ Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp);
+ Float_t yt=yn + tann*zt;
+ zt-=fHlSSD; yt-=fHwSSD;
+ ybest=yt; zbest=zt;
+ qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ());
+ lp[0]=-(-ybest+fYshift[fI]);
+ lp[1]= -zbest+fZshift[fI];
+ lp[2]=0.0025*0.0025; //SigmaY2
+ lp[3]=0.110*0.110; //SigmaZ2
+ lp[4]=qbest; //Q
+ for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2;
+ for (Int_t ilab=0;ilab<3;ilab++){
+ milab[ilab] = pos[ip].GetLabel(ilab);
+ milab[ilab+3] = neg[j].GetLabel(ilab);
+ }
+ //
+ CheckLabels2(milab);
+ ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ());
+ milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det
+ Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fI]};
+ AliITSclusterV2 * cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info);
+ ncl++;
+ cl2->SetChargeRatio(ratio);
+ cl2->SetType(10);
+ pairs[ip][j]=10;
+ if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster
+ cl2->SetType(11);
+ pairs[ip][j]=11;
+ }
+ cused1[ip]++;
+ cused2[j]++;
+ }
+
+ }
+
+ //