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