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