bug fixed
[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
44 #include "AliPMDdigit.h"
45 #include "AliPMDClusterFinder.h"
46 #include "AliPMDClustering.h"
47 #include "AliPMDcluster.h"
48 #include "AliPMDrecpoint1.h"
49
50
51 ClassImp(AliPMDClusterFinder)
52
53 AliPMDClusterFinder::AliPMDClusterFinder()
54 {
55 //
56 // Default Constructor
57 //
58   if (!fRecpoints) fRecpoints = new TClonesArray("AliPMDrecpoint1", 1000);  
59   fNpoint = 0;
60
61   fDebug = 0;
62   fEcut  = 0.;
63
64 }
65 AliPMDClusterFinder::~AliPMDClusterFinder()
66 {
67   // Destructor
68   delete fRecpoints;
69 }
70
71 void AliPMDClusterFinder::OpengAliceFile(Char_t *file, Option_t *option)
72 {
73   // Loads galice.root file and corresponding header, kinematics
74   // hits and sdigits or digits depending on the option
75   //
76   fRunLoader = AliRunLoader::Open(file,AliConfig::fgkDefaultEventFolderName,
77                                   "UPDATE");
78   
79   if (!fRunLoader)
80    {
81      Error("Open","Can not open session for file %s.",file);
82    }
83   
84   fRunLoader->LoadgAlice();
85   fRunLoader->LoadHeader();
86   fRunLoader->LoadKinematics();
87
88   gAlice = fRunLoader->GetAliRun();
89   
90   if (gAlice)
91     {
92       printf("<AliPMDdigitizer::Open> ");
93       printf("AliRun object found on file.\n");
94     }
95   else
96     {
97       printf("<AliPMDdigitizer::Open> ");
98       printf("Could not find AliRun object.\n");
99     }
100   fPMD  = (AliPMD*)gAlice->GetDetector("PMD");
101   fPMDLoader = fRunLoader->GetLoader("PMDLoader");
102   if (fPMDLoader == 0x0)
103     {
104       cerr<<"OpengAlice : Can not find PMD or PMDLoader\n";
105     }
106
107   const char *cDR = strstr(option,"DR");
108
109   if (cDR)
110     {
111       fPMDLoader->LoadDigits("READ");
112       fPMDLoader->LoadRecPoints("recreate");
113     }
114 }
115
116 void AliPMDClusterFinder::Digits2RecPoints(Int_t ievt)
117 {
118   // Converts digits to recpoints after running clustering
119   // algorithm on CPV plane and PREshower plane
120   //
121   Int_t    det = 0,smn = 0;
122   Int_t    cellno;
123   Int_t    xpos,ypos;
124   Float_t  adc;
125   Int_t    isup;
126   Int_t    idet;
127   Float_t  clusdata[7];
128
129   TObjArray *pmdcont = new TObjArray();
130   AliPMDcluster  *pmdcl  = new AliPMDcluster;
131   AliPMDClustering *pmdclust = new AliPMDClustering();
132   pmdclust->SetDebug(fDebug);
133   pmdclust->SetEdepCut(fEcut);
134
135   fRunLoader->GetEvent(ievt);
136   //cout << " ***** Beginning::Digits2RecPoints *****" << endl;
137   fTreeD = fPMDLoader->TreeD();
138   if (fTreeD == 0x0)
139     {
140       cout << " Can not get TreeD" << endl;
141     }
142   AliPMDdigit  *pmddigit;
143   TBranch *branch = fTreeD->GetBranch("PMDDigit");
144   branch->SetAddress(&fDigits);
145
146   ResetRecpoint();
147   fTreeR = fPMDLoader->TreeR();
148   if (fTreeR == 0x0)
149     {
150       fPMDLoader->MakeTree("R");
151       fTreeR = fPMDLoader->TreeR();
152     }
153
154   Int_t bufsize = 16000;
155   fTreeR->Branch("PMDRecpoint", &fRecpoints, bufsize); 
156
157   Int_t nmodules = (Int_t) fTreeD->GetEntries();
158   
159   for (Int_t imodule = 0; imodule < nmodules; imodule++)
160     {
161       ResetCellADC();
162       fTreeD->GetEntry(imodule); 
163       Int_t nentries = fDigits->GetLast();
164       for (Int_t ient = 0; ient < nentries+1; ient++)
165         {
166           pmddigit = (AliPMDdigit*)fDigits->UncheckedAt(ient);
167           
168           det    = pmddigit->GetDetector();
169           smn    = pmddigit->GetSMNumber();
170           cellno = pmddigit->GetCellNumber();
171           adc    = pmddigit->GetADC();
172           //Int_t trno   = pmddigit->GetTrackNumber();
173
174           xpos = cellno/fCol;
175           ypos = cellno - xpos*fCol;
176           fCellADC[xpos][ypos] = (Double_t) adc;
177         }
178
179       idet = det;
180       isup = smn;
181       pmdclust->DoClust(fCellADC,pmdcont);
182       
183       Int_t nentries1 = pmdcont->GetEntries();
184       cout << " nentries1 = " << nentries1 << endl;
185       for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
186         {
187           clusdata[0] = (Float_t) idet;
188           clusdata[1] = (Float_t) isup;
189               
190           pmdcl = (AliPMDcluster*)pmdcont->UncheckedAt(ient1);
191               
192           clusdata[2] = pmdcl->GetClusX();
193           clusdata[3] = pmdcl->GetClusY();
194           clusdata[4] = pmdcl->GetClusADC();
195           clusdata[5] = pmdcl->GetClusCells();
196           clusdata[6] = pmdcl->GetClusRadius();
197           
198           AddRecPoint(clusdata);
199         }
200       pmdcont->Clear();
201       
202       fTreeR->Fill();
203       ResetRecpoint();
204
205     } // modules
206
207   ResetCellADC();
208   
209   fPMDLoader->WriteRecPoints("OVERWRITE");
210
211   //   delete the pointers
212   delete pmdclust;
213   delete pmdcont;
214     
215   //  cout << " ***** End::Digits2RecPoints *****" << endl;
216 }
217
218 void AliPMDClusterFinder::SetCellEdepCut(Float_t ecut)
219 {
220   fEcut = ecut;
221 }
222 void AliPMDClusterFinder::SetDebug(Int_t idebug)
223 {
224   fDebug = idebug;
225 }
226
227 void AliPMDClusterFinder::AddRecPoint(Float_t *clusdata)
228 {
229   // Add Reconstructed points
230   //
231   TClonesArray &lrecpoints = *fRecpoints;
232   AliPMDrecpoint1 *newrecpoint;
233   newrecpoint = new AliPMDrecpoint1(clusdata);
234   new(lrecpoints[fNpoint++]) AliPMDrecpoint1(newrecpoint);
235   delete newrecpoint;
236 }
237 void AliPMDClusterFinder::ResetCellADC()
238 {
239   // Reset the individual cell ADC value to zero
240   //
241   for(Int_t irow = 0; irow < fRow; irow++)
242     {
243       for(Int_t icol = 0; icol < fCol; icol++)
244         {
245           fCellADC[irow][icol] = 0.;
246         }
247     }
248 }
249
250 void AliPMDClusterFinder::ResetRecpoint()
251 {
252   // Clear the list of reconstructed points
253   fNpoint = 0;
254   if (fRecpoints) fRecpoints->Clear();
255 }
256 void AliPMDClusterFinder::UnLoad(Option_t *option)
257 {
258   // Unload all the *.root files
259   //
260   const char *cR = strstr(option,"R");
261
262   fRunLoader->UnloadgAlice();
263   fRunLoader->UnloadHeader();
264   fRunLoader->UnloadKinematics();
265
266   if (cR)
267     {
268       fPMDLoader->UnloadDigits();
269     }
270 }