Adding extra check for GPU_FORCE_64BIT_PTR env var
[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 AliEMCALClusterizerv1::AliEMCALClusterizerv1(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped)
75 : AliEMCALClusterizer(geometry, calib, caloped)
76 {
77   // ctor, geometry and calibration are initialized elsewhere.
78 }
79
80 //____________________________________________________________________________
81   AliEMCALClusterizerv1::~AliEMCALClusterizerv1()
82 {
83   // dtor
84 }
85
86 //____________________________________________________________________________
87 void AliEMCALClusterizerv1::Digits2Clusters(Option_t * option)
88 {
89   // Steering method to perform clusterization for the current event 
90   // in AliEMCALLoader
91   
92   if(strstr(option,"tim"))
93     gBenchmark->Start("EMCALClusterizer"); 
94   
95   if(strstr(option,"print"))
96     Print(""); 
97   
98   //Get calibration parameters from file or digitizer default values.
99   GetCalibrationParameters();
100   
101   //Get dead channel map from file or digitizer default values.
102   GetCaloCalibPedestal();
103         
104   fNumberOfECAClusters = 0;
105   
106   MakeClusters();  //only the real clusters
107   
108   if(fToUnfold){
109     fClusterUnfolding->SetInput(fNumberOfECAClusters,fRecPoints,fDigitsArr);
110     fClusterUnfolding->MakeUnfolding();
111   }
112     
113   //Evaluate position, dispersion and other RecPoint properties for EC section 
114   Int_t index;
115   for(index = 0; index < fRecPoints->GetEntries(); index++) {
116     AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
117     if(rp){
118       rp->EvalAll(fECAW0,fDigitsArr,fJustClusters);
119       //For each rec.point set the distance to the nearest bad crystal
120       if (fCaloPed)
121         rp->EvalDistanceToBadChannels(fCaloPed);
122     }
123     else AliFatal("Null rec point in list!");
124   }
125   
126   fRecPoints->Sort();
127   
128   for(index = 0; index < fRecPoints->GetEntries(); index++) {
129     AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
130     if(rp){
131       rp->SetIndexInList(index);
132       rp->Print();
133     }
134     else AliFatal("Null rec point in list!");
135   }
136   
137   if (fTreeR)
138     fTreeR->Fill();
139   
140   if(strstr(option,"deb") || strstr(option,"all"))  
141     PrintRecPoints(option);
142   
143   AliDebug(1,Form("EMCAL Clusterizer found %d Rec Points",fRecPoints->GetEntriesFast()));
144   
145   if(strstr(option,"tim")){
146     gBenchmark->Stop("EMCALClusterizer");
147     printf("Exec took %f seconds for Clusterizing", 
148            gBenchmark->GetCpuTime("EMCALClusterizer"));
149   }    
150 }
151
152 //____________________________________________________________________________
153 Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2, Bool_t & shared) const
154 {
155   // Gives the neighbourness of two digits = 0 are not neighbour; continue searching 
156   //                                       = 1 are neighbour
157   //                                       = 2 is in different SM; continue searching 
158   // In case it is in different SM, but same phi rack, check if neigbours at eta=0
159   // neighbours are defined as digits having at least a common side 
160   // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster 
161   //                                      which is compared to a digit (d2)  not yet in a cluster  
162   
163   Int_t nSupMod1=0, nModule1=0, nIphi1=0, nIeta1=0, iphi1=0, ieta1=0;
164   Int_t nSupMod2=0, nModule2=0, nIphi2=0, nIeta2=0, iphi2=0, ieta2=0;
165   
166   shared = kFALSE;
167   
168   fGeom->GetCellIndex(d1->GetId(), nSupMod1,nModule1,nIphi1,nIeta1);
169   fGeom->GetCellIndex(d2->GetId(), nSupMod2,nModule2,nIphi2,nIeta2);
170   fGeom->GetCellPhiEtaIndexInSModule(nSupMod1,nModule1,nIphi1,nIeta1, iphi1,ieta1);
171   fGeom->GetCellPhiEtaIndexInSModule(nSupMod2,nModule2,nIphi2,nIeta2, iphi2,ieta2);
172   
173   //If different SM, check if they are in the same phi, then consider cells close to eta=0 as neighbours; May 2010
174   if (nSupMod1 != nSupMod2 ) {
175     //Check if the 2 SM are in the same PHI position (0,1), (2,3), ...
176     Float_t smPhi1 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod1);
177     Float_t smPhi2 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod2);
178     
179     if(!TMath::AreEqualAbs(smPhi1, smPhi2, 1e-3)) return 2; //Not phi rack equal, not neighbours
180     
181     // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
182     // C Side impair SM, nSupMod%2=1; A side pair SM nSupMod%2=0
183     if(nSupMod1%2) ieta1+=AliEMCALGeoParams::fgkEMCALCols;
184     else           ieta2+=AliEMCALGeoParams::fgkEMCALCols;
185     
186     shared = kTRUE; // maybe a shared cluster, we know this later, set it for the moment.
187   } //Different SM, same phi
188   
189   Int_t rowdiff = TMath::Abs(iphi1 - iphi2);  
190   Int_t coldiff = TMath::Abs(ieta1 - ieta2);  
191   
192   // neighbours with at least common side; May 11, 2007
193   if ((coldiff==0 && TMath::Abs(rowdiff)==1) || (rowdiff==0 && TMath::Abs(coldiff)==1)) {  
194     //Diagonal?
195     //if ((coldiff==0 && TMath::Abs(rowdiff==1)) || (rowdiff==0 && TMath::Abs(coldiff==1)) || (TMath::Abs(rowdiff)==1 && TMath::Abs(coldiff==1))) rv = 1;
196     
197     if (gDebug == 2) 
198       printf("AliEMCALClusterizerv1::AreNeighbours(): id1=%d, (row %d, col %d) ; id2=%d, (row %d, col %d), shared %d \n",
199              d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2, shared);   
200     return 1;
201   } //Neighbours
202   else {
203     shared = kFALSE;
204     return 2; 
205   } //Not neighbours
206 }
207
208 //____________________________________________________________________________
209 void AliEMCALClusterizerv1::MakeClusters()
210 {
211   // Steering method to construct the clusters stored in a list of Reconstructed Points
212   // A cluster is defined as a list of neighbour digits
213   // Mar 03, 2007 by PAI
214   
215   if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader");
216   
217   fRecPoints->Delete();
218   
219   // Set up TObjArray with pointers to digits to work on calibrated digits 
220   TObjArray *digitsC = new TObjArray();
221   AliEMCALDigit *digit;
222   Float_t dEnergyCalibrated = 0.0, ehs = 0.0, time = 0.0;
223   TIter nextdigit(fDigitsArr);
224   while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigit())) ) { // calibrate and clean up digits
225     dEnergyCalibrated =  digit->GetAmplitude();
226     time              =  digit->GetTime();
227     Calibrate(dEnergyCalibrated,time,digit->GetId());
228     digit->SetCalibAmp(dEnergyCalibrated);
229     digit->SetTime(time);
230
231     if ( dEnergyCalibrated < fMinECut || time > fTimeMax || time < fTimeMin ){
232       continue;
233     }
234     else if (!fGeom->CheckAbsCellId(digit->GetId()))
235       continue;
236     else{
237       ehs += dEnergyCalibrated;
238       digitsC->AddLast(digit);
239     }
240   } 
241   
242   AliDebug(1,Form("MakeClusters: Number of digits %d  -> (e %f), ehs %f\n",
243                   fDigitsArr->GetEntries(),fMinECut,ehs));
244   
245   TIter nextdigitC(digitsC);
246   while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigitC())) ) { // scan over the list of digitsC
247     TArrayI clusterECAdigitslist(fDigitsArr->GetEntries());
248     dEnergyCalibrated = digit->GetCalibAmp();
249     time              = digit->GetTime();
250     if(fGeom->CheckAbsCellId(digit->GetId()) && ( dEnergyCalibrated > fECAClusteringThreshold  ) ){
251       // start a new Tower RecPoint
252       if(fNumberOfECAClusters >= fRecPoints->GetSize()) fRecPoints->Expand(2*fNumberOfECAClusters+1);
253       
254       AliEMCALRecPoint *recPoint = new  AliEMCALRecPoint(""); 
255       fRecPoints->AddAt(recPoint, fNumberOfECAClusters);
256       recPoint = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(fNumberOfECAClusters)); 
257       if (recPoint) {
258         fNumberOfECAClusters++; 
259         
260         recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
261         recPoint->AddDigit(*digit, digit->GetCalibAmp(), 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(), dEnergyCalibrated, fECAClusteringThreshold));  //Time or TimeR?
267       
268         // Grow cluster by finding neighbours
269         TIter nextClusterDigit(&clusterDigits);
270         
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             //Do not add digits with too different time 
276             Bool_t shared = kFALSE;//cluster shared by 2 SuperModules?
277             if(TMath::Abs(time - digitN->GetTime()) > fTimeCut ) continue; //Time or TimeR?
278             if (AreNeighbours(digit, digitN, shared)==1) {      // call (digit,digitN) in THAT order !!!!! 
279               recPoint->AddDigit(*digitN, digitN->GetCalibAmp(), shared); 
280               clusterDigits.AddLast(digitN); 
281               digitsC->Remove(digitN); 
282             } // if(ineb==1)
283           } // scan over digits
284         } // scan over digits already in cluster
285         
286         AliDebug(2,Form("MakeClusters: %d digitd, energy %f \n", clusterDigits.GetEntries(), recPoint->GetEnergy())); 
287       }//recpoint
288       else AliFatal("Null recpoint in array!");
289     } // If seed found
290   } // while digit 
291   
292   delete digitsC;
293   
294   AliDebug(1,Form("total no of clusters %d from %d digits",fNumberOfECAClusters,fDigitsArr->GetEntriesFast())); 
295 }