41998b0662b5befc2b5ae5254b7d4c5b952e307f
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALClusterizerv1.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 /* $Id$ */
17
18 //-- Author: Yves Schutz (SUBATECH)  & Dmitri Peressounko (SUBATECH & Kurchatov Institute)
19 //--         Gustavo Conesa (LPSC-Grenoble), move common clusterizer functionalities to mother class
20 //////////////////////////////////////////////////////////////////////////////
21 //  Clusterization class. Performs clusterization (collects neighbouring active cells) and 
22 //  unfolds the clusters having several local maxima.  
23 //  Results are stored in TreeR#, branches EMCALTowerRP (EMC recPoints),
24 //  EMCALPreShoRP (CPV RecPoints) and AliEMCALClusterizer (Clusterizer with all 
25 //  parameters including input digits branch title, thresholds etc.)
26 //
27
28 // --- ROOT system ---
29
30 #include <TFile.h> 
31 #include <TMath.h> 
32 #include <TMinuit.h>
33 #include <TTree.h> 
34 #include <TBenchmark.h>
35 #include <TBrowser.h>
36 #include <TROOT.h>
37 #include <TList.h>
38 #include <TClonesArray.h>
39
40 // --- Standard library ---
41 #include <cassert>
42
43 // --- AliRoot header files ---
44 #include "AliLog.h"
45 #include "AliEMCALClusterizerv1.h"
46 #include "AliEMCALRecPoint.h"
47 #include "AliEMCALDigit.h"
48 #include "AliEMCALGeometry.h"
49 #include "AliCaloCalibPedestal.h"
50 #include "AliEMCALCalibData.h"
51 #include "AliESDCaloCluster.h"
52 #include "AliEMCALUnfolding.h"
53
54 ClassImp(AliEMCALClusterizerv1)
55
56 //____________________________________________________________________________
57 AliEMCALClusterizerv1::AliEMCALClusterizerv1(): AliEMCALClusterizer()
58 {
59   // ctor with the indication of the file where header Tree and digits Tree are stored
60   
61   Init() ;
62 }
63
64 //____________________________________________________________________________
65 AliEMCALClusterizerv1::AliEMCALClusterizerv1(AliEMCALGeometry* geometry)
66   : AliEMCALClusterizer(geometry)
67 {
68   // ctor with the indication of the file where header Tree and digits Tree are stored
69   // use this contructor to avoid usage of Init() which uses runloader
70   // change needed by HLT - MP
71
72 }
73
74 //____________________________________________________________________________
75 AliEMCALClusterizerv1::AliEMCALClusterizerv1(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped)
76 : AliEMCALClusterizer(geometry, calib, caloped)
77 {
78         // ctor, geometry and calibration are initialized elsewhere.
79                                 
80 }
81
82
83 //____________________________________________________________________________
84   AliEMCALClusterizerv1::~AliEMCALClusterizerv1()
85 {
86   // dtor
87 }
88
89 //____________________________________________________________________________
90 void AliEMCALClusterizerv1::Digits2Clusters(Option_t * option)
91 {
92   // Steering method to perform clusterization for the current event 
93   // in AliEMCALLoader
94   
95   if(strstr(option,"tim"))
96     gBenchmark->Start("EMCALClusterizer"); 
97   
98   if(strstr(option,"print"))
99     Print("") ; 
100   
101   //Get calibration parameters from file or digitizer default values.
102   GetCalibrationParameters() ;
103   
104   //Get dead channel map from file or digitizer default values.
105   GetCaloCalibPedestal() ;
106         
107   fNumberOfECAClusters = 0;
108   
109   MakeClusters() ;  //only the real clusters
110   
111   if(fToUnfold){
112     fClusterUnfolding->SetInput(fNumberOfECAClusters,fRecPoints,fDigitsArr);
113     fClusterUnfolding->MakeUnfolding();
114   }
115     
116   //Evaluate position, dispersion and other RecPoint properties for EC section 
117   Int_t index ;
118   for(index = 0; index < fRecPoints->GetEntries(); index++) {
119     AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
120     if(rp){
121       rp->EvalAll(fECAW0,fDigitsArr) ;
122       //For each rec.point set the distance to the nearest bad crystal
123             rp->EvalDistanceToBadChannels(fCaloPed);
124     }
125     else AliFatal("Null rec point in list!");
126   }
127   
128   fRecPoints->Sort() ;
129   
130   for(index = 0; index < fRecPoints->GetEntries(); index++) {
131     AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
132     if(rp){
133       rp->SetIndexInList(index) ;
134       rp->Print();
135     }
136     else AliFatal("Null rec point in list!");
137   }
138   
139   fTreeR->Fill();
140   
141   if(strstr(option,"deb") || strstr(option,"all"))  
142     PrintRecPoints(option) ;
143   
144   AliDebug(1,Form("EMCAL Clusterizer found %d Rec Points",fRecPoints->GetEntriesFast()));
145   
146   fRecPoints->Delete();
147   
148   if(strstr(option,"tim")){
149     gBenchmark->Stop("EMCALClusterizer");
150     printf("Exec took %f seconds for Clusterizing", 
151            gBenchmark->GetCpuTime("EMCALClusterizer"));
152   }    
153 }
154
155 //____________________________________________________________________________
156 Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2, Bool_t & shared) const
157 {
158   // Gives the neighbourness of two digits = 0 are not neighbour ; continue searching 
159   //                                       = 1 are neighbour
160   //                                       = 2 is in different SM; continue searching 
161   // In case it is in different SM, but same phi rack, check if neigbours at eta=0
162   // neighbours are defined as digits having at least a common side 
163   // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster 
164   //                                      which is compared to a digit (d2)  not yet in a cluster  
165   
166   static Int_t nSupMod1=0, nModule1=0, nIphi1=0, nIeta1=0, iphi1=0, ieta1=0;
167   static Int_t nSupMod2=0, nModule2=0, nIphi2=0, nIeta2=0, iphi2=0, ieta2=0;
168   
169   shared = kFALSE;
170   
171   fGeom->GetCellIndex(d1->GetId(), nSupMod1,nModule1,nIphi1,nIeta1);
172   fGeom->GetCellIndex(d2->GetId(), nSupMod2,nModule2,nIphi2,nIeta2);
173   fGeom->GetCellPhiEtaIndexInSModule(nSupMod1,nModule1,nIphi1,nIeta1, iphi1,ieta1);
174   fGeom->GetCellPhiEtaIndexInSModule(nSupMod2,nModule2,nIphi2,nIeta2, iphi2,ieta2);
175   
176   //If different SM, check if they are in the same phi, then consider cells close to eta=0 as neighbours; May 2010
177   if(nSupMod1 != nSupMod2 ) {
178     //Check if the 2 SM are in the same PHI position (0,1), (2,3), ...
179     Float_t smPhi1 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod1);
180     Float_t smPhi2 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod2);
181     
182     if(!TMath::AreEqualAbs(smPhi1, smPhi2, 1e-3)) return 2; //Not phi rack equal, not neighbours
183     
184     // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
185     // C Side impair SM, nSupMod%2=1; A side pair SM nSupMod%2=0
186     if(nSupMod1%2) ieta1+=AliEMCALGeoParams::fgkEMCALCols;
187     else           ieta2+=AliEMCALGeoParams::fgkEMCALCols;
188     
189     shared = kTRUE; // maybe a shared cluster, we know this later, set it for the moment.
190     
191   }//Different SM, same phi
192   
193   Int_t rowdiff = TMath::Abs(iphi1 - iphi2);  
194   Int_t coldiff = TMath::Abs(ieta1 - ieta2) ;  
195   
196   // neighbours with at least common side; May 11, 2007
197   if ((coldiff==0 && TMath::Abs(rowdiff)==1) || (rowdiff==0 && TMath::Abs(coldiff)==1)) {  
198     //Diagonal?
199     //if ((coldiff==0 && TMath::Abs(rowdiff==1)) || (rowdiff==0 && TMath::Abs(coldiff==1)) || (TMath::Abs(rowdiff)==1 && TMath::Abs(coldiff==1))) rv = 1;
200     
201     if (gDebug == 2) 
202       printf("AliEMCALClusterizerv1::AreNeighbours(): id1=%d, (row %d, col %d) ; id2=%d, (row %d, col %d), shared %d \n",
203              d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2, shared);   
204     
205     return 1;
206   }//Neighbours
207   else {
208     shared = kFALSE;
209     return 2 ; 
210   }//Not neighbours
211 }
212
213 //____________________________________________________________________________
214 void AliEMCALClusterizerv1::MakeClusters()
215 {
216   // Steering method to construct the clusters stored in a list of Reconstructed Points
217   // A cluster is defined as a list of neighbour digits
218   // Mar 03, 2007 by PAI
219   
220   if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader");
221   
222   fRecPoints->Clear();
223   
224   // Set up TObjArray with pointers to digits to work on 
225   TObjArray *digitsC = new TObjArray();
226   TIter nextdigit(fDigitsArr);
227   AliEMCALDigit *digit;
228   while ( (digit = dynamic_cast<AliEMCALDigit*>(nextdigit())) ) {
229     digitsC->AddLast(digit);
230   }
231   
232   double e = 0.0, ehs = 0.0;
233   TIter nextdigitC(digitsC);
234   while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigitC())) ) { // clean up digits
235     e = Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId());//Time or TimeR?
236     if ( e < fMinECut) //|| digit->GetTimeR() > fTimeCut ) // time window of cell checked in calibrate
237       digitsC->Remove(digit);
238     else    
239       ehs += e;
240   } 
241   AliDebug(1,Form("MakeClusters: Number of digits %d  -> (e %f), ehs %f\n",
242                   fDigitsArr->GetEntries(),fMinECut,ehs));
243   
244   nextdigitC.Reset();
245   
246   while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigitC())) ) { // scan over the list of digitsC
247     TArrayI clusterECAdigitslist(fDigitsArr->GetEntries());
248     
249     if(fGeom->CheckAbsCellId(digit->GetId()) && (Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId()) > fECAClusteringThreshold  ) ){
250       // start a new Tower RecPoint
251       if(fNumberOfECAClusters >= fRecPoints->GetSize()) fRecPoints->Expand(2*fNumberOfECAClusters+1) ;
252       
253       AliEMCALRecPoint *recPoint = new  AliEMCALRecPoint("") ; 
254       fRecPoints->AddAt(recPoint, fNumberOfECAClusters) ;
255       recPoint = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(fNumberOfECAClusters)) ; 
256       if(recPoint){
257         fNumberOfECAClusters++ ; 
258         
259         recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
260         
261         recPoint->AddDigit(*digit, Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId()),kFALSE) ; //Time or TimeR?
262         TObjArray clusterDigits;
263         clusterDigits.AddLast(digit);   
264         digitsC->Remove(digit) ; 
265         
266         AliDebug(1,Form("MakeClusters: OK id = %d, ene = %f , cell.th. = %f \n", digit->GetId(),
267                         Calibrate(digit->GetAmplitude(),digit->GetTime(),digit->GetId()), fECAClusteringThreshold));  //Time or TimeR?
268         Float_t time = digit->GetTime();//Time or TimeR?
269         // Grow cluster by finding neighbours
270         TIter nextClusterDigit(&clusterDigits);
271         while ( (digit = dynamic_cast<AliEMCALDigit*>(nextClusterDigit())) ) { // scan over digits in cluster 
272           TIter nextdigitN(digitsC); 
273           AliEMCALDigit *digitN = 0; // digi neighbor
274           while ( (digitN = (AliEMCALDigit *)nextdigitN()) ) { // scan over all digits to look for neighbours
275             
276             //Do not add digits with too different time 
277             Bool_t shared = kFALSE;//cluster shared by 2 SuperModules?
278             if(TMath::Abs(time - digitN->GetTime()) > fTimeCut ) continue; //Time or TimeR?
279             if (AreNeighbours(digit, digitN, shared)==1) {      // call (digit,digitN) in THAT order !!!!! 
280               recPoint->AddDigit(*digitN, Calibrate(digitN->GetAmplitude(), digitN->GetTime(), digitN->GetId()),shared) ;//Time or TimeR?
281               clusterDigits.AddLast(digitN) ; 
282               digitsC->Remove(digitN) ; 
283             } // if(ineb==1)
284           } // scan over digits
285         } // scan over digits already in cluster
286         
287         AliDebug(2,Form("MakeClusters: %d digitd, energy %f \n", clusterDigits.GetEntries(), recPoint->GetEnergy())); 
288       }//recpoint
289       else AliFatal("Null recpoint in array!");
290     } // If seed found
291   } // while digit 
292   
293   delete digitsC ;
294   
295   AliDebug(1,Form("total no of clusters %d from %d digits",fNumberOfECAClusters,fDigitsArr->GetEntriesFast())); 
296 }
297