]>
Commit | Line | Data |
---|---|---|
d15a28e7 | 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 | ||
b2a60966 | 16 | /* $Id$ */ |
17 | ||
702ab87e | 18 | /* History of cvs commits: |
19 | * | |
20 | * $Log$ | |
fdf65bb5 | 21 | * Revision 1.85 2005/09/27 16:08:08 hristov |
22 | * New version of CDB storage framework (A.Colla) | |
23 | * | |
9e1ceb13 | 24 | * Revision 1.84 2005/09/21 10:02:47 kharlov |
25 | * Reading calibration from CDB (Boris Polichtchouk) | |
26 | * | |
8d6d4fb5 | 27 | * Revision 1.82 2005/09/02 15:43:13 kharlov |
28 | * Add comments in GetCalibrationParameters and Calibrate | |
29 | * | |
5c149aeb | 30 | * Revision 1.81 2005/09/02 14:32:07 kharlov |
31 | * Calibration of raw data | |
32 | * | |
44ae287e | 33 | * Revision 1.80 2005/08/24 15:31:36 kharlov |
34 | * Setting raw digits flag | |
35 | * | |
9735c615 | 36 | * Revision 1.79 2005/07/25 15:53:53 kharlov |
37 | * Read raw data | |
38 | * | |
3cf4f75f | 39 | * Revision 1.78 2005/05/28 14:19:04 schutz |
40 | * Compilation warnings fixed by T.P. | |
41 | * | |
702ab87e | 42 | */ |
43 | ||
9a1398dd | 44 | //*-- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko (SUBATECH & Kurchatov Institute) |
d15a28e7 | 45 | ////////////////////////////////////////////////////////////////////////////// |
9a1398dd | 46 | // Clusterization class. Performs clusterization (collects neighbouring active cells) and |
a4e98857 | 47 | // unfolds the clusters having several local maxima. |
48 | // Results are stored in TreeR#, branches PHOSEmcRP (EMC recPoints), | |
9a1398dd | 49 | // PHOSCpvRP (CPV RecPoints) and AliPHOSClusterizer (Clusterizer with all |
f035f6ce | 50 | // parameters including input digits branch title, thresholds etc.) |
a4e98857 | 51 | // This TTask is normally called from Reconstructioner, but can as well be used in |
52 | // standalone mode. | |
53 | // Use Case: | |
fc12304f | 54 | // root [0] AliPHOSClusterizerv1 * cl = new AliPHOSClusterizerv1("galice.root", "recpointsname", "digitsname") |
a4e98857 | 55 | // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated |
fc12304f | 56 | // // reads gAlice from header file "galice.root", uses digits stored in the branch names "digitsname" (default = "Default") |
57 | // // and saves recpoints in branch named "recpointsname" (default = "digitsname") | |
a4e98857 | 58 | // root [1] cl->ExecuteTask() |
9a1398dd | 59 | // //finds RecPoints in all events stored in galice.root |
a4e98857 | 60 | // root [2] cl->SetDigitsBranch("digits2") |
f035f6ce | 61 | // //sets another title for Digitis (input) branch |
a4e98857 | 62 | // root [3] cl->SetRecPointsBranch("recp2") |
f035f6ce | 63 | // //sets another title four output branches |
a4e98857 | 64 | // root [4] cl->SetEmcLocalMaxCut(0.03) |
9a1398dd | 65 | // //set clusterization parameters |
a4e98857 | 66 | // root [5] cl->ExecuteTask("deb all time") |
9a1398dd | 67 | // //once more finds RecPoints options are |
68 | // // deb - print number of found rec points | |
69 | // // deb all - print number of found RecPoints and some their characteristics | |
70 | // // time - print benchmarking results | |
ed4205d8 | 71 | |
d15a28e7 | 72 | // --- ROOT system --- |
73 | ||
74 | #include "TMath.h" | |
9a1398dd | 75 | #include "TMinuit.h" |
76 | #include "TTree.h" | |
9a1398dd | 77 | #include "TBenchmark.h" |
d15a28e7 | 78 | |
79 | // --- Standard library --- | |
80 | ||
d15a28e7 | 81 | // --- AliRoot header files --- |
351dd634 | 82 | #include "AliLog.h" |
88cb7938 | 83 | #include "AliPHOSGetter.h" |
e957fea8 | 84 | #include "AliPHOSGeometry.h" |
d15a28e7 | 85 | #include "AliPHOSClusterizerv1.h" |
88cb7938 | 86 | #include "AliPHOSEmcRecPoint.h" |
9a1398dd | 87 | #include "AliPHOSCpvRecPoint.h" |
d15a28e7 | 88 | #include "AliPHOSDigit.h" |
9a1398dd | 89 | #include "AliPHOSDigitizer.h" |
ba898748 | 90 | #include "AliPHOSCalibrationDB.h" |
9e1ceb13 | 91 | #include "AliCDBManager.h" |
eaa5012b | 92 | #include "AliCDBStorage.h" |
fdf65bb5 | 93 | #include "AliCDBEntry.h" |
d15a28e7 | 94 | |
95 | ClassImp(AliPHOSClusterizerv1) | |
f035f6ce | 96 | |
d15a28e7 | 97 | //____________________________________________________________________________ |
7b7c1533 | 98 | AliPHOSClusterizerv1::AliPHOSClusterizerv1() : AliPHOSClusterizer() |
d15a28e7 | 99 | { |
f035f6ce | 100 | // default ctor (to be used mainly by Streamer) |
f035f6ce | 101 | |
8d0f3f77 | 102 | InitParameters() ; |
92f521a9 | 103 | fDefaultInit = kTRUE ; |
9a1398dd | 104 | } |
7b7c1533 | 105 | |
9a1398dd | 106 | //____________________________________________________________________________ |
88cb7938 | 107 | AliPHOSClusterizerv1::AliPHOSClusterizerv1(const TString alirunFileName, const TString eventFolderName) |
108 | :AliPHOSClusterizer(alirunFileName, eventFolderName) | |
9a1398dd | 109 | { |
a4e98857 | 110 | // ctor with the indication of the file where header Tree and digits Tree are stored |
9a1398dd | 111 | |
8d0f3f77 | 112 | InitParameters() ; |
2bd5457f | 113 | Init() ; |
92f521a9 | 114 | fDefaultInit = kFALSE ; |
9a1398dd | 115 | } |
fc12304f | 116 | |
9688c1dd | 117 | //____________________________________________________________________________ |
118 | AliPHOSClusterizerv1::~AliPHOSClusterizerv1() | |
119 | { | |
0bc3b8ed | 120 | // dtor |
121 | ||
9688c1dd | 122 | } |
fc12304f | 123 | //____________________________________________________________________________ |
124 | const TString AliPHOSClusterizerv1::BranchName() const | |
125 | { | |
88cb7938 | 126 | return GetName(); |
fc12304f | 127 | } |
128 | ||
9a1398dd | 129 | //____________________________________________________________________________ |
44ae287e | 130 | Float_t AliPHOSClusterizerv1::Calibrate(Int_t amp, Int_t absId) |
88cb7938 | 131 | { |
5c149aeb | 132 | // Convert digitized amplitude into energy. |
133 | // Calibration parameters are taken from calibration data base for raw data, | |
134 | // or from digitizer parameters for simulated data. | |
135 | ||
44ae287e | 136 | if(fCalibData){ |
137 | Int_t relId[4]; | |
138 | AliPHOSGetter *gime = AliPHOSGetter::Instance(); | |
139 | gime->PHOSGeometry()->AbsToRelNumbering(absId,relId) ; | |
140 | Int_t module = relId[0]; | |
8d6d4fb5 | 141 | Int_t column = relId[3]; |
eaa5012b | 142 | Int_t row = relId[2]; |
44ae287e | 143 | if(absId <= fEmcCrystals) { //calibrate as EMC |
144 | fADCchanelEmc = fCalibData->GetADCchannelEmc (module,column,row); | |
145 | fADCpedestalEmc = fCalibData->GetADCpedestalEmc(module,column,row); | |
146 | return fADCpedestalEmc + amp*fADCchanelEmc ; | |
147 | } | |
148 | else //calibrate as CPV, not implemented yet | |
149 | return 0; | |
ba898748 | 150 | } |
151 | else{ //simulation | |
152 | if(absId <= fEmcCrystals) //calibrate as EMC | |
153 | return fADCpedestalEmc + amp*fADCchanelEmc ; | |
154 | else //calibrate as CPV | |
155 | return fADCpedestalCpv+ amp*fADCchanelCpv ; | |
156 | } | |
3758d9fc | 157 | } |
548f0134 | 158 | |
3758d9fc | 159 | //____________________________________________________________________________ |
eabde521 | 160 | void AliPHOSClusterizerv1::Exec(Option_t *option) |
a4e98857 | 161 | { |
eabde521 | 162 | // Steering method to perform clusterization for events |
163 | // in the range from fFirstEvent to fLastEvent. | |
164 | // This range is optionally set by SetEventRange(). | |
165 | // if fLastEvent=-1 (by default), then process events until the end. | |
9a1398dd | 166 | |
167 | if(strstr(option,"tim")) | |
7d493c2c | 168 | gBenchmark->Start("PHOSClusterizer"); |
9a1398dd | 169 | |
eabde521 | 170 | if(strstr(option,"print")) { |
88cb7938 | 171 | Print() ; |
eabde521 | 172 | return ; |
173 | } | |
bed9e3fb | 174 | |
eaa5012b | 175 | GetCalibrationParameters() ; |
176 | ||
ba898748 | 177 | AliPHOSGetter * gime = AliPHOSGetter::Instance() ; |
9735c615 | 178 | if (fRawReader == 0) |
179 | gime->SetRawDigits(kFALSE); | |
180 | else | |
181 | gime->SetRawDigits(kTRUE); | |
88cb7938 | 182 | |
9cae45af | 183 | if (fLastEvent == -1) |
184 | fLastEvent = gime->MaxEvent() - 1 ; | |
185 | else | |
fb43ada4 | 186 | fLastEvent = TMath::Min(fFirstEvent, gime->MaxEvent()); // one event at the time |
eabde521 | 187 | Int_t nEvents = fLastEvent - fFirstEvent + 1; |
ba898748 | 188 | |
9cae45af | 189 | Int_t ievent ; |
9cae45af | 190 | for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) { |
3cf4f75f | 191 | if (fRawReader == 0) |
192 | gime->Event(ievent ,"D"); // Read digits from simulated data | |
193 | else | |
194 | gime->Event(fRawReader,"W"); // Read digits from raw data | |
04f0bda3 | 195 | |
548f0134 | 196 | fNumberOfEmcClusters = fNumberOfCpvClusters = 0 ; |
9a1398dd | 197 | |
88cb7938 | 198 | MakeClusters() ; |
199 | ||
fc12304f | 200 | if(fToUnfold) |
201 | MakeUnfolding() ; | |
7b7c1533 | 202 | |
88cb7938 | 203 | WriteRecPoints(); |
7b7c1533 | 204 | |
94de8339 | 205 | if(strstr(option,"deb")) |
206 | PrintRecPoints(option) ; | |
207 | ||
88cb7938 | 208 | //increment the total number of recpoints per run |
fbf811ec | 209 | fRecPointsInRun += gime->EmcRecPoints()->GetEntriesFast() ; |
210 | fRecPointsInRun += gime->CpvRecPoints()->GetEntriesFast() ; | |
88cb7938 | 211 | } |
212 | ||
ba898748 | 213 | if(fWrite) //do not unload in "on flight" mode |
214 | Unload(); | |
9a1398dd | 215 | |
216 | if(strstr(option,"tim")){ | |
217 | gBenchmark->Stop("PHOSClusterizer"); | |
351dd634 | 218 | AliInfo(Form("took %f seconds for Clusterizing %f seconds per event \n", |
21cd0c07 | 219 | gBenchmark->GetCpuTime("PHOSClusterizer"), |
351dd634 | 220 | gBenchmark->GetCpuTime("PHOSClusterizer")/nEvents )) ; |
88cb7938 | 221 | } |
9a1398dd | 222 | } |
223 | ||
224 | //____________________________________________________________________________ | |
a0636361 | 225 | Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, AliPHOSDigit ** maxAt, Float_t * maxAtEnergy, |
d1de15f5 | 226 | Int_t nPar, Float_t * fitparameters) const |
9a1398dd | 227 | { |
228 | // Calls TMinuit to fit the energy distribution of a cluster with several maxima | |
a4e98857 | 229 | // The initial values for fitting procedure are set equal to the positions of local maxima. |
f035f6ce | 230 | // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers |
9a1398dd | 231 | |
88cb7938 | 232 | |
233 | AliPHOSGetter * gime = AliPHOSGetter::Instance(); | |
234 | TClonesArray * digits = gime->Digits(); | |
7b7c1533 | 235 | |
236 | ||
9a1398dd | 237 | gMinuit->mncler(); // Reset Minuit's list of paramters |
238 | gMinuit->SetPrintLevel(-1) ; // No Printout | |
239 | gMinuit->SetFCN(AliPHOSClusterizerv1::UnfoldingChiSquare) ; | |
240 | // To set the address of the minimization function | |
241 | ||
242 | TList * toMinuit = new TList(); | |
243 | toMinuit->AddAt(emcRP,0) ; | |
7b7c1533 | 244 | toMinuit->AddAt(digits,1) ; |
9a1398dd | 245 | |
246 | gMinuit->SetObjectFit(toMinuit) ; // To tranfer pointer to UnfoldingChiSquare | |
247 | ||
248 | // filling initial values for fit parameters | |
249 | AliPHOSDigit * digit ; | |
250 | ||
251 | Int_t ierflg = 0; | |
252 | Int_t index = 0 ; | |
253 | Int_t nDigits = (Int_t) nPar / 3 ; | |
254 | ||
255 | Int_t iDigit ; | |
256 | ||
7b7c1533 | 257 | const AliPHOSGeometry * geom = gime->PHOSGeometry() ; |
9a1398dd | 258 | |
259 | for(iDigit = 0; iDigit < nDigits; iDigit++){ | |
a0636361 | 260 | digit = maxAt[iDigit]; |
9a1398dd | 261 | |
262 | Int_t relid[4] ; | |
7b7c1533 | 263 | Float_t x = 0.; |
264 | Float_t z = 0.; | |
265 | geom->AbsToRelNumbering(digit->GetId(), relid) ; | |
266 | geom->RelPosInModule(relid, x, z) ; | |
9a1398dd | 267 | |
268 | Float_t energy = maxAtEnergy[iDigit] ; | |
269 | ||
270 | gMinuit->mnparm(index, "x", x, 0.1, 0, 0, ierflg) ; | |
271 | index++ ; | |
272 | if(ierflg != 0){ | |
21cd0c07 | 273 | Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : x = %f\n", x ) ; |
9a1398dd | 274 | return kFALSE; |
275 | } | |
276 | gMinuit->mnparm(index, "z", z, 0.1, 0, 0, ierflg) ; | |
277 | index++ ; | |
278 | if(ierflg != 0){ | |
21cd0c07 | 279 | Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : z =%f\n", z ) ; |
9a1398dd | 280 | return kFALSE; |
281 | } | |
282 | gMinuit->mnparm(index, "Energy", energy , 0.05*energy, 0., 4.*energy, ierflg) ; | |
283 | index++ ; | |
284 | if(ierflg != 0){ | |
21cd0c07 | 285 | Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : energy = %f\n", energy ) ; |
9a1398dd | 286 | return kFALSE; |
287 | } | |
288 | } | |
289 | ||
290 | Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly | |
291 | // depends on it. | |
292 | Double_t p1 = 1.0 ; | |
293 | Double_t p2 = 0.0 ; | |
294 | ||
295 | gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ; // force TMinuit to reduce function calls | |
296 | gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient | |
297 | gMinuit->SetMaxIterations(5); | |
298 | gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ; // No Warnings | |
299 | ||
300 | gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize | |
301 | ||
302 | if(ierflg == 4){ // Minimum not found | |
21cd0c07 | 303 | Warning("FindFit", "PHOS Unfolding fit not converged, cluster abandoned\n" ); |
9a1398dd | 304 | return kFALSE ; |
305 | } | |
306 | for(index = 0; index < nPar; index++){ | |
307 | Double_t err ; | |
308 | Double_t val ; | |
309 | gMinuit->GetParameter(index, val, err) ; // Returns value and error of parameter index | |
310 | fitparameters[index] = val ; | |
311 | } | |
312 | ||
313 | delete toMinuit ; | |
314 | return kTRUE; | |
c3c187e5 | 315 | |
d15a28e7 | 316 | } |
317 | ||
3758d9fc | 318 | //____________________________________________________________________________ |
319 | void AliPHOSClusterizerv1::GetCalibrationParameters() | |
320 | { | |
5c149aeb | 321 | // Set calibration parameters: |
eaa5012b | 322 | // if calibration database exists, they are read from database, |
323 | // otherwise, they are taken from digitizer. | |
5c149aeb | 324 | // |
eaa5012b | 325 | // It is a user responsilibity to open CDB before reconstruction: |
326 | // AliCDBLocal *loc = new AliCDBLocal("CalibDB"); | |
fbf811ec | 327 | |
88cb7938 | 328 | AliPHOSGetter * gime = AliPHOSGetter::Instance(); |
eaa5012b | 329 | |
9e1ceb13 | 330 | if(AliCDBManager::Instance()->IsDefaultStorageSet()){ |
331 | AliCDBEntry *entry = (AliCDBEntry*) AliCDBManager::Instance()->GetDefaultStorage() | |
eaa5012b | 332 | ->Get("PHOS/Calib/GainFactors_and_Pedestals",gAlice->GetRunNumber()); |
9e1ceb13 | 333 | fCalibData = (AliPHOSCalibData*) entry->GetObject(); |
334 | } | |
335 | ||
336 | ||
eaa5012b | 337 | if(!fCalibData) |
338 | { | |
339 | if ( !gime->Digitizer() ) | |
340 | gime->LoadDigitizer(); | |
341 | AliPHOSDigitizer * dig = gime->Digitizer(); | |
342 | fADCchanelEmc = dig->GetEMCchannel() ; | |
343 | fADCpedestalEmc = dig->GetEMCpedestal(); | |
ba898748 | 344 | |
eaa5012b | 345 | fADCchanelCpv = dig->GetCPVchannel() ; |
346 | fADCpedestalCpv = dig->GetCPVpedestal() ; | |
347 | } | |
3758d9fc | 348 | } |
548f0134 | 349 | |
d15a28e7 | 350 | //____________________________________________________________________________ |
a4e98857 | 351 | void AliPHOSClusterizerv1::Init() |
352 | { | |
2bd5457f | 353 | // Make all memory allocations which can not be done in default constructor. |
354 | // Attach the Clusterizer task to the list of PHOS tasks | |
ba898748 | 355 | |
356 | AliPHOSGetter* gime = AliPHOSGetter::Instance() ; | |
357 | if(!gime) | |
358 | gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName.Data()); | |
fbf811ec | 359 | |
88cb7938 | 360 | AliPHOSGeometry * geom = gime->PHOSGeometry(); |
ba898748 | 361 | |
3758d9fc | 362 | fEmcCrystals = geom->GetNModules() * geom->GetNCristalsInModule() ; |
363 | ||
7b7c1533 | 364 | if(!gMinuit) |
88cb7938 | 365 | gMinuit = new TMinuit(100); |
fbf811ec | 366 | |
ba898748 | 367 | if ( !gime->Clusterizer() ) { |
88cb7938 | 368 | gime->PostClusterizer(this); |
ba898748 | 369 | } |
9a1398dd | 370 | } |
7b7c1533 | 371 | |
8d0f3f77 | 372 | //____________________________________________________________________________ |
373 | void AliPHOSClusterizerv1::InitParameters() | |
374 | { | |
375 | ||
376 | fNumberOfCpvClusters = 0 ; | |
377 | fNumberOfEmcClusters = 0 ; | |
378 | ||
379 | fCpvClusteringThreshold = 0.0; | |
380 | fEmcClusteringThreshold = 0.2; | |
381 | ||
382 | fEmcLocMaxCut = 0.03 ; | |
383 | fCpvLocMaxCut = 0.03 ; | |
ba898748 | 384 | |
385 | fEmcMinE = 0.01 ; | |
386 | fCpvMinE = 0.0 ; | |
387 | ||
8d0f3f77 | 388 | fW0 = 4.5 ; |
389 | fW0CPV = 4.0 ; | |
390 | ||
391 | fEmcTimeGate = 1.e-8 ; | |
392 | ||
393 | fToUnfold = kTRUE ; | |
88cb7938 | 394 | |
8d0f3f77 | 395 | fRecPointsInRun = 0 ; |
c2cd5471 | 396 | |
ba898748 | 397 | fWrite = kTRUE ; |
398 | ||
44ae287e | 399 | fCalibData = 0 ; |
ae911977 | 400 | |
eabde521 | 401 | SetEventRange(0,-1) ; |
8d0f3f77 | 402 | } |
403 | ||
9a1398dd | 404 | //____________________________________________________________________________ |
405 | Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)const | |
d15a28e7 | 406 | { |
b2a60966 | 407 | // Gives the neighbourness of two digits = 0 are not neighbour but continue searching |
408 | // = 1 are neighbour | |
409 | // = 2 are not neighbour but do not continue searching | |
a4e98857 | 410 | // neighbours are defined as digits having at least a common vertex |
b2a60966 | 411 | // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster |
412 | // which is compared to a digit (d2) not yet in a cluster | |
413 | ||
88cb7938 | 414 | AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ; |
7b7c1533 | 415 | |
d15a28e7 | 416 | Int_t rv = 0 ; |
417 | ||
d15a28e7 | 418 | Int_t relid1[4] ; |
7b7c1533 | 419 | geom->AbsToRelNumbering(d1->GetId(), relid1) ; |
d15a28e7 | 420 | |
421 | Int_t relid2[4] ; | |
7b7c1533 | 422 | geom->AbsToRelNumbering(d2->GetId(), relid2) ; |
d15a28e7 | 423 | |
9688c1dd | 424 | if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module |
92862013 | 425 | Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ; |
426 | Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ; | |
d15a28e7 | 427 | |
92862013 | 428 | if (( coldiff <= 1 ) && ( rowdiff <= 1 )){ |
9688c1dd | 429 | if((relid1[1] != 0) || (TMath::Abs(d1->GetTime() - d2->GetTime() ) < fEmcTimeGate)) |
d15a28e7 | 430 | rv = 1 ; |
431 | } | |
432 | else { | |
433 | if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1)) | |
88cb7938 | 434 | rv = 2; // Difference in row numbers is too large to look further |
d15a28e7 | 435 | } |
436 | ||
437 | } | |
438 | else { | |
439 | ||
9688c1dd | 440 | if( (relid1[0] < relid2[0]) || (relid1[1] != relid2[1]) ) |
d15a28e7 | 441 | rv=2 ; |
442 | ||
443 | } | |
d72dfbc3 | 444 | |
d15a28e7 | 445 | return rv ; |
446 | } | |
ba898748 | 447 | //____________________________________________________________________________ |
448 | void AliPHOSClusterizerv1::CleanDigits(TClonesArray * digits){ | |
449 | for(Int_t i=0; i<digits->GetEntriesFast(); i++){ | |
450 | AliPHOSDigit * digit = static_cast<AliPHOSDigit*>(digits->At(i)) ; | |
451 | Float_t cut = IsInEmc(digit) ? fEmcMinE : fCpvMinE ; | |
452 | if(Calibrate(digit->GetAmp(),digit->GetId()) < cut) | |
453 | digits->RemoveAt(i) ; | |
454 | } | |
455 | digits->Compress() ; | |
456 | for (Int_t i = 0 ; i < digits->GetEntriesFast() ; i++) { | |
457 | AliPHOSDigit *digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ; | |
458 | digit->SetIndexInList(i) ; | |
459 | } | |
d15a28e7 | 460 | |
ba898748 | 461 | } |
d15a28e7 | 462 | //____________________________________________________________________________ |
9a1398dd | 463 | Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit) const |
d15a28e7 | 464 | { |
b2a60966 | 465 | // Tells if (true) or not (false) the digit is in a PHOS-EMC module |
466 | ||
9f616d61 | 467 | Bool_t rv = kFALSE ; |
88cb7938 | 468 | AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ; |
d15a28e7 | 469 | |
2b629790 | 470 | Int_t nEMC = geom->GetNModules()*geom->GetNPhi()*geom->GetNZ(); |
d15a28e7 | 471 | |
2b629790 | 472 | if(digit->GetId() <= nEMC ) rv = kTRUE; |
ed4205d8 | 473 | |
474 | return rv ; | |
475 | } | |
476 | ||
ed4205d8 | 477 | //____________________________________________________________________________ |
9a1398dd | 478 | Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit) const |
ed4205d8 | 479 | { |
fad3e5b9 | 480 | // Tells if (true) or not (false) the digit is in a PHOS-CPV module |
ed4205d8 | 481 | |
482 | Bool_t rv = kFALSE ; | |
88cb7938 | 483 | |
484 | AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ; | |
ed4205d8 | 485 | |
2b629790 | 486 | Int_t nEMC = geom->GetNModules()*geom->GetNPhi()*geom->GetNZ(); |
ed4205d8 | 487 | |
2b629790 | 488 | if(digit->GetId() > nEMC ) rv = kTRUE; |
d15a28e7 | 489 | |
490 | return rv ; | |
491 | } | |
9a1398dd | 492 | |
493 | //____________________________________________________________________________ | |
88cb7938 | 494 | void AliPHOSClusterizerv1::Unload() |
495 | { | |
496 | AliPHOSGetter * gime = AliPHOSGetter::Instance() ; | |
497 | gime->PhosLoader()->UnloadDigits() ; | |
498 | gime->PhosLoader()->UnloadRecPoints() ; | |
499 | } | |
500 | ||
501 | //____________________________________________________________________________ | |
502 | void AliPHOSClusterizerv1::WriteRecPoints() | |
a4e98857 | 503 | { |
7b7c1533 | 504 | |
505 | // Creates new branches with given title | |
a4e98857 | 506 | // fills and writes into TreeR. |
88cb7938 | 507 | |
508 | AliPHOSGetter * gime = AliPHOSGetter::Instance(); | |
509 | ||
fbf811ec | 510 | TObjArray * emcRecPoints = gime->EmcRecPoints() ; |
511 | TObjArray * cpvRecPoints = gime->CpvRecPoints() ; | |
512 | TClonesArray * digits = gime->Digits() ; | |
ba898748 | 513 | |
514 | ||
9a1398dd | 515 | Int_t index ; |
ba898748 | 516 | //Evaluate position, dispersion and other RecPoint properties.. |
517 | Int_t nEmc = emcRecPoints->GetEntriesFast(); | |
518 | for(index = 0; index < nEmc; index++){ | |
519 | AliPHOSEmcRecPoint * rp = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) ); | |
520 | rp->Purify(fEmcMinE) ; | |
521 | if(rp->GetMultiplicity()>0.) //If this RP is not empty | |
522 | rp->EvalAll(fW0,digits) ; | |
523 | else{ | |
524 | emcRecPoints->RemoveAt(index) ; | |
525 | delete rp ; | |
526 | } | |
527 | } | |
528 | emcRecPoints->Compress() ; | |
7b7c1533 | 529 | emcRecPoints->Sort() ; |
ba898748 | 530 | // emcRecPoints->Expand(emcRecPoints->GetEntriesFast()) ; |
531 | for(index = 0; index < emcRecPoints->GetEntries(); index++){ | |
88cb7938 | 532 | dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) )->SetIndexInList(index) ; |
ba898748 | 533 | } |
fbf811ec | 534 | |
9a1398dd | 535 | //Now the same for CPV |
ba898748 | 536 | for(index = 0; index < cpvRecPoints->GetEntries(); index++){ |
537 | AliPHOSCpvRecPoint * rp = dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(index) ); | |
538 | rp->EvalAll(fW0CPV,digits) ; | |
539 | } | |
7b7c1533 | 540 | cpvRecPoints->Sort() ; |
fbf811ec | 541 | |
7b7c1533 | 542 | for(index = 0; index < cpvRecPoints->GetEntries(); index++) |
2bb500e5 | 543 | dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(index) )->SetIndexInList(index) ; |
fbf811ec | 544 | |
7b7c1533 | 545 | cpvRecPoints->Expand(cpvRecPoints->GetEntriesFast()) ; |
9a1398dd | 546 | |
ba898748 | 547 | if(fWrite){ //We write TreeR |
548 | TTree * treeR = gime->TreeR(); | |
9a1398dd | 549 | |
ba898748 | 550 | Int_t bufferSize = 32000 ; |
551 | Int_t splitlevel = 0 ; | |
9a1398dd | 552 | |
ba898748 | 553 | //First EMC |
554 | TBranch * emcBranch = treeR->Branch("PHOSEmcRP","TObjArray",&emcRecPoints,bufferSize,splitlevel); | |
555 | emcBranch->SetTitle(BranchName()); | |
556 | ||
557 | //Now CPV branch | |
558 | TBranch * cpvBranch = treeR->Branch("PHOSCpvRP","TObjArray",&cpvRecPoints,bufferSize,splitlevel); | |
559 | cpvBranch->SetTitle(BranchName()); | |
560 | ||
561 | emcBranch ->Fill() ; | |
562 | cpvBranch ->Fill() ; | |
563 | ||
564 | gime->WriteRecPoints("OVERWRITE"); | |
565 | gime->WriteClusterizer("OVERWRITE"); | |
566 | } | |
9a1398dd | 567 | } |
568 | ||
569 | //____________________________________________________________________________ | |
570 | void AliPHOSClusterizerv1::MakeClusters() | |
571 | { | |
572 | // Steering method to construct the clusters stored in a list of Reconstructed Points | |
573 | // A cluster is defined as a list of neighbour digits | |
d15a28e7 | 574 | |
88cb7938 | 575 | |
576 | AliPHOSGetter * gime = AliPHOSGetter::Instance(); | |
577 | ||
578 | TObjArray * emcRecPoints = gime->EmcRecPoints() ; | |
579 | TObjArray * cpvRecPoints = gime->CpvRecPoints() ; | |
7b7c1533 | 580 | |
7d493c2c | 581 | emcRecPoints->Delete() ; |
582 | cpvRecPoints->Delete() ; | |
7b7c1533 | 583 | |
fbf811ec | 584 | TClonesArray * digits = gime->Digits() ; |
88cb7938 | 585 | |
ba898748 | 586 | //Remove digits below threshold |
587 | CleanDigits(digits) ; | |
588 | ||
589 | ||
88cb7938 | 590 | TClonesArray * digitsC = static_cast<TClonesArray*>( digits->Clone() ) ; |
7b7c1533 | 591 | |
592 | ||
d15a28e7 | 593 | // Clusterization starts |
9a1398dd | 594 | |
7b7c1533 | 595 | TIter nextdigit(digitsC) ; |
d15a28e7 | 596 | AliPHOSDigit * digit ; |
92862013 | 597 | Bool_t notremoved = kTRUE ; |
6ad0bfa0 | 598 | |
88cb7938 | 599 | while ( (digit = dynamic_cast<AliPHOSDigit *>( nextdigit()) ) ) { // scan over the list of digitsC |
600 | ||
601 | ||
9a1398dd | 602 | AliPHOSRecPoint * clu = 0 ; |
c161df70 | 603 | |
d1de15f5 | 604 | TArrayI clusterdigitslist(1500) ; |
d15a28e7 | 605 | Int_t index ; |
f2bc1b87 | 606 | |
3758d9fc | 607 | if (( IsInEmc (digit) && Calibrate(digit->GetAmp(),digit->GetId()) > fEmcClusteringThreshold ) || |
608 | ( IsInCpv (digit) && Calibrate(digit->GetAmp(),digit->GetId()) > fCpvClusteringThreshold ) ) { | |
d15a28e7 | 609 | Int_t iDigitInCluster = 0 ; |
7b7c1533 | 610 | |
6ad0bfa0 | 611 | if ( IsInEmc(digit) ) { |
88cb7938 | 612 | // start a new EMC RecPoint |
613 | if(fNumberOfEmcClusters >= emcRecPoints->GetSize()) | |
614 | emcRecPoints->Expand(2*fNumberOfEmcClusters+1) ; | |
615 | ||
616 | emcRecPoints->AddAt(new AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ; | |
617 | clu = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(fNumberOfEmcClusters) ) ; | |
618 | fNumberOfEmcClusters++ ; | |
619 | clu->AddDigit(*digit, Calibrate(digit->GetAmp(),digit->GetId())) ; | |
620 | clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ; | |
621 | iDigitInCluster++ ; | |
622 | digitsC->Remove(digit) ; | |
ca77b032 | 623 | |
d72dfbc3 | 624 | } else { |
88cb7938 | 625 | |
626 | // start a new CPV cluster | |
627 | if(fNumberOfCpvClusters >= cpvRecPoints->GetSize()) | |
628 | cpvRecPoints->Expand(2*fNumberOfCpvClusters+1); | |
629 | ||
630 | cpvRecPoints->AddAt(new AliPHOSCpvRecPoint(""), fNumberOfCpvClusters) ; | |
631 | ||
632 | clu = dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(fNumberOfCpvClusters) ) ; | |
633 | fNumberOfCpvClusters++ ; | |
634 | clu->AddDigit(*digit, Calibrate(digit->GetAmp(),digit->GetId()) ) ; | |
635 | clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ; | |
636 | iDigitInCluster++ ; | |
637 | digitsC->Remove(digit) ; | |
d15a28e7 | 638 | nextdigit.Reset() ; |
88cb7938 | 639 | |
640 | // Here we remove remaining EMC digits, which cannot make a cluster | |
641 | ||
92862013 | 642 | if( notremoved ) { |
88cb7938 | 643 | while( ( digit = dynamic_cast<AliPHOSDigit *>( nextdigit() ) ) ) { |
9f616d61 | 644 | if( IsInEmc(digit) ) |
88cb7938 | 645 | digitsC->Remove(digit) ; |
d15a28e7 | 646 | else |
88cb7938 | 647 | break ; |
648 | } | |
649 | notremoved = kFALSE ; | |
650 | } | |
651 | ||
d15a28e7 | 652 | } // else |
653 | ||
654 | nextdigit.Reset() ; | |
fad3e5b9 | 655 | |
d15a28e7 | 656 | AliPHOSDigit * digitN ; |
657 | index = 0 ; | |
9f616d61 | 658 | while (index < iDigitInCluster){ // scan over digits already in cluster |
88cb7938 | 659 | digit = dynamic_cast<AliPHOSDigit*>( digits->At(clusterdigitslist[index]) ) ; |
660 | index++ ; | |
661 | while ( (digitN = dynamic_cast<AliPHOSDigit *>( nextdigit() ) ) ) { // scan over the reduced list of digits | |
662 | Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!! | |
d15a28e7 | 663 | switch (ineb ) { |
9f616d61 | 664 | case 0 : // not a neighbour |
88cb7938 | 665 | break ; |
666 | case 1 : // are neighbours | |
667 | clu->AddDigit(*digitN, Calibrate( digitN->GetAmp(), digitN->GetId() ) ) ; | |
668 | clusterdigitslist[iDigitInCluster] = digitN->GetIndexInList() ; | |
669 | iDigitInCluster++ ; | |
670 | digitsC->Remove(digitN) ; | |
671 | break ; | |
9f616d61 | 672 | case 2 : // too far from each other |
88cb7938 | 673 | goto endofloop; |
674 | } // switch | |
675 | ||
676 | } // while digitN | |
677 | ||
d15a28e7 | 678 | endofloop: ; |
88cb7938 | 679 | nextdigit.Reset() ; |
680 | ||
d15a28e7 | 681 | } // loop over cluster |
ad8cfaf4 | 682 | |
ad8cfaf4 | 683 | } // energy theshold |
ad8cfaf4 | 684 | |
31aa6d6c | 685 | |
d15a28e7 | 686 | } // while digit |
687 | ||
7b7c1533 | 688 | delete digitsC ; |
9688c1dd | 689 | |
9a1398dd | 690 | } |
691 | ||
692 | //____________________________________________________________________________ | |
a4e98857 | 693 | void AliPHOSClusterizerv1::MakeUnfolding() |
694 | { | |
695 | // Unfolds clusters using the shape of an ElectroMagnetic shower | |
9688c1dd | 696 | // Performs unfolding of all EMC/CPV clusters |
9a1398dd | 697 | |
88cb7938 | 698 | AliPHOSGetter * gime = AliPHOSGetter::Instance(); |
fc12304f | 699 | |
7b7c1533 | 700 | const AliPHOSGeometry * geom = gime->PHOSGeometry() ; |
88cb7938 | 701 | |
702 | TObjArray * emcRecPoints = gime->EmcRecPoints() ; | |
703 | TObjArray * cpvRecPoints = gime->CpvRecPoints() ; | |
fbf811ec | 704 | TClonesArray * digits = gime->Digits() ; |
7b7c1533 | 705 | |
a4e98857 | 706 | // Unfold first EMC clusters |
9a1398dd | 707 | if(fNumberOfEmcClusters > 0){ |
708 | ||
7b7c1533 | 709 | Int_t nModulesToUnfold = geom->GetNModules() ; |
9a1398dd | 710 | |
711 | Int_t numberofNotUnfolded = fNumberOfEmcClusters ; | |
712 | Int_t index ; | |
713 | for(index = 0 ; index < numberofNotUnfolded ; index++){ | |
714 | ||
88cb7938 | 715 | AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) ) ; |
9a1398dd | 716 | if(emcRecPoint->GetPHOSMod()> nModulesToUnfold) |
88cb7938 | 717 | break ; |
9a1398dd | 718 | |
719 | Int_t nMultipl = emcRecPoint->GetMultiplicity() ; | |
a0636361 | 720 | AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ; |
9a1398dd | 721 | Float_t * maxAtEnergy = new Float_t[nMultipl] ; |
7b7c1533 | 722 | Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fEmcLocMaxCut,digits) ; |
9a1398dd | 723 | |
724 | if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0 | |
88cb7938 | 725 | UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ; |
726 | emcRecPoints->Remove(emcRecPoint); | |
727 | emcRecPoints->Compress() ; | |
728 | index-- ; | |
729 | fNumberOfEmcClusters -- ; | |
730 | numberofNotUnfolded-- ; | |
9a1398dd | 731 | } |
c6bf27f2 | 732 | else{ |
733 | emcRecPoint->SetNExMax(1) ; //Only one local maximum | |
734 | } | |
9a1398dd | 735 | |
736 | delete[] maxAt ; | |
737 | delete[] maxAtEnergy ; | |
738 | } | |
739 | } | |
a4e98857 | 740 | // Unfolding of EMC clusters finished |
9a1398dd | 741 | |
742 | ||
a4e98857 | 743 | // Unfold now CPV clusters |
9a1398dd | 744 | if(fNumberOfCpvClusters > 0){ |
745 | ||
9688c1dd | 746 | Int_t nModulesToUnfold = geom->GetNModules() ; |
9a1398dd | 747 | |
748 | Int_t numberofCpvNotUnfolded = fNumberOfCpvClusters ; | |
749 | Int_t index ; | |
750 | for(index = 0 ; index < numberofCpvNotUnfolded ; index++){ | |
751 | ||
88cb7938 | 752 | AliPHOSRecPoint * recPoint = dynamic_cast<AliPHOSRecPoint *>( cpvRecPoints->At(index) ) ; |
9a1398dd | 753 | |
754 | if(recPoint->GetPHOSMod()> nModulesToUnfold) | |
88cb7938 | 755 | break ; |
9a1398dd | 756 | |
88cb7938 | 757 | AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint*>(recPoint) ; |
9a1398dd | 758 | |
759 | Int_t nMultipl = emcRecPoint->GetMultiplicity() ; | |
a0636361 | 760 | AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ; |
9a1398dd | 761 | Float_t * maxAtEnergy = new Float_t[nMultipl] ; |
7b7c1533 | 762 | Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fCpvLocMaxCut,digits) ; |
9a1398dd | 763 | |
764 | if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0 | |
88cb7938 | 765 | UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ; |
766 | cpvRecPoints->Remove(emcRecPoint); | |
767 | cpvRecPoints->Compress() ; | |
768 | index-- ; | |
769 | numberofCpvNotUnfolded-- ; | |
770 | fNumberOfCpvClusters-- ; | |
9a1398dd | 771 | } |
772 | ||
773 | delete[] maxAt ; | |
774 | delete[] maxAtEnergy ; | |
775 | } | |
776 | } | |
777 | //Unfolding of Cpv clusters finished | |
778 | ||
779 | } | |
780 | ||
9a1398dd | 781 | //____________________________________________________________________________ |
782 | Double_t AliPHOSClusterizerv1::ShowerShape(Double_t r) | |
783 | { | |
784 | // Shape of the shower (see PHOS TDR) | |
a4e98857 | 785 | // If you change this function, change also the gradient evaluation in ChiSquare() |
9a1398dd | 786 | |
787 | Double_t r4 = r*r*r*r ; | |
788 | Double_t r295 = TMath::Power(r, 2.95) ; | |
789 | Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ; | |
790 | return shape ; | |
791 | } | |
792 | ||
793 | //____________________________________________________________________________ | |
794 | void AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc, | |
88cb7938 | 795 | Int_t nMax, |
796 | AliPHOSDigit ** maxAt, | |
797 | Float_t * maxAtEnergy) | |
9a1398dd | 798 | { |
799 | // Performs the unfolding of a cluster with nMax overlapping showers | |
800 | ||
88cb7938 | 801 | AliPHOSGetter * gime = AliPHOSGetter::Instance(); |
fc12304f | 802 | |
7b7c1533 | 803 | const AliPHOSGeometry * geom = gime->PHOSGeometry() ; |
88cb7938 | 804 | |
fbf811ec | 805 | const TClonesArray * digits = gime->Digits() ; |
806 | TObjArray * emcRecPoints = gime->EmcRecPoints() ; | |
807 | TObjArray * cpvRecPoints = gime->CpvRecPoints() ; | |
7b7c1533 | 808 | |
9a1398dd | 809 | Int_t nPar = 3 * nMax ; |
810 | Float_t * fitparameters = new Float_t[nPar] ; | |
811 | ||
812 | Bool_t rv = FindFit(iniEmc, maxAt, maxAtEnergy, nPar, fitparameters) ; | |
813 | if( !rv ) { | |
88cb7938 | 814 | // Fit failed, return and remove cluster |
c6bf27f2 | 815 | iniEmc->SetNExMax(-1) ; |
9a1398dd | 816 | delete[] fitparameters ; |
817 | return ; | |
818 | } | |
819 | ||
820 | // create ufolded rec points and fill them with new energy lists | |
821 | // First calculate energy deposited in each sell in accordance with fit (without fluctuations): efit[] | |
822 | // and later correct this number in acordance with actual energy deposition | |
823 | ||
824 | Int_t nDigits = iniEmc->GetMultiplicity() ; | |
825 | Float_t * efit = new Float_t[nDigits] ; | |
7b7c1533 | 826 | Float_t xDigit=0.,zDigit=0.,distance=0. ; |
827 | Float_t xpar=0.,zpar=0.,epar=0. ; | |
9a1398dd | 828 | Int_t relid[4] ; |
7b7c1533 | 829 | AliPHOSDigit * digit = 0 ; |
9a1398dd | 830 | Int_t * emcDigits = iniEmc->GetDigitsList() ; |
831 | ||
832 | Int_t iparam ; | |
833 | Int_t iDigit ; | |
834 | for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){ | |
88cb7938 | 835 | digit = dynamic_cast<AliPHOSDigit*>( digits->At(emcDigits[iDigit] ) ) ; |
7b7c1533 | 836 | geom->AbsToRelNumbering(digit->GetId(), relid) ; |
837 | geom->RelPosInModule(relid, xDigit, zDigit) ; | |
9a1398dd | 838 | efit[iDigit] = 0; |
839 | ||
840 | iparam = 0 ; | |
841 | while(iparam < nPar ){ | |
842 | xpar = fitparameters[iparam] ; | |
843 | zpar = fitparameters[iparam+1] ; | |
844 | epar = fitparameters[iparam+2] ; | |
845 | iparam += 3 ; | |
846 | distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ; | |
847 | distance = TMath::Sqrt(distance) ; | |
848 | efit[iDigit] += epar * ShowerShape(distance) ; | |
849 | } | |
850 | } | |
851 | ||
852 | ||
853 | // Now create new RecPoints and fill energy lists with efit corrected to fluctuations | |
854 | // so that energy deposited in each cell is distributed betwin new clusters proportionally | |
855 | // to its contribution to efit | |
856 | ||
857 | Float_t * emcEnergies = iniEmc->GetEnergiesList() ; | |
858 | Float_t ratio ; | |
859 | ||
860 | iparam = 0 ; | |
861 | while(iparam < nPar ){ | |
862 | xpar = fitparameters[iparam] ; | |
863 | zpar = fitparameters[iparam+1] ; | |
864 | epar = fitparameters[iparam+2] ; | |
865 | iparam += 3 ; | |
866 | ||
7b7c1533 | 867 | AliPHOSEmcRecPoint * emcRP = 0 ; |
9a1398dd | 868 | |
869 | if(iniEmc->IsEmc()){ //create new entries in fEmcRecPoints... | |
870 | ||
7b7c1533 | 871 | if(fNumberOfEmcClusters >= emcRecPoints->GetSize()) |
88cb7938 | 872 | emcRecPoints->Expand(2*fNumberOfEmcClusters) ; |
9a1398dd | 873 | |
73a68ccb | 874 | (*emcRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint("") ; |
88cb7938 | 875 | emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(fNumberOfEmcClusters) ) ; |
9a1398dd | 876 | fNumberOfEmcClusters++ ; |
c6bf27f2 | 877 | emcRP->SetNExMax((Int_t)nPar/3) ; |
9a1398dd | 878 | } |
879 | else{//create new entries in fCpvRecPoints | |
7b7c1533 | 880 | if(fNumberOfCpvClusters >= cpvRecPoints->GetSize()) |
88cb7938 | 881 | cpvRecPoints->Expand(2*fNumberOfCpvClusters) ; |
9a1398dd | 882 | |
73a68ccb | 883 | (*cpvRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint("") ; |
88cb7938 | 884 | emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( cpvRecPoints->At(fNumberOfCpvClusters) ) ; |
9a1398dd | 885 | fNumberOfCpvClusters++ ; |
886 | } | |
887 | ||
888 | Float_t eDigit ; | |
889 | for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){ | |
88cb7938 | 890 | digit = dynamic_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) ) ; |
7b7c1533 | 891 | geom->AbsToRelNumbering(digit->GetId(), relid) ; |
892 | geom->RelPosInModule(relid, xDigit, zDigit) ; | |
9a1398dd | 893 | distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ; |
894 | distance = TMath::Sqrt(distance) ; | |
895 | ratio = epar * ShowerShape(distance) / efit[iDigit] ; | |
896 | eDigit = emcEnergies[iDigit] * ratio ; | |
897 | emcRP->AddDigit( *digit, eDigit ) ; | |
88cb7938 | 898 | } |
9a1398dd | 899 | } |
900 | ||
901 | delete[] fitparameters ; | |
902 | delete[] efit ; | |
903 | ||
904 | } | |
905 | ||
906 | //_____________________________________________________________________________ | |
907 | void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag) | |
908 | { | |
a4e98857 | 909 | // Calculates the Chi square for the cluster unfolding minimization |
9a1398dd | 910 | // Number of parameters, Gradient, Chi squared, parameters, what to do |
911 | ||
88cb7938 | 912 | TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ; |
9a1398dd | 913 | |
88cb7938 | 914 | AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint*>( toMinuit->At(0) ) ; |
915 | TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) ) ; | |
9a1398dd | 916 | |
917 | ||
918 | ||
88cb7938 | 919 | // AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( gMinuit->GetObjectFit() ) ; // EmcRecPoint to fit |
9a1398dd | 920 | |
921 | Int_t * emcDigits = emcRP->GetDigitsList() ; | |
922 | ||
7b7c1533 | 923 | Int_t nOdigits = emcRP->GetDigitsMultiplicity() ; |
9a1398dd | 924 | |
925 | Float_t * emcEnergies = emcRP->GetEnergiesList() ; | |
88cb7938 | 926 | |
927 | const AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ; | |
9a1398dd | 928 | fret = 0. ; |
929 | Int_t iparam ; | |
930 | ||
931 | if(iflag == 2) | |
932 | for(iparam = 0 ; iparam < nPar ; iparam++) | |
933 | Grad[iparam] = 0 ; // Will evaluate gradient | |
934 | ||
935 | Double_t efit ; | |
936 | ||
937 | AliPHOSDigit * digit ; | |
938 | Int_t iDigit ; | |
939 | ||
7b7c1533 | 940 | for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) { |
9a1398dd | 941 | |
88cb7938 | 942 | digit = dynamic_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) ); |
9a1398dd | 943 | |
944 | Int_t relid[4] ; | |
945 | Float_t xDigit ; | |
946 | Float_t zDigit ; | |
947 | ||
948 | geom->AbsToRelNumbering(digit->GetId(), relid) ; | |
949 | ||
950 | geom->RelPosInModule(relid, xDigit, zDigit) ; | |
951 | ||
952 | if(iflag == 2){ // calculate gradient | |
953 | Int_t iParam = 0 ; | |
954 | efit = 0 ; | |
955 | while(iParam < nPar ){ | |
88cb7938 | 956 | Double_t distance = (xDigit - x[iParam]) * (xDigit - x[iParam]) ; |
957 | iParam++ ; | |
958 | distance += (zDigit - x[iParam]) * (zDigit - x[iParam]) ; | |
959 | distance = TMath::Sqrt( distance ) ; | |
960 | iParam++ ; | |
961 | efit += x[iParam] * ShowerShape(distance) ; | |
962 | iParam++ ; | |
9a1398dd | 963 | } |
88cb7938 | 964 | Double_t sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E) |
9a1398dd | 965 | iParam = 0 ; |
966 | while(iParam < nPar ){ | |
88cb7938 | 967 | Double_t xpar = x[iParam] ; |
968 | Double_t zpar = x[iParam+1] ; | |
969 | Double_t epar = x[iParam+2] ; | |
970 | Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ); | |
971 | Double_t shape = sum * ShowerShape(dr) ; | |
972 | Double_t r4 = dr*dr*dr*dr ; | |
973 | Double_t r295 = TMath::Power(dr,2.95) ; | |
974 | Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) + | |
975 | 0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ; | |
976 | ||
977 | Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ; // Derivative over x | |
978 | iParam++ ; | |
979 | Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ; // Derivative over z | |
980 | iParam++ ; | |
981 | Grad[iParam] += shape ; // Derivative over energy | |
982 | iParam++ ; | |
9a1398dd | 983 | } |
984 | } | |
985 | efit = 0; | |
986 | iparam = 0 ; | |
987 | ||
988 | while(iparam < nPar ){ | |
989 | Double_t xpar = x[iparam] ; | |
990 | Double_t zpar = x[iparam+1] ; | |
991 | Double_t epar = x[iparam+2] ; | |
992 | iparam += 3 ; | |
993 | Double_t distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ; | |
994 | distance = TMath::Sqrt(distance) ; | |
995 | efit += epar * ShowerShape(distance) ; | |
996 | } | |
997 | ||
88cb7938 | 998 | fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/emcEnergies[iDigit] ; |
9a1398dd | 999 | // Here we assume, that sigma = sqrt(E) |
1000 | } | |
83974468 | 1001 | |
d15a28e7 | 1002 | } |
1003 | ||
1004 | //____________________________________________________________________________ | |
702ab87e | 1005 | void AliPHOSClusterizerv1::Print(const Option_t *)const |
d15a28e7 | 1006 | { |
d1de15f5 | 1007 | // Print clusterizer parameters |
1008 | ||
21cd0c07 | 1009 | TString message ; |
1010 | TString taskName(GetName()) ; | |
1011 | taskName.ReplaceAll(Version(), "") ; | |
1012 | ||
1013 | if( strcmp(GetName(), "") !=0 ) { | |
9a1398dd | 1014 | // Print parameters |
21cd0c07 | 1015 | message = "\n--------------- %s %s -----------\n" ; |
1016 | message += "Clusterizing digits from the file: %s\n" ; | |
1017 | message += " Branch: %s\n" ; | |
1018 | message += " EMC Clustering threshold = %f\n" ; | |
1019 | message += " EMC Local Maximum cut = %f\n" ; | |
1020 | message += " EMC Logarothmic weight = %f\n" ; | |
1021 | message += " CPV Clustering threshold = %f\n" ; | |
1022 | message += " CPV Local Maximum cut = %f\n" ; | |
1023 | message += " CPV Logarothmic weight = %f\n" ; | |
9a1398dd | 1024 | if(fToUnfold) |
21cd0c07 | 1025 | message += " Unfolding on\n" ; |
9a1398dd | 1026 | else |
21cd0c07 | 1027 | message += " Unfolding off\n" ; |
9a1398dd | 1028 | |
21cd0c07 | 1029 | message += "------------------------------------------------------------------" ; |
9a1398dd | 1030 | } |
21cd0c07 | 1031 | else |
1032 | message = " AliPHOSClusterizerv1 not initialized " ; | |
1033 | ||
351dd634 | 1034 | AliInfo(Form("%s, %s %s %s %s %s %s %s %s %s %s", message.Data(), |
21cd0c07 | 1035 | taskName.Data(), |
1036 | GetTitle(), | |
1037 | taskName.Data(), | |
1038 | GetName(), | |
1039 | fEmcClusteringThreshold, | |
1040 | fEmcLocMaxCut, | |
1041 | fW0, | |
1042 | fCpvClusteringThreshold, | |
1043 | fCpvLocMaxCut, | |
351dd634 | 1044 | fW0CPV )) ; |
9a1398dd | 1045 | } |
21cd0c07 | 1046 | |
1047 | ||
9a1398dd | 1048 | //____________________________________________________________________________ |
a4e98857 | 1049 | void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option) |
1050 | { | |
1051 | // Prints list of RecPoints produced at the current pass of AliPHOSClusterizer | |
9a1398dd | 1052 | |
88cb7938 | 1053 | AliPHOSGetter * gime = AliPHOSGetter::Instance(); |
1054 | ||
1055 | TObjArray * emcRecPoints = gime->EmcRecPoints() ; | |
1056 | TObjArray * cpvRecPoints = gime->CpvRecPoints() ; | |
7b7c1533 | 1057 | |
351dd634 | 1058 | AliInfo(Form("\nevent %d \n Found %d EMC RecPoints and %d CPV RecPoints", |
1059 | gAlice->GetEvNumber(), | |
1060 | emcRecPoints->GetEntriesFast(), | |
1061 | cpvRecPoints->GetEntriesFast() )) ; | |
88cb7938 | 1062 | |
1063 | fRecPointsInRun += emcRecPoints->GetEntriesFast() ; | |
1064 | fRecPointsInRun += cpvRecPoints->GetEntriesFast() ; | |
1065 | ||
11f9c5ff | 1066 | |
9a1398dd | 1067 | if(strstr(option,"all")) { |
709e117a | 1068 | printf("\n EMC clusters \n") ; |
1069 | printf("Index Ene(MeV) Multi Module X Y Z Lambdas_1 Lambda_2 # of prim Primaries list\n") ; | |
9a1398dd | 1070 | Int_t index ; |
7b7c1533 | 1071 | for (index = 0 ; index < emcRecPoints->GetEntries() ; index++) { |
1072 | AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )emcRecPoints->At(index) ; | |
9a1398dd | 1073 | TVector3 locpos; |
1074 | rp->GetLocalPosition(locpos); | |
9a1398dd | 1075 | Float_t lambda[2]; |
1076 | rp->GetElipsAxis(lambda); | |
9a1398dd | 1077 | Int_t * primaries; |
1078 | Int_t nprimaries; | |
1079 | primaries = rp->GetPrimaries(nprimaries); | |
709e117a | 1080 | printf("\n%6d %8.2f %3d %2d %4.1f %4.1f %4.1f %4f %4f %2d : ", |
11f9c5ff | 1081 | rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetPHOSMod(), |
1082 | locpos.X(), locpos.Y(), locpos.Z(), lambda[0], lambda[1], nprimaries) ; | |
3bf72d32 | 1083 | |
21cd0c07 | 1084 | for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) { |
709e117a | 1085 | printf("%d ", primaries[iprimary] ) ; |
21cd0c07 | 1086 | } |
709e117a | 1087 | printf("\n") ; |
1088 | } | |
1089 | ||
1090 | //Now plot CPV recPoints | |
1091 | printf("\n CPV clusters \n") ; | |
1092 | printf("Index Ene(MeV) Module X Y Z \n") ; | |
1093 | for (index = 0 ; index < cpvRecPoints->GetEntries() ; index++) { | |
2bb500e5 | 1094 | AliPHOSCpvRecPoint * rp = (AliPHOSCpvRecPoint * )cpvRecPoints->At(index) ; |
1095 | ||
1096 | TVector3 locpos; | |
1097 | rp->GetLocalPosition(locpos); | |
1098 | ||
1099 | printf("\n%6d %8.2f %2d %4.1f %4.1f %4.1f \n", | |
1100 | rp->GetIndexInList(), rp->GetEnergy(), rp->GetPHOSMod(), | |
1101 | locpos.X(), locpos.Y(), locpos.Z()) ; | |
88cb7938 | 1102 | } |
9a1398dd | 1103 | } |
d15a28e7 | 1104 | } |
9a1398dd | 1105 |