]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PMD/AliPMDtracker.cxx
defects fixed
[u/mrichter/AliRoot.git] / PMD / AliPMDtracker.cxx
index ca7c6ce100483db70fa2234a4e49dd290e43357f..27ef93aa3b5ec88786230a8372124cc28cd23109 100644 (file)
 #include <TNtuple.h>
 #include <TParticle.h>
 
+#include <TGeoMatrix.h>
+
+#include "AliGeomManager.h"
+
 #include "AliPMDcluster.h"
 #include "AliPMDclupid.h"
 #include "AliPMDrecpoint1.h"
+#include "AliPMDrecdata.h"
+#include "AliPMDrechit.h"
 #include "AliPMDUtility.h"
 #include "AliPMDDiscriminator.h"
 #include "AliPMDEmpDiscriminator.h"
@@ -49,6 +55,7 @@ ClassImp(AliPMDtracker)
 AliPMDtracker::AliPMDtracker():
   fTreeR(0),
   fRecpoints(new TClonesArray("AliPMDrecpoint1", 10)),
+  fRechits(new TClonesArray("AliPMDrechit", 10)),
   fPMDcontin(new TObjArray()),
   fPMDcontout(new TObjArray()),
   fPMDutil(new AliPMDUtility()),
@@ -71,6 +78,7 @@ AliPMDtracker:: AliPMDtracker(const AliPMDtracker & /* tracker */):
   TObject(/* tracker */),
   fTreeR(0),
   fRecpoints(NULL),
+  fRechits(NULL),
   fPMDcontin(NULL),
   fPMDcontout(NULL),
   fPMDutil(NULL),
@@ -104,6 +112,11 @@ AliPMDtracker::~AliPMDtracker()
     {
       fRecpoints->Clear();
     }
+  if (fRechits)
+    {
+      fRechits->Clear();
+    }
+
   if (fPMDcontin)
     {
       fPMDcontin->Delete();
@@ -133,9 +146,18 @@ void AliPMDtracker::Clusters2Tracks(AliESDEvent *event)
   // algorithm on CPV plane and PREshower plane
   //
 
-  Int_t   idet;
-  Int_t   ismn;
-  Float_t clusdata[6];
+  Int_t   idet = 0;
+  Int_t   ismn = 0;
+  Int_t   trackno = 1, trackpid = 0;
+  Float_t clusdata[6] = {0.,0.,0.,0.,0.,0.};
+  
+  Int_t *irow;
+  Int_t *icol;
+  Int_t *itra;
+  Int_t *ipid;
+  Float_t *cadc;
+
+  AliPMDrechit *rechit = 0x0;
 
   TBranch *branch = fTreeR->GetBranch("PMDRecpoint");
   if (!branch)
@@ -144,7 +166,16 @@ void AliPMDtracker::Clusters2Tracks(AliESDEvent *event)
       return;
     }
   branch->SetAddress(&fRecpoints);  
-  
+
+  TBranch *branch1 = fTreeR->GetBranch("PMDRechit");
+  if (!branch1)
+    {
+      AliError("PMDRechit branch not found");
+      return;
+    }
+  branch1->SetAddress(&fRechits);  
+
+  Int_t ncrhit = 0;
   Int_t   nmodules = (Int_t) branch->GetEntries();
   
   AliDebug(1,Form("Number of modules filled in treeR = %d",nmodules));
@@ -154,6 +185,7 @@ void AliPMDtracker::Clusters2Tracks(AliESDEvent *event)
       Int_t nentries = fRecpoints->GetLast();
       AliDebug(2,Form("Number of clusters per modules filled in treeR = %d"
                      ,nentries));
+
       for(Int_t ient = 0; ient < nentries+1; ient++)
        {
          fPMDrecpoint = (AliPMDrecpoint1*)fRecpoints->UncheckedAt(ient);
@@ -166,22 +198,96 @@ void AliPMDtracker::Clusters2Tracks(AliESDEvent *event)
          clusdata[4] = fPMDrecpoint->GetClusSigmaX();
          clusdata[5] = fPMDrecpoint->GetClusSigmaY();
 
-         fPMDclin = new AliPMDrecpoint1(idet,ismn,clusdata);
-         fPMDcontin->Add(fPMDclin);
+         if (clusdata[4] >= 0. && clusdata[5] >= 0.)
+           { 
+             // extract the associated cell information
+             branch1->GetEntry(ncrhit); 
+             Int_t nenbr1 = fRechits->GetLast() + 1;
+
+             irow = new Int_t[nenbr1];
+             icol = new Int_t[nenbr1];
+             itra = new Int_t[nenbr1];
+             ipid = new Int_t[nenbr1];
+             cadc = new Float_t[nenbr1];
+             
+             for (Int_t ient1 = 0; ient1 < nenbr1; ient1++)
+               {
+                 irow[ient1] = -99;
+                 icol[ient1] = -99;
+                 itra[ient1] = -99;
+                 ipid[ient1] = -99;
+                 cadc[ient1] = 0.;
+               }
+             for (Int_t ient1 = 0; ient1 < nenbr1; ient1++)
+               {
+                 rechit = (AliPMDrechit*)fRechits->UncheckedAt(ient1);
+                 //irow[ient1] = rechit->GetCellX();
+                 //icol[ient1] = rechit->GetCellY();
+                 itra[ient1] = rechit->GetCellTrack();
+                 ipid[ient1] = rechit->GetCellPid();
+                 cadc[ient1] = rechit->GetCellAdc();
+               }
+             if (idet == 0)
+               {
+                 AssignTrPidToCluster(nenbr1, itra, ipid, cadc,
+                                      trackno, trackpid);
+               }
+             else if (idet == 1)
+               {
+                 trackno  = itra[0];
+                 trackpid = ipid[0];
+               }
+
+             delete [] irow;
+             delete [] icol;
+             delete [] itra;
+             delete [] ipid;
+             delete [] cadc;
+
+             fPMDclin = new AliPMDrecdata(idet,ismn,trackno,trackpid,clusdata);
+             fPMDcontin->Add(fPMDclin);
+
+             ncrhit++;
+           }
        }
     }
 
-  AliPMDDiscriminator *pmddiscriminator = new AliPMDEmpDiscriminator();
-  pmddiscriminator->Discrimination(fPMDcontin,fPMDcontout);
+  AliPMDEmpDiscriminator pmddiscriminator;
+  pmddiscriminator.Discrimination(fPMDcontin,fPMDcontout);
+
+  // alignment implemention
+
+  Double_t sectr[4][3] = { {0.,0.,0.},{0.,0.,0.},{0.,0.,0.},{0.,0.,0.}};
+  TString snsector="PMD/Sector";
+  TString symname;
+  TGeoHMatrix gpmdor;
+  
+  for(Int_t isector=1; isector<=4; isector++)
+    {
+      symname = snsector;
+      symname += isector;
+      TGeoHMatrix *gpmdal = AliGeomManager::GetMatrix(symname);
+      Double_t *tral = gpmdal->GetTranslation();
+
+      AliGeomManager::GetOrigGlobalMatrix(symname, gpmdor);
+      Double_t *tror = gpmdor.GetTranslation();
+      
+      for(Int_t ixyz=0; ixyz<3; ixyz++)
+       {
+         sectr[isector-1][ixyz] = tral[ixyz] - tror[ixyz];
+       }
+    }
 
   const Float_t kzpos = 361.5;    // middle of the PMD
 
-  Int_t   det,smn;
-  Float_t xpos,ypos;
-  Float_t adc, ncell, rad;
+  Int_t   ix = -1, iy = -1;
+  Int_t   det = 0, smn = 0, trno = 1, trpid = 0, mstat = 0;
+  Float_t xpos = 0., ypos = 0.;
+  Float_t adc = 0., ncell = 0., radx = 0., rady = 0.;
   Float_t xglobal = 0., yglobal = 0., zglobal = 0;
-  Float_t pid;
+  Float_t pid = 0.;
 
+  fPMDutil->ApplyAlignment(sectr);
 
   Int_t nentries2 = fPMDcontout->GetEntries();
   AliDebug(1,Form("Number of clusters coming after discrimination = %d"
@@ -192,11 +298,36 @@ void AliPMDtracker::Clusters2Tracks(AliESDEvent *event)
       
       det   = fPMDclout->GetDetector();
       smn   = fPMDclout->GetSMN();
+      trno  = fPMDclout->GetClusTrackNo();
+      trpid = fPMDclout->GetClusTrackPid();
+      mstat = fPMDclout->GetClusMatching();
       xpos  = fPMDclout->GetClusX();
       ypos  = fPMDclout->GetClusY();
       adc   = fPMDclout->GetClusADC();
       ncell = fPMDclout->GetClusCells();
-      rad   = fPMDclout->GetClusRadius();
+      radx  = fPMDclout->GetClusSigmaX();
+      // Here in the variable "rady" we are keeping the row and col
+      // of the single isolated cluster having ncell = 1 for offline
+      // calibration
+      
+      if ((radx > 999. && radx < 1000.) && ncell == 1)
+       {
+         if (smn < 12)
+           {
+             ix = (Int_t) (ypos +0.5);
+             iy = (Int_t) xpos;
+           }
+         else if (smn > 12 && smn < 24)
+           {
+             ix = (Int_t) xpos;
+             iy = (Int_t) (ypos +0.5);
+           }
+         rady = (Float_t) (ix*100 + iy);
+       }
+      else
+       {
+         rady  = fPMDclout->GetClusSigmaY();
+       }
       pid   = fPMDclout->GetClusPID();
       
       //
@@ -213,35 +344,147 @@ void AliPMDtracker::Clusters2Tracks(AliESDEvent *event)
        **********************************************************************/
       //
 
-      fPMDutil->RectGeomCellPos(smn,xpos,ypos,xglobal,yglobal);
 
       if (det == 0)
        {
-         zglobal = kzpos + 1.6; // PREshower plane
+         zglobal = kzpos + 1.65; // PREshower plane
        }
       else if (det == 1)
        {
-         zglobal = kzpos - 1.7; // CPV plane
+         zglobal = kzpos - 1.65; // CPV plane
        }
 
+      fPMDutil->RectGeomCellPos(smn,xpos,ypos,xglobal,yglobal,zglobal);
+
       // Fill ESD
 
       AliESDPmdTrack *esdpmdtr = new  AliESDPmdTrack();
 
       esdpmdtr->SetDetector(det);
+      esdpmdtr->SetSmn(smn);
+      esdpmdtr->SetClusterTrackNo(trno);
+      esdpmdtr->SetClusterTrackPid(trpid);
+      esdpmdtr->SetClusterMatching(mstat);
+      
       esdpmdtr->SetClusterX(xglobal);
       esdpmdtr->SetClusterY(yglobal);
       esdpmdtr->SetClusterZ(zglobal);
       esdpmdtr->SetClusterADC(adc);
       esdpmdtr->SetClusterCells(ncell);
       esdpmdtr->SetClusterPID(pid);
+      esdpmdtr->SetClusterSigmaX(radx);
+      esdpmdtr->SetClusterSigmaY(rady);
 
       event->AddPmdTrack(esdpmdtr);
+      delete esdpmdtr;
     }
 
   fPMDcontin->Delete();
   fPMDcontout->Delete();
 
+}
+//--------------------------------------------------------------------//
+void AliPMDtracker::AssignTrPidToCluster(Int_t nentry, Int_t *itra,
+                                        Int_t *ipid, Float_t *cadc,
+                                        Int_t &trackno, Int_t &trackpid)
+{
+  // assign the track number and the corresponding pid to a cluster
+  // split cluster part will be done at the time of calculating eff/pur
+
+  Int_t *phentry = new Int_t [nentry];
+  Int_t *hadentry = new Int_t [nentry];
+  Int_t *trenergy = 0x0;
+  Int_t *trpid    = 0x0;
+  Int_t *sortcoord = 0x0;
+
+  Int_t ngtrack = 0;
+  Int_t nhtrack = 0;
+  for (Int_t i = 0; i < nentry; i++)
+    {
+      phentry[i] = -1;
+      hadentry[i] = -1;
+
+      if (ipid[i] == 22)
+       {
+         phentry[ngtrack] = i;
+         ngtrack++;
+       }
+      else if (ipid[i] != 22)
+       {
+         hadentry[nhtrack] = i;
+         nhtrack++;
+       }
+    }
+  
+  Int_t nghadtrack = ngtrack + nhtrack;
+
+  if (ngtrack == 0)
+    {
+      // hadron track
+      // no need of track number, set to -1
+      trackpid = 8;
+      trackno  = -1;
+    }
+  else if (ngtrack >= 1)
+    {
+      // one or more than one photon track + charged track
+      // find out which track deposits maximum energy and
+      // assign that track number and track pid
+
+      trenergy  = new Int_t [nghadtrack];
+      trpid     = new Int_t [nghadtrack];
+      sortcoord = new Int_t [nghadtrack];
+      for (Int_t i = 0; i < ngtrack; i++)
+       {
+         trenergy[i] = 0.;
+         trpid[i]    = -1;
+         for (Int_t j = 0; j < nentry; j++)
+           {
+             if (ipid[j] == 22 && itra[j] == itra[phentry[i]])
+               {
+                 trenergy[i] += cadc[j];
+                 trpid[i]     = 22;
+               }
+           }
+       }
+      for (Int_t i = ngtrack; i < nghadtrack; i++)
+       {
+         trenergy[i] = 0.;
+         trpid[i]    = -1;
+         for (Int_t j = 0; j < nentry; j++)
+           {
+             if (ipid[j] != 22 && itra[j] == itra[hadentry[i-ngtrack]])
+               {
+                 trenergy[i] += cadc[j];
+                 trpid[i]     = ipid[j];
+               }
+           }
+       }
+      
+      Bool_t jsort = true;
+      TMath::Sort(nghadtrack,trenergy,sortcoord,jsort);
+      
+      Int_t gtr = sortcoord[0];   
+      if (trpid[gtr] == 22)
+       {
+         trackpid = 22;
+         trackno  = itra[phentry[gtr]];   // highest adc track
+       }
+      else
+       {
+         trackpid = 8;
+         trackno = -1;
+       }
+      
+      delete [] trenergy;
+      delete [] trpid;
+      delete [] sortcoord;
+      
+    }   // end of ngtrack >= 1
+
+  delete [] phentry;
+  delete [] hadentry;
+  
 }
 //--------------------------------------------------------------------//
 void AliPMDtracker::SetVertex(Double_t vtx[3], Double_t evtx[3])