correct mask for V0 charge decoding in STU payload
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALClusterizerv2.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 "AliEMCALClusterizerv2.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(AliEMCALClusterizerv2)
55
56 //____________________________________________________________________________
57 AliEMCALClusterizerv2::AliEMCALClusterizerv2() 
58   : AliEMCALClusterizerv1(), fDoEnGradCut(1)
59 {
60   // ctor with the indication of the file where header Tree and digits Tree are stored
61 }
62
63 //____________________________________________________________________________
64 AliEMCALClusterizerv2::AliEMCALClusterizerv2(AliEMCALGeometry* geometry)
65   : AliEMCALClusterizerv1(geometry), fDoEnGradCut(1)
66 {
67   // ctor with the indication of the file where header Tree and digits Tree are stored
68   // use this contructor to avoid usage of Init() which uses runloader
69   // change needed by HLT - MP
70 }
71
72 //____________________________________________________________________________
73 AliEMCALClusterizerv2::AliEMCALClusterizerv2(AliEMCALGeometry* geometry, AliEMCALCalibData* calib, AliCaloCalibPedestal* caloped)
74   : AliEMCALClusterizerv1(geometry, calib, caloped), fDoEnGradCut(1)
75 {
76   // ctor, geometry and calibration are initialized elsewhere.
77 }
78
79 //____________________________________________________________________________
80 AliEMCALClusterizerv2::~AliEMCALClusterizerv2()
81 {
82   // dtor
83 }
84
85 //____________________________________________________________________________
86 Int_t AliEMCALClusterizerv2::AreNeighbours(AliEMCALDigit* d1, AliEMCALDigit* d2, Bool_t& shared) const
87
88   // Gives the neighbourness of two digits = 0 are not neighbour; continue searching 
89   //                                       = 1 are neighbour
90   //                                       = 2 is in different SM; continue searching 
91   // In case it is in different SM, but same phi rack, check if neigbours at eta=0
92   // neighbours are defined as digits having at least a common side 
93   // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster 
94   //                                      which is compared to a digit (d2)  not yet in a cluster  
95   
96   if (fDoEnGradCut) {
97     if (d2->GetCalibAmp()>d1->GetCalibAmp())
98       return 3; // energy of neighboring cell should be smaller in order to become a neighbor
99   }
100   return AliEMCALClusterizerv1::AreNeighbours(d1,d2,shared);
101 }
102
103 //____________________________________________________________________________
104 void AliEMCALClusterizerv2::MakeClusters()
105 {
106   // Make list of clusters. Start from highest energy cell.
107   
108   if (fGeom==0) 
109     AliFatal("Did not get geometry from EMCALLoader");
110   
111   fRecPoints->Delete();
112
113   // set up TObjArray with pointers to digits to work on, calibrate digits 
114   TObjArray digitsC;
115   Double_t ehs = 0.0;
116   AliEMCALDigit *digit = 0;
117   TIter nextdigit(fDigitsArr);
118   while ( (digit = static_cast<AliEMCALDigit*>(nextdigit())) ) {
119     Float_t dEnergyCalibrated = digit->GetAmplitude();
120     Float_t time              = digit->GetTime();
121     Calibrate(dEnergyCalibrated, time ,digit->GetId());
122     digit->SetCalibAmp(dEnergyCalibrated);
123     digit->SetTime(time);
124     if (dEnergyCalibrated < fMinECut) 
125       continue;
126     if (!fGeom->CheckAbsCellId(digit->GetId()))
127       continue;
128     ehs += dEnergyCalibrated;
129     digitsC.AddLast(digit);
130   }
131
132   AliDebug(1,Form("MakeClusters: Number of digits %d  -> ehs %f (minE %f)\n",
133                   fDigitsArr->GetEntries(),ehs,fMinECut));
134
135   TIter nextdigitC(&digitsC);
136   while (1) {
137     Int_t   iMaxEnergyDigit = -1;
138     Float_t dMaxEnergyDigit = -1;
139     AliEMCALDigit *pMaxEnergyDigit = 0;
140     nextdigitC.Reset();
141     while ( (digit = static_cast<AliEMCALDigit *>(nextdigitC())) ) { 
142       Float_t dEnergyCalibrated = digit->GetCalibAmp();
143       if (dEnergyCalibrated>fECAClusteringThreshold && dEnergyCalibrated>dMaxEnergyDigit) {
144         dMaxEnergyDigit = dEnergyCalibrated;
145         iMaxEnergyDigit = digit->GetId();
146         pMaxEnergyDigit = digit;
147       }
148     }
149     if (iMaxEnergyDigit<0 || digitsC.GetEntries() <= 0) {
150       break;
151     }
152
153     if (fNumberOfECAClusters>=fRecPoints->GetSize()) 
154       fRecPoints->Expand(2*fNumberOfECAClusters+1);
155
156     AliEMCALRecPoint *recPoint = new  AliEMCALRecPoint(""); 
157     recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
158     recPoint->AddDigit(*pMaxEnergyDigit, pMaxEnergyDigit->GetCalibAmp(), kFALSE);
159     fRecPoints->AddAt(recPoint, fNumberOfECAClusters++);
160     digitsC.Remove(pMaxEnergyDigit); 
161     TObjArray clusterDigits;
162     clusterDigits.AddLast(pMaxEnergyDigit);     
163     TIter nextClusterDigit(&clusterDigits);
164     Float_t time = pMaxEnergyDigit->GetTime(); 
165
166     AliDebug(1,Form("MakeClusters: Max digit found id = %d, ene = %f , clus.th. = %f \n", 
167                     iMaxEnergyDigit, dMaxEnergyDigit, fECAClusteringThreshold));
168
169     while ( (digit = static_cast<AliEMCALDigit*>(nextClusterDigit())) ) { // scan over digits in cluster 
170       TIter nextdigitN(&digitsC); 
171       AliEMCALDigit *digitN = 0; // digi neighbor
172       while ( (digitN = static_cast<AliEMCALDigit*>(nextdigitN())) ) { // scan over all digits to look for neighbours
173         //Do not add digits with too different time 
174         if (TMath::Abs(time - digitN->GetTime()) > fTimeCut ) 
175           continue;
176         Bool_t shared = kFALSE; //cluster shared by 2 SuperModules?
177         if (AreNeighbours(digit, digitN, shared)==1) {
178           recPoint->AddDigit(*digitN, digitN->GetCalibAmp(), shared);
179           clusterDigits.AddLast(digitN); 
180           digitsC.Remove(digitN); 
181         } 
182       }
183     }
184     AliDebug(2,Form("MakeClusters: %d digitd, energy %f \n", clusterDigits.GetEntries(), recPoint->GetEnergy())); 
185   }
186   AliDebug(1,Form("total no of clusters %d from %d digits",fNumberOfECAClusters,fDigitsArr->GetEntriesFast())); 
187 }
188