method Raw2SDigits added
[u/mrichter/AliRoot.git] / PMD / AliPMDDiscriminator.cxx
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
37 #include "AliLog.h"
38 #include "AliPMDcluster.h"
39 #include "AliPMDclupid.h"
40 #include "AliPMDDiscriminator.h"
41
42 ClassImp(AliPMDDiscriminator)
43
44 AliPMDDiscriminator::AliPMDDiscriminator() :
45   fDiscrim(0)
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();
97
98   AliDebug(1,Form("Number of total clusters from CPV PRE = %d",nentries1));
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];
220           //PH    clusdata[5] = (Float_t) clupidpre[jjsort];
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