discrimination based on PRE-CPV matching
[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 <TTree.h>
27 #include <TObjArray.h>
28 #include <TClonesArray.h>
29
30 #include "AliLog.h"
31 #include "AliRunLoader.h"
32 #include "AliLoader.h"
33 #include "AliRawReader.h"
34
35 #include "AliPMDdigit.h"
36 #include "AliPMDClusterFinder.h"
37 #include "AliPMDClustering.h"
38 #include "AliPMDcluster.h"
39 #include "AliPMDrecpoint1.h"
40 #include "AliPMDRawStream.h"
41
42
43
44 ClassImp(AliPMDClusterFinder)
45
46 AliPMDClusterFinder::AliPMDClusterFinder():
47   fRunLoader(0),
48   fPMDLoader(0),
49   fTreeD(0),
50   fTreeR(0),
51   fDigits(new TClonesArray("AliPMDdigit", 1000)),
52   fRecpoints(new TClonesArray("AliPMDrecpoint1", 1000)),
53   fNpoint(0),
54   fEcut(0.)
55 {
56 //
57 // Constructor
58 //
59 }
60 // ------------------------------------------------------------------------- //
61 AliPMDClusterFinder::AliPMDClusterFinder(AliRunLoader* runLoader):
62   fRunLoader(runLoader),
63   fPMDLoader(runLoader->GetLoader("PMDLoader")),
64   fTreeD(0),
65   fTreeR(0),
66   fDigits(new TClonesArray("AliPMDdigit", 1000)),
67   fRecpoints(new TClonesArray("AliPMDrecpoint1", 1000)),
68   fNpoint(0),
69   fEcut(0.)
70 {
71 //
72 // Constructor
73 //
74 }
75 // ------------------------------------------------------------------------- //
76 AliPMDClusterFinder::~AliPMDClusterFinder()
77 {
78   // Destructor
79   if (fDigits)
80     {
81       fDigits->Delete();
82       delete fDigits;
83       fDigits=0;
84     }
85   if (fRecpoints)
86     {
87       fRecpoints->Delete();
88       delete fRecpoints;
89       fRecpoints=0;
90     }
91 }
92 // ------------------------------------------------------------------------- //
93
94 void AliPMDClusterFinder::Digits2RecPoints(Int_t ievt)
95 {
96   // Converts digits to recpoints after running clustering
97   // algorithm on CPV plane and PREshower plane
98   //
99   Int_t    det = 0,smn = 0;
100   Int_t    xpos,ypos;
101   Float_t  adc;
102   Int_t    ismn;
103   Int_t    idet;
104   Float_t  clusdata[5];
105
106   TObjArray *pmdcont = new TObjArray();
107   AliPMDClustering *pmdclust = new AliPMDClustering();
108
109   pmdclust->SetEdepCut(fEcut);
110
111   fRunLoader->GetEvent(ievt);
112
113
114   fTreeD = fPMDLoader->TreeD();
115   if (fTreeD == 0x0)
116     {
117       AliFatal("AliPMDClusterFinder: Can not get TreeD");
118
119     }
120   AliPMDdigit  *pmddigit;
121   TBranch *branch = fTreeD->GetBranch("PMDDigit");
122   branch->SetAddress(&fDigits);
123
124   ResetRecpoint();
125
126   fTreeR = fPMDLoader->TreeR();
127   if (fTreeR == 0x0)
128     {
129       fPMDLoader->MakeTree("R");
130       fTreeR = fPMDLoader->TreeR();
131     }
132
133   Int_t bufsize = 16000;
134   fTreeR->Branch("PMDRecpoint", &fRecpoints, bufsize); 
135
136   Int_t nmodules = (Int_t) fTreeD->GetEntries();
137
138   for (Int_t imodule = 0; imodule < nmodules; imodule++)
139     {
140       ResetCellADC();
141       fTreeD->GetEntry(imodule); 
142       Int_t nentries = fDigits->GetLast();
143       for (Int_t ient = 0; ient < nentries+1; ient++)
144         {
145           pmddigit = (AliPMDdigit*)fDigits->UncheckedAt(ient);
146           
147           det    = pmddigit->GetDetector();
148           smn    = pmddigit->GetSMNumber();
149           xpos   = pmddigit->GetRow();
150           ypos   = pmddigit->GetColumn();
151           adc    = pmddigit->GetADC();
152           //Int_t trno   = pmddigit->GetTrackNumber();
153           fCellADC[xpos][ypos] = (Double_t) adc;
154         }
155
156       idet = det;
157       ismn = smn;
158       pmdclust->DoClust(idet,ismn,fCellADC,pmdcont);
159       
160       Int_t nentries1 = pmdcont->GetEntries();
161
162       AliDebug(1,Form("Total number of clusters/module = %d",nentries1));
163
164       for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
165         {
166           AliPMDcluster *pmdcl = (AliPMDcluster*)pmdcont->UncheckedAt(ient1);
167           idet        = pmdcl->GetDetector();
168           ismn        = pmdcl->GetSMN();
169           clusdata[0] = pmdcl->GetClusX();
170           clusdata[1] = pmdcl->GetClusY();
171           clusdata[2] = pmdcl->GetClusADC();
172           clusdata[3] = pmdcl->GetClusCells();
173           clusdata[4] = pmdcl->GetClusRadius();
174
175           AddRecPoint(idet,ismn,clusdata);
176         }
177       pmdcont->Clear();
178       
179       fTreeR->Fill();
180       ResetRecpoint();
181
182     } // modules
183
184   ResetCellADC();
185   fPMDLoader = fRunLoader->GetLoader("PMDLoader");  
186   fPMDLoader->WriteRecPoints("OVERWRITE");
187
188   //   delete the pointers
189   delete pmdclust;
190   delete pmdcont;
191     
192 }
193 // ------------------------------------------------------------------------- //
194
195 void AliPMDClusterFinder::Digits2RecPoints(AliRawReader *rawReader,
196                                            TTree *clustersTree)
197 {
198   // Converts RAW data to recpoints after running clustering
199   // algorithm on CPV and PREshower plane
200   //
201
202   Float_t  clusdata[5];
203
204   TObjArray *pmdcont = new TObjArray();
205   AliPMDClustering *pmdclust = new AliPMDClustering();
206
207   pmdclust->SetEdepCut(fEcut);
208
209   ResetRecpoint();
210
211   Int_t bufsize = 16000;
212   clustersTree->Branch("PMDRecpoint", &fRecpoints, bufsize); 
213
214   const Int_t kDDL = 6;
215   const Int_t kRow = 48;
216   const Int_t kCol = 96;
217
218   Int_t idet = 0;
219   Int_t iSMN = 0;
220   
221   for (Int_t indexDDL = 0; indexDDL < kDDL; indexDDL++)
222     {
223       if (indexDDL < 4)
224         {
225           iSMN = 6;
226         }
227       else if (indexDDL >= 4)
228         {
229           iSMN = 12;
230         }
231       Int_t ***precpvADC;
232       precpvADC = new int **[iSMN];
233       for (Int_t i=0; i<iSMN; i++) precpvADC[i] = new int *[kRow];
234       for (Int_t i=0; i<iSMN;i++)
235         {
236           for (Int_t j=0; j<kRow; j++) precpvADC[i][j] = new int [kCol];
237         }
238       for (Int_t i = 0; i < iSMN; i++)
239         {
240           for (Int_t j = 0; j < kRow; j++)
241             {
242               for (Int_t k = 0; k < kCol; k++)
243                 {
244                   precpvADC[i][j][k] = 0;
245                 }
246             }
247         }
248       ResetCellADC();
249       rawReader->Reset();
250       AliPMDRawStream pmdinput(rawReader);
251       rawReader->Select(12, indexDDL, indexDDL);
252       while(pmdinput.Next())
253         {
254           Int_t det = pmdinput.GetDetector();
255           Int_t smn = pmdinput.GetSMN();
256           //Int_t mcm = pmdinput.GetMCM();
257           //Int_t chno = pmdinput.GetChannel();
258           Int_t row = pmdinput.GetRow();
259           Int_t col = pmdinput.GetColumn();
260           Int_t sig = pmdinput.GetSignal();
261           
262           Int_t indexsmn = 0;
263
264           if (indexDDL < 4)
265             {
266               if (det != 0)
267                 AliError(Form("*DDL %d and Detector NUMBER %d NOT MATCHING *",
268                               indexDDL, det));
269               indexsmn = smn - indexDDL * 6;
270             }
271           else if (indexDDL == 4)
272             {
273               if (det != 1)
274                 AliError(Form("*DDL %d and Detector NUMBER %d NOT MATCHING *",
275                               indexDDL, det));
276               if (smn < 6)
277                 {
278                   indexsmn = smn;
279                 }
280               else if (smn >= 12 && smn < 18)
281                 {
282                   indexsmn = smn - 6;
283                 }
284             }
285           else if (indexDDL == 5)
286             {
287               if (det != 1)
288                 AliError(Form("*DDL %d and Detector NUMBER %d NOT MATCHING *",
289                               indexDDL, det));
290               if (smn >= 6 && smn < 12)
291                 {
292                   indexsmn = smn - 6;
293                 }
294               else if (smn >= 18 && smn < 24)
295                 {
296                   indexsmn = smn - 12;
297                 }
298             }         
299           precpvADC[indexsmn][row][col] = sig;
300         } // while loop
301
302       Int_t ismn = 0;
303       for (Int_t indexsmn = 0; indexsmn < iSMN; indexsmn++)
304         {
305           ResetCellADC();
306           for (Int_t irow = 0; irow < kRow; irow++)
307             {
308               for (Int_t icol = 0; icol < kCol; icol++)
309                 {
310                   fCellADC[irow][icol] = 
311                     (Double_t) precpvADC[indexsmn][irow][icol];
312                 } // row
313             }     // col
314           if (indexDDL < 4)
315             {
316               ismn = indexsmn + indexDDL * 6;
317               idet = 0;
318             }
319           else if (indexDDL == 4)
320             {
321               if (indexsmn < 6)
322                 {
323                   ismn = indexsmn;
324                 }
325               else if (indexsmn >= 6 && indexsmn < 12)
326                 {
327                   ismn = indexsmn + 6;
328                 }
329               idet = 1;
330             }
331           else if (indexDDL == 5)
332             {
333               if (indexsmn < 6)
334                 {
335                   ismn = indexsmn + 6;
336                 }
337               else if (indexsmn >= 6 && indexsmn < 12)
338                 {
339                   ismn = indexsmn + 12;
340                 }
341               idet = 1;
342             }
343
344
345           pmdclust->DoClust(idet,ismn,fCellADC,pmdcont);
346           Int_t nentries1 = pmdcont->GetEntries();
347
348           AliDebug(1,Form("Total number of clusters/module = %d",nentries1));
349
350           for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
351             {
352               AliPMDcluster *pmdcl = 
353                 (AliPMDcluster*)pmdcont->UncheckedAt(ient1);
354               idet        = pmdcl->GetDetector();
355               ismn        = pmdcl->GetSMN();
356               clusdata[0] = pmdcl->GetClusX();
357               clusdata[1] = pmdcl->GetClusY();
358               clusdata[2] = pmdcl->GetClusADC();
359               clusdata[3] = pmdcl->GetClusCells();
360               clusdata[4] = pmdcl->GetClusRadius();
361
362               AddRecPoint(idet,ismn,clusdata);
363             }
364           pmdcont->Clear();
365           
366           //fTreeR->Fill();
367           clustersTree->Fill();
368           ResetRecpoint();
369
370
371         } // smn
372
373       for (Int_t i=0; i<iSMN; i++)
374         {
375           for (Int_t j=0; j<kRow; j++) delete [] precpvADC[i][j];
376         }
377       for (Int_t i=0; i<iSMN; i++) delete [] precpvADC[i];
378       delete precpvADC;
379     } // DDL Loop
380   
381   ResetCellADC();
382   
383   //   delete the pointers
384   delete pmdclust;
385   delete pmdcont;
386
387 }
388 // ------------------------------------------------------------------------- //
389
390 void AliPMDClusterFinder::Digits2RecPoints(Int_t ievt, AliRawReader *rawReader)
391 {
392   // Converts RAW data to recpoints after running clustering
393   // algorithm on CPV and PREshower plane
394   //
395
396   Float_t  clusdata[5];
397
398   TObjArray *pmdcont = new TObjArray();
399   AliPMDClustering *pmdclust = new AliPMDClustering();
400
401   pmdclust->SetEdepCut(fEcut);
402
403   fRunLoader->GetEvent(ievt);
404
405   ResetRecpoint();
406
407   fTreeR = fPMDLoader->TreeR();
408   if (fTreeR == 0x0)
409     {
410       fPMDLoader->MakeTree("R");
411       fTreeR = fPMDLoader->TreeR();
412     }
413   Int_t bufsize = 16000;
414   fTreeR->Branch("PMDRecpoint", &fRecpoints, bufsize); 
415
416   const Int_t kDDL = 6;
417   const Int_t kRow = 48;
418   const Int_t kCol = 96;
419
420   Int_t idet = 0;
421   Int_t iSMN = 0;
422   
423   for (Int_t indexDDL = 0; indexDDL < kDDL; indexDDL++)
424     {
425       if (indexDDL < 4)
426         {
427           iSMN = 6;
428         }
429       else if (indexDDL >= 4)
430         {
431           iSMN = 12;
432         }
433       Int_t ***precpvADC;
434       precpvADC = new int **[iSMN];
435       for (Int_t i=0; i<iSMN; i++) precpvADC[i] = new int *[kRow];
436       for (Int_t i=0; i<iSMN;i++)
437         {
438           for (Int_t j=0; j<kRow; j++) precpvADC[i][j] = new int [kCol];
439         }
440       for (Int_t i = 0; i < iSMN; i++)
441         {
442           for (Int_t j = 0; j < kRow; j++)
443             {
444               for (Int_t k = 0; k < kCol; k++)
445                 {
446                   precpvADC[i][j][k] = 0;
447                 }
448             }
449         }
450       ResetCellADC();
451       rawReader->Reset();
452       AliPMDRawStream pmdinput(rawReader);
453       rawReader->Select(12, indexDDL, indexDDL);
454       while(pmdinput.Next())
455         {
456           Int_t det = pmdinput.GetDetector();
457           Int_t smn = pmdinput.GetSMN();
458           //Int_t mcm = pmdinput.GetMCM();
459           //Int_t chno = pmdinput.GetChannel();
460           Int_t row = pmdinput.GetRow();
461           Int_t col = pmdinput.GetColumn();
462           Int_t sig = pmdinput.GetSignal();
463           
464           Int_t indexsmn = 0;
465
466           if (indexDDL < 4)
467             {
468               if (det != 0)
469                 AliError(Form("*DDL %d and Detector NUMBER %d NOT MATCHING *",
470                               indexDDL, det));
471               indexsmn = smn - indexDDL * 6;
472             }
473           else if (indexDDL == 4)
474             {
475               if (det != 1)
476                 AliError(Form("*DDL %d and Detector NUMBER %d NOT MATCHING *",
477                               indexDDL, det));
478               if (smn < 6)
479                 {
480                   indexsmn = smn;
481                 }
482               else if (smn >= 12 && smn < 18)
483                 {
484                   indexsmn = smn - 6;
485                 }
486             }
487           else if (indexDDL == 5)
488             {
489               if (det != 1)
490                 AliError(Form("*DDL %d and Detector NUMBER %d NOT MATCHING *",
491                               indexDDL, det));
492               if (smn >= 6 && smn < 12)
493                 {
494                   indexsmn = smn - 6;
495                 }
496               else if (smn >= 18 && smn < 24)
497                 {
498                   indexsmn = smn - 12;
499                 }
500             }         
501           precpvADC[indexsmn][row][col] = sig;
502         } // while loop
503
504       Int_t ismn = 0;
505       for (Int_t indexsmn = 0; indexsmn < iSMN; indexsmn++)
506         {
507           ResetCellADC();
508           for (Int_t irow = 0; irow < kRow; irow++)
509             {
510               for (Int_t icol = 0; icol < kCol; icol++)
511                 {
512                   fCellADC[irow][icol] = 
513                     (Double_t) precpvADC[indexsmn][irow][icol];
514                 } // row
515             }     // col
516           if (indexDDL < 4)
517             {
518               ismn = indexsmn + indexDDL * 6;
519               idet = 0;
520             }
521           else if (indexDDL == 4)
522             {
523               if (indexsmn < 6)
524                 {
525                   ismn = indexsmn;
526                 }
527               else if (indexsmn >= 6 && indexsmn < 12)
528                 {
529                   ismn = indexsmn + 6;
530                 }
531               idet = 1;
532             }
533           else if (indexDDL == 5)
534             {
535               if (indexsmn < 6)
536                 {
537                   ismn = indexsmn + 6;
538                 }
539               else if (indexsmn >= 6 && indexsmn < 12)
540                 {
541                   ismn = indexsmn + 12;
542                 }
543               idet = 1;
544             }
545
546
547           pmdclust->DoClust(idet,ismn,fCellADC,pmdcont);
548           Int_t nentries1 = pmdcont->GetEntries();
549
550           AliDebug(1,Form("Total number of clusters/module = %d",nentries1));
551
552           for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
553             {
554               AliPMDcluster *pmdcl = 
555                 (AliPMDcluster*)pmdcont->UncheckedAt(ient1);
556               idet        = pmdcl->GetDetector();
557               ismn        = pmdcl->GetSMN();
558               clusdata[0] = pmdcl->GetClusX();
559               clusdata[1] = pmdcl->GetClusY();
560               clusdata[2] = pmdcl->GetClusADC();
561               clusdata[3] = pmdcl->GetClusCells();
562               clusdata[4] = pmdcl->GetClusRadius();
563
564               AddRecPoint(idet,ismn,clusdata);
565             }
566           pmdcont->Clear();
567           
568           fTreeR->Fill();
569           ResetRecpoint();
570
571
572         } // smn
573
574       for (Int_t i=0; i<iSMN; i++)
575         {
576           for (Int_t j=0; j<kRow; j++) delete [] precpvADC[i][j];
577         }
578       for (Int_t i=0; i<iSMN; i++) delete [] precpvADC[i];
579       delete precpvADC;
580     } // DDL Loop
581   
582   ResetCellADC();
583   
584   fPMDLoader = fRunLoader->GetLoader("PMDLoader");  
585   fPMDLoader->WriteRecPoints("OVERWRITE");
586
587   //   delete the pointers
588   delete pmdclust;
589   delete pmdcont;
590
591 }
592 // ------------------------------------------------------------------------- //
593 void AliPMDClusterFinder::SetCellEdepCut(Float_t ecut)
594 {
595   fEcut = ecut;
596 }
597 // ------------------------------------------------------------------------- //
598 void AliPMDClusterFinder::AddRecPoint(Int_t idet,Int_t ismn,Float_t *clusdata)
599 {
600   // Add Reconstructed points
601   //
602   TClonesArray &lrecpoints = *fRecpoints;
603   AliPMDrecpoint1 *newrecpoint;
604   newrecpoint = new AliPMDrecpoint1(idet, ismn, clusdata);
605   new(lrecpoints[fNpoint++]) AliPMDrecpoint1(newrecpoint);
606   delete newrecpoint;
607 }
608 // ------------------------------------------------------------------------- //
609 void AliPMDClusterFinder::ResetCellADC()
610 {
611   // Reset the individual cell ADC value to zero
612   //
613   for(Int_t irow = 0; irow < fgkRow; irow++)
614     {
615       for(Int_t icol = 0; icol < fgkCol; icol++)
616         {
617           fCellADC[irow][icol] = 0.;
618         }
619     }
620 }
621 // ------------------------------------------------------------------------- //
622
623 void AliPMDClusterFinder::ResetRecpoint()
624 {
625   // Clear the list of reconstructed points
626   fNpoint = 0;
627   if (fRecpoints) fRecpoints->Clear();
628 }
629 // ------------------------------------------------------------------------- //
630 void AliPMDClusterFinder::Load()
631 {
632   // Load all the *.root files
633   //
634   fPMDLoader->LoadDigits("READ");
635   fPMDLoader->LoadRecPoints("recreate");
636 }
637 // ------------------------------------------------------------------------- //
638 void AliPMDClusterFinder::LoadClusters()
639 {
640   // Load all the *.root files
641   //
642   fPMDLoader->LoadRecPoints("recreate");
643 }
644 // ------------------------------------------------------------------------- //
645 void AliPMDClusterFinder::UnLoad()
646 {
647   // Unload all the *.root files
648   //
649   fPMDLoader->UnloadDigits();
650   fPMDLoader->UnloadRecPoints();
651 }
652 // ------------------------------------------------------------------------- //
653 void AliPMDClusterFinder::UnLoadClusters()
654 {
655   // Unload all the *.root files
656   //
657   fPMDLoader->UnloadRecPoints();
658 }
659 // ------------------------------------------------------------------------- //