]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
New version from Boris and Enrico with improved ghost removal
authorbarbera <barbera@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 15 Jun 2001 16:04:50 +0000 (16:04 +0000)
committerbarbera <barbera@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 15 Jun 2001 16:04:50 +0000 (16:04 +0000)
ITS/AliITSClusterFinderSSD.cxx
ITS/SSD_ntuple.C
ITS/SSDrecpointTest.C

index 3798b7623b189501e7aaebb4945a84185baebffd..3086be4cfd8aa4ef123fe018da600c0662de3a72 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
+/**************************************************************************
+ * The package was revised and changed by Boris Batiounia in the time     *
+ * period of March - June 2001                                            *
+ **************************************************************************/
+
+
 #include <iostream.h>
 #include <TArrayI.h>
 #include "AliRun.h"
@@ -130,11 +136,14 @@ void AliITSClusterFinderSSD::FindRawClusters(Int_t module)
 // 3. If necesery, resolves for each group of neighbouring digits 
 //    how many clusters creates it.
 // 4. Colculate the x and z coordinate  
-
   Int_t lay, lad, detect;
   AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
   AliITSgeom *geom = aliITS->GetITSgeom();
+
            geom->GetModuleId(module,lay, lad, detect);
+
+           //cout<<"FindRawCl: layer,module ="<<lay<<","<<module<<endl; 
+
             if ( lay == 6 )
              ((AliITSsegmentationSSD*)fSegmentation)->SetLayer(6);
            if ( lay == 5 )
@@ -147,6 +156,7 @@ void AliITSClusterFinderSSD::FindRawClusters(Int_t module)
   FindNeighbouringDigits(); //ad. 2
   //SeparateOverlappedClusters();  //ad. 3
   ClustersToPackages();  //ad. 4
+  AliITSRecPoint rnew;
   fMap->ClearMap();
 }
 
@@ -592,7 +602,6 @@ void AliITSClusterFinderSSD::ClustersToPackages()
                                           // MB: well, that's not true that one
                                           //cannot sort objects in TClonesArray
   
-  //AliITSpackageSSD *currentpkg;
   AliITSclusterSSD *currentP;
   AliITSclusterSSD *currentN;
   Int_t j1, j2;    
@@ -600,31 +609,32 @@ void AliITSClusterFinderSSD::ClustersToPackages()
 //Fills in One Side Clusters Index Array
   FillClIndexArrays(oneSclP,oneSclN); 
 //Sorts filled Arrays
-  SortClusters(oneSclP,oneSclN);                   
+  //SortClusters(oneSclP,oneSclN);                   
 
   fNPackages=1;      
   new ((*fPackages)[0]) AliITSpackageSSD(fClusterP,fClusterN);
-  //currentpkg = (AliITSpackageSSD*)((*fPackages)[0]);
 
-  // Take all pairs of recpoints x coordinates in both sides  
+
+  //This part was includede by Boris Batiounia in March 2001.
+
+  // Take all recpoint pairs (x coordinates) in both P and N sides  
   // to calculate z coordinates of the recpoints
+
   for (j1=0;j1<fNClusterP;j1++) {  
       currentP = GetPSideCluster(oneSclP[j1]);
       Double_t xP = currentP->GetPosition();
-      //Int_t NxP = currentP->GetNumOfDigits();
       Float_t signalP = currentP->GetTotalSignal();
     for (j2=0;j2<fNClusterN;j2++) {  
       currentN = GetNSideCluster(oneSclN[j2]);
       Double_t xN = currentN->GetPosition();
-      //Int_t NxN = currentN->GetNumOfDigits();
       Float_t signalN = currentN->GetTotalSignal();
       CreateNewRecPoint(xP,1,xN,1,signalP,signalN,currentP, currentN, 0.75);
 
     }
   }
 
-   delete [] oneSclP;
-   delete [] oneSclN;
+   delete oneSclP;
+   delete oneSclN;
 
 }
 
@@ -658,6 +668,7 @@ CreateNewRecPoint(Float_t P, Float_t dP, Float_t N, Float_t dN,
      Float_t signal = 0;
      Float_t dedx = 0;
 
+     
      if(SigP>SigN) {
        signal = SigP;
        dedx = SigP*kADCtoKeV;
@@ -666,20 +677,14 @@ CreateNewRecPoint(Float_t P, Float_t dP, Float_t N, Float_t dN,
        dedx = SigN*kADCtoKeV;
      }        
 
-     Float_t signalCut = TMath::Abs((SigP-SigN)/signal);
-
-
      tr = (Int_t*) clusterP->GetTracks(n);
      NTracks = clusterP->GetNTracks();
-     //cout<<"NewRec: NTracks,tr0,tr1,tr2 ="<<NTracks<<","<<tr[0]<<","<<tr[1]<<","<<tr[2]<<endl;
-     if(NTracks>0 && signalCut<0.18) {
-     // the cut of the signals difference in P and N sides to
-     // remove the "ghosts"
+
        cnew.fSignalP=SigP;
        cnew.fSignalN=SigN;
        cnew.fMultiplicity=nstripsP;
        cnew.fMultiplicityN=nstripsN;
-       cnew.fQErr=signalCut;
+       cnew.fQErr=TMath::Abs(SigP-SigN);
        cnew.fNtracks=NTracks;
 
        fITS->AddCluster(2,&cnew);
@@ -698,7 +703,7 @@ CreateNewRecPoint(Float_t P, Float_t dP, Float_t N, Float_t dN,
        fITS->AddRecPoint(rnew);
 
        return kTRUE;
-     }
+       //}
   }
      return kFALSE;  
 }
@@ -770,6 +775,10 @@ Bool_t AliITSClusterFinderSSD::GetCrossing (Float_t &P, Float_t &N)
 { 
 // get crossing 
 
+  // This function was rivised and changed by Boris Batiounia in March 2001
+
+
+
   Float_t Dx = fSegmentation->Dx(); // detector size in x direction, microns
   Float_t Dz = fSegmentation->Dz(); // detector size in z direction, microns
 
@@ -786,11 +795,13 @@ Bool_t AliITSClusterFinderSSD::GetCrossing (Float_t &P, Float_t &N)
   fTanN=TMath::Tan(StereoN);
   Float_t kP = fTanP; // Tangent of 0.0075 mrad
   Float_t kN = fTanN; // Tangent of 0.0275 mrad
-
+  
     P *= fPitch;
     N *= fPitch; 
+
     xP = N;      // change the mistake for the P/N
     xN = P;      // coordinates correspondence in this function
+
     x = xP + kP*(Dz*kN-xP+xN)/(kP+kN);
     z = (Dz*kN-xP+xN)/(kP+kN); 
     xL = x - Dx/2;
@@ -799,6 +810,7 @@ Bool_t AliITSClusterFinderSSD::GetCrossing (Float_t &P, Float_t &N)
     N = zL;  
 
     if(TMath::Abs(xL) > Dx/2 || TMath::Abs(zL) > Dz/2) return kFALSE;
+    
     // Check that xL and zL are inside the detector for the 
     // correspondent xP and xN coordinates
 
index 1ef7d97efbfa82da0a82f3ab00f9669bc6edd6e1..f89e196c26a312c909b2bdb2b3ffbf827420f5be 100644 (file)
@@ -6,7 +6,7 @@ TFile *f = new TFile("SSD_his.root");
 //gStyle->SetOptStat(1111111);
 //gStyle->SetOptLogy();
 //TCanvas *c1 = new TCanvas("c1","SPD clusters",400,10,600,700);
-TCanvas *c2 = new TCanvas("c2","SPD clusters",400,10,600,700);
+TCanvas *c2 = new TCanvas("c2","SSD clusters",400,10,600,700);
 //c1->Divide(2,2);
 c2->Divide(2,2);
 
@@ -41,7 +41,46 @@ c2->Divide(2,2);
 // -------------------------------------------------------------------------
 
 
+/*
+c2->cd(1);
+gPad->SetFillColor(33);
+      ntuple->SetFillColor(42);
+      ntuple->Draw("dx","lay == 5 && hitprim == 1");
+c2->cd(2);
+gPad->SetFillColor(33);
+      ntuple->SetFillColor(46);
+      ntuple->Draw("dz","lay == 5 && hitprim == 1");
+c2->cd(3);
+gPad->SetFillColor(33);
+      ntuple->SetFillColor(42);
+      ntuple->Draw("dx","lay == 6 && hitprim == 1");
+c2->cd(4);
+gPad->SetFillColor(33);
+      ntuple->SetFillColor(46);
+      ntuple->Draw("dz","lay == 6 && hitprim == 1");
+*/
 
+/*
+c2->cd(1);
+gPad->SetFillColor(33);
+      ntuple->SetFillColor(42);
+      ntuple->Draw("dx","lay == 5 && hitprim == 1&&ntrover==1");
+c2->cd(2);
+gPad->SetFillColor(33);
+      ntuple->SetFillColor(46);
+      ntuple->Draw("dz","lay == 5 && hitprim == 1&&ntrover==1");
+c2->cd(3);
+gPad->SetFillColor(33);
+      ntuple->SetFillColor(42);
+      ntuple->Draw("dx","lay == 6 && hitprim == 1&&ntrover==1");
+c2->cd(4);
+gPad->SetFillColor(33);
+      ntuple->SetFillColor(46);
+      ntuple->Draw("dz","lay == 6 && hitprim == 1&&ntrover==1");
+*/
+
+
+/*
 c2->cd(1);
 gPad->SetFillColor(33);
       ntuple->SetFillColor(42);
@@ -58,7 +97,26 @@ c2->cd(4);
 gPad->SetFillColor(33);
       ntuple->SetFillColor(46);
       ntuple->Draw("dz","lay == 6 && hitprim == 1&&abs(dz)<5000");
+*/
 
+/*
+c2->cd(1);
+gPad->SetFillColor(33);
+      ntuple->SetFillColor(42);
+      ntuple->Draw("dx","lay == 5 && hitprim == 1&&abs(dx)<200&&ntrover==1");
+c2->cd(2);
+gPad->SetFillColor(33);
+      ntuple->SetFillColor(46);
+      ntuple->Draw("dz","lay == 5 && hitprim == 1&&abs(dz)<5000&&ntrover==1");
+c2->cd(3);
+gPad->SetFillColor(33);
+      ntuple->SetFillColor(42);
+      ntuple->Draw("dx","lay == 6 && hitprim == 1&&abs(dx)<200&&ntrover==1");
+c2->cd(4);
+gPad->SetFillColor(33);
+      ntuple->SetFillColor(46);
+      ntuple->Draw("dz","lay == 6 && hitprim == 1&&abs(dz)<5000&&ntrover==1");
+*/
 
 /*
 c2->cd(1);
@@ -79,6 +137,101 @@ gPad->SetFillColor(33);
       ntuple1->Draw("nxN","lay == 6 && noverprim>=0");
 */
 
+/*
+c2->cd(1);
+gPad->SetFillColor(33);
+      ntuple1->SetFillColor(42);
+      ntuple1->Draw("dx","lay == 5&&noverprim>=0&&abs(dx)<50&&nxP==2&&nxN==2");
+c2->cd(2);
+gPad->SetFillColor(33);
+      ntuple1->SetFillColor(46);
+      ntuple1->Draw("dz","lay == 5&&noverprim>=0&&abs(dz)<2000&&nxP==2&&nxN==2");
+c2->cd(3);
+gPad->SetFillColor(33);
+      ntuple1->SetFillColor(42);
+      ntuple1->Draw("dx","lay == 6 && noverprim>=0&&abs(dx)<200");
+c2->cd(4);
+gPad->SetFillColor(33);
+      ntuple1->SetFillColor(46);
+      ntuple1->Draw("dz","lay == 6 && noverprim>=0&&abs(dz)<5000");
+*/
+
+/*      
+c2->cd(1);
+gPad->SetFillColor(33);
+      ntuple1->SetFillColor(42);
+      ntuple1->Draw("dx","lay == 5&&noverprim>=0&&noverlaps==0&&qcut<0.18");
+c2->cd(2);
+gPad->SetFillColor(33);
+      ntuple1->SetFillColor(46);
+      ntuple1->Draw("dz","lay == 5&&noverprim>=0&&noverlaps==0&&qcut<0.18");
+c2->cd(3);
+gPad->SetFillColor(33);
+      ntuple1->SetFillColor(42);
+      ntuple1->Draw("dx","lay == 6 && noverprim>=0&&noverlaps==0&&qcut<0.18");
+c2->cd(4);
+gPad->SetFillColor(33);
+      ntuple1->SetFillColor(46);
+      ntuple1->Draw("dz","lay == 6 && noverprim>=0&&noverlaps==0&&qcut<0.18");
+*/      
+
+/*      
+c2->cd(1);
+gPad->SetFillColor(33);
+      ntuple1->SetFillColor(42);
+      ntuple1->Draw("dx","lay == 5&&noverprim>=0&&noverlaps==0&&abs(dx)<100&&qcut<0.18");
+c2->cd(2);
+gPad->SetFillColor(33);
+      ntuple1->SetFillColor(46);
+      ntuple1->Draw("dz","lay == 5&&noverprim>=0&&noverlaps==0&&abs(dz)<5000&&qcut<0.18");
+c2->cd(3);
+gPad->SetFillColor(33);
+      ntuple1->SetFillColor(42);
+      ntuple1->Draw("dx","lay == 6&&noverprim>=0&&noverlaps==0&&abs(dx)<100&&qcut<0.18");
+c2->cd(4);
+gPad->SetFillColor(33);
+      ntuple1->SetFillColor(46);
+      ntuple1->Draw("dz","lay == 6&&noverprim>=0&&noverlaps==0&&abs(dz)<5000&&qcut<0.18");
+*/      
+
+/*
+c2->cd(1);
+gPad->SetFillColor(33);
+      ntuple1->SetFillColor(42);
+      ntuple1->Draw("noverprim","lay == 5&&noverprim>=-1");
+c2->cd(2);
+gPad->SetFillColor(33);
+      ntuple1->SetFillColor(46);
+      ntuple1->Draw("noverlaps","lay == 5&&noverprim>=-1");
+c2->cd(3);
+gPad->SetFillColor(33);
+      ntuple1->SetFillColor(42);
+      ntuple1->Draw("noverprim","lay == 6 && noverprim>=-1");
+c2->cd(4);
+gPad->SetFillColor(33);
+      ntuple1->SetFillColor(46);
+      ntuple1->Draw("noverlaps","lay == 6 && noverprim>=-1");
+*/
+
+      /*
+c2->cd(1);
+gPad->SetFillColor(33);
+      ntuple1->SetFillColor(42);
+      ntuple1->Draw("ntrover","lay == 5&&noverprim>=0");
+c2->cd(2);
+gPad->SetFillColor(33);
+      ntuple1->SetFillColor(46);
+      ntuple1->Draw("noverlaps","lay == 5&&noverprim>=0");
+c2->cd(3);
+gPad->SetFillColor(33);
+      ntuple1->SetFillColor(42);
+      ntuple1->Draw("ntrover","lay == 6 && noverprim>=0");
+c2->cd(4);
+gPad->SetFillColor(33);
+      ntuple1->SetFillColor(46);
+      ntuple1->Draw("noverlaps","lay == 6 && noverprim>=0");
+      */
+
 /*
 c2->cd(1);
 gPad->SetFillColor(33);
@@ -101,15 +254,34 @@ gPad->SetFillColor(33);
 
 /////////////////////   Histogramm/ntuple  analysis  ////////////////////////
 
+                      
+c2->cd(1);
+gPad->SetFillColor(33);
+      adcPadcN5all->SetFillColor(42);
+      adcPadcN5all->Draw();
+c2->cd(2);
+gPad->SetFillColor(33);
+      adcPadcN6all->SetFillColor(42);
+      adcPadcN6all->Draw();
+c2->cd(3);
+gPad->SetFillColor(33);
+      adcPadcN5cut->SetFillColor(42);
+      adcPadcN5cut->Draw();
+c2->cd(4);
+gPad->SetFillColor(33);
+      adcPadcN6cut->SetFillColor(46);
+      adcPadcN6cut->Draw();
+                       
+
       /*
 c2->cd(1);
 gPad->SetFillColor(33);
       ntuple1->SetFillColor(42);
-      ntuple1->Draw("qclP","noverprim>=0");
+      ntuple1->Draw("qclP","noverprim>=0&&qclP<500");
 c2->cd(2);
 gPad->SetFillColor(33);
       ntuple1->SetFillColor(42);
-      ntuple1->Draw("qclN","noverprim>=0");
+      ntuple1->Draw("qclN","noverprim>=0&&qclN<500");
 c2->cd(3);
 gPad->SetFillColor(33);
       adcPadcN5cut->SetFillColor(42);
index f3e9294e420970b8103b30230f9b36d53f9cd370..4cc9bc091a4beb4196d989e0a0bbf5e4c823e155 100644 (file)
@@ -70,6 +70,7 @@ void SSDrecpointTest (Int_t evNumber1=0,Int_t evNumber2=0)
             Float_t qclP;
             Float_t qclN;
             Float_t qrec;
+           Float_t qcut;
             Float_t dx;
             Float_t dz;
           } ntuple1_st;
@@ -103,6 +104,7 @@ void SSDrecpointTest (Int_t evNumber1=0,Int_t evNumber2=0)
          ntuple1->Branch("qclP",&ntuple1_st.qclP,"qclP/F");
           ntuple1->Branch("qclN",&ntuple1_st.qclN,"qclN/F");
           ntuple1->Branch("qrec",&ntuple1_st.qrec,"qrec/F");
+          ntuple1->Branch("qcut",&ntuple1_st.qcut,"qcut/F");
           ntuple1->Branch("dx",&ntuple1_st.dx,"dx/F");
           ntuple1->Branch("dz",&ntuple1_st.dz,"dz/F");
          ntuple1->Branch("noverlaps",&ntuple1_st.noverlaps,"noverlaps/I");
@@ -237,14 +239,14 @@ void SSDrecpointTest (Int_t evNumber1=0,Int_t evNumber2=0)
      printf("det type %d first2, last2 %d %d \n",2,first2,last2);
 
      // module loop for the SSD
-     for (mod=first2; mod<last2+1; mod++) {  // for the "ALL" option
-       //for (mod=0; mod<last2-first2+1; mod++) { //for the "SSD" option
-
+     //for (mod=first2; mod<last2+1; mod++) {  // for the "ALL" option
+     for (mod=0; mod<last2-first2+1; mod++) { //for the "SSD" option
+       
        TTree *TR = gAlice->TreeR();
        Int_t nentrec=TR->GetEntries();
        //printf("Found %d entries in the RecPoints tree\n",nentrec);
-      
-              //cout << "CLUSTERS: reset" << endl;
+       
+       //cout << "CLUSTERS: reset" << endl;
        ITS->ResetClusters();
        //cout << "CLUSTERS: get" << endl;
        TC->GetEvent(mod);
@@ -262,24 +264,25 @@ void SSDrecpointTest (Int_t evNumber1=0,Int_t evNumber2=0)
        totclust += nclusters;
        //if (nclusters) printf("Found %d clusters for module %d\n",nrecc,mod);
        
-       //AliITSmodule *Mod = (AliITSmodule *)fITSmodules->At(mod+first2);
+       AliITSmodule *Mod = (AliITSmodule *)fITSmodules->At(mod+first2);
        // for the "SSD" option
 
-       AliITSmodule *Mod = (AliITSmodule *)fITSmodules->At(mod);
+       //AliITSmodule *Mod = (AliITSmodule *)fITSmodules->At(mod);
        // for the "ALL" option
 
        //       printf("Mod: %X\n",Mod);
        Int_t nhits = Mod->GetNhits();
        Float_t epart = 0;
-       cout <<" module,nrecp,nclusters,nhits ="<<mod<<","<<nrecp<<","<<nclusters<<","<<nhits<< endl;
+       //cout <<" module,nrecp,nclusters,nhits ="<<mod<<","<<nrecp<<","<<nclusters<<","<<nhits<< endl;
 
        // ---------------- cluster/hit analysis ---------------------
 
 
-     Float_t pathInSSD = 300.;
+       Float_t pathInSSD = 300.;
 
        // ---- Recpoint loop
        for (Int_t pnt=0;pnt<nrecp;pnt++) {
+
         itsPnt  = (AliITSRecPoint*)ITSrec->At(pnt);
         if(!itsPnt) continue;
         itsClu  = (AliITSRawClusterSSD*)ITSclu->At(pnt);
@@ -290,13 +293,11 @@ void SSDrecpointTest (Int_t evNumber1=0,Int_t evNumber2=0)
         Int_t ntrover = itsClu->fNtracks;
         Float_t qclP = itsClu->fSignalP;     // in ADC
         Float_t qclN = itsClu->fSignalN;     // in ADC
-        //Float_t dq = qclP - qclN;
-        Float_t qcut = itsClu->fQErr;        // abs(dq)/signal,
-                                             // where signal is
-                                             // max of qclP,qclN        
+        Float_t dq = TMath::Abs(qclP - qclN);    
         Float_t xrec = 10000*itsPnt->GetX();
         Float_t zrec = 10000*itsPnt->GetZ();
         Float_t qrec = itsPnt->GetQ();      // in ADC, maximum from fSignalP/N
+        Float_t qcut = dq/qrec;
         //Float_t dedx = itsPnt->GetdEdX();   // in KeV (ADC * 2.16)
         Float_t dedx = itsPnt->fdEdX;   // in KeV (ADC * 2.16)
          Int_t ii = 0;
@@ -305,44 +306,46 @@ void SSDrecpointTest (Int_t evNumber1=0,Int_t evNumber2=0)
         Int_t tr2 = itsPnt->GetLabel(ii);
          Int_t ii = 2;
         Int_t tr3 = itsPnt->GetLabel(ii);
-        // fill ntuple2
-            ntuple2_st.nxP = nxP;
-             ntuple2_st.nxN = nxN;
-            ntuple2_st.x = xrec/1000;
-             ntuple2_st.z = zrec/1000;
-
-             if(qcut < 0.18) ntuple2->Fill();
-
-
-         Int_t noverlaps = 0;
-         Int_t noverprim = 0;
-         Int_t flaghit = 0;
-          Float_t xhit0 = 1e+7;
-          Float_t yhit0 = 1e+7;
-          Float_t zhit0 = 1e+7;
-
-       // Hit loop
-        for (Int_t hit=0;hit<nhits;hit++) {
-
-        itsHit   = (AliITShit*)Mod->GetHit(hit);
-
-        Int_t flagtrack = 0;
-        Int_t hitlayer = itsHit->GetLayer();
-        Int_t hitladder= itsHit->GetLadder();
-        Int_t hitdet= itsHit->GetDetector();
-
-        Int_t track = itsHit->GetTrack();
-        Int_t dray = 0;
-        Int_t hitstat = itsHit->GetTrackStatus();
-
-         Float_t zhit = 10000*itsHit->GetZL();
-         Float_t xhit = 10000*itsHit->GetXL();
-         Float_t yhit = 10000*itsHit->GetYL();
-         Float_t ehit = 1.0e+6*itsHit->GetIonization(); // hit energy, KeV 
 
+        // fill ntuple2
+        ntuple2_st.nxP = nxP;
+        ntuple2_st.nxN = nxN;
+        ntuple2_st.x = xrec/1000;
+        ntuple2_st.z = zrec/1000;
+        
+        //if(qcut < 0.18) ntuple2->Fill();
+        
+        Int_t noverlaps = 0;
+        Int_t noverprim = 0;
+        Int_t flaghit = 0;
+        Float_t xhit0 = 1e+7;
+        Float_t yhit0 = 1e+7;
+        Float_t zhit0 = 1e+7;
+        Float_t dxprimbest = 1e+7;
+        Float_t dzprimbest = 1e+7;
+        
+        // Hit loop
+        for (Int_t hit=0;hit<nhits;hit++) {
+
+          itsHit   = (AliITShit*)Mod->GetHit(hit);
+          
+          Int_t flagtrack = 0;
+          Int_t hitlayer = itsHit->GetLayer();
+          Int_t hitladder= itsHit->GetLadder();
+          Int_t hitdet= itsHit->GetDetector();
+          
+          Int_t track = itsHit->GetTrack();
+          Int_t dray = 0;
+          Int_t hitstat = itsHit->GetTrackStatus();
+          
+          Float_t zhit = 10000*itsHit->GetZL();
+          Float_t xhit = 10000*itsHit->GetXL();
+          Float_t yhit = 10000*itsHit->GetYL();
+          Float_t ehit = 1.0e+6*itsHit->GetIonization(); // hit energy, KeV 
+          
           Int_t parent = itsHit->GetParticle()->GetFirstMother();
           Int_t partcode = itsHit->GetParticle()->GetPdgCode();
-
+          
    //  partcode (pdgCode): 11 - e-, 13 - mu-, 22 - gamma, 111 - pi0, 211 - i+
    //  310 - K0s, 321 - K+, 2112 - n, 2212 - p, 3122 - lambda
 
@@ -350,143 +353,156 @@ void SSDrecpointTest (Int_t evNumber1=0,Int_t evNumber2=0)
                                                      // vertex
           pmod *= 1.0e+3;
 
-         if(hitstat == 66 && yhit < -ylim) {  // entering hit
+         if(hitstat == 66 && yhit < -ylim) {
            xhit0 = xhit;
            yhit0 = yhit;
            zhit0 = zhit;
          }
 
          if(hitstat == 66) continue; // Take the not entering hits only 
-
+         
          if(xhit0 > 9e+6 || zhit0 > 9e+6 || yhit0 > 9e+6) {
-           //cout<<"default xhit0,zhit0,yhit0 ="<<xhit0<<","<<zhit0<<","<<yhit0<<endl;
            continue;
          }
-
-
-
+         
          // Consider the hits only with the track number equaled to one
          // of the recpoint
          if((track == tr1) || (track == tr2) || (track == tr3)) flagtrack = 1;
+         
+         if(flagtrack == 1) {     // the hit corresponds to the recpoint
+           
+           flaghit = 1;
+           
+           //Float_t px = itsHit->GetPXL(); // the momenta at this GEANT point
+           //Float_t py = itsHit->GetPYL();
+           //Float_t pz = itsHit->GetPZL();
+           
+           Int_t hitprim = 0;
+           
+           if(partcode == 11 && pmod < 6) dray = 1; // delta ray is e-
+           // at p < 6 MeV/c
+           
+           if((hitstat == 68 || hitstat == 33) && dray == 0)  noverlaps=noverlaps + 1;
+           // overlapps for all hits but
+           // not for delta ray which
+           // also went out from the
+           // detector and returned
+           // again
+           
+           
+           // x,z resolution colculation
+           if((hitstat == 68 || hitsat == 33) && dray == 0) {
+             Float_t xmed = (xhit + xhit0)/2;
+             Float_t zmed = (zhit + zhit0)/2;
+             Float_t xdif = xmed - xrec;
+             Float_t zdif = zmed - zrec;
+             
+             if(parent < 0)  {
+               hitprim = 1; // hitprim=1 for the primery particles
+               noverprim += 1;
+             }
+             pathInSSD = TMath::Sqrt((xhit0-xhit)*(xhit0-xhit)+(yhit0-yhit)*(yhit0-yhit)+(zhit0-zhit)*(zhit0-zhit));
+             
+             // Find the best xdif and zdif from any ones for the primery
+             // particles (to remove the wronge xdif and zdif if the hit
+             // belongs to the other package containing the same P/N cluster
+             
+             if(hitprim > 0) {
+               if(TMath::Abs(dxprimbest)>TMath::Abs(xdif)) dxprimbest = xdif;          
+               if(TMath::Abs(dzprimbest)>TMath::Abs(zdif)) dzprimbest = zdif;         
+             }
+             
+             // fill ntuple
+             ntuple_st.lay = hitlayer;
+             ntuple_st.nxP = nxP;
+             ntuple_st.nxN = nxN;
+             ntuple_st.hitprim = hitprim;
+             ntuple_st.partcode = partcode;
+             ntuple_st.ntrover = ntrover;
+             ntuple_st.x = xrec/1000;
+             ntuple_st.z = zrec/1000;
+             ntuple_st.dx = xdif;
+             ntuple_st.dz = zdif;
+             ntuple_st.pmod = pmod;
+             
+             //if(qcut < 0.18) ntuple->Fill();
+             ntuple->Fill();
+             
+             //if(hitlayer == 5 && qcut < 0.18) {
+             
+             if(hitlayer == 5 ) {
+               Xres5->Fill(xdif);
+               Zres5->Fill(zdif);
+               Path5->Fill(pathInSSD);
+             }
+             //if(hitlayer == 6 && qcut < 0.18) {
+             if(hitlayer == 6) {
+               Xres6->Fill(xdif);
+               Zres6->Fill(zdif);
+               Path6->Fill(pathInSSD);
+             }
+           } // hitstat 68/33
+         } else {       // non correspondent hit
+           xhit0 = 1e+7;
+           zhit0 = 1e+7;
+         } // end of hit-recpoint correspondence
+        } // hit loop       
+        
+        if(flaghit == 1) {
+          
+          if(noverlaps == 0) noverlaps = 1; // cluster contains one or more
+          // delta rays only
+   
+          
+          // fill ntuple1
+          ntuple1_st.lay = hitlayer;
+          ntuple1_st.lad = hitladder;
+          ntuple1_st.det = hitdet;
+          ntuple1_st.nxP = nxP;
+          ntuple1_st.nxN = nxN;
+          ntuple1_st.qclP = qclP*300/pathInSSD; 
+          ntuple1_st.qclN = qclN*300/pathInSSD; 
+          ntuple1_st.qrec = qrec*300/pathInSSD; 
+           ntuple1_st.qcut = qcut;
+          ntuple1_st.dx = dxprimbest;
+          ntuple1_st.dz = dzprimbest;
+          noverlaps -= 1;
+          noverprim -= 1;
+          ntuple1_st.noverlaps = noverlaps;
+          ntuple1_st.noverprim = noverprim;
+          ntuple1_st.ntrover = ntrover;
+          
+          //if(qcut < 0.18) ntuple1->Fill();
+          ntuple1->Fill();
+          
+          Float_t de = dedx*300./pathInSSD;
+          dEdX->Fill(de);
+          if(noverprim >=0) {
+            if(hitlayer == 5 ) {
+              adcPadcN5all->Fill(qclP,qclN);
+            }
+            if(hitlayer == 6 ) {
+              adcPadcN6all->Fill(qclP,qclN);
+            }
+            if(hitlayer == 5 && qcut < 0.18) {
+              //if(hitlayer == 5 && noverlaps == 0) {
+              adcPadcN5cut->Fill(qclP,qclN);
+              NxP5->Fill(nxP);
+              NxN5->Fill(nxN);
+            }
+            if(hitlayer == 6 && qcut < 0.18) {
+              //if(hitlayer == 6 && noverlaps == 0) {
+              adcPadcN6cut->Fill(qclP,qclN);
+              NxP6->Fill(nxP);
+              NxN6->Fill(nxN);
+            }
+          }
+        } // flaghit = 1
 
-         if(flagtrack == 1) {     // the hit corresponds to the recpoint
-
-          flaghit = 1;
-
-          //Float_t px = itsHit->GetPXL(); // the momenta at this GEANT point
-          //Float_t py = itsHit->GetPYL();
-          //Float_t pz = itsHit->GetPZL();
-
-         Int_t hitprim = 0;
-
-        if(partcode == 11 && pmod < 6) dray = 1; // delta ray is e-
-                                                 // at p < 6 MeV/c
-
-         if((hitstat == 68 || hitstat == 33) && dray == 0)  noverlaps=noverlaps + 1;
-                                                  // overlapps for all hits but
-                                                 // not for delta ray which
-                                                 // also went out from the
-                                                 // detector and returned
-                                                 // again
-
-
-         // x,z resolution colculation
-          if((hitstat == 68 || hitsat == 33) && dray == 0) {
-            Float_t xmed = (xhit + xhit0)/2;
-            Float_t zmed = (zhit + zhit0)/2;
-            Float_t xdif = xmed - xrec;
-            Float_t zdif = zmed - zrec;
-
-            if(parent < 0)  {
-             hitprim = 1; // hitprim=1 for the primery particles
-             noverprim += 1;
-           }
-            pathInSSD = TMath::Sqrt((xhit0-xhit)*(xhit0-xhit)+(yhit0-yhit)*(yhit0-yhit)+(zhit0-zhit)*(zhit0-zhit));
-
-            // fill ntuple
-             ntuple_st.lay = hitlayer;
-            ntuple_st.nxP = nxP;
-             ntuple_st.nxN = nxN;
-            ntuple_st.hitprim = hitprim;
-             ntuple_st.partcode = partcode;
-             ntuple_st.ntrover = ntrover;
-            ntuple_st.x = xrec/1000;
-             ntuple_st.z = zrec/1000;
-            ntuple_st.dx = xdif;
-             ntuple_st.dz = zdif;
-             ntuple_st.pmod = pmod;
-
-             //if(qcut < 0.18) ntuple->Fill();
-             ntuple->Fill();
-
-            //if(hitlayer == 5 && qcut < 0.18) {
-             
-           if(hitlayer == 5 ) {
-             Xres5->Fill(xdif);
-             Zres5->Fill(zdif);
-             Path5->Fill(pathInSSD);
-           }
-            //if(hitlayer == 6 && qcut < 0.18) {
-            if(hitlayer == 6) {
-             Xres6->Fill(xdif);
-             Zres6->Fill(zdif);
-             Path6->Fill(pathInSSD);
-            }
-         } // hitstat 68/33
-        } else {       // non correspondent hit
-         xhit0 = 1e+7;
-         zhit0 = 1e+7;
-        } // end of hit-recpoint correspondence
-       } // hit loop       
-
-       if(flaghit == 1) {
-
-         if(noverlaps == 0) noverlaps = 1; // cluster contains one or more
-         // delta rays only
-
-         // fill ntuple1
-         ntuple1_st.lay = hitlayer;
-         ntuple1_st.lad = hitladder;
-         ntuple1_st.det = hitdet;
-         ntuple1_st.nxP = nxP;
-         ntuple1_st.nxN = nxN;
-         ntuple1_st.qclP = qclP*300/pathInSSD; 
-         ntuple1_st.qclN = qclN*300/pathInSSD; 
-         ntuple1_st.qrec = qrec*300/pathInSSD; 
-         ntuple1_st.dx = xdif;
-         ntuple1_st.dz = zdif;
-         noverlaps -= 1;
-         noverprim -= 1;
-         ntuple1_st.noverlaps = noverlaps;
-         ntuple1_st.noverprim = noverprim;
-          ntuple1_st.ntrover = ntrover;
-
-         //if(qcut < 0.18) ntuple1->Fill();
-         ntuple1->Fill();
-
-          Float_t de = dedx*300./pathInSSD;
-          dEdX->Fill(de);
-           if(hitlayer == 5 ) {
-             adcPadcN5all->Fill(qclP,qclN);
-            }
-           if(hitlayer == 6 ) {
-             adcPadcN6all->Fill(qclP,qclN);
-            }
-           if(hitlayer == 5 && qcut < 0.18) {
-             adcPadcN5cut->Fill(qclP,qclN);
-             NxP5->Fill(nxP);
-             NxN5->Fill(nxN);
-            }
-           if(hitlayer == 6 && qcut < 0.18) {
-             adcPadcN6cut->Fill(qclP,qclN);
-             NxP6->Fill(nxP);
-             NxN6->Fill(nxN);
-            }
-       } // flaghit = 1
        } //b.b. recpoint loop
      } //b.b. module loop
    } //b.b. evnt loop
-
+   
    TFile fhistos("SSD_his.root","RECREATE");
 
    ntuple->Write();