new file to test SDigits
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALClusterizerNxN.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 // This class derives from AliEMCALClustrerizer but also keeps the API of AliEMCALClusterizerv1
17 // Algorithm:
18 // 1. peek the most energetic cell
19 // 2. assign it as a center of the cluster and add cells surrounding it: 3x3, 5x5...
20 // 3. remove the cells contributing to the cluster
21 // 4. start from 1 for the remaining clusters
22 // 5. cluster splitting (not implemented yet) - use the shape analysis to resolve the energy sharing
23 // - for high energy clusters check the surrounding of the 3x3 clusters for extra energy 
24 // (merge 3x3 clusters and resolve the internal energy sharing - case for 2 clusters merged)
25 // Use Case:
26 //  root [0] AliEMCALClusterizerNxN * cl = new AliEMCALClusterizerNxN("galice.root")  
27 //  Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
28 //               //reads gAlice from header file "..."                      
29 //  root [1] cl->ExecuteTask()  
30 //               //finds RecPoints in all events stored in galice.root
31 //  root [2] cl->SetDigitsBranch("digits2") 
32 //               //sets another title for Digitis (input) branch
33 //  root [3] cl->SetRecPointsBranch("recp2")  
34 //               //sets another title four output branches
35 //  root [4] cl->SetTowerLocalMaxCut(0.03)  
36 //               //set clusterization parameters
37 //  root [5] cl->ExecuteTask("deb all time")  
38 //               //once more finds RecPoints options are 
39 //               // deb - print number of found rec points
40 //               // deb all - print number of found RecPoints and some their characteristics 
41 //               // time - print benchmarking results
42
43 // --- ROOT system ---
44 #include <TMath.h> 
45 #include <TMinuit.h>
46 #include <TTree.h> 
47 #include <TBenchmark.h>
48 #include <TBrowser.h>
49 #include <TROOT.h>
50 #include <TClonesArray.h>
51
52 // --- Standard library ---
53 #include <cassert>
54
55 // --- AliRoot header files ---
56 #include "AliLog.h"
57 #include "AliEMCALClusterizerNxN.h"
58 #include "AliEMCALRecPoint.h"
59 #include "AliEMCALDigit.h"
60 #include "AliEMCALGeometry.h"
61 #include "AliCaloCalibPedestal.h"
62 #include "AliEMCALCalibData.h"
63 #include "AliESDCaloCluster.h"
64 #include "AliEMCALUnfolding.h"
65
66 ClassImp(AliEMCALClusterizerNxN)
67
68 //____________________________________________________________________________
69 AliEMCALClusterizerNxN::AliEMCALClusterizerNxN()
70   : AliEMCALClusterizer(), fNRowDiff(1), fNColDiff(1), fEnergyGrad(0)
71 {
72   // ctor with the indication of the file where header Tree and digits Tree are stored
73 }
74
75 //____________________________________________________________________________
76 AliEMCALClusterizerNxN::AliEMCALClusterizerNxN(AliEMCALGeometry* geometry)
77   : AliEMCALClusterizer(geometry), fNRowDiff(1), fNColDiff(1), fEnergyGrad(0)
78 {
79   // ctor with the indication of the file where header Tree and digits Tree are stored
80   // use this contructor to avoid usage of Init() which uses runloader
81   // change needed by HLT - MP
82 }
83
84 //____________________________________________________________________________
85 AliEMCALClusterizerNxN::AliEMCALClusterizerNxN(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped)
86 : AliEMCALClusterizer(geometry, calib, caloped), fNRowDiff(1), fNColDiff(1), fEnergyGrad(0)
87 {
88   // ctor, geometry and calibration are initialized elsewhere.
89 }
90
91 //____________________________________________________________________________
92 AliEMCALClusterizerNxN::~AliEMCALClusterizerNxN()
93 {
94   // dtor
95 }
96
97 //____________________________________________________________________________
98 void AliEMCALClusterizerNxN::Digits2Clusters(Option_t * option)
99 {
100   // Steering method to perform clusterization for the current event 
101   // in AliEMCALLoader
102   
103   if (strstr(option,"tim"))
104     gBenchmark->Start("EMCALClusterizer"); 
105   
106   if (strstr(option,"print"))
107     Print(""); 
108   
109   //Get calibration parameters from file or digitizer default values.
110   GetCalibrationParameters();
111   
112   //Get dead channel map from file or digitizer default values.
113   GetCaloCalibPedestal();
114         
115   MakeClusters();  //only the real clusters
116   
117   if (fToUnfold) {
118     fClusterUnfolding->SetInput(fNumberOfECAClusters,fRecPoints,fDigitsArr);
119     fClusterUnfolding->MakeUnfolding();
120   }
121   
122   //Evaluate position, dispersion and other RecPoint properties for EC section 
123   for (Int_t index = 0; index < fRecPoints->GetEntries(); index++) { 
124     AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
125     if (rp) {
126       rp->EvalAll(fECAW0,fDigitsArr,fJustClusters);
127       AliDebug(5, Form("MAX INDEX %d ", rp->GetMaximalEnergyIndex()));
128       //For each rec.point set the distance to the nearest bad crystal
129       if (fCaloPed)
130         rp->EvalDistanceToBadChannels(fCaloPed);
131     }
132   }
133   
134   fRecPoints->Sort();
135   
136   for (Int_t index = 0; index < fRecPoints->GetEntries(); index++) {
137     AliEMCALRecPoint *rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
138     if (rp) {
139       rp->SetIndexInList(index);
140     }
141     else AliFatal("RecPoint NULL!!");
142   }
143   
144   if (fTreeR)
145     fTreeR->Fill();
146   
147   if (strstr(option,"deb") || strstr(option,"all"))  
148     PrintRecPoints(option);
149   
150   AliDebug(1,Form("EMCAL Clusterizer found %d Rec Points",fRecPoints->GetEntriesFast()));
151   
152   if (strstr(option,"tim")) {
153     gBenchmark->Stop("EMCALClusterizer");
154     printf("Exec took %f seconds for Clusterizing", 
155            gBenchmark->GetCpuTime("EMCALClusterizer"));
156   }    
157 }
158
159 //____________________________________________________________________________
160 Int_t AliEMCALClusterizerNxN::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2, Bool_t & shared) const
161 {
162   // Gives the neighbourness of two digits = 0 are not neighbour ; continue searching 
163   //                                       = 1 are neighbour
164   //                                       = 2 is in different SM; continue searching 
165   // In case it is in different SM, but same phi rack, check if neigbours at eta=0
166   // neighbours are defined as digits having at least a common side 
167   // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster 
168   //                                      which is compared to a digit (d2)  not yet in a cluster  
169   
170   if (fEnergyGrad) { //false by default
171     if (d2->GetCalibAmp()>d1->GetCalibAmp())
172       return 3; // energy of neighboring cell should be smaller in order to become a neighbor
173   }
174
175   Int_t nSupMod1=0, nModule1=0, nIphi1=0, nIeta1=0, iphi1=0, ieta1=0;
176   Int_t nSupMod2=0, nModule2=0, nIphi2=0, nIeta2=0, iphi2=0, ieta2=0;
177   Int_t rowdiff=0, coldiff=0;
178   
179   shared = kFALSE;
180   
181   fGeom->GetCellIndex(d1->GetId(), nSupMod1,nModule1,nIphi1,nIeta1);
182   fGeom->GetCellIndex(d2->GetId(), nSupMod2,nModule2,nIphi2,nIeta2);
183   fGeom->GetCellPhiEtaIndexInSModule(nSupMod1,nModule1,nIphi1,nIeta1, iphi1,ieta1);
184   fGeom->GetCellPhiEtaIndexInSModule(nSupMod2,nModule2,nIphi2,nIeta2, iphi2,ieta2);
185   
186   //If different SM, check if they are in the same phi, then consider cells close to eta=0 as neighbours; May 2010
187   if (nSupMod1 != nSupMod2 ) 
188     {
189       //Check if the 2 SM are in the same PHI position (0,1), (2,3), ...
190       Float_t smPhi1 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod1);
191       Float_t smPhi2 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod2);
192       
193       if (!TMath::AreEqualAbs(smPhi1, smPhi2, 1e-3)) return 2; //Not phi rack equal, not neighbours
194       
195       // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
196       // C Side impair SM, nSupMod%2=1; A side pair SM nSupMod%2=0
197       if (nSupMod1%2) ieta1+=AliEMCALGeoParams::fgkEMCALCols;
198       else           ieta2+=AliEMCALGeoParams::fgkEMCALCols;
199       
200       shared = kTRUE; // maybe a shared cluster, we know this later, set it for the moment.
201       
202     }//Different SM, same phi
203   
204   rowdiff = TMath::Abs(iphi1 - iphi2);  
205   coldiff = TMath::Abs(ieta1 - ieta2);  
206   
207   // neighbours +-1 in col and row
208   if ( TMath::Abs(coldiff) <= fNColDiff && TMath::Abs(rowdiff) <= fNRowDiff)
209     {
210       
211       AliDebug(9, Form("AliEMCALClusterizerNxN::AreNeighbours(): id1=%d, (row %d, col %d) ; id2=%d, (row %d, col %d), shared %d \n",
212                        d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2, shared));
213       
214       return 1;
215     }//Neighbours
216   else 
217     {
218       AliDebug(9, Form("NOT AliEMCALClusterizerNxN::AreNeighbours(): id1=%d, (row %d, col %d) ; id2=%d, (row %d, col %d), shared %d \n",
219                        d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2, shared));
220       shared = kFALSE;
221       return 2; 
222     }//Not neighbours
223 }
224
225 //____________________________________________________________________________
226 void AliEMCALClusterizerNxN::MakeClusters()
227 {
228   // Make clusters
229   
230   if (fGeom==0) 
231     AliFatal("Did not get geometry from EMCALLoader");
232   
233   fNumberOfECAClusters = 0;
234   fRecPoints->Delete();
235   
236   // Set up TObjArray with pointers to digits to work on, calibrate digits 
237   TObjArray digitsC;
238   TIter nextdigit(fDigitsArr);
239   AliEMCALDigit *digit = 0;
240   while ( (digit = static_cast<AliEMCALDigit*>(nextdigit())) ) {
241     Float_t dEnergyCalibrated = digit->GetAmplitude();
242     Float_t time              = digit->GetTime();
243     Calibrate(dEnergyCalibrated,time ,digit->GetId());
244     digit->SetCalibAmp(dEnergyCalibrated);
245     digit->SetTime(time);    
246     digitsC.AddLast(digit);
247   }
248   
249   TIter nextdigitC(&digitsC);
250   
251   AliDebug(1,Form("MakeClusters: Number of digits %d  -> (e %f)\n",
252                   fDigitsArr->GetEntries(),fMinECut));
253   
254   Bool_t bDone = kFALSE;
255   while ( bDone != kTRUE )
256   {
257     //first sort the digits:
258     Int_t iMaxEnergyDigit = -1;
259     Float_t dMaxEnergyDigit = -1;
260     AliEMCALDigit *pMaxEnergyDigit = 0;
261     nextdigitC.Reset();
262     while ( (digit = static_cast<AliEMCALDigit *>(nextdigitC())) ) 
263     { // scan over the list of digitsC
264       Float_t dEnergyCalibrated = digit->GetCalibAmp();
265       Float_t time              = digit->GetTime();
266       if (fGeom->CheckAbsCellId(digit->GetId()) &&
267           dEnergyCalibrated > fMinECut          &&
268           time              < fTimeMax          &&
269           time              > fTimeMin             ) // no threshold by default!
270       {                                              // needs to be set in OCDB!
271         if (dEnergyCalibrated > dMaxEnergyDigit) 
272         {
273           dMaxEnergyDigit = dEnergyCalibrated;
274           iMaxEnergyDigit = digit->GetId();
275           pMaxEnergyDigit = digit;
276         }
277       }
278     }
279     if (iMaxEnergyDigit < 0 || digitsC.GetEntries() <= 0) 
280     {
281       bDone = kTRUE;
282       continue;
283     }
284     
285     AliDebug (2, Form("Max digit found: %1.5f AbsId: %d", dMaxEnergyDigit, iMaxEnergyDigit));
286     
287     // keep the candidate digits in a list
288     TList clusterDigitList;
289     clusterDigitList.SetOwner(kFALSE);
290     clusterDigitList.AddLast(pMaxEnergyDigit);   
291     
292     Double_t clusterCandidateEnergy = dMaxEnergyDigit;
293     
294     // now loop over the rest of the digits and cluster into NxN cluster 
295     // we do not actually cluster yet: we keep them in the list clusterDigitList
296     nextdigitC.Reset();
297     while ( (digit = static_cast<AliEMCALDigit *>(nextdigitC())) ) 
298     { // scan over the list of digitsC
299       if (digit == pMaxEnergyDigit) continue;
300       Float_t dEnergyCalibrated = digit->GetCalibAmp();
301       AliDebug(5, Form("-> Digit ENERGY: %1.5f", dEnergyCalibrated));
302       if (fGeom->CheckAbsCellId(digit->GetId()) && dEnergyCalibrated > 0.0  )
303       {
304         Float_t time = pMaxEnergyDigit->GetTime(); //Time or TimeR?
305         if (TMath::Abs(time - digit->GetTime()) > fTimeCut ) continue; //Time or TimeR?
306         Bool_t shared = kFALSE; //cluster shared by 2 SuperModules?
307         if (AreNeighbours(pMaxEnergyDigit, digit, shared) == 1) // call (digit,digitN) in THAT order !!!!! 
308         {      
309           clusterDigitList.AddLast(digit);
310           clusterCandidateEnergy += dEnergyCalibrated;
311         }
312       }
313     }// loop over the next digits
314     
315     // start a cluster here only if a cluster energy is larger than clustering threshold
316     AliDebug(5, Form("Clusterization threshold is %f MeV", fECAClusteringThreshold));
317     if (clusterCandidateEnergy > fECAClusteringThreshold)
318     {
319       if (fNumberOfECAClusters >= fRecPoints->GetSize()) 
320         fRecPoints->Expand(2*fNumberOfECAClusters+1);
321       
322       (*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("") ;
323       AliEMCALRecPoint *recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(fNumberOfECAClusters) ) ;
324       // AliEMCALRecPoint *recPoint = new  AliEMCALRecPoint(""); 
325       // fRecPoints->AddAt(recPoint, fNumberOfECAClusters);
326       // recPoint = static_cast<AliEMCALRecPoint *>(fRecPoints->At(fNumberOfECAClusters)); 
327       if (recPoint) {
328         fNumberOfECAClusters++;       
329         recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
330         
331         AliDebug(9, Form("Number of cells per cluster (max is 9!): %d", clusterDigitList.GetEntries()));
332         for (Int_t idig = 0; idig < clusterDigitList.GetEntries(); idig++)
333         {
334           digit = (AliEMCALDigit*)clusterDigitList.At(idig);
335           AliDebug(5, Form(" Adding digit %d", digit->GetId()));
336           // note: this way the sharing info is lost!
337           recPoint->AddDigit(*digit, digit->GetCalibAmp(), kFALSE); //Time or TimeR?
338           digitsC.Remove(digit);                  
339         }
340       }// recpoint
341     }
342     else
343     {
344       // we do not want to start clustering in the same spot!
345       // but in this case we may NOT reuse this seed for another cluster!
346       // need a better bookeeping?
347       digitsC.Remove(pMaxEnergyDigit);
348     }
349     
350     AliDebug (2, Form("Number of digits left: %d", digitsC.GetEntries()));      
351   } // while ! done 
352   
353   AliDebug(1,Form("total no of clusters %d from %d digits",fNumberOfECAClusters,fDigitsArr->GetEntriesFast())); 
354 }