[u/mrichter/AliRoot.git] / PMD / AliPMDtracker.cxx

/***************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

//-----------------------------------------------------//
//                                                     //
//           Date   : March 25 2004                    //
//  This reads the file PMD.RecPoints.root(TreeR),     //
//  calls the Clustering algorithm and stores the      //
//  clustering output in PMD.RecPoints.root(TreeR)     // 
//                                                     //
//-----------------------------------------------------//

#include <Riostream.h>
#include <TMath.h>
#include <TTree.h>
#include <TObjArray.h>
#include <TClonesArray.h>
#include <TFile.h>
#include <TBranch.h>
#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"
#include "AliPMDtracker.h"

#include "AliESDPmdTrack.h"
#include "AliESDEvent.h"
#include "AliLog.h"

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()),
  fPMDrecpoint(0),
  fPMDclin(0),
  fPMDclout(0),
  fXvertex(0.),
  fYvertex(0.),
  fZvertex(0.),
  fSigmaX(0.),
  fSigmaY(0.),
  fSigmaZ(0.)
{
  //
  // Default Constructor
  //
}
//--------------------------------------------------------------------//
AliPMDtracker:: AliPMDtracker(const AliPMDtracker & /* tracker */):
  TObject(/* tracker */),
  fTreeR(0),
  fRecpoints(NULL),
  fRechits(NULL),
  fPMDcontin(NULL),
  fPMDcontout(NULL),
  fPMDutil(NULL),
  fPMDrecpoint(0),
  fPMDclin(0),
  fPMDclout(0),
  fXvertex(0.),
  fYvertex(0.),
  fZvertex(0.),
  fSigmaX(0.),
  fSigmaY(0.),
  fSigmaZ(0.)
{
  // copy constructor
  AliError("Copy constructor not allowed");
}

//--------------------------------------------------------------------//
AliPMDtracker& AliPMDtracker::operator=(const AliPMDtracker & /* tracker */)
{
 // assignment operator
  AliError("Assignment operator not allowed");
  return *this;
}

//--------------------------------------------------------------------//
AliPMDtracker::~AliPMDtracker()
{
  // Destructor
  if (fRecpoints)
    {
      fRecpoints->Clear();
    }
  if (fRechits)
    {
      fRechits->Clear();
    }

  if (fPMDcontin)
    {
      fPMDcontin->Delete();
      delete fPMDcontin;
      fPMDcontin=0;
      
    }
  if (fPMDcontout)
  {
      fPMDcontout->Delete();
      delete fPMDcontout;
      fPMDcontout=0;

    }
  delete fPMDutil;
}
//--------------------------------------------------------------------//
void AliPMDtracker::LoadClusters(TTree *treein)
{
  // Load the Reconstructed tree
  fTreeR = treein;
}
//--------------------------------------------------------------------//
void AliPMDtracker::Clusters2Tracks(AliESDEvent *event)
{
  // Converts digits to recpoints after running clustering
  // algorithm on CPV plane and PREshower plane
  //

  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)
    {
      AliError("PMDRecpoint branch not found");
      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));
  for (Int_t imodule = 0; imodule < nmodules; imodule++)
    {
      branch->GetEntry(imodule); 
      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);
	  idet        = fPMDrecpoint->GetDetector();
	  ismn        = fPMDrecpoint->GetSMNumber();
	  clusdata[0] = fPMDrecpoint->GetClusX();
	  clusdata[1] = fPMDrecpoint->GetClusY();
	  clusdata[2] = fPMDrecpoint->GetClusADC();
	  clusdata[3] = fPMDrecpoint->GetClusCells();
	  clusdata[4] = fPMDrecpoint->GetClusSigmaX();
	  clusdata[5] = fPMDrecpoint->GetClusSigmaY();

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

  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   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 = 0.;

  fPMDutil->ApplyAlignment(sectr);

  Int_t nentries2 = fPMDcontout->GetEntries();
  AliDebug(1,Form("Number of clusters coming after discrimination = %d"
		  ,nentries2));
  for (Int_t ient1 = 0; ient1 < nentries2; ient1++)
    {
      fPMDclout = (AliPMDclupid*)fPMDcontout->UncheckedAt(ient1);
      
      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();
      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();
      
      //
      /**********************************************************************
       *    det   : Detector, 0: PRE & 1:CPV                                *
       *    smn   : Serial Module Number 0 to 23 for each plane             *
       *    xpos  : x-position of the cluster                               *
       *    ypos  : y-position of the cluster                               *
       *            THESE xpos & ypos are not the true xpos and ypos        *
       *            for some of the unit modules. They are rotated.         *
       *    adc   : ADC contained in the cluster                            *
       *    ncell : Number of cells contained in the cluster                *
       *    rad   : radius of the cluster (1d fit)                          *
       **********************************************************************/
      //


      if (det == 0)
	{
	  zglobal = kzpos + 1.65; // PREshower plane
	}
      else if (det == 1)
	{
	  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 *trpid    = 0x0;
  Int_t *sortcoord = 0x0;
  Float_t *trenergy = 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 Float_t [nghadtrack];
      trpid     = new Int_t [nghadtrack];
      sortcoord = new Int_t [2*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])
{
  fXvertex = vtx[0];
  fYvertex = vtx[1];
  fZvertex = vtx[2];
  fSigmaX  = evtx[0];
  fSigmaY  = evtx[1];
  fSigmaZ  = evtx[2];
}
//--------------------------------------------------------------------//
void AliPMDtracker::ResetClusters()
{
  if (fRecpoints) fRecpoints->Clear();
}
//--------------------------------------------------------------------//
Commit	Line	Data
	1	/***************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	//-----------------------------------------------------//
	17	// //
	18	// Date : March 25 2004 //
	19	// This reads the file PMD.RecPoints.root(TreeR), //
	20	// calls the Clustering algorithm and stores the //
	21	// clustering output in PMD.RecPoints.root(TreeR) //
	22	// //
	23	//-----------------------------------------------------//
	24
	25	#include <Riostream.h>
	26	#include <TMath.h>
	27	#include <TTree.h>
	28	#include <TObjArray.h>
	29	#include <TClonesArray.h>
	30	#include <TFile.h>
	31	#include <TBranch.h>
	32	#include <TNtuple.h>
	33	#include <TParticle.h>
	34
	35	#include <TGeoMatrix.h>
	36
	37	#include "AliGeomManager.h"
	38
	39	#include "AliPMDcluster.h"
	40	#include "AliPMDclupid.h"
	41	#include "AliPMDrecpoint1.h"
	42	#include "AliPMDrecdata.h"
	43	#include "AliPMDrechit.h"
	44	#include "AliPMDUtility.h"
	45	#include "AliPMDDiscriminator.h"
	46	#include "AliPMDEmpDiscriminator.h"
	47	#include "AliPMDtracker.h"
	48
	49	#include "AliESDPmdTrack.h"
	50	#include "AliESDEvent.h"
	51	#include "AliLog.h"
	52
	53	ClassImp(AliPMDtracker)
	54
	55	AliPMDtracker::AliPMDtracker():
	56	fTreeR(0),
	57	fRecpoints(new TClonesArray("AliPMDrecpoint1", 10)),
	58	fRechits(new TClonesArray("AliPMDrechit", 10)),
	59	fPMDcontin(new TObjArray()),
	60	fPMDcontout(new TObjArray()),
	61	fPMDutil(new AliPMDUtility()),
	62	fPMDrecpoint(0),
	63	fPMDclin(0),
	64	fPMDclout(0),
	65	fXvertex(0.),
	66	fYvertex(0.),
	67	fZvertex(0.),
	68	fSigmaX(0.),
	69	fSigmaY(0.),
	70	fSigmaZ(0.)
	71	{
	72	//
	73	// Default Constructor
	74	//
	75	}
	76	//--------------------------------------------------------------------//
	77	AliPMDtracker:: AliPMDtracker(const AliPMDtracker & /* tracker */):
	78	TObject(/* tracker */),
	79	fTreeR(0),
	80	fRecpoints(NULL),
	81	fRechits(NULL),
	82	fPMDcontin(NULL),
	83	fPMDcontout(NULL),
	84	fPMDutil(NULL),
	85	fPMDrecpoint(0),
	86	fPMDclin(0),
	87	fPMDclout(0),
	88	fXvertex(0.),
	89	fYvertex(0.),
	90	fZvertex(0.),
	91	fSigmaX(0.),
	92	fSigmaY(0.),
	93	fSigmaZ(0.)
	94	{
	95	// copy constructor
	96	AliError("Copy constructor not allowed");
	97	}
	98
	99	//--------------------------------------------------------------------//
	100	AliPMDtracker& AliPMDtracker::operator=(const AliPMDtracker & /* tracker */)
	101	{
	102	// assignment operator
	103	AliError("Assignment operator not allowed");
	104	return *this;
	105	}
	106
	107	//--------------------------------------------------------------------//
	108	AliPMDtracker::~AliPMDtracker()
	109	{
	110	// Destructor
	111	if (fRecpoints)
	112	{
	113	fRecpoints->Clear();
	114	}
	115	if (fRechits)
	116	{
	117	fRechits->Clear();
	118	}
	119
	120	if (fPMDcontin)
	121	{
	122	fPMDcontin->Delete();
	123	delete fPMDcontin;
	124	fPMDcontin=0;
	125
	126	}
	127	if (fPMDcontout)
	128	{
	129	fPMDcontout->Delete();
	130	delete fPMDcontout;
	131	fPMDcontout=0;
	132
	133	}
	134	delete fPMDutil;
	135	}
	136	//--------------------------------------------------------------------//
	137	void AliPMDtracker::LoadClusters(TTree *treein)
	138	{
	139	// Load the Reconstructed tree
	140	fTreeR = treein;
	141	}
	142	//--------------------------------------------------------------------//
	143	void AliPMDtracker::Clusters2Tracks(AliESDEvent *event)
	144	{
	145	// Converts digits to recpoints after running clustering
	146	// algorithm on CPV plane and PREshower plane
	147	//
	148
	149	Int_t idet = 0;
	150	Int_t ismn = 0;
	151	Int_t trackno = 1, trackpid = 0;
	152	Float_t clusdata[6] = {0.,0.,0.,0.,0.,0.};
	153
	154	Int_t *irow;
	155	Int_t *icol;
	156	Int_t *itra;
	157	Int_t *ipid;
	158	Float_t *cadc;
	159
	160	AliPMDrechit *rechit = 0x0;
	161
	162	TBranch *branch = fTreeR->GetBranch("PMDRecpoint");
	163	if (!branch)
	164	{
	165	AliError("PMDRecpoint branch not found");
	166	return;
	167	}
	168	branch->SetAddress(&fRecpoints);
	169
	170	TBranch *branch1 = fTreeR->GetBranch("PMDRechit");
	171	if (!branch1)
	172	{
	173	AliError("PMDRechit branch not found");
	174	return;
	175	}
	176	branch1->SetAddress(&fRechits);
	177
	178	Int_t ncrhit = 0;
	179	Int_t nmodules = (Int_t) branch->GetEntries();
	180
	181	AliDebug(1,Form("Number of modules filled in treeR = %d",nmodules));
	182	for (Int_t imodule = 0; imodule < nmodules; imodule++)
	183	{
	184	branch->GetEntry(imodule);
	185	Int_t nentries = fRecpoints->GetLast();
	186	AliDebug(2,Form("Number of clusters per modules filled in treeR = %d"
	187	,nentries));
	188
	189	for(Int_t ient = 0; ient < nentries+1; ient++)
	190	{
	191	fPMDrecpoint = (AliPMDrecpoint1*)fRecpoints->UncheckedAt(ient);
	192	idet = fPMDrecpoint->GetDetector();
	193	ismn = fPMDrecpoint->GetSMNumber();
	194	clusdata[0] = fPMDrecpoint->GetClusX();
	195	clusdata[1] = fPMDrecpoint->GetClusY();
	196	clusdata[2] = fPMDrecpoint->GetClusADC();
	197	clusdata[3] = fPMDrecpoint->GetClusCells();
	198	clusdata[4] = fPMDrecpoint->GetClusSigmaX();
	199	clusdata[5] = fPMDrecpoint->GetClusSigmaY();
	200
	201	if (clusdata[4] >= 0. && clusdata[5] >= 0.)
	202	{
	203	// extract the associated cell information
	204	branch1->GetEntry(ncrhit);
	205	Int_t nenbr1 = fRechits->GetLast() + 1;
	206
	207	irow = new Int_t[nenbr1];
	208	icol = new Int_t[nenbr1];
	209	itra = new Int_t[nenbr1];
	210	ipid = new Int_t[nenbr1];
	211	cadc = new Float_t[nenbr1];
	212
	213	for (Int_t ient1 = 0; ient1 < nenbr1; ient1++)
	214	{
	215	irow[ient1] = -99;
	216	icol[ient1] = -99;
	217	itra[ient1] = -99;
	218	ipid[ient1] = -99;
	219	cadc[ient1] = 0.;
	220	}
	221	for (Int_t ient1 = 0; ient1 < nenbr1; ient1++)
	222	{
	223	rechit = (AliPMDrechit*)fRechits->UncheckedAt(ient1);
	224	//irow[ient1] = rechit->GetCellX();
	225	//icol[ient1] = rechit->GetCellY();
	226	itra[ient1] = rechit->GetCellTrack();
	227	ipid[ient1] = rechit->GetCellPid();
	228	cadc[ient1] = rechit->GetCellAdc();
	229	}
	230	if (idet == 0)
	231	{
	232	AssignTrPidToCluster(nenbr1, itra, ipid, cadc,
	233	trackno, trackpid);
	234	}
	235	else if (idet == 1)
	236	{
	237	trackno = itra[0];
	238	trackpid = ipid[0];
	239	}
	240
	241	delete [] irow;
	242	delete [] icol;
	243	delete [] itra;
	244	delete [] ipid;
	245	delete [] cadc;
	246
	247	fPMDclin = new AliPMDrecdata(idet,ismn,trackno,trackpid,clusdata);
	248	fPMDcontin->Add(fPMDclin);
	249
	250	ncrhit++;
	251	}
	252	}
	253	}
	254
	255	AliPMDEmpDiscriminator pmddiscriminator;
	256	pmddiscriminator.Discrimination(fPMDcontin,fPMDcontout);
	257
	258	// alignment implemention
	259
	260	Double_t sectr[4][3] = { {0.,0.,0.},{0.,0.,0.},{0.,0.,0.},{0.,0.,0.}};
	261	TString snsector="PMD/Sector";
	262	TString symname;
	263	TGeoHMatrix gpmdor;
	264
	265	for(Int_t isector=1; isector<=4; isector++)
	266	{
	267	symname = snsector;
	268	symname += isector;
	269	TGeoHMatrix *gpmdal = AliGeomManager::GetMatrix(symname);
	270	Double_t *tral = gpmdal->GetTranslation();
	271
	272	AliGeomManager::GetOrigGlobalMatrix(symname, gpmdor);
	273	Double_t *tror = gpmdor.GetTranslation();
	274
	275	for(Int_t ixyz=0; ixyz<3; ixyz++)
	276	{
	277	sectr[isector-1][ixyz] = tral[ixyz] - tror[ixyz];
	278	}
	279	}
	280
	281	const Float_t kzpos = 361.5; // middle of the PMD
	282
	283	Int_t ix = -1, iy = -1;
	284	Int_t det = 0, smn = 0, trno = 1, trpid = 0, mstat = 0;
	285	Float_t xpos = 0., ypos = 0.;
	286	Float_t adc = 0., ncell = 0., radx = 0., rady = 0.;
	287	Float_t xglobal = 0., yglobal = 0., zglobal = 0;
	288	Float_t pid = 0.;
	289
	290	fPMDutil->ApplyAlignment(sectr);
	291
	292	Int_t nentries2 = fPMDcontout->GetEntries();
	293	AliDebug(1,Form("Number of clusters coming after discrimination = %d"
	294	,nentries2));
	295	for (Int_t ient1 = 0; ient1 < nentries2; ient1++)
	296	{
	297	fPMDclout = (AliPMDclupid*)fPMDcontout->UncheckedAt(ient1);
	298
	299	det = fPMDclout->GetDetector();
	300	smn = fPMDclout->GetSMN();
	301	trno = fPMDclout->GetClusTrackNo();
	302	trpid = fPMDclout->GetClusTrackPid();
	303	mstat = fPMDclout->GetClusMatching();
	304	xpos = fPMDclout->GetClusX();
	305	ypos = fPMDclout->GetClusY();
	306	adc = fPMDclout->GetClusADC();
	307	ncell = fPMDclout->GetClusCells();
	308	radx = fPMDclout->GetClusSigmaX();
	309	// Here in the variable "rady" we are keeping the row and col
	310	// of the single isolated cluster having ncell = 1 for offline
	311	// calibration
	312
	313	if ((radx > 999. && radx < 1000.) && ncell == 1)
	314	{
	315	if (smn < 12)
	316	{
	317	ix = (Int_t) (ypos +0.5);
	318	iy = (Int_t) xpos;
	319	}
	320	else if (smn > 12 && smn < 24)
	321	{
	322	ix = (Int_t) xpos;
	323	iy = (Int_t) (ypos +0.5);
	324	}
	325	rady = (Float_t) (ix*100 + iy);
	326	}
	327	else
	328	{
	329	rady = fPMDclout->GetClusSigmaY();
	330	}
	331	pid = fPMDclout->GetClusPID();
	332
	333	//
	334	/**********************************************************************
	335	* det : Detector, 0: PRE & 1:CPV *
	336	* smn : Serial Module Number 0 to 23 for each plane *
	337	* xpos : x-position of the cluster *
	338	* ypos : y-position of the cluster *
	339	* THESE xpos & ypos are not the true xpos and ypos *
	340	* for some of the unit modules. They are rotated. *
	341	* adc : ADC contained in the cluster *
	342	* ncell : Number of cells contained in the cluster *
	343	* rad : radius of the cluster (1d fit) *
	344	**********************************************************************/
	345	//
	346
	347
	348	if (det == 0)
	349	{
	350	zglobal = kzpos + 1.65; // PREshower plane
	351	}
	352	else if (det == 1)
	353	{
	354	zglobal = kzpos - 1.65; // CPV plane
	355	}
	356
	357	fPMDutil->RectGeomCellPos(smn,xpos,ypos,xglobal,yglobal,zglobal);
	358
	359	// Fill ESD
	360
	361	AliESDPmdTrack *esdpmdtr = new AliESDPmdTrack();
	362
	363	esdpmdtr->SetDetector(det);
	364	esdpmdtr->SetSmn(smn);
	365	esdpmdtr->SetClusterTrackNo(trno);
	366	esdpmdtr->SetClusterTrackPid(trpid);
	367	esdpmdtr->SetClusterMatching(mstat);
	368
	369	esdpmdtr->SetClusterX(xglobal);
	370	esdpmdtr->SetClusterY(yglobal);
	371	esdpmdtr->SetClusterZ(zglobal);
	372	esdpmdtr->SetClusterADC(adc);
	373	esdpmdtr->SetClusterCells(ncell);
	374	esdpmdtr->SetClusterPID(pid);
	375	esdpmdtr->SetClusterSigmaX(radx);
	376	esdpmdtr->SetClusterSigmaY(rady);
	377
	378	event->AddPmdTrack(esdpmdtr);
	379	delete esdpmdtr;
	380	}
	381
	382	fPMDcontin->Delete();
	383	fPMDcontout->Delete();
	384
	385	}
	386	//--------------------------------------------------------------------//
	387	void AliPMDtracker::AssignTrPidToCluster(Int_t nentry, Int_t *itra,
	388	Int_t ipid, Float_t cadc,
	389	Int_t &trackno, Int_t &trackpid)
	390	{
	391	// assign the track number and the corresponding pid to a cluster
	392	// split cluster part will be done at the time of calculating eff/pur
	393
	394	Int_t *phentry = new Int_t [nentry];
	395	Int_t *hadentry = new Int_t [nentry];
	396	Int_t *trpid = 0x0;
	397	Int_t *sortcoord = 0x0;
	398	Float_t *trenergy = 0x0;
	399
	400	Int_t ngtrack = 0;
	401	Int_t nhtrack = 0;
	402	for (Int_t i = 0; i < nentry; i++)
	403	{
	404	phentry[i] = -1;
	405	hadentry[i] = -1;
	406
	407	if (ipid[i] == 22)
	408	{
	409	phentry[ngtrack] = i;
	410	ngtrack++;
	411	}
	412	else if (ipid[i] != 22)
	413	{
	414	hadentry[nhtrack] = i;
	415	nhtrack++;
	416	}
	417	}
	418
	419	Int_t nghadtrack = ngtrack + nhtrack;
	420
	421	if (ngtrack == 0)
	422	{
	423	// hadron track
	424	// no need of track number, set to -1
	425	trackpid = 8;
	426	trackno = -1;
	427	}
	428	else if (ngtrack >= 1)
	429	{
	430	// one or more than one photon track + charged track
	431	// find out which track deposits maximum energy and
	432	// assign that track number and track pid
	433
	434	trenergy = new Float_t [nghadtrack];
	435	trpid = new Int_t [nghadtrack];
	436	sortcoord = new Int_t [2*nghadtrack];
	437	for (Int_t i = 0; i < ngtrack; i++)
	438	{
	439	trenergy[i] = 0.;
	440	trpid[i] = -1;
	441	for (Int_t j = 0; j < nentry; j++)
	442	{
	443	if (ipid[j] == 22 && itra[j] == itra[phentry[i]])
	444	{
	445	trenergy[i] += cadc[j];
	446	trpid[i] = 22;
	447	}
	448	}
	449	}
	450	for (Int_t i = ngtrack; i < nghadtrack; i++)
	451	{
	452	trenergy[i] = 0.;
	453	trpid[i] = -1;
	454	for (Int_t j = 0; j < nentry; j++)
	455	{
	456	if (ipid[j] != 22 && itra[j] == itra[hadentry[i-ngtrack]])
	457	{
	458	trenergy[i] += cadc[j];
	459	trpid[i] = ipid[j];
	460	}
	461	}
	462	}
	463
	464	Bool_t jsort = true;
	465	TMath::Sort(nghadtrack,trenergy,sortcoord,jsort);
	466
	467	Int_t gtr = sortcoord[0];
	468	if (trpid[gtr] == 22)
	469	{
	470	trackpid = 22;
	471	trackno = itra[phentry[gtr]]; // highest adc track
	472	}
	473	else
	474	{
	475	trackpid = 8;
	476	trackno = -1;
	477	}
	478
	479	delete [] trenergy;
	480	delete [] trpid;
	481	delete [] sortcoord;
	482
	483	} // end of ngtrack >= 1
	484
	485	delete [] phentry;
	486	delete [] hadentry;
	487
	488	}
	489	//--------------------------------------------------------------------//
	490	void AliPMDtracker::SetVertex(Double_t vtx[3], Double_t evtx[3])
	491	{
	492	fXvertex = vtx[0];
	493	fYvertex = vtx[1];
	494	fZvertex = vtx[2];
	495	fSigmaX = evtx[0];
	496	fSigmaY = evtx[1];
	497	fSigmaZ = evtx[2];
	498	}
	499	//--------------------------------------------------------------------//
	500	void AliPMDtracker::ResetClusters()
	501	{
	502	if (fRecpoints) fRecpoints->Clear();
	503	}
	504	//--------------------------------------------------------------------//