[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 "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;
  Int_t   ismn;
  Int_t   trackno, trackpid;
  Float_t clusdata[6];
  
  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++)
		{
		  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);

  const Float_t kzpos = 361.5;    // middle of the PMD

  Int_t   ix = -1, iy = -1;
  Int_t   det,smn,trno,trpid,mstat;
  Float_t xpos,ypos;
  Float_t adc, ncell, radx, rady;
  Float_t xglobal = 0., yglobal = 0., zglobal = 0;
  Float_t pid;

  fPMDutil->ApplyAlignment();

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