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