724ac2b0c9b49eb5f000f258714d773f771979cc
[u/mrichter/AliRoot.git] / PMD / AliPMDClusterFinder.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 "AliRun.h"
38 #include "AliPMD.h"
39 #include "AliDetector.h"
40 #include "AliRunLoader.h"
41 #include "AliLoader.h"
42 #include "AliHeader.h"
43 #include "AliRawReader.h"
44
45 #include "AliPMDdigit.h"
46 #include "AliPMDClusterFinder.h"
47 #include "AliPMDClustering.h"
48 #include "AliPMDcluster.h"
49 #include "AliPMDrecpoint1.h"
50 #include "AliPMDRawStream.h"
51
52 ClassImp(AliPMDClusterFinder)
53
54 AliPMDClusterFinder::AliPMDClusterFinder(AliRunLoader* runLoader):
55   fRunLoader(runLoader),
56   fPMDLoader(runLoader->GetLoader("PMDLoader")),
57   fTreeD(0),
58   fTreeR(0),
59   fDigits(new TClonesArray("AliPMDdigit", 1000)),
60   fRecpoints(new TClonesArray("AliPMDrecpoint1", 1000)),
61   fNpoint(0),
62   fDebug(0),
63   fEcut(0.)
64 {
65 //
66 // Constructor
67 //
68 }
69 // ------------------------------------------------------------------------- //
70 AliPMDClusterFinder::~AliPMDClusterFinder()
71 {
72   // Destructor
73   if (fDigits)
74     {
75       fDigits->Delete();
76       delete fDigits;
77       fDigits=0;
78     }
79   if (fRecpoints)
80     {
81       fRecpoints->Delete();
82       delete fRecpoints;
83       fRecpoints=0;
84     }
85 }
86 // ------------------------------------------------------------------------- //
87
88 void AliPMDClusterFinder::Digits2RecPoints(Int_t ievt)
89 {
90   // Converts digits to recpoints after running clustering
91   // algorithm on CPV plane and PREshower plane
92   //
93   Int_t    det = 0,smn = 0;
94   Int_t    xpos,ypos;
95   Float_t  adc;
96   Int_t    ismn;
97   Int_t    idet;
98   Float_t  clusdata[5];
99
100   TObjArray *pmdcont = new TObjArray();
101   AliPMDClustering *pmdclust = new AliPMDClustering();
102   pmdclust->SetDebug(fDebug);
103   pmdclust->SetEdepCut(fEcut);
104
105   fRunLoader->GetEvent(ievt);
106   //cout << " ***** Beginning::Digits2RecPoints *****" << endl;
107   fTreeD = fPMDLoader->TreeD();
108   if (fTreeD == 0x0)
109     {
110       cout << " Can not get TreeD" << endl;
111     }
112   AliPMDdigit  *pmddigit;
113   TBranch *branch = fTreeD->GetBranch("PMDDigit");
114   branch->SetAddress(&fDigits);
115
116   ResetRecpoint();
117   fTreeR = fPMDLoader->TreeR();
118   if (fTreeR == 0x0)
119     {
120       fPMDLoader->MakeTree("R");
121       fTreeR = fPMDLoader->TreeR();
122     }
123
124   Int_t bufsize = 16000;
125   fTreeR->Branch("PMDRecpoint", &fRecpoints, bufsize); 
126
127   Int_t nmodules = (Int_t) fTreeD->GetEntries();
128   
129   for (Int_t imodule = 0; imodule < nmodules; imodule++)
130     {
131       ResetCellADC();
132       fTreeD->GetEntry(imodule); 
133       Int_t nentries = fDigits->GetLast();
134       for (Int_t ient = 0; ient < nentries+1; ient++)
135         {
136           pmddigit = (AliPMDdigit*)fDigits->UncheckedAt(ient);
137           
138           det    = pmddigit->GetDetector();
139           smn    = pmddigit->GetSMNumber();
140           xpos   = pmddigit->GetRow();
141           ypos   = pmddigit->GetColumn();
142           adc    = pmddigit->GetADC();
143           //Int_t trno   = pmddigit->GetTrackNumber();
144           fCellADC[xpos][ypos] = (Double_t) adc;
145         }
146
147       idet = det;
148       ismn = smn;
149       //      pmdclust->DoClust(idet,ismn,fCellADC,pmdcont);
150       
151       Int_t nentries1 = pmdcont->GetEntries();
152 //      cout << " nentries1 = " << nentries1 << endl;
153       for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
154         {
155           AliPMDcluster *pmdcl = (AliPMDcluster*)pmdcont->UncheckedAt(ient1);
156           idet        = pmdcl->GetDetector();
157           ismn        = pmdcl->GetSMN();
158           clusdata[0] = pmdcl->GetClusX();
159           clusdata[1] = pmdcl->GetClusY();
160           clusdata[2] = pmdcl->GetClusADC();
161           clusdata[3] = pmdcl->GetClusCells();
162           clusdata[4] = pmdcl->GetClusRadius();
163           
164           AddRecPoint(idet,ismn,clusdata);
165         }
166       pmdcont->Clear();
167       
168       fTreeR->Fill();
169       ResetRecpoint();
170
171     } // modules
172
173   ResetCellADC();
174   fPMDLoader = fRunLoader->GetLoader("PMDLoader");  
175   fPMDLoader->WriteRecPoints("OVERWRITE");
176
177   //   delete the pointers
178   delete pmdclust;
179   delete pmdcont;
180     
181   //  cout << " ***** End::Digits2RecPoints *****" << endl;
182 }
183 // ------------------------------------------------------------------------- //
184
185 void AliPMDClusterFinder::Digits2RecPoints(Int_t ievt, AliRawReader *rawReader)
186 {
187   // Converts digits to recpoints after running clustering
188   // algorithm on CPV plane and PREshower plane
189   //
190
191   Float_t  clusdata[5];
192
193
194   rawReader->Reset();
195   AliPMDRawStream pmdinput(rawReader);
196
197   TObjArray *pmdcont = new TObjArray();
198   AliPMDClustering *pmdclust = new AliPMDClustering();
199   pmdclust->SetDebug(fDebug);
200   pmdclust->SetEdepCut(fEcut);
201
202   fRunLoader->GetEvent(ievt);
203
204   ResetRecpoint();
205   fTreeR = fPMDLoader->TreeR();
206   if (fTreeR == 0x0)
207     {
208       fPMDLoader->MakeTree("R");
209       fTreeR = fPMDLoader->TreeR();
210     }
211
212   Int_t bufsize = 16000;
213   fTreeR->Branch("PMDRecpoint", &fRecpoints, bufsize); 
214
215   Int_t irownew  = 0;
216   Int_t icolnew  = 0;
217   Int_t irownew1 = 0;
218   Int_t icolnew1 = 0;
219   const Int_t kDet = 2;
220   const Int_t kSMN = 24;
221   const Int_t kRow = 48;
222   const Int_t kCol = 96;
223   Int_t preADC[kSMN][kRow][kCol];
224   Int_t cpvADC[kSMN][kRow][kCol];
225
226   for (Int_t i = 0; i < kSMN; i++)
227     {
228       for (Int_t j = 0; j < kRow; j++)
229         {
230           for (Int_t k = 0; k < kCol; k++)
231             {
232               preADC[i][j][k] = 0;
233               cpvADC[i][j][k] = 0;
234             }
235         }
236     }
237   ResetCellADC();
238
239   while(pmdinput.Next())
240     {
241       Int_t det = pmdinput.GetDetector();
242       Int_t smn = pmdinput.GetSMN();
243       //Int_t mcm = pmdinput.GetMCM();
244       //Int_t chno = pmdinput.GetChannel();
245       Int_t row = pmdinput.GetRow();
246       Int_t col = pmdinput.GetColumn();
247       Int_t sig = pmdinput.GetSignal();
248
249       //cout << " det = " << det << " smn = " << smn 
250       //   << " row = " << row << " col = " << col
251       //   << " sig = " << sig << endl;
252
253       // Transform all the (0,0) coordinates to the geant frame
254       if(smn < 6)
255         {
256           irownew1 = 95 - row;
257           icolnew1 = col;
258         }
259       else if(smn >= 6 && smn < 12)
260         {
261           irownew1 = row;
262           icolnew1 = 47 - col;
263         }
264       else if(smn >= 12 && smn < 18)
265         {
266           irownew1 = 47 - row;
267           icolnew1 = col;
268         }
269       else if(smn >= 18 && smn < 24)
270         {
271           irownew1 = row;
272           icolnew1 = 95 - col;
273         }
274
275       // for smn < 12          : row = 96, column = 48
276       // for smn>= 12 and < 24 : row = 48, column = 96
277       // In order to make it uniform dimension, smn < 12 are inverted
278       // i.i., row becomes column and column becomes row
279       // for others it remains same
280       // This is further inverted back while calculating eta and phi
281       if(smn < 12)
282         {
283           // SupeModule 1 and 2 : Rows are inverted to columns and vice versa
284           // and at the time of calculating the eta,phi it is again reverted
285           // back
286           irownew = icolnew1;
287           icolnew = irownew1;
288         }
289       else if( smn >= 12 && smn < 24)
290         {
291           irownew = irownew1;
292           icolnew = icolnew1;
293         }
294
295       if (det == 0)
296         {
297           preADC[smn][irownew][icolnew] = sig;
298         }
299       else if (det == 1)
300         {
301           cpvADC[smn][irownew][icolnew] = sig;
302         }
303
304     } // while loop
305
306   for (Int_t idet = 0; idet < kDet; idet++)
307     {
308       for (Int_t ismn = 0; ismn < kSMN; ismn++)
309         {
310           for (Int_t irow = 0; irow < kRow; irow++)
311             {
312               for (Int_t icol = 0; icol < kCol; icol++)
313                 {
314                   if (idet == 0)
315                     {
316                       fCellADC[irow][icol] = 
317                         (Double_t) preADC[ismn][irow][icol];
318                     }
319                   else if (idet == 1)
320                     {
321                       fCellADC[irow][icol] = 
322                         (Double_t) cpvADC[ismn][irow][icol];
323                     }
324                 } // row
325             }     // col
326           pmdclust->DoClust(idet,ismn,fCellADC,pmdcont);
327           Int_t nentries1 = pmdcont->GetEntries();
328           //      cout << " nentries1 = " << nentries1 << endl;
329           for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
330             {
331               AliPMDcluster *pmdcl = 
332                 (AliPMDcluster*)pmdcont->UncheckedAt(ient1);
333               idet        = pmdcl->GetDetector();
334               ismn        = pmdcl->GetSMN();
335               clusdata[0] = pmdcl->GetClusX();
336               clusdata[1] = pmdcl->GetClusY();
337               clusdata[2] = pmdcl->GetClusADC();
338               clusdata[3] = pmdcl->GetClusCells();
339               clusdata[4] = pmdcl->GetClusRadius();
340               
341               AddRecPoint(idet,ismn,clusdata);
342             }
343           pmdcont->Clear();
344           
345           fTreeR->Fill();
346           ResetRecpoint();
347           
348         } // smn
349     }     // det (pre, cpv)
350
351   ResetCellADC();
352   fPMDLoader = fRunLoader->GetLoader("PMDLoader");  
353   fPMDLoader->WriteRecPoints("OVERWRITE");
354
355   //   delete the pointers
356   delete pmdclust;
357   delete pmdcont;
358
359   //  cout << " ***** End::Digits2RecPoints :: Raw *****" << endl;
360 }
361 // ------------------------------------------------------------------------- //
362 void AliPMDClusterFinder::SetCellEdepCut(Float_t ecut)
363 {
364   fEcut = ecut;
365 }
366 // ------------------------------------------------------------------------- //
367 void AliPMDClusterFinder::SetDebug(Int_t idebug)
368 {
369   fDebug = idebug;
370 }
371 // ------------------------------------------------------------------------- //
372 void AliPMDClusterFinder::AddRecPoint(Int_t idet,Int_t ismn,Float_t *clusdata)
373 {
374   // Add Reconstructed points
375   //
376   TClonesArray &lrecpoints = *fRecpoints;
377   AliPMDrecpoint1 *newrecpoint;
378   newrecpoint = new AliPMDrecpoint1(idet, ismn, clusdata);
379   new(lrecpoints[fNpoint++]) AliPMDrecpoint1(newrecpoint);
380   delete newrecpoint;
381 }
382 // ------------------------------------------------------------------------- //
383 void AliPMDClusterFinder::ResetCellADC()
384 {
385   // Reset the individual cell ADC value to zero
386   //
387   for(Int_t irow = 0; irow < fgkRow; irow++)
388     {
389       for(Int_t icol = 0; icol < fgkCol; icol++)
390         {
391           fCellADC[irow][icol] = 0.;
392         }
393     }
394 }
395 // ------------------------------------------------------------------------- //
396
397 void AliPMDClusterFinder::ResetRecpoint()
398 {
399   // Clear the list of reconstructed points
400   fNpoint = 0;
401   if (fRecpoints) fRecpoints->Clear();
402 }
403 // ------------------------------------------------------------------------- //
404 void AliPMDClusterFinder::Load()
405 {
406   // Load all the *.root files
407   //
408   fPMDLoader->LoadDigits("READ");
409   fPMDLoader->LoadRecPoints("recreate");
410 }
411 // ------------------------------------------------------------------------- //
412 void AliPMDClusterFinder::UnLoad()
413 {
414   // Unload all the *.root files
415   //
416   fPMDLoader->UnloadDigits();
417   fPMDLoader->UnloadRecPoints();
418 }
419 // ------------------------------------------------------------------------- //