]>
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 | ||
d72dfbc3 | 18 | /* $Log: |
19 | 1 October 2000. Yuri Kharlov: | |
20 | AreNeighbours() | |
21 | PPSD upper layer is considered if number of layers>1 | |
22 | ||
23 | 18 October 2000. Yuri Kharlov: | |
24 | AliPHOSClusterizerv1() | |
25 | CPV clusterizing parameters added | |
26 | ||
27 | MakeClusters() | |
28 | After first PPSD digit remove EMC digits only once | |
29 | */ | |
9a1398dd | 30 | //*-- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko (SUBATECH & Kurchatov Institute) |
d15a28e7 | 31 | ////////////////////////////////////////////////////////////////////////////// |
9a1398dd | 32 | // Clusterization class. Performs clusterization (collects neighbouring active cells) and |
a4e98857 | 33 | // unfolds the clusters having several local maxima. |
34 | // Results are stored in TreeR#, branches PHOSEmcRP (EMC recPoints), | |
9a1398dd | 35 | // PHOSCpvRP (CPV RecPoints) and AliPHOSClusterizer (Clusterizer with all |
f035f6ce | 36 | // parameters including input digits branch title, thresholds etc.) |
a4e98857 | 37 | // This TTask is normally called from Reconstructioner, but can as well be used in |
38 | // standalone mode. | |
39 | // Use Case: | |
40 | // root [0] AliPHOSClusterizerv1 * cl = new AliPHOSClusterizerv1("galice.root") | |
41 | // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated | |
9a1398dd | 42 | // //reads gAlice from header file "..." |
a4e98857 | 43 | // root [1] cl->ExecuteTask() |
9a1398dd | 44 | // //finds RecPoints in all events stored in galice.root |
a4e98857 | 45 | // root [2] cl->SetDigitsBranch("digits2") |
f035f6ce | 46 | // //sets another title for Digitis (input) branch |
a4e98857 | 47 | // root [3] cl->SetRecPointsBranch("recp2") |
f035f6ce | 48 | // //sets another title four output branches |
a4e98857 | 49 | // root [4] cl->SetEmcLocalMaxCut(0.03) |
9a1398dd | 50 | // //set clusterization parameters |
a4e98857 | 51 | // root [5] cl->ExecuteTask("deb all time") |
9a1398dd | 52 | // //once more finds RecPoints options are |
53 | // // deb - print number of found rec points | |
54 | // // deb all - print number of found RecPoints and some their characteristics | |
55 | // // time - print benchmarking results | |
ed4205d8 | 56 | |
d15a28e7 | 57 | // --- ROOT system --- |
58 | ||
9a1398dd | 59 | #include "TROOT.h" |
60 | #include "TFile.h" | |
61 | #include "TFolder.h" | |
d15a28e7 | 62 | #include "TMath.h" |
9a1398dd | 63 | #include "TMinuit.h" |
64 | #include "TTree.h" | |
65 | #include "TSystem.h" | |
66 | #include "TBenchmark.h" | |
d15a28e7 | 67 | |
68 | // --- Standard library --- | |
69 | ||
de9ec31b | 70 | #include <iostream.h> |
9a1398dd | 71 | #include <iomanip.h> |
d15a28e7 | 72 | |
73 | // --- AliRoot header files --- | |
74 | ||
75 | #include "AliPHOSClusterizerv1.h" | |
9a1398dd | 76 | #include "AliPHOSCpvRecPoint.h" |
d15a28e7 | 77 | #include "AliPHOSDigit.h" |
9a1398dd | 78 | #include "AliPHOSDigitizer.h" |
d15a28e7 | 79 | #include "AliPHOSEmcRecPoint.h" |
9a1398dd | 80 | #include "AliPHOS.h" |
d15a28e7 | 81 | #include "AliPHOSPpsdRecPoint.h" |
9a1398dd | 82 | #include "AliRun.h" |
d15a28e7 | 83 | |
84 | ClassImp(AliPHOSClusterizerv1) | |
f035f6ce | 85 | |
d15a28e7 | 86 | //____________________________________________________________________________ |
9a1398dd | 87 | AliPHOSClusterizerv1::AliPHOSClusterizerv1():AliPHOSClusterizer() |
d15a28e7 | 88 | { |
f035f6ce | 89 | // default ctor (to be used mainly by Streamer) |
90 | SetName("AliPHOSClusterizer"); | |
9a1398dd | 91 | SetTitle("Version 1") ; |
f035f6ce | 92 | |
9a1398dd | 93 | fNumberOfCpvClusters = 0 ; |
d15a28e7 | 94 | fNumberOfEmcClusters = 0 ; |
f035f6ce | 95 | |
9a1398dd | 96 | fCpvClusteringThreshold = 0.0; |
c3c187e5 | 97 | fEmcClusteringThreshold = 0.2; |
c3c187e5 | 98 | fPpsdClusteringThreshold = 0.0000002 ; |
9a1398dd | 99 | |
100 | fEmcLocMaxCut = 0.03 ; | |
101 | fCpvLocMaxCut = 0.03 ; | |
102 | ||
103 | fW0 = 4.5 ; | |
104 | fW0CPV = 4.0 ; | |
f035f6ce | 105 | |
9a1398dd | 106 | fGeom = 0 ; |
107 | ||
108 | fDigits = 0 ; | |
109 | fDigitizer = 0 ; | |
110 | fEmcRecPoints = 0 ; | |
111 | fCpvRecPoints = 0 ; | |
112 | ||
113 | fIsInitialized = kFALSE ; | |
114 | ||
115 | } | |
116 | //____________________________________________________________________________ | |
f035f6ce | 117 | AliPHOSClusterizerv1::AliPHOSClusterizerv1(const char* headerFile,const char* digitsFile):AliPHOSClusterizer() |
9a1398dd | 118 | { |
a4e98857 | 119 | // ctor with the indication of the file where header Tree and digits Tree are stored |
9a1398dd | 120 | SetName("AliPHOSClusterizer"); |
121 | SetTitle("Version 1") ; | |
122 | ||
123 | fNumberOfCpvClusters = 0 ; | |
124 | fNumberOfEmcClusters = 0 ; | |
f035f6ce | 125 | |
ed4205d8 | 126 | fCpvClusteringThreshold = 0.0; |
9a1398dd | 127 | fEmcClusteringThreshold = 0.2; |
128 | fPpsdClusteringThreshold = 0.0000002 ; | |
129 | ||
130 | fEmcLocMaxCut = 0.03 ; | |
131 | fCpvLocMaxCut = 0.03 ; | |
132 | ||
d15a28e7 | 133 | fW0 = 4.5 ; |
c3c187e5 | 134 | fW0CPV = 4.0 ; |
9a1398dd | 135 | |
136 | fToUnfold = kTRUE ; | |
137 | ||
138 | fHeaderFileName = headerFile ; | |
139 | fDigitsBranchTitle = digitsFile ; | |
140 | ||
141 | TFile * file = (TFile*) gROOT->GetFile(fHeaderFileName.Data() ) ; | |
142 | ||
143 | if(file == 0){ | |
f035f6ce | 144 | if(fHeaderFileName.Contains("rfio")) // if we read file using HPSS |
145 | file = TFile::Open(fHeaderFileName.Data(),"update") ; | |
146 | else | |
147 | file = new TFile(fHeaderFileName.Data(),"update") ; | |
9a1398dd | 148 | gAlice = (AliRun *) file->Get("gAlice") ; |
149 | } | |
150 | ||
151 | AliPHOS * phos = (AliPHOS *) gAlice->GetDetector("PHOS") ; | |
152 | fGeom = AliPHOSGeometry::GetInstance(phos->GetGeometry()->GetName(),phos->GetGeometry()->GetTitle() ); | |
153 | ||
154 | fDigits = new TClonesArray("AliPHOSDigit",10) ; | |
155 | fDigitizer = new AliPHOSDigitizer() ; | |
156 | fEmcRecPoints = new TObjArray(200) ; | |
157 | fCpvRecPoints = new TObjArray(200) ; | |
7f5d510c | 158 | fEmcRecPoints->SetOwner(); // This lets Clear() really detete rec.points in array |
159 | fCpvRecPoints->SetOwner(); | |
9a1398dd | 160 | |
161 | if(!gMinuit) gMinuit = new TMinuit(100) ; | |
162 | ||
163 | // add Task to //root/Tasks folder | |
164 | TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ; | |
165 | roottasks->Add(this) ; | |
166 | ||
167 | ||
168 | fIsInitialized = kTRUE ; | |
169 | ||
170 | } | |
171 | //____________________________________________________________________________ | |
a4e98857 | 172 | void AliPHOSClusterizerv1::Exec(Option_t * option) |
173 | { | |
174 | // Steering function | |
9a1398dd | 175 | |
176 | if(!fIsInitialized) Init() ; | |
177 | ||
178 | if(strstr(option,"tim")) | |
179 | gBenchmark->Start("PHOSClusterizer"); | |
180 | ||
181 | Int_t nEvents = (Int_t) gAlice->TreeE()->GetEntries() ; | |
182 | ||
183 | for(fEvent = 0;fEvent< nEvents; fEvent++){ | |
184 | if(!ReadDigits()) //reads digits for event fEvent | |
185 | return; | |
186 | MakeClusters() ; | |
187 | ||
188 | if(fToUnfold) MakeUnfolding() ; | |
189 | WriteRecPoints() ; | |
190 | if(strstr(option,"deb")) | |
191 | PrintRecPoints(option) ; | |
192 | } | |
193 | ||
194 | if(strstr(option,"tim")){ | |
195 | gBenchmark->Stop("PHOSClusterizer"); | |
196 | cout << "AliPHOSClusterizer:" << endl ; | |
197 | cout << " took " << gBenchmark->GetCpuTime("PHOSClusterizer") << " seconds for Clusterizing " | |
198 | << gBenchmark->GetCpuTime("PHOSClusterizer")/nEvents << " seconds per event " << endl ; | |
199 | cout << endl ; | |
200 | } | |
201 | ||
202 | ||
203 | } | |
204 | ||
205 | //____________________________________________________________________________ | |
206 | Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, int * maxAt, Float_t * maxAtEnergy, | |
207 | Int_t nPar, Float_t * fitparameters) | |
208 | { | |
209 | // Calls TMinuit to fit the energy distribution of a cluster with several maxima | |
a4e98857 | 210 | // The initial values for fitting procedure are set equal to the positions of local maxima. |
f035f6ce | 211 | // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers |
9a1398dd | 212 | |
213 | gMinuit->mncler(); // Reset Minuit's list of paramters | |
214 | gMinuit->SetPrintLevel(-1) ; // No Printout | |
215 | gMinuit->SetFCN(AliPHOSClusterizerv1::UnfoldingChiSquare) ; | |
216 | // To set the address of the minimization function | |
217 | ||
218 | TList * toMinuit = new TList(); | |
219 | toMinuit->AddAt(emcRP,0) ; | |
220 | toMinuit->AddAt(fDigits,1) ; | |
221 | ||
222 | gMinuit->SetObjectFit(toMinuit) ; // To tranfer pointer to UnfoldingChiSquare | |
223 | ||
224 | // filling initial values for fit parameters | |
225 | AliPHOSDigit * digit ; | |
226 | ||
227 | Int_t ierflg = 0; | |
228 | Int_t index = 0 ; | |
229 | Int_t nDigits = (Int_t) nPar / 3 ; | |
230 | ||
231 | Int_t iDigit ; | |
232 | ||
233 | ||
234 | for(iDigit = 0; iDigit < nDigits; iDigit++){ | |
235 | digit = (AliPHOSDigit *) maxAt[iDigit]; | |
236 | ||
237 | Int_t relid[4] ; | |
238 | Float_t x ; | |
239 | Float_t z ; | |
240 | fGeom->AbsToRelNumbering(digit->GetId(), relid) ; | |
241 | fGeom->RelPosInModule(relid, x, z) ; | |
242 | ||
243 | Float_t energy = maxAtEnergy[iDigit] ; | |
244 | ||
245 | gMinuit->mnparm(index, "x", x, 0.1, 0, 0, ierflg) ; | |
246 | index++ ; | |
247 | if(ierflg != 0){ | |
248 | cout << "PHOS Unfolding> Unable to set initial value for fit procedure : x = " << x << endl ; | |
249 | return kFALSE; | |
250 | } | |
251 | gMinuit->mnparm(index, "z", z, 0.1, 0, 0, ierflg) ; | |
252 | index++ ; | |
253 | if(ierflg != 0){ | |
254 | cout << "PHOS Unfolding> Unable to set initial value for fit procedure : z = " << z << endl ; | |
255 | return kFALSE; | |
256 | } | |
257 | gMinuit->mnparm(index, "Energy", energy , 0.05*energy, 0., 4.*energy, ierflg) ; | |
258 | index++ ; | |
259 | if(ierflg != 0){ | |
260 | cout << "PHOS Unfolding> Unable to set initial value for fit procedure : energy = " << energy << endl ; | |
261 | return kFALSE; | |
262 | } | |
263 | } | |
264 | ||
265 | Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly | |
266 | // depends on it. | |
267 | Double_t p1 = 1.0 ; | |
268 | Double_t p2 = 0.0 ; | |
269 | ||
270 | gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ; // force TMinuit to reduce function calls | |
271 | gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient | |
272 | gMinuit->SetMaxIterations(5); | |
273 | gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ; // No Warnings | |
274 | ||
275 | gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize | |
276 | ||
277 | if(ierflg == 4){ // Minimum not found | |
278 | cout << "PHOS Unfolding> Fit not converged, cluster abandoned "<< endl ; | |
279 | return kFALSE ; | |
280 | } | |
281 | for(index = 0; index < nPar; index++){ | |
282 | Double_t err ; | |
283 | Double_t val ; | |
284 | gMinuit->GetParameter(index, val, err) ; // Returns value and error of parameter index | |
285 | fitparameters[index] = val ; | |
286 | } | |
287 | ||
288 | delete toMinuit ; | |
289 | return kTRUE; | |
c3c187e5 | 290 | |
d15a28e7 | 291 | } |
292 | ||
293 | //____________________________________________________________________________ | |
a4e98857 | 294 | void AliPHOSClusterizerv1::Init() |
295 | { | |
f035f6ce | 296 | //Make all memory allocations which can not be done in default constructor. |
9a1398dd | 297 | if(!fIsInitialized){ |
298 | if(fHeaderFileName.IsNull()) | |
299 | fHeaderFileName = "galice.root" ; | |
300 | ||
a4e98857 | 301 | |
9a1398dd | 302 | TFile * file = (TFile*) gROOT->GetFile(fHeaderFileName.Data() ) ; |
a4e98857 | 303 | |
9a1398dd | 304 | if(file == 0){ |
305 | file = new TFile(fHeaderFileName.Data(),"update") ; | |
306 | gAlice = (AliRun *) file->Get("gAlice") ; | |
307 | } | |
a4e98857 | 308 | |
9a1398dd | 309 | AliPHOS * phos = (AliPHOS *) gAlice->GetDetector("PHOS") ; |
310 | fGeom = AliPHOSGeometry::GetInstance(phos->GetGeometry()->GetName(),phos->GetGeometry()->GetTitle() ); | |
a4e98857 | 311 | |
9a1398dd | 312 | fDigits = new TClonesArray("AliPHOSDigit",10) ; |
313 | fDigitizer = new AliPHOSDigitizer() ; | |
314 | fEmcRecPoints = new TObjArray(200) ; | |
315 | fCpvRecPoints = new TObjArray(200) ; | |
7f5d510c | 316 | fEmcRecPoints->SetOwner(); // This lets Clear() really detete rec.points in array |
317 | fCpvRecPoints->SetOwner(); | |
a4e98857 | 318 | |
9a1398dd | 319 | if(!gMinuit) gMinuit = new TMinuit(100) ; |
320 | ||
321 | // add Task to //root/Tasks folder | |
322 | TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ; | |
323 | roottasks->Add(this) ; | |
324 | ||
325 | fIsInitialized = kTRUE ; | |
326 | } | |
327 | } | |
328 | //____________________________________________________________________________ | |
329 | Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)const | |
d15a28e7 | 330 | { |
b2a60966 | 331 | // Gives the neighbourness of two digits = 0 are not neighbour but continue searching |
332 | // = 1 are neighbour | |
333 | // = 2 are not neighbour but do not continue searching | |
a4e98857 | 334 | // neighbours are defined as digits having at least a common vertex |
b2a60966 | 335 | // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster |
336 | // which is compared to a digit (d2) not yet in a cluster | |
337 | ||
d15a28e7 | 338 | Int_t rv = 0 ; |
339 | ||
d15a28e7 | 340 | Int_t relid1[4] ; |
6b9b768e | 341 | fGeom->AbsToRelNumbering(d1->GetId(), relid1) ; |
d15a28e7 | 342 | |
343 | Int_t relid2[4] ; | |
6b9b768e | 344 | fGeom->AbsToRelNumbering(d2->GetId(), relid2) ; |
d15a28e7 | 345 | |
346 | if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module and the same PPSD Module | |
92862013 | 347 | Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ; |
348 | Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ; | |
d15a28e7 | 349 | |
92862013 | 350 | if (( coldiff <= 1 ) && ( rowdiff <= 1 )){ |
d15a28e7 | 351 | rv = 1 ; |
352 | } | |
353 | else { | |
354 | if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1)) | |
355 | rv = 2; // Difference in row numbers is too large to look further | |
356 | } | |
357 | ||
358 | } | |
359 | else { | |
360 | ||
361 | if( (relid1[0] < relid2[0]) || (relid1[1] < relid2[1]) ) | |
362 | rv=2 ; | |
363 | ||
364 | } | |
d72dfbc3 | 365 | |
366 | //Do NOT clusterize upper PPSD | |
ed4205d8 | 367 | if( IsInPpsd(d1) && IsInPpsd(d2) && |
368 | relid1[1] > 0 && | |
d72dfbc3 | 369 | relid1[1] < fGeom->GetNumberOfPadsPhi()*fGeom->GetNumberOfPadsPhi() ) rv = 2 ; |
370 | ||
d15a28e7 | 371 | return rv ; |
372 | } | |
373 | ||
d15a28e7 | 374 | |
375 | //____________________________________________________________________________ | |
9a1398dd | 376 | Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit) const |
d15a28e7 | 377 | { |
b2a60966 | 378 | // Tells if (true) or not (false) the digit is in a PHOS-EMC module |
379 | ||
9f616d61 | 380 | Bool_t rv = kFALSE ; |
d15a28e7 | 381 | |
d15a28e7 | 382 | Int_t relid[4] ; |
6b9b768e | 383 | fGeom->AbsToRelNumbering(digit->GetId(), relid) ; |
d15a28e7 | 384 | |
ed4205d8 | 385 | if ( relid[1] == 0 ) rv = kTRUE; |
386 | ||
387 | return rv ; | |
388 | } | |
389 | ||
390 | //____________________________________________________________________________ | |
9a1398dd | 391 | Bool_t AliPHOSClusterizerv1::IsInPpsd(AliPHOSDigit * digit) const |
ed4205d8 | 392 | { |
fad3e5b9 | 393 | // Tells if (true) or not (false) the digit is in a PHOS-PPSD module |
ed4205d8 | 394 | |
395 | Bool_t rv = kFALSE ; | |
396 | ||
397 | Int_t relid[4] ; | |
398 | fGeom->AbsToRelNumbering(digit->GetId(), relid) ; | |
399 | ||
400 | if ( relid[1] > 0 && relid[0] > fGeom->GetNCPVModules() ) rv = kTRUE; | |
401 | ||
402 | return rv ; | |
403 | } | |
404 | ||
405 | //____________________________________________________________________________ | |
9a1398dd | 406 | Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit) const |
ed4205d8 | 407 | { |
fad3e5b9 | 408 | // Tells if (true) or not (false) the digit is in a PHOS-CPV module |
ed4205d8 | 409 | |
410 | Bool_t rv = kFALSE ; | |
411 | ||
412 | Int_t relid[4] ; | |
413 | fGeom->AbsToRelNumbering(digit->GetId(), relid) ; | |
414 | ||
415 | if ( relid[1] > 0 && relid[0] <= fGeom->GetNCPVModules() ) rv = kTRUE; | |
d15a28e7 | 416 | |
417 | return rv ; | |
418 | } | |
d15a28e7 | 419 | //____________________________________________________________________________ |
a4e98857 | 420 | Bool_t AliPHOSClusterizerv1::ReadDigits() |
421 | { | |
422 | // reads digitis with specified title from TreeD | |
c161df70 | 423 | |
b73f246d | 424 | fNumberOfEmcClusters = 0 ; |
b73f246d | 425 | fNumberOfCpvClusters = 0 ; |
c161df70 | 426 | |
9a1398dd | 427 | // Get Digits Tree header from file |
9a1398dd | 428 | gAlice->GetEvent(fEvent) ; |
429 | gAlice->SetEvent(fEvent) ; | |
430 | ||
f035f6ce | 431 | TTree * treeD = gAlice->TreeD() ; |
9a1398dd | 432 | |
433 | if(treeD==0){ | |
f035f6ce | 434 | char treeName[20]; |
435 | sprintf(treeName,"TreeD%d",fEvent); | |
9a1398dd | 436 | cout << "Error in AliPHOSClusterizerv1 : no "<<treeName << endl ; |
437 | cout << " Do nothing " << endl ; | |
438 | return kFALSE ; | |
439 | } | |
440 | ||
441 | TBranch * digitsBranch = 0; | |
442 | TBranch * digitizerBranch = 0; | |
443 | ||
444 | TObjArray * branches = treeD->GetListOfBranches() ; | |
445 | Int_t ibranch; | |
446 | Bool_t phosNotFound = kTRUE ; | |
447 | Bool_t digitizerNotFound = kTRUE ; | |
448 | ||
449 | for(ibranch = 0;ibranch <branches->GetEntries();ibranch++){ | |
450 | ||
451 | if(phosNotFound){ | |
452 | digitsBranch=(TBranch *) branches->At(ibranch) ; | |
453 | if( fDigitsBranchTitle.CompareTo(digitsBranch->GetTitle())==0 ) | |
454 | if( strcmp(digitsBranch->GetName(),"PHOS") == 0) | |
455 | phosNotFound = kFALSE ; | |
456 | } | |
457 | ||
458 | if(digitizerNotFound){ | |
459 | digitizerBranch = (TBranch *) branches->At(ibranch) ; | |
460 | if( fDigitsBranchTitle.CompareTo(digitizerBranch->GetTitle()) == 0) | |
461 | if( strcmp(digitizerBranch->GetName(),"AliPHOSDigitizer") == 0) | |
462 | digitizerNotFound = kFALSE ; | |
463 | } | |
464 | ||
465 | } | |
466 | ||
467 | if(digitizerNotFound || phosNotFound){ | |
468 | cout << "ERROR in AliPHOSClusterizerv1: " << endl ; | |
469 | cout << " Can't find Branch with digits or Digitizer "<< endl ; ; | |
470 | cout << " Do nothing" <<endl ; | |
471 | return kFALSE ; | |
472 | } | |
473 | ||
474 | digitsBranch->SetAddress(&fDigits) ; | |
475 | digitizerBranch->SetAddress(&fDigitizer) ; | |
761e34c0 | 476 | |
477 | digitsBranch->GetEntry(0) ; | |
478 | digitizerBranch->GetEntry(0) ; | |
9a1398dd | 479 | |
480 | fPedestal = fDigitizer->GetPedestal() ; | |
481 | fSlope = fDigitizer->GetSlope() ; | |
482 | return kTRUE ; | |
483 | } | |
484 | ||
485 | //____________________________________________________________________________ | |
a4e98857 | 486 | void AliPHOSClusterizerv1::WriteRecPoints() |
487 | { | |
488 | // Checks, if PHOSEmcRP etc. branches with given title already exist, | |
489 | // if yes exits without writing, | |
490 | // else creates new branches with given title | |
491 | // fills and writes into TreeR. | |
9a1398dd | 492 | |
493 | Int_t index ; | |
494 | //Evaluate poisition, dispersion and other RecPoint properties... | |
495 | for(index = 0; index < fEmcRecPoints->GetEntries(); index++) | |
496 | ((AliPHOSEmcRecPoint *)fEmcRecPoints->At(index))->EvalAll(fW0,fDigits) ; | |
497 | ||
498 | fEmcRecPoints->Sort() ; | |
499 | ||
500 | for(index = 0; index < fEmcRecPoints->GetEntries(); index++) | |
501 | ((AliPHOSEmcRecPoint *)fEmcRecPoints->At(index))->SetIndexInList(index) ; | |
502 | ||
761e34c0 | 503 | fEmcRecPoints->Expand(fEmcRecPoints->GetEntriesFast()) ; |
504 | ||
9a1398dd | 505 | //Now the same for CPV |
506 | for(index = 0; index < fCpvRecPoints->GetEntries(); index++) | |
507 | ((AliPHOSRecPoint *)fCpvRecPoints->At(index))->EvalAll(fW0CPV,fDigits) ; | |
508 | ||
509 | fCpvRecPoints->Sort() ; | |
510 | ||
511 | for(index = 0; index < fCpvRecPoints->GetEntries(); index++) | |
512 | ((AliPHOSRecPoint *)fCpvRecPoints->At(index))->SetIndexInList(index) ; | |
513 | ||
761e34c0 | 514 | fCpvRecPoints->Expand(fCpvRecPoints->GetEntriesFast()) ; |
515 | ||
9a1398dd | 516 | if(gAlice->TreeR()==0) |
517 | gAlice->MakeTree("R") ; | |
518 | ||
519 | ||
520 | //Check, if branches already exist | |
521 | TBranch * emcBranch = 0; | |
522 | TBranch * cpvBranch = 0; | |
523 | TBranch * clusterizerBranch = 0; | |
524 | ||
525 | TObjArray * branches = gAlice->TreeR()->GetListOfBranches() ; | |
526 | Int_t ibranch; | |
527 | Bool_t emcNotFound = kTRUE ; | |
528 | Bool_t cpvNotFound = kTRUE ; | |
529 | Bool_t clusterizerNotFound = kTRUE ; | |
530 | ||
531 | for(ibranch = 0;ibranch <branches->GetEntries();ibranch++){ | |
532 | ||
533 | if(emcNotFound){ | |
534 | emcBranch=(TBranch *) branches->At(ibranch) ; | |
535 | if( fRecPointsBranchTitle.CompareTo(emcBranch->GetTitle())==0 ){ | |
536 | if( strcmp(emcBranch->GetName(),"PHOSEmcRP") == 0) { | |
537 | emcNotFound = kFALSE ; | |
538 | } | |
539 | } | |
540 | } | |
541 | if(cpvNotFound){ | |
542 | cpvBranch=(TBranch *) branches->At(ibranch) ; | |
543 | if( fRecPointsBranchTitle.CompareTo(cpvBranch->GetTitle())==0 ){ | |
544 | if( strcmp(cpvBranch->GetName(),"PHOSCpvRP") == 0) { | |
545 | cpvNotFound = kFALSE ; | |
546 | } | |
547 | } | |
548 | } | |
549 | ||
550 | if(clusterizerNotFound){ | |
551 | clusterizerBranch = (TBranch *) branches->At(ibranch) ; | |
552 | if( fRecPointsBranchTitle.CompareTo(clusterizerBranch->GetTitle()) == 0){ | |
553 | if( strcmp(clusterizerBranch->GetName(),"AliPHOSClusterizer") == 0) { | |
554 | clusterizerNotFound = kFALSE ; | |
555 | } | |
556 | } | |
557 | } | |
558 | ||
559 | } | |
560 | ||
561 | if(!(clusterizerNotFound && emcNotFound && cpvNotFound)){ | |
562 | cout << "AliPHOSClusterizer error" << endl; | |
563 | cout << " Branches PHOSEmcRP, PHOSCpvRP and AliPHOSClusterizer " << endl ; | |
564 | cout << " with title '" << fRecPointsBranchTitle.Data() <<"' already exist" << endl ; | |
565 | cout << " can not overwrite " << endl ; | |
566 | return ; | |
567 | } | |
568 | ||
569 | //Make branches in TreeR for RecPoints and Clusterizer | |
570 | char * filename = 0; | |
571 | if(gSystem->Getenv("CONFIG_SPLIT_FILE")!=0){ //generating file name | |
572 | filename = new char[strlen(gAlice->GetBaseFile())+20] ; | |
573 | sprintf(filename,"%s/PHOS.Reco.root",gAlice->GetBaseFile()) ; | |
574 | } | |
575 | ||
576 | //Make new branches | |
577 | TDirectory *cwd = gDirectory; | |
578 | ||
579 | //First EMC | |
580 | Int_t bufferSize = 32000 ; | |
581 | Int_t splitlevel = 0 ; | |
582 | emcBranch = gAlice->TreeR()->Branch("PHOSEmcRP","TObjArray",&fEmcRecPoints,bufferSize,splitlevel); | |
583 | emcBranch->SetTitle(fRecPointsBranchTitle.Data()); | |
584 | if (filename) { | |
585 | emcBranch->SetFile(filename); | |
586 | TIter next( emcBranch->GetListOfBranches()); | |
761e34c0 | 587 | TBranch * sb ; |
588 | while ((sb=(TBranch*)next())) { | |
589 | sb->SetFile(filename); | |
9a1398dd | 590 | } |
591 | cwd->cd(); | |
592 | } | |
593 | ||
594 | //Now CPV branch | |
595 | cpvBranch = gAlice->TreeR()->Branch("PHOSCpvRP","TObjArray",&fCpvRecPoints,bufferSize,splitlevel); | |
596 | cpvBranch->SetTitle(fRecPointsBranchTitle.Data()); | |
597 | if (filename) { | |
598 | cpvBranch->SetFile(filename); | |
599 | TIter next( cpvBranch->GetListOfBranches()); | |
761e34c0 | 600 | TBranch * sb; |
601 | while ((sb=(TBranch*)next())) { | |
602 | sb->SetFile(filename); | |
9a1398dd | 603 | } |
604 | cwd->cd(); | |
605 | } | |
606 | ||
607 | //And Finally clusterizer branch | |
608 | AliPHOSClusterizerv1 * cl = this ; | |
609 | clusterizerBranch = gAlice->TreeR()->Branch("AliPHOSClusterizer","AliPHOSClusterizerv1", | |
610 | &cl,bufferSize,splitlevel); | |
611 | clusterizerBranch->SetTitle(fRecPointsBranchTitle.Data()); | |
612 | if (filename) { | |
613 | clusterizerBranch->SetFile(filename); | |
614 | TIter next( clusterizerBranch->GetListOfBranches()); | |
761e34c0 | 615 | TBranch * sb ; |
616 | while ((sb=(TBranch*)next())) { | |
617 | sb->SetFile(filename); | |
9a1398dd | 618 | } |
619 | cwd->cd(); | |
620 | } | |
621 | ||
9a1398dd | 622 | |
761e34c0 | 623 | emcBranch->Fill() ; |
624 | cpvBranch->Fill() ; | |
625 | clusterizerBranch->Fill() ; | |
9a1398dd | 626 | gAlice->TreeR()->Write(0,kOverwrite) ; |
627 | ||
628 | } | |
629 | ||
630 | //____________________________________________________________________________ | |
631 | void AliPHOSClusterizerv1::MakeClusters() | |
632 | { | |
633 | // Steering method to construct the clusters stored in a list of Reconstructed Points | |
634 | // A cluster is defined as a list of neighbour digits | |
635 | fEmcRecPoints->Clear() ; | |
636 | fCpvRecPoints->Clear() ; | |
637 | ||
d15a28e7 | 638 | |
639 | // Clusterization starts | |
9a1398dd | 640 | TClonesArray * digits = (TClonesArray*)fDigits->Clone() ; |
641 | ||
642 | TIter nextdigit(digits) ; | |
d15a28e7 | 643 | AliPHOSDigit * digit ; |
92862013 | 644 | Bool_t notremoved = kTRUE ; |
6ad0bfa0 | 645 | |
d15a28e7 | 646 | while ( (digit = (AliPHOSDigit *)nextdigit()) ) { // scan over the list of digits |
9a1398dd | 647 | AliPHOSRecPoint * clu = 0 ; |
c161df70 | 648 | |
9a1398dd | 649 | TArrayI clusterdigitslist(1000) ; |
d15a28e7 | 650 | Int_t index ; |
f2bc1b87 | 651 | |
b73f246d | 652 | if (( IsInEmc (digit) && Calibrate(digit->GetAmp()) > fEmcClusteringThreshold ) || |
653 | ( IsInPpsd(digit) && Calibrate(digit->GetAmp()) > fPpsdClusteringThreshold ) || | |
654 | ( IsInCpv (digit) && Calibrate(digit->GetAmp()) > fCpvClusteringThreshold ) ) { | |
f2bc1b87 | 655 | |
d15a28e7 | 656 | Int_t iDigitInCluster = 0 ; |
657 | ||
6ad0bfa0 | 658 | if ( IsInEmc(digit) ) { |
83974468 | 659 | // start a new EMC RecPoint |
9a1398dd | 660 | if(fNumberOfEmcClusters >= fEmcRecPoints->GetSize()) fEmcRecPoints->Expand(2*fNumberOfEmcClusters+1) ; |
7f5d510c | 661 | fEmcRecPoints->AddAt(new AliPHOSEmcRecPoint(), fNumberOfEmcClusters) ; |
9a1398dd | 662 | clu = (AliPHOSEmcRecPoint *) fEmcRecPoints->At(fNumberOfEmcClusters) ; |
7f5d510c | 663 | fNumberOfEmcClusters++ ; |
9f616d61 | 664 | clu->AddDigit(*digit, Calibrate(digit->GetAmp())) ; |
9a1398dd | 665 | clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ; |
9f616d61 | 666 | iDigitInCluster++ ; |
9a1398dd | 667 | digits->Remove(digit) ; |
ca77b032 | 668 | |
d72dfbc3 | 669 | } else { |
83974468 | 670 | |
fad3e5b9 | 671 | // start a new PPSD/CPV cluster |
9a1398dd | 672 | if(fNumberOfCpvClusters >= fCpvRecPoints->GetSize()) fCpvRecPoints->Expand(2*fNumberOfCpvClusters+1); |
673 | if(IsInPpsd(digit)) | |
674 | fCpvRecPoints->AddAt(new AliPHOSPpsdRecPoint(),fNumberOfCpvClusters) ; | |
fad3e5b9 | 675 | else |
9a1398dd | 676 | fCpvRecPoints->AddAt(new AliPHOSCpvRecPoint(), fNumberOfCpvClusters) ; |
677 | clu = (AliPHOSPpsdRecPoint *) fCpvRecPoints->At(fNumberOfCpvClusters) ; | |
678 | fNumberOfCpvClusters++ ; | |
fad3e5b9 | 679 | |
6ad0bfa0 | 680 | clu->AddDigit(*digit, Calibrate(digit->GetAmp()) ) ; |
9a1398dd | 681 | clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ; |
9f616d61 | 682 | iDigitInCluster++ ; |
9a1398dd | 683 | digits->Remove(digit) ; |
d15a28e7 | 684 | nextdigit.Reset() ; |
685 | ||
9f616d61 | 686 | // Here we remove resting EMC digits, which cannot make cluster |
fad3e5b9 | 687 | |
92862013 | 688 | if( notremoved ) { |
9f616d61 | 689 | while( ( digit = (AliPHOSDigit *)nextdigit() ) ) { |
9f616d61 | 690 | if( IsInEmc(digit) ) |
9a1398dd | 691 | digits->Remove(digit) ; |
d15a28e7 | 692 | else |
693 | break ; | |
d72dfbc3 | 694 | } |
695 | notremoved = kFALSE ; | |
696 | } | |
d15a28e7 | 697 | |
698 | } // else | |
699 | ||
700 | nextdigit.Reset() ; | |
fad3e5b9 | 701 | |
d15a28e7 | 702 | AliPHOSDigit * digitN ; |
703 | index = 0 ; | |
9f616d61 | 704 | while (index < iDigitInCluster){ // scan over digits already in cluster |
9a1398dd | 705 | digit = (AliPHOSDigit*)fDigits->At(clusterdigitslist[index]) ; |
9f616d61 | 706 | index++ ; |
d15a28e7 | 707 | while ( (digitN = (AliPHOSDigit *)nextdigit()) ) { // scan over the reduced list of digits |
d72dfbc3 | 708 | Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!! |
d15a28e7 | 709 | switch (ineb ) { |
9f616d61 | 710 | case 0 : // not a neighbour |
d72dfbc3 | 711 | break ; |
9f616d61 | 712 | case 1 : // are neighbours |
83974468 | 713 | clu->AddDigit(*digitN, Calibrate( digitN->GetAmp() ) ) ; |
9a1398dd | 714 | clusterdigitslist[iDigitInCluster] = digitN->GetIndexInList() ; |
9f616d61 | 715 | iDigitInCluster++ ; |
9a1398dd | 716 | digits->Remove(digitN) ; |
d15a28e7 | 717 | break ; |
9f616d61 | 718 | case 2 : // too far from each other |
d15a28e7 | 719 | goto endofloop; |
720 | } // switch | |
721 | ||
722 | } // while digitN | |
fad3e5b9 | 723 | |
d15a28e7 | 724 | endofloop: ; |
725 | nextdigit.Reset() ; | |
726 | ||
727 | } // loop over cluster | |
ad8cfaf4 | 728 | |
ad8cfaf4 | 729 | } // energy theshold |
ad8cfaf4 | 730 | |
31aa6d6c | 731 | |
d15a28e7 | 732 | } // while digit |
733 | ||
9a1398dd | 734 | delete digits ; |
735 | ||
736 | } | |
737 | ||
738 | //____________________________________________________________________________ | |
a4e98857 | 739 | void AliPHOSClusterizerv1::MakeUnfolding() |
740 | { | |
741 | // Unfolds clusters using the shape of an ElectroMagnetic shower | |
9a1398dd | 742 | // Performs unfolding of all EMC/CPV but NOT ppsd clusters |
743 | ||
a4e98857 | 744 | // Unfold first EMC clusters |
9a1398dd | 745 | if(fNumberOfEmcClusters > 0){ |
746 | ||
747 | Int_t nModulesToUnfold = fGeom->GetNModules() ; | |
748 | ||
749 | Int_t numberofNotUnfolded = fNumberOfEmcClusters ; | |
750 | Int_t index ; | |
751 | for(index = 0 ; index < numberofNotUnfolded ; index++){ | |
752 | ||
753 | AliPHOSEmcRecPoint * emcRecPoint = (AliPHOSEmcRecPoint *) fEmcRecPoints->At(index) ; | |
754 | if(emcRecPoint->GetPHOSMod()> nModulesToUnfold) | |
755 | break ; | |
756 | ||
757 | Int_t nMultipl = emcRecPoint->GetMultiplicity() ; | |
758 | Int_t * maxAt = new Int_t[nMultipl] ; | |
759 | Float_t * maxAtEnergy = new Float_t[nMultipl] ; | |
760 | Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fEmcLocMaxCut,fDigits) ; | |
761 | ||
762 | if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0 | |
763 | UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ; | |
764 | fEmcRecPoints->Remove(emcRecPoint); | |
765 | fEmcRecPoints->Compress() ; | |
766 | index-- ; | |
767 | fNumberOfEmcClusters -- ; | |
768 | numberofNotUnfolded-- ; | |
769 | } | |
770 | ||
771 | delete[] maxAt ; | |
772 | delete[] maxAtEnergy ; | |
773 | } | |
774 | } | |
a4e98857 | 775 | // Unfolding of EMC clusters finished |
9a1398dd | 776 | |
777 | ||
a4e98857 | 778 | // Unfold now CPV clusters |
9a1398dd | 779 | if(fNumberOfCpvClusters > 0){ |
780 | ||
781 | Int_t nModulesToUnfold = fGeom->GetNCPVModules() ; | |
782 | ||
783 | Int_t numberofCpvNotUnfolded = fNumberOfCpvClusters ; | |
784 | Int_t index ; | |
785 | for(index = 0 ; index < numberofCpvNotUnfolded ; index++){ | |
786 | ||
787 | AliPHOSRecPoint * recPoint = (AliPHOSRecPoint *) fCpvRecPoints->At(index) ; | |
788 | ||
789 | if(recPoint->GetPHOSMod()> nModulesToUnfold) | |
790 | break ; | |
791 | ||
792 | AliPHOSEmcRecPoint * emcRecPoint = (AliPHOSEmcRecPoint*) recPoint ; | |
793 | ||
794 | Int_t nMultipl = emcRecPoint->GetMultiplicity() ; | |
795 | Int_t * maxAt = new Int_t[nMultipl] ; | |
796 | Float_t * maxAtEnergy = new Float_t[nMultipl] ; | |
797 | Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fCpvLocMaxCut,fDigits) ; | |
798 | ||
799 | if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0 | |
800 | UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ; | |
801 | fCpvRecPoints->Remove(emcRecPoint); | |
802 | fCpvRecPoints->Compress() ; | |
803 | index-- ; | |
804 | numberofCpvNotUnfolded-- ; | |
805 | fNumberOfCpvClusters-- ; | |
806 | } | |
807 | ||
808 | delete[] maxAt ; | |
809 | delete[] maxAtEnergy ; | |
810 | } | |
811 | } | |
812 | //Unfolding of Cpv clusters finished | |
813 | ||
814 | } | |
815 | ||
9a1398dd | 816 | //____________________________________________________________________________ |
817 | Double_t AliPHOSClusterizerv1::ShowerShape(Double_t r) | |
818 | { | |
819 | // Shape of the shower (see PHOS TDR) | |
a4e98857 | 820 | // If you change this function, change also the gradient evaluation in ChiSquare() |
9a1398dd | 821 | |
822 | Double_t r4 = r*r*r*r ; | |
823 | Double_t r295 = TMath::Power(r, 2.95) ; | |
824 | Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ; | |
825 | return shape ; | |
826 | } | |
827 | ||
828 | //____________________________________________________________________________ | |
829 | void AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc, | |
830 | Int_t nMax, | |
831 | int * maxAt, | |
832 | Float_t * maxAtEnergy) | |
833 | { | |
834 | // Performs the unfolding of a cluster with nMax overlapping showers | |
835 | ||
836 | Int_t nPar = 3 * nMax ; | |
837 | Float_t * fitparameters = new Float_t[nPar] ; | |
838 | ||
839 | Bool_t rv = FindFit(iniEmc, maxAt, maxAtEnergy, nPar, fitparameters) ; | |
840 | if( !rv ) { | |
841 | // Fit failed, return and remove cluster | |
842 | delete[] fitparameters ; | |
843 | return ; | |
844 | } | |
845 | ||
846 | // create ufolded rec points and fill them with new energy lists | |
847 | // First calculate energy deposited in each sell in accordance with fit (without fluctuations): efit[] | |
848 | // and later correct this number in acordance with actual energy deposition | |
849 | ||
850 | Int_t nDigits = iniEmc->GetMultiplicity() ; | |
851 | Float_t * efit = new Float_t[nDigits] ; | |
852 | Float_t xDigit,zDigit,distance ; | |
853 | Float_t xpar,zpar,epar ; | |
854 | Int_t relid[4] ; | |
855 | AliPHOSDigit * digit ; | |
856 | Int_t * emcDigits = iniEmc->GetDigitsList() ; | |
857 | ||
858 | Int_t iparam ; | |
859 | Int_t iDigit ; | |
860 | for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){ | |
861 | digit = (AliPHOSDigit*) fDigits->At(emcDigits[iDigit] ) ; | |
862 | fGeom->AbsToRelNumbering(digit->GetId(), relid) ; | |
863 | fGeom->RelPosInModule(relid, xDigit, zDigit) ; | |
864 | efit[iDigit] = 0; | |
865 | ||
866 | iparam = 0 ; | |
867 | while(iparam < nPar ){ | |
868 | xpar = fitparameters[iparam] ; | |
869 | zpar = fitparameters[iparam+1] ; | |
870 | epar = fitparameters[iparam+2] ; | |
871 | iparam += 3 ; | |
872 | distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ; | |
873 | distance = TMath::Sqrt(distance) ; | |
874 | efit[iDigit] += epar * ShowerShape(distance) ; | |
875 | } | |
876 | } | |
877 | ||
878 | ||
879 | // Now create new RecPoints and fill energy lists with efit corrected to fluctuations | |
880 | // so that energy deposited in each cell is distributed betwin new clusters proportionally | |
881 | // to its contribution to efit | |
882 | ||
883 | Float_t * emcEnergies = iniEmc->GetEnergiesList() ; | |
884 | Float_t ratio ; | |
885 | ||
886 | iparam = 0 ; | |
887 | while(iparam < nPar ){ | |
888 | xpar = fitparameters[iparam] ; | |
889 | zpar = fitparameters[iparam+1] ; | |
890 | epar = fitparameters[iparam+2] ; | |
891 | iparam += 3 ; | |
892 | ||
893 | AliPHOSEmcRecPoint * emcRP ; | |
894 | ||
895 | if(iniEmc->IsEmc()){ //create new entries in fEmcRecPoints... | |
896 | ||
897 | if(fNumberOfEmcClusters >= fEmcRecPoints->GetSize()) | |
898 | fEmcRecPoints->Expand(2*fNumberOfEmcClusters) ; | |
899 | ||
900 | (*fEmcRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint() ; | |
901 | emcRP = (AliPHOSEmcRecPoint *) fEmcRecPoints->At(fNumberOfEmcClusters); | |
902 | fNumberOfEmcClusters++ ; | |
903 | } | |
904 | else{//create new entries in fCpvRecPoints | |
905 | if(fNumberOfCpvClusters >= fCpvRecPoints->GetSize()) | |
906 | fCpvRecPoints->Expand(2*fNumberOfCpvClusters) ; | |
907 | ||
908 | (*fCpvRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint() ; | |
909 | emcRP = (AliPHOSEmcRecPoint *) fCpvRecPoints->At(fNumberOfCpvClusters); | |
910 | fNumberOfCpvClusters++ ; | |
911 | } | |
912 | ||
913 | Float_t eDigit ; | |
914 | for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){ | |
915 | digit = (AliPHOSDigit*) fDigits->At( emcDigits[iDigit] ) ; | |
916 | fGeom->AbsToRelNumbering(digit->GetId(), relid) ; | |
917 | fGeom->RelPosInModule(relid, xDigit, zDigit) ; | |
918 | distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ; | |
919 | distance = TMath::Sqrt(distance) ; | |
920 | ratio = epar * ShowerShape(distance) / efit[iDigit] ; | |
921 | eDigit = emcEnergies[iDigit] * ratio ; | |
922 | emcRP->AddDigit( *digit, eDigit ) ; | |
923 | } | |
924 | } | |
925 | ||
926 | delete[] fitparameters ; | |
927 | delete[] efit ; | |
928 | ||
929 | } | |
930 | ||
931 | //_____________________________________________________________________________ | |
932 | void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag) | |
933 | { | |
a4e98857 | 934 | // Calculates the Chi square for the cluster unfolding minimization |
9a1398dd | 935 | // Number of parameters, Gradient, Chi squared, parameters, what to do |
936 | ||
937 | ||
938 | TList * toMinuit = (TList*) gMinuit->GetObjectFit() ; | |
939 | ||
940 | AliPHOSEmcRecPoint * emcRP = (AliPHOSEmcRecPoint*) toMinuit->At(0) ; | |
941 | TClonesArray * digits = (TClonesArray*)toMinuit->At(1) ; | |
942 | ||
943 | ||
944 | ||
945 | // AliPHOSEmcRecPoint * emcRP = (AliPHOSEmcRecPoint *) gMinuit->GetObjectFit() ; // EmcRecPoint to fit | |
946 | ||
947 | Int_t * emcDigits = emcRP->GetDigitsList() ; | |
948 | ||
949 | Int_t nOfDigits = emcRP->GetDigitsMultiplicity() ; | |
950 | ||
951 | Float_t * emcEnergies = emcRP->GetEnergiesList() ; | |
952 | ||
953 | AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ; | |
954 | ||
955 | fret = 0. ; | |
956 | Int_t iparam ; | |
957 | ||
958 | if(iflag == 2) | |
959 | for(iparam = 0 ; iparam < nPar ; iparam++) | |
960 | Grad[iparam] = 0 ; // Will evaluate gradient | |
961 | ||
962 | Double_t efit ; | |
963 | ||
964 | AliPHOSDigit * digit ; | |
965 | Int_t iDigit ; | |
966 | ||
967 | for( iDigit = 0 ; iDigit < nOfDigits ; iDigit++) { | |
968 | ||
969 | digit = (AliPHOSDigit*) digits->At( emcDigits[iDigit] ) ; | |
970 | ||
971 | Int_t relid[4] ; | |
972 | Float_t xDigit ; | |
973 | Float_t zDigit ; | |
974 | ||
975 | geom->AbsToRelNumbering(digit->GetId(), relid) ; | |
976 | ||
977 | geom->RelPosInModule(relid, xDigit, zDigit) ; | |
978 | ||
979 | if(iflag == 2){ // calculate gradient | |
980 | Int_t iParam = 0 ; | |
981 | efit = 0 ; | |
982 | while(iParam < nPar ){ | |
983 | Double_t distance = (xDigit - x[iParam]) * (xDigit - x[iParam]) ; | |
984 | iParam++ ; | |
985 | distance += (zDigit - x[iParam]) * (zDigit - x[iParam]) ; | |
986 | distance = TMath::Sqrt( distance ) ; | |
987 | iParam++ ; | |
988 | efit += x[iParam] * ShowerShape(distance) ; | |
989 | iParam++ ; | |
990 | } | |
991 | Double_t sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E) | |
992 | iParam = 0 ; | |
993 | while(iParam < nPar ){ | |
994 | Double_t xpar = x[iParam] ; | |
995 | Double_t zpar = x[iParam+1] ; | |
996 | Double_t epar = x[iParam+2] ; | |
997 | Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ); | |
998 | Double_t shape = sum * ShowerShape(dr) ; | |
999 | Double_t r4 = dr*dr*dr*dr ; | |
1000 | Double_t r295 = TMath::Power(dr,2.95) ; | |
1001 | Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) + | |
1002 | 0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ; | |
1003 | ||
1004 | Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ; // Derivative over x | |
1005 | iParam++ ; | |
1006 | Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ; // Derivative over z | |
1007 | iParam++ ; | |
1008 | Grad[iParam] += shape ; // Derivative over energy | |
1009 | iParam++ ; | |
1010 | } | |
1011 | } | |
1012 | efit = 0; | |
1013 | iparam = 0 ; | |
1014 | ||
1015 | while(iparam < nPar ){ | |
1016 | Double_t xpar = x[iparam] ; | |
1017 | Double_t zpar = x[iparam+1] ; | |
1018 | Double_t epar = x[iparam+2] ; | |
1019 | iparam += 3 ; | |
1020 | Double_t distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ; | |
1021 | distance = TMath::Sqrt(distance) ; | |
1022 | efit += epar * ShowerShape(distance) ; | |
1023 | } | |
1024 | ||
1025 | fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/emcEnergies[iDigit] ; | |
1026 | // Here we assume, that sigma = sqrt(E) | |
1027 | } | |
83974468 | 1028 | |
d15a28e7 | 1029 | } |
1030 | ||
1031 | //____________________________________________________________________________ | |
9a1398dd | 1032 | void AliPHOSClusterizerv1::Print(Option_t * option)const |
d15a28e7 | 1033 | { |
a4e98857 | 1034 | // Prints the parameters of the clusterizer |
9a1398dd | 1035 | if(fIsInitialized){ |
1036 | ||
1037 | // Print parameters | |
1038 | ||
1039 | cout << "---------------"<< GetName() << " " << GetTitle()<< "-----------" << endl | |
1040 | << "Clusterizing digits from the file: " << fHeaderFileName.Data() << endl | |
1041 | << " Branch: " << fDigitsBranchTitle.Data() << endl | |
1042 | << endl | |
1043 | << " EMC Clustering threshold = " << fEmcClusteringThreshold << endl | |
1044 | << " EMC Local Maximum cut = " << fEmcLocMaxCut << endl | |
1045 | << " EMC Logarothmic weight = " << fW0 << endl | |
1046 | << endl | |
1047 | << " CPV Clustering threshold = " << fCpvClusteringThreshold << endl | |
1048 | << " CPV Local Maximum cut = " << fCpvLocMaxCut << endl | |
1049 | << " CPV Logarothmic weight = " << fW0CPV << endl | |
1050 | << endl | |
1051 | << " PPSD Clustering threshold = " << fPpsdClusteringThreshold << endl; | |
1052 | if(fToUnfold) | |
1053 | cout << " Unfolding on " << endl ; | |
1054 | else | |
1055 | cout << " Unfolding off " << endl ; | |
1056 | ||
1057 | cout << "------------------------------------------------------------------" <<endl ; | |
1058 | } | |
1059 | else | |
1060 | cout << " AliPHOSClusterizerv1 not initialized " << endl ; | |
1061 | } | |
1062 | //____________________________________________________________________________ | |
a4e98857 | 1063 | void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option) |
1064 | { | |
1065 | // Prints list of RecPoints produced at the current pass of AliPHOSClusterizer | |
9a1398dd | 1066 | |
1067 | cout << "AliPHOSClusterizerv1: " << endl ; | |
1068 | cout << " Found "<< fEmcRecPoints->GetEntriesFast() << " EMC Rec Points and " | |
1069 | << fCpvRecPoints->GetEntriesFast() << " CPV RecPoints" << endl ; | |
1070 | ||
1071 | if(strstr(option,"all")) { | |
1072 | cout << "EMC clusters " << endl ; | |
1073 | cout << " Index " | |
1074 | << " Ene(MeV) " | |
1075 | << " Multi " | |
1076 | << " Module " | |
1077 | << " X " | |
1078 | << " Y " | |
1079 | << " Z " | |
1080 | << " Lambda 1 " | |
1081 | << " Lambda 2 " | |
1082 | << " MaxEnergy " | |
1083 | << " # of prim " | |
1084 | << " Primaries list " << endl; | |
1085 | ||
1086 | Int_t index ; | |
1087 | for (index = 0 ; index < fEmcRecPoints->GetEntries() ; index++) { | |
1088 | AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )fEmcRecPoints->At(index) ; | |
1089 | ||
1090 | cout << setw(6) << rp->GetIndexInList() << " "; | |
1091 | cout << setw(6) << rp->GetEnergy() << " "; | |
1092 | cout << setw(6) << rp->GetMultiplicity()<< " "; | |
1093 | cout << setw(6) << rp->GetPHOSMod() << " "; | |
1094 | ||
1095 | TVector3 locpos; | |
1096 | rp->GetLocalPosition(locpos); | |
1097 | cout << setw(8) << locpos.X() << " "; | |
1098 | cout << setw(8) << locpos.Y() << " "; | |
1099 | cout << setw(8) << locpos.Z() << " "; | |
1100 | ||
1101 | Float_t lambda[2]; | |
1102 | rp->GetElipsAxis(lambda); | |
1103 | cout << setw(10)<< lambda[0] << " "; | |
1104 | cout << setw(10)<< lambda[1] << " "; | |
1105 | ||
1106 | ||
1107 | Int_t * primaries; | |
1108 | Int_t nprimaries; | |
1109 | primaries = rp->GetPrimaries(nprimaries); | |
1110 | cout << setw(8) << primaries << " "; | |
1111 | ||
1112 | for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) | |
1113 | cout << setw(4) << primaries[iprimary] << " "; | |
1114 | cout << endl; | |
1115 | } | |
1116 | cout << endl ; | |
1117 | ||
1118 | //Now plot PCV/PPSD recPoints | |
1119 | cout << "EMC clusters " << endl ; | |
1120 | cout << " Index " | |
1121 | << " Multi " | |
1122 | << " Module " | |
1123 | << " Layer " | |
1124 | << " X " | |
1125 | << " Y " | |
1126 | << " Z " | |
1127 | << " # of prim " | |
1128 | << " Primaries list " << endl; | |
1129 | ||
1130 | for (index = 0 ; index < fEmcRecPoints->GetEntries() ; index++) { | |
1131 | AliPHOSRecPoint * rp = (AliPHOSRecPoint * )fCpvRecPoints->At(index) ; | |
1132 | cout << setw(6) << rp->GetIndexInList() << " "; | |
1133 | cout << setw(6) << rp->GetPHOSMod() << " "; | |
1134 | ||
1135 | if( (strcmp(rp->ClassName() , "AliPHOSPpsdRecPoint" )) == 0){ | |
1136 | AliPHOSPpsdRecPoint * ppsd = (AliPHOSPpsdRecPoint*) rp ; | |
1137 | if(ppsd->GetUp()) | |
1138 | cout <<" CPV "; | |
1139 | else | |
1140 | cout <<" PPSD "; | |
1141 | } | |
1142 | else | |
1143 | cout <<" CPV "; | |
1144 | ||
1145 | TVector3 locpos; | |
1146 | rp->GetLocalPosition(locpos); | |
1147 | cout << setw(8) << locpos.X() << " "; | |
1148 | cout << setw(8) << locpos.Y() << " "; | |
1149 | cout << setw(8) << locpos.Z() << " "; | |
1150 | ||
1151 | Int_t * primaries; | |
1152 | Int_t nprimaries; | |
1153 | primaries = rp->GetPrimaries(nprimaries); | |
1154 | cout << setw(8) << primaries << " "; | |
1155 | ||
1156 | for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) | |
1157 | cout << setw(4) << primaries[iprimary] << " "; | |
1158 | cout << endl; | |
1159 | } | |
1160 | ||
1161 | ||
1162 | cout << "-------------------------------------------------"<<endl ; | |
1163 | } | |
d15a28e7 | 1164 | } |
9a1398dd | 1165 |