- 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;
-}