+//__________________________________________________________________________
+void AliITSClusterFinderSDD::ResolveClusters(){
+ // The function to resolve clusters if the clusters overlapping exists
+/* AliITS *iTS=(AliITS*)gAlice->GetModule("ITS");
+ // get number of clusters for this module
+ Int_t nofClusters = fClusters->GetEntriesFast();
+ nofClusters -= fNclusters;
+ //cout<<"Resolve Cl: nofClusters, fNclusters ="<<nofClusters<<","
+ // <<fNclusters<<endl;
+ Int_t fNofMaps = fSegmentation->Npz();
+ Int_t fNofAnodes = fNofMaps/2;
+ Int_t dummy=0;
+ Double_t fTimeStep = fSegmentation->Dpx(dummy);
+ Double_t fSddLength = fSegmentation->Dx();
+ Double_t fDriftSpeed = fResponse->DriftSpeed();
+ Double_t anodePitch = fSegmentation->Dpz(dummy);
+ Float_t n, baseline;
+ fResponse->GetNoiseParam(n,baseline);
+ Float_t dzz_1A = anodePitch * anodePitch / 12;
+ // fill Map of signals
+ fMap->FillMap();
+ Int_t j,i,ii,ianode,anode,itime;
+ Int_t wing,astart,astop,tstart,tstop,nanode;
+ Double_t fadc,ClusterTime;
+ Double_t q[400],x[400],z[400]; // digit charges and coordinates
+ for(j=0; j<nofClusters; j++) {
+ AliITSRawClusterSDD *clusterJ=(AliITSRawClusterSDD*) fClusters->At(j);
+ Int_t ndigits = 0;
+ astart=clusterJ->Astart();
+ astop=clusterJ->Astop();
+ tstart=clusterJ->Tstartf();
+ tstop=clusterJ->Tstopf();
+ nanode=clusterJ->Anodes(); // <- Ernesto
+ wing=(Int_t)clusterJ->W();
+ if(wing == 2) {
+ astart += fNofAnodes;
+ astop += fNofAnodes;
+ } // end if
+ // cout<<"astart,astop,tstart,tstop ="<<astart<<","<<astop<<","
+ // <<tstart<<","<<tstop<<endl;
+ // clear the digit arrays
+ for(ii=0; ii<400; ii++) {
+ q[ii] = 0.;
+ x[ii] = 0.;
+ z[ii] = 0.;
+ } // end for ii
+
+ for(ianode=astart; ianode<=astop; ianode++) {
+ for(itime=tstart; itime<=tstop; itime++) {
+ fadc=fMap->GetSignal(ianode,itime);
+ if(fadc>baseline) {
+ fadc-=(Double_t)baseline;
+ q[ndigits] = fadc*(fTimeStep/160); // KeV
+ anode = ianode;
+ if(wing == 2) anode -= fNofAnodes;
+ z[ndigits] = (anode + 0.5 - fNofAnodes/2)*anodePitch;
+ ClusterTime = itime*fTimeStep;
+ if(ClusterTime > fTimeCorr) ClusterTime -= fTimeCorr;// ns
+ x[ndigits] = fSddLength - ClusterTime*fDriftSpeed;
+ if(wing == 1) x[ndigits] *= (-1);
+ // cout<<"ianode,itime,fadc ="<<ianode<<","<<itime<<","
+ // <<fadc<<endl;
+ // cout<<"wing,anode,ndigits,charge ="<<wing<<","
+ // <<anode<<","<<ndigits<<","<<q[ndigits]<<endl;
+ ndigits++;
+ continue;
+ } // end if
+ fadc=0;
+ // cout<<"fadc=0, ndigits ="<<ndigits<<endl;
+ } // time loop
+ } // anode loop
+ // cout<<"for new cluster ndigits ="<<ndigits<<endl;
+ // Fit cluster to resolve for two separate ones --------------------
+ Double_t qq=0., xm=0., zm=0., xx=0., zz=0., xz=0.;
+ Double_t dxx=0., dzz=0., dxz=0.;
+ Double_t scl = 0., tmp, tga, elps = -1.;
+ Double_t xfit[2], zfit[2], qfit[2];
+ Double_t pitchz = anodePitch*1.e-4; // cm
+ Double_t pitchx = fTimeStep*fDriftSpeed*1.e-4; // cm
+ Double_t sigma2;
+ Int_t nfhits;
+ Int_t nbins = ndigits;
+ Int_t separate = 0;
+ // now, all lengths are in microns
+ for (ii=0; ii<nbins; ii++) {
+ qq += q[ii];
+ xm += x[ii]*q[ii];
+ zm += z[ii]*q[ii];
+ xx += x[ii]*x[ii]*q[ii];
+ zz += z[ii]*z[ii]*q[ii];
+ xz += x[ii]*z[ii]*q[ii];
+ } // end for ii
+ xm /= qq;
+ zm /= qq;
+ xx /= qq;
+ zz /= qq;
+ xz /= qq;
+ dxx = xx - xm*xm;
+ dzz = zz - zm*zm;
+ dxz = xz - xm*zm;
+
+ // shrink the cluster in the time direction proportionaly to the
+ // dxx/dzz, which lineary depends from the drift path
+ // new Ernesto........
+ if( nanode == 1 ){
+ dzz = dzz_1A; // for one anode cluster dzz = anode**2/12
+ scl = TMath::Sqrt( 7.2/(-0.57*xm*1.e-3+71.8) );
+ } // end if
+ if( nanode == 2 ){
+ scl = TMath::Sqrt( (-0.18*xm*1.e-3+21.3)/(-0.57*xm*1.e-3+71.8) );
+ } // end if
+ if( nanode == 3 ){
+ scl = TMath::Sqrt( (-0.5*xm*1.e-3+34.5)/(-0.57*xm*1.e-3+71.8) );
+ } // end if
+ if( nanode > 3 ){
+ scl = TMath::Sqrt( (1.3*xm*1.e-3+49.)/(-0.57*xm*1.e-3+71.8) );
+ } // end if
+ // cout<<"1 microns: zm,dzz,xm,dxx,dxz,qq ="<<zm<<","<<dzz<<","
+ // <<xm<<","<<dxx<<","<<dxz<<","<<qq<<endl;
+ // old Boris.........
+ // tmp=29730. - 585.*fabs(xm/1000.);
+ // scl=TMath::Sqrt(tmp/130000.);
+
+ xm *= scl;
+ xx *= scl*scl;
+ xz *= scl;
+
+ dxx = xx - xm*xm;
+ // dzz = zz - zm*zm;
+ dxz = xz - xm*zm;
+ // cout<<"microns: zm,dzz,xm,dxx,xz,dxz,qq ="<<zm<<","<<dzz<<","
+ // <<xm<<","<<dxx<<","<<xz<<","<<dxz<<","<<qq<<endl;
+ // if(dzz < 7200.) dzz=7200.;//for one anode cluster dzz = anode**2/12
+
+ if (dxx < 0.) dxx=0.;
+ // the data if no cluster overlapping (the coordunates are in cm)
+ nfhits = 1;
+ xfit[0] = xm*1.e-4;
+ zfit[0] = zm*1.e-4;
+ qfit[0] = qq;
+ // if(nbins < 7) cout<<"**** nbins ="<<nbins<<endl;
+
+ if (nbins >= 7) {
+ if (dxz==0.) tga=0.;
+ else {
+ tmp=0.5*(dzz-dxx)/dxz;
+ tga = (dxz<0.) ? tmp-TMath::Sqrt(tmp*tmp+1) :
+ tmp+TMath::Sqrt(tmp*tmp+1);
+ } // end if dxz
+ elps=(tga*tga*dxx-2*tga*dxz+dzz)/(dxx+2*tga*dxz+tga*tga*dzz);
+ // change from microns to cm
+ xm *= 1.e-4;
+ zm *= 1.e-4;
+ zz *= 1.e-8;
+ xx *= 1.e-8;
+ xz *= 1.e-8;
+ dxz *= 1.e-8;
+ dxx *= 1.e-8;
+ dzz *= 1.e-8;
+ // cout<<"cm: zm,dzz,xm,dxx,xz,dxz,qq ="<<zm<<","<<dzz<<","
+ // <<xm<<","<<dxx<<","<<xz<<","<<dxz<<","<<qq<<endl;
+ for (i=0; i<nbins; i++) {
+ x[i] = x[i] *= scl;
+ x[i] = x[i] *= 1.e-4;
+ z[i] = z[i] *= 1.e-4;
+ } // end for i
+ // cout<<"!!! elps ="<<elps<<endl;
+ if (elps < 0.3) { // try to separate hits
+ separate = 1;
+ tmp=atan(tga);
+ Double_t cosa=cos(tmp),sina=sin(tmp);
+ Double_t a1=0., x1=0., xxx=0.;
+ for (i=0; i<nbins; i++) {
+ tmp=x[i]*cosa + z[i]*sina;
+ if (q[i] > a1) {
+ a1=q[i];
+ x1=tmp;
+ } // end if
+ xxx += tmp*tmp*tmp*q[i];
+ } // end for i
+ xxx /= qq;
+ Double_t z12=-sina*xm + cosa*zm;
+ sigma2=(sina*sina*xx-2*cosa*sina*xz+cosa*cosa*zz) - z12*z12;
+ xm=cosa*xm + sina*zm;
+ xx=cosa*cosa*xx + 2*cosa*sina*xz + sina*sina*zz;
+ Double_t x2=(xx - xm*x1 - sigma2)/(xm - x1);
+ Double_t r=a1*2*TMath::ACos(-1.)*sigma2/(qq*pitchx*pitchz);
+ for (i=0; i<33; i++) { // solve a system of equations
+ Double_t x1_old=x1, x2_old=x2, r_old=r;
+ Double_t c11=x1-x2;
+ Double_t c12=r;
+ Double_t c13=1-r;
+ Double_t c21=x1*x1 - x2*x2;
+ Double_t c22=2*r*x1;
+ Double_t c23=2*(1-r)*x2;
+ Double_t c31=3*sigma2*(x1-x2) + x1*x1*x1 - x2*x2*x2;
+ Double_t c32=3*r*(sigma2 + x1*x1);
+ Double_t c33=3*(1-r)*(sigma2 + x2*x2);
+ Double_t f1=-(r*x1 + (1-r)*x2 - xm);
+ Double_t f2=-(r*(sigma2+x1*x1)+(1-r)*(sigma2+x2*x2)- xx);
+ Double_t f3=-(r*x1*(3*sigma2+x1*x1)+(1-r)*x2*
+ (3*sigma2+x2*x2)-xxx);
+ Double_t d=c11*c22*c33+c21*c32*c13+c12*c23*c31-
+ c31*c22*c13 - c21*c12*c33 - c32*c23*c11;
+ if (d==0.) {
+ cout<<"*********** d=0 ***********\n";
+ break;
+ } // end if
+ Double_t dr=f1*c22*c33 + f2*c32*c13 + c12*c23*f3 -
+ f3*c22*c13 - f2*c12*c33 - c32*c23*f1;
+ Double_t d1=c11*f2*c33 + c21*f3*c13 + f1*c23*c31 -
+ c31*f2*c13 - c21*f1*c33 - f3*c23*c11;
+ Double_t d2=c11*c22*f3 + c21*c32*f1 + c12*f2*c31 -
+ c31*c22*f1 - c21*c12*f3 - c32*f2*c11;
+ r += dr/d;
+ x1 += d1/d;
+ x2 += d2/d;
+ if (fabs(x1-x1_old) > 0.0001) continue;
+ if (fabs(x2-x2_old) > 0.0001) continue;
+ if (fabs(r-r_old)/5 > 0.001) continue;
+ a1=r*qq*pitchx*pitchz/(2*TMath::ACos(-1.)*sigma2);
+ Double_t a2=a1*(1-r)/r;
+ qfit[0]=a1; xfit[0]=x1*cosa - z12*sina; zfit[0]=x1*sina +
+ z12*cosa;
+ qfit[1]=a2; xfit[1]=x2*cosa - z12*sina; zfit[1]=x2*sina +
+ z12*cosa;
+ nfhits=2;
+ break; // Ok !
+ } // end for i
+ if (i==33) cerr<<"No more iterations ! "<<endl;
+ } // end of attempt to separate overlapped clusters
+ } // end of nbins cut
+ if(elps < 0.) cout<<" elps=-1 ="<<elps<<endl;
+ if(elps >0. && elps< 0.3 && nfhits == 1) cout<<" small elps, nfh=1 ="
+ <<elps<<","<<nfhits<<endl;
+ if(nfhits == 2) cout<<" nfhits=2 ="<<nfhits<<endl;
+ for (i=0; i<nfhits; i++) {
+ xfit[i] *= (1.e+4/scl);
+ if(wing == 1) xfit[i] *= (-1);
+ zfit[i] *= 1.e+4;
+ // cout<<" --------- i,xfiti,zfiti,qfiti ="<<i<<","
+ // <<xfit[i]<<","<<zfit[i]<<","<<qfit[i]<<endl;
+ } // end for i
+ Int_t ncl = nfhits;
+ if(nfhits == 1 && separate == 1) {
+ cout<<"!!!!! no separate"<<endl;
+ ncl = -2;
+ } // end if
+ if(nfhits == 2) {
+ cout << "Split cluster: " << endl;
+ clusterJ->PrintInfo();
+ cout << " in: " << endl;
+ for (i=0; i<nfhits; i++) {
+ // AliITSRawClusterSDD *clust = new AliITSRawClusterSDD(wing,
+ -1,-1,(Float_t)qfit[i],ncl,0,0,
+ (Float_t)xfit[i],
+ (Float_t)zfit[i],0,0,0,0,
+ tstart,tstop,astart,astop);
+ // AliITSRawClusterSDD *clust = new AliITSRawClusterSDD(wing,-1,
+ // -1,(Float_t)qfit[i],0,0,0,
+ // (Float_t)xfit[i],
+ // (Float_t)zfit[i],0,0,0,0,
+ // tstart,tstop,astart,astop,ncl);
+ // ???????????
+ // if(wing == 1) xfit[i] *= (-1);
+ Float_t Anode = (zfit[i]/anodePitch+fNofAnodes/2-0.5);
+ Float_t Time = (fSddLength - xfit[i])/fDriftSpeed;
+ Float_t clusterPeakAmplitude = clusterJ->PeakAmpl();
+ Float_t peakpos = clusterJ->PeakPos();
+ Float_t clusteranodePath = (Anode - fNofAnodes/2)*anodePitch;
+ Float_t clusterDriftPath = Time*fDriftSpeed;
+ clusterDriftPath = fSddLength-clusterDriftPath;
+ AliITSRawClusterSDD *clust = new AliITSRawClusterSDD(wing,Anode,
+ Time,qfit[i],
+ clusterPeakAmplitude,peakpos,
+ 0.,0.,clusterDriftPath,
+ clusteranodePath,clusterJ->Samples()/2
+ ,tstart,tstop,0,0,0,astart,astop);
+ clust->PrintInfo();
+ iTS->AddCluster(1,clust);
+ // cout<<"new cluster added: tstart,tstop,astart,astop,x,ncl ="
+ // <<tstart<<","<<tstop<<","<<astart<<","<<astop<<","<<xfit[i]
+ // <<","<<ncl<<endl;
+ delete clust;
+ }// nfhits loop
+ fClusters->RemoveAt(j);
+ } // if nfhits = 2
+} // cluster loop
+fClusters->Compress();
+fMap->ClearMap();
+*/
+ return;
+}
+//______________________________________________________________________
+void AliITSClusterFinderSDD::GetRecPoints(){
+ // get rec points
+ static AliITS *iTS=(AliITS*)gAlice->GetModule("ITS");
+ // get number of clusters for this module
+ Int_t nofClusters = fClusters->GetEntriesFast();
+ nofClusters -= fNclusters;
+ const Float_t kconvGeV = 1.e-6; // GeV -> KeV
+ const Float_t kconv = 1.0e-4;
+ const Float_t kRMSx = 38.0*kconv; // microns->cm ITS TDR Table 1.3
+ const Float_t kRMSz = 28.0*kconv; // microns->cm ITS TDR Table 1.3
+ Int_t i;
+ Int_t ix, iz, idx=-1;
+ AliITSdigitSDD *dig=0;
+ Int_t ndigits=fDigits->GetEntriesFast();
+ for(i=0; i<nofClusters; i++) {
+ AliITSRawClusterSDD *clusterI = (AliITSRawClusterSDD*)fClusters->At(i);
+ if(!clusterI) Error("SDD: GetRecPoints","i clusterI ",i,clusterI);
+ if(clusterI) idx=clusterI->PeakPos();
+ if(idx>ndigits) Error("SDD: GetRecPoints","idx ndigits",idx,ndigits);
+ // try peak neighbours - to be done
+ if(idx&&idx<= ndigits) dig =(AliITSdigitSDD*)fDigits->UncheckedAt(idx);
+ if(!dig) {
+ // try cog
+ fSegmentation->GetPadIxz(clusterI->X(),clusterI->Z(),ix,iz);
+ dig = (AliITSdigitSDD*)fMap->GetHit(iz-1,ix-1);
+ // if null try neighbours
+ if (!dig) dig = (AliITSdigitSDD*)fMap->GetHit(iz-1,ix);
+ if (!dig) dig = (AliITSdigitSDD*)fMap->GetHit(iz-1,ix+1);
+ if (!dig) printf("SDD: cannot assign the track number!\n");
+ } // end if !dig
+ AliITSRecPoint rnew;
+ rnew.SetX(clusterI->X());
+ rnew.SetZ(clusterI->Z());
+ rnew.SetQ(clusterI->Q()); // in KeV - should be ADC
+ rnew.SetdEdX(kconvGeV*clusterI->Q());
+ rnew.SetSigmaX2(kRMSx*kRMSx);
+ rnew.SetSigmaZ2(kRMSz*kRMSz);