]>
Commit | Line | Data |
---|---|---|
483b0559 | 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 | **************************************************************************/ | |
173558f2 | 15 | |
483b0559 | 16 | /* $Id$ */ |
803d1ab0 | 17 | |
3a8be91c | 18 | //-- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko (SUBATECH & Kurchatov Institute) |
ee08edde | 19 | //-- Gustavo Conesa (LPSC-Grenoble), move common clusterizer functionalities to mother class |
483b0559 | 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.) | |
ee08edde | 26 | // |
483b0559 | 27 | |
28 | // --- ROOT system --- | |
29 | ||
e52475ed | 30 | #include <TFile.h> |
e52475ed | 31 | #include <TMath.h> |
32 | #include <TMinuit.h> | |
33 | #include <TTree.h> | |
e52475ed | 34 | #include <TBenchmark.h> |
35 | #include <TBrowser.h> | |
a1e17193 | 36 | #include <TROOT.h> |
ee08edde | 37 | #include <TList.h> |
38 | #include <TClonesArray.h> | |
1d59832c | 39 | |
483b0559 | 40 | // --- Standard library --- |
ee08edde | 41 | #include <cassert> |
173558f2 | 42 | |
483b0559 | 43 | // --- AliRoot header files --- |
ee08edde | 44 | #include "AliLog.h" |
483b0559 | 45 | #include "AliEMCALClusterizerv1.h" |
70a93198 | 46 | #include "AliEMCALRecPoint.h" |
483b0559 | 47 | #include "AliEMCALDigit.h" |
05a92d59 | 48 | #include "AliEMCALGeometry.h" |
40164976 | 49 | #include "AliCaloCalibPedestal.h" |
a435f763 | 50 | #include "AliEMCALCalibData.h" |
ee08edde | 51 | #include "AliESDCaloCluster.h" |
483b0559 | 52 | |
53 | ClassImp(AliEMCALClusterizerv1) | |
1963b290 | 54 | |
483b0559 | 55 | //____________________________________________________________________________ |
ee08edde | 56 | AliEMCALClusterizerv1::AliEMCALClusterizerv1(): AliEMCALClusterizer() |
483b0559 | 57 | { |
58 | // ctor with the indication of the file where header Tree and digits Tree are stored | |
59 | ||
483b0559 | 60 | Init() ; |
483b0559 | 61 | } |
05a92d59 | 62 | |
5544799a | 63 | //____________________________________________________________________________ |
64 | AliEMCALClusterizerv1::AliEMCALClusterizerv1(AliEMCALGeometry* geometry) | |
ee08edde | 65 | : AliEMCALClusterizer(geometry) |
5544799a | 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 | ||
5544799a | 71 | } |
72 | ||
0c5b726e | 73 | //____________________________________________________________________________ |
40164976 | 74 | AliEMCALClusterizerv1::AliEMCALClusterizerv1(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped) |
ee08edde | 75 | : AliEMCALClusterizer(geometry, calib, caloped) |
0c5b726e | 76 | { |
77 | // ctor, geometry and calibration are initialized elsewhere. | |
ee08edde | 78 | |
0c5b726e | 79 | } |
80 | ||
81 | ||
483b0559 | 82 | //____________________________________________________________________________ |
83 | AliEMCALClusterizerv1::~AliEMCALClusterizerv1() | |
84 | { | |
ef305168 | 85 | // dtor |
ef305168 | 86 | } |
87 | ||
483b0559 | 88 | //____________________________________________________________________________ |
c47157cd | 89 | void AliEMCALClusterizerv1::Digits2Clusters(Option_t * option) |
483b0559 | 90 | { |
98e9578e | 91 | // Steering method to perform clusterization for the current event |
92 | // in AliEMCALLoader | |
51e4198e | 93 | |
483b0559 | 94 | if(strstr(option,"tim")) |
95 | gBenchmark->Start("EMCALClusterizer"); | |
96 | ||
97 | if(strstr(option,"print")) | |
98 | Print("") ; | |
51e4198e | 99 | |
1bd98442 | 100 | //Get calibration parameters from file or digitizer default values. |
101 | GetCalibrationParameters() ; | |
51e4198e | 102 | |
40164976 | 103 | //Get dead channel map from file or digitizer default values. |
104 | GetCaloCalibPedestal() ; | |
105 | ||
98e9578e | 106 | fNumberOfECAClusters = 0; |
51e4198e | 107 | |
0e7c6655 | 108 | MakeClusters() ; //only the real clusters |
51e4198e | 109 | |
98e9578e | 110 | if(fToUnfold) |
111 | MakeUnfolding() ; | |
51e4198e | 112 | |
bbac7e6f | 113 | Int_t index ; |
51e4198e | 114 | |
bbac7e6f | 115 | //Evaluate position, dispersion and other RecPoint properties for EC section |
116 | for(index = 0; index < fRecPoints->GetEntries(); index++) { | |
51e4198e | 117 | AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index)); |
118 | if(rp){ | |
119 | rp->EvalAll(fECAW0,fDigitsArr) ; | |
120 | //For each rec.point set the distance to the nearest bad crystal | |
121 | rp->EvalDistanceToBadChannels(fCaloPed); | |
122 | } | |
123 | else AliFatal("Null rec point in list!"); | |
bbac7e6f | 124 | } |
51e4198e | 125 | |
bbac7e6f | 126 | fRecPoints->Sort() ; |
51e4198e | 127 | |
bbac7e6f | 128 | for(index = 0; index < fRecPoints->GetEntries(); index++) { |
51e4198e | 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!"); | |
bbac7e6f | 135 | } |
51e4198e | 136 | |
c47157cd | 137 | fTreeR->Fill(); |
138 | ||
98e9578e | 139 | if(strstr(option,"deb") || strstr(option,"all")) |
140 | PrintRecPoints(option) ; | |
51e4198e | 141 | |
c47157cd | 142 | AliDebug(1,Form("EMCAL Clusterizer found %d Rec Points",fRecPoints->GetEntriesFast())); |
51e4198e | 143 | |
3f9ad99f | 144 | fRecPoints->Delete(); |
51e4198e | 145 | |
483b0559 | 146 | if(strstr(option,"tim")){ |
147 | gBenchmark->Stop("EMCALClusterizer"); | |
98e9578e | 148 | printf("Exec took %f seconds for Clusterizing", |
51e4198e | 149 | gBenchmark->GetCpuTime("EMCALClusterizer")); |
98e9578e | 150 | } |
483b0559 | 151 | } |
152 | ||
153 | //____________________________________________________________________________ | |
d804d556 | 154 | Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALRecPoint * recPoint, AliEMCALDigit ** maxAt, |
155 | const Float_t* maxAtEnergy, | |
225cd96d | 156 | Int_t nPar, Float_t * fitparameters) const |
157 | { | |
158 | // Calls TMinuit to fit the energy distribution of a cluster with several maxima | |
159 | // The initial values for fitting procedure are set equal to the | |
160 | // positions of local maxima. | |
161 | // Cluster will be fitted as a superposition of nPar/3 | |
162 | // electromagnetic showers | |
163 | ||
164 | if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader"); | |
8de94e0f | 165 | |
166 | if(!gMinuit) | |
167 | gMinuit = new TMinuit(100) ; | |
483b0559 | 168 | |
483b0559 | 169 | gMinuit->mncler(); // Reset Minuit's list of paramters |
170 | gMinuit->SetPrintLevel(-1) ; // No Printout | |
225cd96d | 171 | gMinuit->SetFCN(AliEMCALClusterizerv1::UnfoldingChiSquare) ; |
172 | // To set the address of the minimization function | |
483b0559 | 173 | TList * toMinuit = new TList(); |
d804d556 | 174 | toMinuit->AddAt(recPoint,0) ; |
c47157cd | 175 | toMinuit->AddAt(fDigitsArr,1) ; |
225cd96d | 176 | toMinuit->AddAt(fGeom,2) ; |
177 | ||
483b0559 | 178 | gMinuit->SetObjectFit(toMinuit) ; // To tranfer pointer to UnfoldingChiSquare |
179 | ||
180 | // filling initial values for fit parameters | |
181 | AliEMCALDigit * digit ; | |
182 | ||
225cd96d | 183 | Int_t ierflg = 0; |
483b0559 | 184 | Int_t index = 0 ; |
185 | Int_t nDigits = (Int_t) nPar / 3 ; | |
186 | ||
187 | Int_t iDigit ; | |
188 | ||
483b0559 | 189 | for(iDigit = 0; iDigit < nDigits; iDigit++){ |
225cd96d | 190 | digit = maxAt[iDigit]; |
191 | Double_t x = 0.; | |
192 | Double_t y = 0.; | |
193 | Double_t z = 0.; | |
483b0559 | 194 | |
225cd96d | 195 | fGeom->RelPosCellInSModule(digit->GetId(), y, x, z); |
483b0559 | 196 | |
197 | Float_t energy = maxAtEnergy[iDigit] ; | |
198 | ||
199 | gMinuit->mnparm(index, "x", x, 0.1, 0, 0, ierflg) ; | |
225cd96d | 200 | index++ ; |
201 | if(ierflg != 0){ | |
202 | Error("FindFit", "EMCAL Unfolding Unable to set initial value for fit procedure : x = %f", x ) ; | |
483b0559 | 203 | return kFALSE; |
204 | } | |
205 | gMinuit->mnparm(index, "z", z, 0.1, 0, 0, ierflg) ; | |
225cd96d | 206 | index++ ; |
483b0559 | 207 | if(ierflg != 0){ |
225cd96d | 208 | Error("FindFit", "EMCAL Unfolding Unable to set initial value for fit procedure : z = %f", z) ; |
483b0559 | 209 | return kFALSE; |
210 | } | |
211 | gMinuit->mnparm(index, "Energy", energy , 0.05*energy, 0., 4.*energy, ierflg) ; | |
225cd96d | 212 | index++ ; |
483b0559 | 213 | if(ierflg != 0){ |
225cd96d | 214 | Error("FindFit", "EMCAL Unfolding Unable to set initial value for fit procedure : energy = %f", energy) ; |
483b0559 | 215 | return kFALSE; |
216 | } | |
217 | } | |
218 | ||
225cd96d | 219 | Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; |
220 | // The number of function call slightly depends on it. | |
cb1efdae | 221 | //Double_t p1 = 1.0 ; |
483b0559 | 222 | Double_t p2 = 0.0 ; |
223 | ||
225cd96d | 224 | gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ; // force TMinuit to reduce function calls |
0561e246 | 225 | // gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient |
483b0559 | 226 | gMinuit->SetMaxIterations(5); |
227 | gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ; // No Warnings | |
225cd96d | 228 | gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize |
483b0559 | 229 | |
225cd96d | 230 | if(ierflg == 4){ // Minimum not found |
231 | Error("FindFit", "EMCAL Unfolding Fit not converged, cluster abandoned " ) ; | |
483b0559 | 232 | return kFALSE ; |
225cd96d | 233 | } |
483b0559 | 234 | for(index = 0; index < nPar; index++){ |
ee08edde | 235 | Double_t err = 0. ; |
236 | Double_t val = 0. ; | |
483b0559 | 237 | gMinuit->GetParameter(index, val, err) ; // Returns value and error of parameter index |
238 | fitparameters[index] = val ; | |
225cd96d | 239 | } |
483b0559 | 240 | |
241 | delete toMinuit ; | |
242 | return kTRUE; | |
243 | ||
244 | } | |
245 | ||
483b0559 | 246 | //____________________________________________________________________________ |
25bb3dcb | 247 | Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2, Bool_t & shared) const |
483b0559 | 248 | { |
25bb3dcb | 249 | // Gives the neighbourness of two digits = 0 are not neighbour ; continue searching |
250 | // = 1 are neighbour | |
251 | // = 2 is in different SM; continue searching | |
252 | // In case it is in different SM, but same phi rack, check if neigbours at eta=0 | |
253 | // neighbours are defined as digits having at least a common side | |
254 | // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster | |
255 | // which is compared to a digit (d2) not yet in a cluster | |
256 | ||
257 | static Int_t nSupMod1=0, nModule1=0, nIphi1=0, nIeta1=0, iphi1=0, ieta1=0; | |
258 | static Int_t nSupMod2=0, nModule2=0, nIphi2=0, nIeta2=0, iphi2=0, ieta2=0; | |
829ba234 | 259 | |
25bb3dcb | 260 | shared = kFALSE; |
829ba234 | 261 | |
25bb3dcb | 262 | fGeom->GetCellIndex(d1->GetId(), nSupMod1,nModule1,nIphi1,nIeta1); |
263 | fGeom->GetCellIndex(d2->GetId(), nSupMod2,nModule2,nIphi2,nIeta2); | |
264 | fGeom->GetCellPhiEtaIndexInSModule(nSupMod1,nModule1,nIphi1,nIeta1, iphi1,ieta1); | |
265 | fGeom->GetCellPhiEtaIndexInSModule(nSupMod2,nModule2,nIphi2,nIeta2, iphi2,ieta2); | |
266 | ||
267 | //If different SM, check if they are in the same phi, then consider cells close to eta=0 as neighbours; May 2010 | |
268 | if(nSupMod1 != nSupMod2 ) { | |
269 | //Check if the 2 SM are in the same PHI position (0,1), (2,3), ... | |
270 | Float_t smPhi1 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod1); | |
271 | Float_t smPhi2 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod2); | |
272 | ||
273 | if(!TMath::AreEqualAbs(smPhi1, smPhi2, 1e-3)) return 2; //Not phi rack equal, not neighbours | |
274 | ||
275 | // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2 | |
276 | // C Side impair SM, nSupMod%2=1; A side pair SM nSupMod%2=0 | |
277 | if(nSupMod1%2) ieta1+=AliEMCALGeoParams::fgkEMCALCols; | |
278 | else ieta2+=AliEMCALGeoParams::fgkEMCALCols; | |
279 | ||
280 | shared = kTRUE; // maybe a shared cluster, we know this later, set it for the moment. | |
281 | ||
282 | }//Different SM, same phi | |
283 | ||
ee08edde | 284 | Int_t rowdiff = TMath::Abs(iphi1 - iphi2); |
285 | Int_t coldiff = TMath::Abs(ieta1 - ieta2) ; | |
25bb3dcb | 286 | |
287 | // neighbours with at least common side; May 11, 2007 | |
288 | if ((coldiff==0 && TMath::Abs(rowdiff)==1) || (rowdiff==0 && TMath::Abs(coldiff)==1)) { | |
289 | //Diagonal? | |
290 | //if ((coldiff==0 && TMath::Abs(rowdiff==1)) || (rowdiff==0 && TMath::Abs(coldiff==1)) || (TMath::Abs(rowdiff)==1 && TMath::Abs(coldiff==1))) rv = 1; | |
291 | ||
292 | if (gDebug == 2) | |
293 | printf("AliEMCALClusterizerv1::AreNeighbours(): id1=%d, (row %d, col %d) ; id2=%d, (row %d, col %d), shared %d \n", | |
294 | d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2, shared); | |
295 | ||
296 | return 1; | |
297 | }//Neighbours | |
298 | else { | |
299 | shared = kFALSE; | |
300 | return 2 ; | |
301 | }//Not neighbours | |
483b0559 | 302 | } |
303 | ||
a5c60732 | 304 | //____________________________________________________________________________ |
0e7c6655 | 305 | void AliEMCALClusterizerv1::MakeClusters() |
483b0559 | 306 | { |
307 | // Steering method to construct the clusters stored in a list of Reconstructed Points | |
308 | // A cluster is defined as a list of neighbour digits | |
f1487f22 | 309 | // Mar 03, 2007 by PAI |
51e4198e | 310 | |
f1487f22 | 311 | if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader"); |
51e4198e | 312 | |
c47157cd | 313 | fRecPoints->Clear(); |
51e4198e | 314 | |
98e9578e | 315 | // Set up TObjArray with pointers to digits to work on |
316 | TObjArray *digitsC = new TObjArray(); | |
c47157cd | 317 | TIter nextdigit(fDigitsArr); |
98e9578e | 318 | AliEMCALDigit *digit; |
319 | while ( (digit = dynamic_cast<AliEMCALDigit*>(nextdigit())) ) { | |
320 | digitsC->AddLast(digit); | |
321 | } | |
51e4198e | 322 | |
1d46d1f6 | 323 | double e = 0.0, ehs = 0.0; |
98e9578e | 324 | TIter nextdigitC(digitsC); |
1d46d1f6 | 325 | while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigitC())) ) { // clean up digits |
829ba234 | 326 | e = Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId());//Time or TimeR? |
327 | if ( e < fMinECut) //|| digit->GetTimeR() > fTimeCut ) // time window of cell checked in calibrate | |
98e9578e | 328 | digitsC->Remove(digit); |
329 | else | |
330 | ehs += e; | |
1d59832c | 331 | } |
29b7e56e | 332 | AliDebug(1,Form("MakeClusters: Number of digits %d -> (e %f), ehs %f\n", |
51e4198e | 333 | fDigitsArr->GetEntries(),fMinECut,ehs)); |
334 | ||
1d46d1f6 | 335 | nextdigitC.Reset(); |
51e4198e | 336 | |
1d46d1f6 | 337 | while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigitC())) ) { // scan over the list of digitsC |
c47157cd | 338 | TArrayI clusterECAdigitslist(fDigitsArr->GetEntries()); |
51e4198e | 339 | |
829ba234 | 340 | if(fGeom->CheckAbsCellId(digit->GetId()) && (Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId()) > fECAClusteringThreshold ) ){ |
98e9578e | 341 | // start a new Tower RecPoint |
c47157cd | 342 | if(fNumberOfECAClusters >= fRecPoints->GetSize()) fRecPoints->Expand(2*fNumberOfECAClusters+1) ; |
51e4198e | 343 | |
98e9578e | 344 | AliEMCALRecPoint *recPoint = new AliEMCALRecPoint("") ; |
c47157cd | 345 | fRecPoints->AddAt(recPoint, fNumberOfECAClusters) ; |
346 | recPoint = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(fNumberOfECAClusters)) ; | |
51e4198e | 347 | if(recPoint){ |
348 | fNumberOfECAClusters++ ; | |
349 | ||
350 | recPoint->SetClusterType(AliVCluster::kEMCALClusterv1); | |
351 | ||
352 | recPoint->AddDigit(*digit, Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId()),kFALSE) ; //Time or TimeR? | |
353 | TObjArray clusterDigits; | |
354 | clusterDigits.AddLast(digit); | |
355 | digitsC->Remove(digit) ; | |
356 | ||
357 | AliDebug(1,Form("MakeClusters: OK id = %d, ene = %f , cell.th. = %f \n", digit->GetId(), | |
358 | Calibrate(digit->GetAmplitude(),digit->GetTime(),digit->GetId()), fECAClusteringThreshold)); //Time or TimeR? | |
359 | Float_t time = digit->GetTime();//Time or TimeR? | |
360 | // Grow cluster by finding neighbours | |
361 | TIter nextClusterDigit(&clusterDigits); | |
362 | while ( (digit = dynamic_cast<AliEMCALDigit*>(nextClusterDigit())) ) { // scan over digits in cluster | |
363 | TIter nextdigitN(digitsC); | |
364 | AliEMCALDigit *digitN = 0; // digi neighbor | |
365 | while ( (digitN = (AliEMCALDigit *)nextdigitN()) ) { // scan over all digits to look for neighbours | |
366 | ||
367 | //Do not add digits with too different time | |
368 | Bool_t shared = kFALSE;//cluster shared by 2 SuperModules? | |
369 | if(TMath::Abs(time - digitN->GetTime()) > fTimeCut ) continue; //Time or TimeR? | |
370 | if (AreNeighbours(digit, digitN, shared)==1) { // call (digit,digitN) in THAT order !!!!! | |
371 | recPoint->AddDigit(*digitN, Calibrate(digitN->GetAmplitude(), digitN->GetTime(), digitN->GetId()),shared) ;//Time or TimeR? | |
372 | clusterDigits.AddLast(digitN) ; | |
373 | digitsC->Remove(digitN) ; | |
374 | } // if(ineb==1) | |
375 | } // scan over digits | |
376 | } // scan over digits already in cluster | |
377 | ||
98e9578e | 378 | AliDebug(2,Form("MakeClusters: %d digitd, energy %f \n", clusterDigits.GetEntries(), recPoint->GetEnergy())); |
51e4198e | 379 | }//recpoint |
380 | else AliFatal("Null recpoint in array!"); | |
98e9578e | 381 | } // If seed found |
1d59832c | 382 | } // while digit |
51e4198e | 383 | |
483b0559 | 384 | delete digitsC ; |
ab6a174f | 385 | |
c47157cd | 386 | AliDebug(1,Form("total no of clusters %d from %d digits",fNumberOfECAClusters,fDigitsArr->GetEntriesFast())); |
483b0559 | 387 | } |
388 | ||
7ea6391b | 389 | //____________________________________________________________________________ |
225cd96d | 390 | void AliEMCALClusterizerv1::MakeUnfolding() |
483b0559 | 391 | { |
225cd96d | 392 | // Unfolds clusters using the shape of an ElectroMagnetic shower |
393 | // Performs unfolding of all clusters | |
51e4198e | 394 | |
225cd96d | 395 | if(fNumberOfECAClusters > 0){ |
396 | if (fGeom==0) | |
397 | AliFatal("Did not get geometry from EMCALLoader") ; | |
398 | Int_t nModulesToUnfold = fGeom->GetNCells(); | |
51e4198e | 399 | |
225cd96d | 400 | Int_t numberofNotUnfolded = fNumberOfECAClusters ; |
401 | Int_t index ; | |
402 | for(index = 0 ; index < numberofNotUnfolded ; index++){ | |
51e4198e | 403 | |
d804d556 | 404 | AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(index) ) ; |
51e4198e | 405 | if(recPoint){ |
406 | TVector3 gpos; | |
407 | Int_t absId = -1; | |
408 | recPoint->GetGlobalPosition(gpos); | |
409 | fGeom->GetAbsCellIdFromEtaPhi(gpos.Eta(),gpos.Phi(),absId); | |
410 | if(absId > nModulesToUnfold) | |
411 | break ; | |
412 | ||
413 | Int_t nMultipl = recPoint->GetMultiplicity() ; | |
414 | AliEMCALDigit ** maxAt = new AliEMCALDigit*[nMultipl] ; | |
415 | Float_t * maxAtEnergy = new Float_t[nMultipl] ; | |
416 | Int_t nMax = recPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fECALocMaxCut,fDigitsArr) ; | |
417 | ||
418 | if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0 | |
419 | UnfoldCluster(recPoint, nMax, maxAt, maxAtEnergy) ; | |
420 | fRecPoints->Remove(recPoint); | |
421 | fRecPoints->Compress() ; | |
422 | index-- ; | |
423 | fNumberOfECAClusters-- ; | |
424 | numberofNotUnfolded-- ; | |
425 | } | |
426 | else{ | |
427 | recPoint->SetNExMax(1) ; //Only one local maximum | |
428 | } | |
429 | ||
430 | delete[] maxAt ; | |
431 | delete[] maxAtEnergy ; | |
225cd96d | 432 | } |
51e4198e | 433 | else AliFatal("Null recpoint in Array!"); |
225cd96d | 434 | } |
435 | } | |
436 | // End of Unfolding of clusters | |
9859bfc0 | 437 | } |
173558f2 | 438 | |
9859bfc0 | 439 | //____________________________________________________________________________ |
225cd96d | 440 | Double_t AliEMCALClusterizerv1::ShowerShape(Double_t x, Double_t y) |
9859bfc0 | 441 | { |
225cd96d | 442 | // Shape of the shower |
9859bfc0 | 443 | // If you change this function, change also the gradient evaluation in ChiSquare() |
173558f2 | 444 | |
225cd96d | 445 | Double_t r = sqrt(x*x+y*y); |
446 | Double_t r133 = TMath::Power(r, 1.33) ; | |
447 | Double_t r669 = TMath::Power(r, 6.69) ; | |
448 | Double_t shape = TMath::Exp( -r133 * (1. / (1.57 + 0.0860 * r133) - 0.55 / (1 + 0.000563 * r669) ) ) ; | |
9859bfc0 | 449 | return shape ; |
450 | } | |
173558f2 | 451 | |
9859bfc0 | 452 | //____________________________________________________________________________ |
225cd96d | 453 | void AliEMCALClusterizerv1::UnfoldCluster(AliEMCALRecPoint * iniTower, |
454 | Int_t nMax, | |
455 | AliEMCALDigit ** maxAt, | |
456 | Float_t * maxAtEnergy) | |
9859bfc0 | 457 | { |
458 | // Performs the unfolding of a cluster with nMax overlapping showers | |
225cd96d | 459 | Int_t nPar = 3 * nMax ; |
460 | Float_t * fitparameters = new Float_t[nPar] ; | |
51e4198e | 461 | |
225cd96d | 462 | if (fGeom==0) |
463 | AliFatal("Did not get geometry from EMCALLoader") ; | |
51e4198e | 464 | |
225cd96d | 465 | Bool_t rv = FindFit(iniTower, maxAt, maxAtEnergy, nPar, fitparameters) ; |
466 | if( !rv ) { | |
467 | // Fit failed, return and remove cluster | |
468 | iniTower->SetNExMax(-1) ; | |
469 | delete[] fitparameters ; | |
470 | return ; | |
471 | } | |
51e4198e | 472 | |
225cd96d | 473 | // create unfolded rec points and fill them with new energy lists |
474 | // First calculate energy deposited in each sell in accordance with | |
475 | // fit (without fluctuations): efit[] | |
476 | // and later correct this number in acordance with actual energy | |
477 | // deposition | |
51e4198e | 478 | |
225cd96d | 479 | Int_t nDigits = iniTower->GetMultiplicity() ; |
480 | Float_t * efit = new Float_t[nDigits] ; | |
481 | Double_t xDigit=0.,yDigit=0.,zDigit=0. ; | |
482 | Float_t xpar=0.,zpar=0.,epar=0. ; | |
51e4198e | 483 | |
225cd96d | 484 | AliEMCALDigit * digit = 0 ; |
d804d556 | 485 | Int_t * digitsList = iniTower->GetDigitsList() ; |
ee08edde | 486 | |
487 | Int_t iparam = 0 ; | |
225cd96d | 488 | Int_t iDigit ; |
489 | for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){ | |
d804d556 | 490 | digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At(digitsList[iDigit] ) ) ; |
51e4198e | 491 | if(digit){ |
492 | fGeom->RelPosCellInSModule(digit->GetId(), yDigit, xDigit, zDigit); | |
493 | efit[iDigit] = 0; | |
494 | ||
495 | while(iparam < nPar ){ | |
496 | xpar = fitparameters[iparam] ; | |
497 | zpar = fitparameters[iparam+1] ; | |
498 | epar = fitparameters[iparam+2] ; | |
499 | iparam += 3 ; | |
500 | efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ; | |
501 | } | |
225cd96d | 502 | } |
51e4198e | 503 | else AliFatal("Null digit in array!"); |
225cd96d | 504 | } |
51e4198e | 505 | |
225cd96d | 506 | // Now create new RecPoints and fill energy lists with efit corrected to fluctuations |
507 | // so that energy deposited in each cell is distributed between new clusters proportionally | |
508 | // to its contribution to efit | |
51e4198e | 509 | |
d804d556 | 510 | Float_t * energiesList = iniTower->GetEnergiesList() ; |
ee08edde | 511 | Float_t ratio = 0 ; |
51e4198e | 512 | |
225cd96d | 513 | iparam = 0 ; |
514 | while(iparam < nPar ){ | |
515 | xpar = fitparameters[iparam] ; | |
516 | zpar = fitparameters[iparam+1] ; | |
517 | epar = fitparameters[iparam+2] ; | |
518 | iparam += 3 ; | |
51e4198e | 519 | |
d804d556 | 520 | AliEMCALRecPoint * recPoint = 0 ; |
51e4198e | 521 | |
225cd96d | 522 | if(fNumberOfECAClusters >= fRecPoints->GetSize()) |
523 | fRecPoints->Expand(2*fNumberOfECAClusters) ; | |
51e4198e | 524 | |
225cd96d | 525 | (*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("") ; |
d804d556 | 526 | recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(fNumberOfECAClusters) ) ; |
51e4198e | 527 | if(recPoint){ |
528 | fNumberOfECAClusters++ ; | |
529 | recPoint->SetNExMax((Int_t)nPar/3) ; | |
530 | ||
531 | Float_t eDigit = 0. ; | |
532 | for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){ | |
533 | digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ; | |
534 | if(digit){ | |
535 | fGeom->RelPosCellInSModule(digit->GetId(), yDigit, xDigit, zDigit); | |
536 | ||
537 | ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar) / efit[iDigit] ; | |
538 | eDigit = energiesList[iDigit] * ratio ; | |
539 | recPoint->AddDigit( *digit, eDigit, kFALSE ) ; //FIXME, need to study the shared case | |
540 | } | |
541 | else AliFatal("Null digit in array!"); | |
542 | } | |
225cd96d | 543 | } |
51e4198e | 544 | else AliFatal("Null recpoint in array!"); |
225cd96d | 545 | } |
51e4198e | 546 | |
225cd96d | 547 | delete[] fitparameters ; |
548 | delete[] efit ; | |
51e4198e | 549 | |
9859bfc0 | 550 | } |
173558f2 | 551 | |
9859bfc0 | 552 | //_____________________________________________________________________________ |
225cd96d | 553 | void AliEMCALClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, |
554 | Double_t & fret, | |
555 | Double_t * x, Int_t iflag) | |
9859bfc0 | 556 | { |
557 | // Calculates the Chi square for the cluster unfolding minimization | |
558 | // Number of parameters, Gradient, Chi squared, parameters, what to do | |
51e4198e | 559 | |
225cd96d | 560 | TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ; |
51e4198e | 561 | if(!toMinuit){ |
562 | printf("Unfolding not possible!\n"); | |
563 | return; | |
564 | } | |
565 | ||
d804d556 | 566 | AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint*>( toMinuit->At(0) ) ; |
225cd96d | 567 | TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) ) ; |
568 | // A bit buggy way to get an access to the geometry | |
569 | // To be revised! | |
570 | AliEMCALGeometry *geom = dynamic_cast<AliEMCALGeometry *>(toMinuit->At(2)); | |
51e4198e | 571 | |
572 | if(!recPoint || !digits || !geom){ | |
573 | printf("Unfolding not possible!\n"); | |
574 | return; | |
575 | } | |
576 | ||
d804d556 | 577 | Int_t * digitsList = recPoint->GetDigitsList() ; |
51e4198e | 578 | |
d804d556 | 579 | Int_t nOdigits = recPoint->GetDigitsMultiplicity() ; |
51e4198e | 580 | |
d804d556 | 581 | Float_t * energiesList = recPoint->GetEnergiesList() ; |
51e4198e | 582 | |
225cd96d | 583 | fret = 0. ; |
584 | Int_t iparam ; | |
51e4198e | 585 | |
225cd96d | 586 | if(iflag == 2) |
587 | for(iparam = 0 ; iparam < nPar ; iparam++) | |
588 | Grad[iparam] = 0 ; // Will evaluate gradient | |
51e4198e | 589 | |
ee08edde | 590 | Double_t efit = 0. ; |
51e4198e | 591 | |
225cd96d | 592 | AliEMCALDigit * digit ; |
593 | Int_t iDigit ; | |
51e4198e | 594 | |
225cd96d | 595 | for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) { |
51e4198e | 596 | |
d804d556 | 597 | digit = dynamic_cast<AliEMCALDigit*>( digits->At( digitsList[iDigit] ) ); |
51e4198e | 598 | if (digit) { |
599 | ||
600 | Double_t xDigit=0 ; | |
601 | Double_t zDigit=0 ; | |
602 | Double_t yDigit=0 ;//not used yet, assumed to be 0 | |
603 | ||
604 | geom->RelPosCellInSModule(digit->GetId(), yDigit, xDigit, zDigit); | |
605 | ||
606 | if(iflag == 2){ // calculate gradient | |
607 | Int_t iParam = 0 ; | |
608 | efit = 0. ; | |
609 | while(iParam < nPar ){ | |
610 | Double_t dx = (xDigit - x[iParam]) ; | |
611 | iParam++ ; | |
612 | Double_t dz = (zDigit - x[iParam]) ; | |
613 | iParam++ ; | |
614 | efit += x[iParam] * ShowerShape(dx,dz) ; | |
615 | iParam++ ; | |
616 | } | |
617 | Double_t sum = 2. * (efit - energiesList[iDigit]) / energiesList[iDigit] ; // Here we assume, that sigma = sqrt(E) | |
618 | iParam = 0 ; | |
619 | while(iParam < nPar ){ | |
620 | Double_t xpar = x[iParam] ; | |
621 | Double_t zpar = x[iParam+1] ; | |
622 | Double_t epar = x[iParam+2] ; | |
623 | Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ); | |
624 | Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar) ; | |
625 | Double_t r133 = TMath::Power(dr, 1.33); | |
626 | Double_t r669 = TMath::Power(dr,6.69); | |
627 | Double_t deriv =-1.33 * TMath::Power(dr,0.33)*dr * ( 1.57 / ( (1.57 + 0.0860 * r133) * (1.57 + 0.0860 * r133) ) | |
628 | - 0.55 / (1 + 0.000563 * r669) / ( (1 + 0.000563 * r669) * (1 + 0.000563 * r669) ) ) ; | |
629 | ||
630 | Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ; // Derivative over x | |
631 | iParam++ ; | |
632 | Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ; // Derivative over z | |
633 | iParam++ ; | |
634 | Grad[iParam] += shape ; // Derivative over energy | |
635 | iParam++ ; | |
636 | } | |
225cd96d | 637 | } |
51e4198e | 638 | efit = 0; |
639 | iparam = 0 ; | |
640 | ||
641 | ||
642 | while(iparam < nPar ){ | |
643 | Double_t xpar = x[iparam] ; | |
644 | Double_t zpar = x[iparam+1] ; | |
645 | Double_t epar = x[iparam+2] ; | |
646 | iparam += 3 ; | |
647 | efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ; | |
225cd96d | 648 | } |
51e4198e | 649 | |
650 | fret += (efit-energiesList[iDigit])*(efit-energiesList[iDigit])/energiesList[iDigit] ; | |
651 | // Here we assume, that sigma = sqrt(E) | |
225cd96d | 652 | } |
225cd96d | 653 | } |
9859bfc0 | 654 | } |