[u/mrichter/AliRoot.git] / PMD / AliPMDDiscriminator.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   : August 05 2003                   //
//  This reads the file PMD.digits.root(TreeD),        //
//  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 <TNtuple.h>
#include <TParticle.h>

#include "AliLog.h"
#include "AliPMDcluster.h"
#include "AliPMDclupid.h"
#include "AliPMDDiscriminator.h"

ClassImp(AliPMDDiscriminator)

AliPMDDiscriminator::AliPMDDiscriminator() :
  fDiscrim(0)
{
//
// Default Constructor
//

}
// -----------------------------------------------------------------------
AliPMDDiscriminator::~AliPMDDiscriminator()
{
  // Destructor
}

// -----------------------------------------------------------------------

void AliPMDDiscriminator::Discrimination(TObjArray *pmdcontin, TObjArray *pmdcontout)
{
  // Does Photon/Hadron discrimination

  if(fDiscrim == 0)
    {
      EmpDiscrimination(pmdcontin, pmdcontout);
    }
  else if(fDiscrim == 1)
    {
      NNDiscrimination();
    }
}
// -----------------------------------------------------------------------

void AliPMDDiscriminator::EmpDiscrimination(TObjArray *pmdcontin, TObjArray *pmdcontout)
{
  // Does Photon/Hadron discrimination
  // matching the clusters of CPV and PREshower plane
  //
  const  Int_t kumperdet = 24;
  static Int_t neibx[6]={1,0,-1,-1,0,1}, neiby[6]={0,1,1,0,-1,-1}; 
  Int_t   det,smn;
  Int_t   iprecount[24], icpvcount[24];
  Float_t xpos,ypos;
  Float_t adc, ncell, rad;
  Float_t clusdata[6];

  for(Int_t i = 0; i < kumperdet; i++)
    {
      iprecount[i] = 0;
      icpvcount[i] = 0;
    }
  AliPMDcluster  *pmdcl    = 0;
  AliPMDclupid   *pmdclout = 0;

  Int_t nentries1 = pmdcontin->GetEntries();

  AliDebug(1,Form("Number of total clusters from CPV PRE = %d",nentries1));
  for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
    {
      pmdcl = (AliPMDcluster*)pmdcontin->UncheckedAt(ient1);

      det   = pmdcl->GetDetector();
      smn   = pmdcl->GetSMN();
      if(det == 0) iprecount[smn]++;
      if(det == 1) icpvcount[smn]++;
    } // Entries of TObjArray loop

  Int_t   idet, ismn;
  Float_t edepcpv[48][96];
  Int_t   statuscpv[48][96];

  for(Int_t i = 0; i < kumperdet; i++) // unit module
    {
      // Initialisation of the ADC of CPV (1UM)
      for (Int_t ix = 0; ix < 48;ix++)
	{
	  for (Int_t iy = 0; iy < 96;iy++)
	    {
	      edepcpv[ix][iy]   = 0.;
	      statuscpv[ix][iy] = 0;
	    }
	}
      Int_t precounter   = iprecount[i];
      Int_t cpvcounter   = icpvcount[i];

      Float_t *xpadpre   = new Float_t[precounter];
      Float_t *ypadpre   = new Float_t[precounter];
      Float_t *adcpre    = new Float_t[precounter];
      Float_t *ncellpre  = new Float_t[precounter];
      Float_t *radpre    = new Float_t[precounter];
      Int_t   *sortcoord = new Int_t[precounter];
      Int_t   *clupidpre = new Int_t[precounter];

      Float_t *xpadcpv   = new Float_t[cpvcounter];
      Float_t *ypadcpv   = new Float_t[cpvcounter];
      Float_t *adccpv    = new Float_t[cpvcounter];
      Float_t *ncellcpv  = new Float_t[cpvcounter];
      Float_t *radcpv    = new Float_t[cpvcounter];

      Int_t ii = 0;
      Int_t ij = 0;
      for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
	{
	  pmdcl = (AliPMDcluster*)pmdcontin->UncheckedAt(ient1);
	  
	  det   = pmdcl->GetDetector();
	  smn   = pmdcl->GetSMN();
	  xpos  = pmdcl->GetClusX();
	  ypos  = pmdcl->GetClusY();
	  adc   = pmdcl->GetClusADC();
	  ncell = pmdcl->GetClusCells();
	  rad   = pmdcl->GetClusRadius();
	  
	  if(det == 0 && smn == i)
	    {
	      xpadpre[ii]  = xpos;
	      ypadpre[ii]  = ypos;
	      adcpre[ii]   = adc;
	      ncellpre[ii] = ncell;
	      radpre[ii]   = rad;
	      ii++;
	    }
	  if(det == 1 && smn == i)
	    {
	      Int_t ix = (Int_t) (xpos+0.5);
	      Int_t iy = (Int_t) (ypos+0.5);
	      if(ix > 47) ix = 47;
	      if(iy > 95) iy = 95;
	      edepcpv[ix][iy] = adc;
	      xpadcpv[ij]  = xpos;
	      ypadcpv[ij]  = ypos;
	      adccpv[ij]   = adc;
	      ncellcpv[ij] = ncell;
	      radcpv[ij]   = rad;
	      ij++;
	    }
	} // Entries of TObjArray loop
      // sorts from lowest ADC to highest ADC
      // and returns the coordinates
      Bool_t jsort = false;
      TMath::Sort(precounter,adcpre,sortcoord,jsort);

      Int_t jjsort = 0;
      for(Int_t jj=0; jj<precounter; jj++)
	{
	  // PRE information
	  // PIDs for PRE clusters are 0(photon) and 1(hadron)

	  jjsort = sortcoord[jj];

	  Int_t ix = (Int_t) (xpadpre[jjsort]+0.5);
	  Int_t iy = (Int_t) (ypadpre[jjsort]+0.5);
	  if(ix > 47) ix = 47;
	  if(iy > 95) iy = 95;

	  for(Int_t jk=0; jk<6; jk++)
	    {
	      Int_t jd1 = ix + neibx[jk]; 
	      Int_t jd2 = iy + neiby[jk];
	      if(jd1 <0 ) jd1 = 0;
	      if(jd1 >47) jd1 = 47;
	      if(jd2 <0 ) jd2 = 0;
	      if(jd2 >47) jd2 = 47;
	      if(edepcpv[jd1][jd2] > 0.0 && statuscpv[jd1][jd2] == 0)
		{
		  statuscpv[jd1][jd2] = 1;
		  clupidpre[jjsort]   = 1;
		  break;
		}
	    }

	  idet        = 0;
	  ismn        = i;
	  clusdata[0] = xpadpre[jjsort];
	  clusdata[1] = ypadpre[jjsort];
	  clusdata[2] = adcpre[jjsort];
	  clusdata[3] = ncellpre[jjsort];
	  clusdata[4] = radpre[jjsort];
	  //PH	  clusdata[5] = (Float_t) clupidpre[jjsort];

	  // Temporary the cluster PID is set to 1 if the
	  // adc > 3MIP units which will be changed later on.

	  if (adcpre[jjsort] > 4500.)
	    {
	      clusdata[5] = 1.0;
	    }
	  else
	    {
	      clusdata[5] = 0.0;
	    }
	  pmdclout = new AliPMDclupid(idet,ismn,clusdata);
	  pmdcontout->Add(pmdclout);
	} // jj loop

      for(Int_t jj=0; jj<cpvcounter; jj++)
	{
	  // CPV information
	  // PID for CPV clusters is 1

	  idet        = 1;
	  ismn        = i;
	  clusdata[0] = xpadcpv[jj];
	  clusdata[1] = ypadcpv[jj];
	  clusdata[2] = adccpv[jj];
	  clusdata[3] = ncellcpv[jj];
	  clusdata[4] = radcpv[jj];
	  clusdata[5] = 0.;

	  pmdclout = new AliPMDclupid(idet,ismn,clusdata);
	  pmdcontout->Add(pmdclout);
	}

      // delete all the pointers
      delete [] xpadpre;
      delete [] ypadpre;
      delete [] adcpre;
      delete [] ncellpre;
      delete [] radpre;
      delete [] clupidpre;
      delete [] xpadcpv;
      delete [] ypadcpv;
      delete [] adccpv;
      delete [] ncellcpv;
      delete [] radcpv;
    } // i loop

}
// -----------------------------------------------------------------------
void AliPMDDiscriminator::NNDiscrimination()
{
  // This method does discrimination using Neural Network technique

}
// -----------------------------------------------------------------------
void AliPMDDiscriminator::SetDiscrimination(Int_t idiscrim)
{
  fDiscrim = idiscrim;
}
// -----------------------------------------------------------------------
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 : August 05 2003 //
	19	// This reads the file PMD.digits.root(TreeD), //
	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 <TNtuple.h>
	32	#include <TParticle.h>
	33
ecee2a1a	34	#include "AliLog.h"
7dcaf913	35	#include "AliPMDcluster.h"
	36	#include "AliPMDclupid.h"
	37	#include "AliPMDDiscriminator.h"
	38
	39	ClassImp(AliPMDDiscriminator)
	40
	41	AliPMDDiscriminator::AliPMDDiscriminator() :
ecee2a1a	42	fDiscrim(0)
7dcaf913	43	{
	44	//
	45	// Default Constructor
	46	//
	47
	48	}
	49	// -----------------------------------------------------------------------
	50	AliPMDDiscriminator::~AliPMDDiscriminator()
	51	{
	52	// Destructor
	53	}
	54
	55	// -----------------------------------------------------------------------
	56
	57	void AliPMDDiscriminator::Discrimination(TObjArray pmdcontin, TObjArray pmdcontout)
	58	{
	59	// Does Photon/Hadron discrimination
	60
	61	if(fDiscrim == 0)
	62	{
	63	EmpDiscrimination(pmdcontin, pmdcontout);
	64	}
	65	else if(fDiscrim == 1)
	66	{
	67	NNDiscrimination();
	68	}
	69	}
	70	// -----------------------------------------------------------------------
	71
	72	void AliPMDDiscriminator::EmpDiscrimination(TObjArray pmdcontin, TObjArray pmdcontout)
	73	{
	74	// Does Photon/Hadron discrimination
	75	// matching the clusters of CPV and PREshower plane
	76	//
	77	const Int_t kumperdet = 24;
	78	static Int_t neibx[6]={1,0,-1,-1,0,1}, neiby[6]={0,1,1,0,-1,-1};
	79	Int_t det,smn;
	80	Int_t iprecount[24], icpvcount[24];
	81	Float_t xpos,ypos;
	82	Float_t adc, ncell, rad;
	83	Float_t clusdata[6];
	84
	85	for(Int_t i = 0; i < kumperdet; i++)
	86	{
	87	iprecount[i] = 0;
	88	icpvcount[i] = 0;
	89	}
	90	AliPMDcluster *pmdcl = 0;
	91	AliPMDclupid *pmdclout = 0;
	92
	93	Int_t nentries1 = pmdcontin->GetEntries();
ecee2a1a	94
ecee2a1a	95	AliDebug(1,Form("Number of total clusters from CPV PRE = %d",nentries1));
7dcaf913	96	for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
	97	{
	98	pmdcl = (AliPMDcluster*)pmdcontin->UncheckedAt(ient1);
	99
	100	det = pmdcl->GetDetector();
	101	smn = pmdcl->GetSMN();
	102	if(det == 0) iprecount[smn]++;
	103	if(det == 1) icpvcount[smn]++;
	104	} // Entries of TObjArray loop
	105
	106	Int_t idet, ismn;
	107	Float_t edepcpv[48][96];
	108	Int_t statuscpv[48][96];
	109
	110	for(Int_t i = 0; i < kumperdet; i++) // unit module
	111	{
	112	// Initialisation of the ADC of CPV (1UM)
	113	for (Int_t ix = 0; ix < 48;ix++)
	114	{
	115	for (Int_t iy = 0; iy < 96;iy++)
	116	{
	117	edepcpv[ix][iy] = 0.;
	118	statuscpv[ix][iy] = 0;
	119	}
	120	}
	121	Int_t precounter = iprecount[i];
	122	Int_t cpvcounter = icpvcount[i];
	123
	124	Float_t *xpadpre = new Float_t[precounter];
	125	Float_t *ypadpre = new Float_t[precounter];
	126	Float_t *adcpre = new Float_t[precounter];
	127	Float_t *ncellpre = new Float_t[precounter];
	128	Float_t *radpre = new Float_t[precounter];
	129	Int_t *sortcoord = new Int_t[precounter];
	130	Int_t *clupidpre = new Int_t[precounter];
	131
	132	Float_t *xpadcpv = new Float_t[cpvcounter];
	133	Float_t *ypadcpv = new Float_t[cpvcounter];
	134	Float_t *adccpv = new Float_t[cpvcounter];
	135	Float_t *ncellcpv = new Float_t[cpvcounter];
	136	Float_t *radcpv = new Float_t[cpvcounter];
	137
	138	Int_t ii = 0;
	139	Int_t ij = 0;
	140	for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
	141	{
	142	pmdcl = (AliPMDcluster*)pmdcontin->UncheckedAt(ient1);
	143
	144	det = pmdcl->GetDetector();
	145	smn = pmdcl->GetSMN();
	146	xpos = pmdcl->GetClusX();
	147	ypos = pmdcl->GetClusY();
	148	adc = pmdcl->GetClusADC();
	149	ncell = pmdcl->GetClusCells();
	150	rad = pmdcl->GetClusRadius();
	151
	152	if(det == 0 && smn == i)
	153	{
	154	xpadpre[ii] = xpos;
	155	ypadpre[ii] = ypos;
	156	adcpre[ii] = adc;
	157	ncellpre[ii] = ncell;
	158	radpre[ii] = rad;
	159	ii++;
160	}
161	if(det == 1 && smn == i)
162	{
163	Int_t ix = (Int_t) (xpos+0.5);
164	Int_t iy = (Int_t) (ypos+0.5);
165	if(ix > 47) ix = 47;
166	if(iy > 95) iy = 95;
167	edepcpv[ix][iy] = adc;
168	xpadcpv[ij] = xpos;
169	ypadcpv[ij] = ypos;
170	adccpv[ij] = adc;
171	ncellcpv[ij] = ncell;
172	radcpv[ij] = rad;
173	ij++;
174	}
175	} // Entries of TObjArray loop
176	// sorts from lowest ADC to highest ADC
177	// and returns the coordinates
178	Bool_t jsort = false;
179	TMath::Sort(precounter,adcpre,sortcoord,jsort);
180
181	Int_t jjsort = 0;
182	for(Int_t jj=0; jj<precounter; jj++)
183	{
184	// PRE information
185	// PIDs for PRE clusters are 0(photon) and 1(hadron)
186
187	jjsort = sortcoord[jj];
188
189	Int_t ix = (Int_t) (xpadpre[jjsort]+0.5);
190	Int_t iy = (Int_t) (ypadpre[jjsort]+0.5);
191	if(ix > 47) ix = 47;
192	if(iy > 95) iy = 95;
193
194	for(Int_t jk=0; jk<6; jk++)
195	{
196	Int_t jd1 = ix + neibx[jk];
197	Int_t jd2 = iy + neiby[jk];
198	if(jd1 <0 ) jd1 = 0;
199	if(jd1 >47) jd1 = 47;
200	if(jd2 <0 ) jd2 = 0;
201	if(jd2 >47) jd2 = 47;
202	if(edepcpv[jd1][jd2] > 0.0 && statuscpv[jd1][jd2] == 0)
203	{
204	statuscpv[jd1][jd2] = 1;
205	clupidpre[jjsort] = 1;
206	break;
207	}
208	}
209
210	idet = 0;
211	ismn = i;
212	clusdata[0] = xpadpre[jjsort];
213	clusdata[1] = ypadpre[jjsort];
214	clusdata[2] = adcpre[jjsort];
215	clusdata[3] = ncellpre[jjsort];
216	clusdata[4] = radpre[jjsort];
278747fb	217	//PH clusdata[5] = (Float_t) clupidpre[jjsort];
7dcaf913	218
	219	// Temporary the cluster PID is set to 1 if the
	220	// adc > 3MIP units which will be changed later on.
	221
	222	if (adcpre[jjsort] > 4500.)
	223	{
	224	clusdata[5] = 1.0;
	225	}
	226	else
	227	{
	228	clusdata[5] = 0.0;
	229	}
	230	pmdclout = new AliPMDclupid(idet,ismn,clusdata);
	231	pmdcontout->Add(pmdclout);
	232	} // jj loop
	233
	234	for(Int_t jj=0; jj<cpvcounter; jj++)
	235	{
	236	// CPV information
	237	// PID for CPV clusters is 1
	238
	239	idet = 1;
	240	ismn = i;
	241	clusdata[0] = xpadcpv[jj];
	242	clusdata[1] = ypadcpv[jj];
	243	clusdata[2] = adccpv[jj];
	244	clusdata[3] = ncellcpv[jj];
	245	clusdata[4] = radcpv[jj];
	246	clusdata[5] = 0.;
	247
	248	pmdclout = new AliPMDclupid(idet,ismn,clusdata);
	249	pmdcontout->Add(pmdclout);
	250	}
	251
	252	// delete all the pointers
	253	delete [] xpadpre;
	254	delete [] ypadpre;
	255	delete [] adcpre;
	256	delete [] ncellpre;
	257	delete [] radpre;
	258	delete [] clupidpre;
	259	delete [] xpadcpv;
	260	delete [] ypadcpv;
	261	delete [] adccpv;
	262	delete [] ncellcpv;
	263	delete [] radcpv;
	264	} // i loop
	265
	266	}
	267	// -----------------------------------------------------------------------
	268	void AliPMDDiscriminator::NNDiscrimination()
	269	{
	270	// This method does discrimination using Neural Network technique
	271
	272	}
	273	// -----------------------------------------------------------------------
	274	void AliPMDDiscriminator::SetDiscrimination(Int_t idiscrim)
	275	{
	276	fDiscrim = idiscrim;
	277	}
	278	// -----------------------------------------------------------------------
	279