Version where the process for HLT and vdrift on DAQ are off(Raphaelle)
[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 <TTree.h>
28 #include <TObjArray.h>
29 #include <TClonesArray.h>
30 #include <TFile.h>
31 #include <TNtuple.h>
32 #include <TParticle.h>
33
34 #include "AliLog.h"
35 #include "AliPMDcluster.h"
36 #include "AliPMDclupid.h"
37 #include "AliPMDDiscriminator.h"
38
39 ClassImp(AliPMDDiscriminator)
40
41 AliPMDDiscriminator::AliPMDDiscriminator() :
42   fDiscrim(0)
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 = 0, smn = 0;
80   Int_t   iprecount[24], icpvcount[24];
81   Float_t xpos = 0., ypos = 0.;
82   Float_t adc = 0., ncell = 0., rad = 0.;
83   Float_t clusdata[6] = {0.,0.,0.,0.,0.,0.};
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();
94
95   AliDebug(1,Form("Number of total clusters from CPV PRE = %d",nentries1));
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 = 0, ismn = 0;
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];
217           //PH    clusdata[5] = (Float_t) clupidpre[jjsort];
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