1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 /* History of cvs commits:
21 * Revision 1.100 2007/01/10 11:07:26 kharlov
22 * Raw digits writing to file (B.Polichtchouk)
24 * Revision 1.99 2006/11/07 16:49:51 kharlov
25 * Corrections for next event switch in case of raw data (B.Polichtchouk)
27 * Revision 1.98 2006/10/27 17:14:27 kharlov
28 * Introduce AliDebug and AliLog (B.Polichtchouk)
30 * Revision 1.97 2006/08/29 11:41:19 kharlov
31 * Missing implementation of ctors and = operator are added
33 * Revision 1.96 2006/08/25 16:56:30 kharlov
34 * Compliance with Effective C++
36 * Revision 1.95 2006/08/11 12:36:26 cvetan
37 * Update of the PHOS code needed in order to read and reconstruct the beam test raw data (i.e. without an existing galice.root)
39 * Revision 1.94 2006/08/07 12:27:49 hristov
40 * Removing obsolete code which affected the event numbering scheme
42 * Revision 1.93 2006/08/01 12:20:17 cvetan
43 * 1. Adding a possibility to read and reconstruct an old rcu formatted raw data. This is controlled by an option of AliReconstruction and AliPHOSReconstructor. 2. In case of raw data processing (without galice.root) create the default AliPHOSGeometry object. Most likely this should be moved to the CDB
45 * Revision 1.92 2006/04/29 20:26:46 hristov
46 * Separate EMC and CPV calibration (Yu.Kharlov)
48 * Revision 1.91 2006/04/22 10:30:17 hristov
49 * Add fEnergy to AliPHOSDigit and operate with EMC amplitude in energy units (Yu.Kharlov)
51 * Revision 1.90 2006/04/11 15:22:59 hristov
52 * run number in query set to -1: forces AliCDBManager to use its run number (A.Colla)
54 * Revision 1.89 2006/03/13 14:05:42 kharlov
55 * Calibration objects for EMC and CPV
57 * Revision 1.88 2006/01/11 08:54:52 hristov
58 * Additional protection in case no calibration entry was found
60 * Revision 1.87 2005/11/22 08:46:43 kharlov
61 * Updated to new CDB (Boris Polichtchouk)
63 * Revision 1.86 2005/11/14 21:52:43 hristov
66 * Revision 1.85 2005/09/27 16:08:08 hristov
67 * New version of CDB storage framework (A.Colla)
69 * Revision 1.84 2005/09/21 10:02:47 kharlov
70 * Reading calibration from CDB (Boris Polichtchouk)
72 * Revision 1.82 2005/09/02 15:43:13 kharlov
73 * Add comments in GetCalibrationParameters and Calibrate
75 * Revision 1.81 2005/09/02 14:32:07 kharlov
76 * Calibration of raw data
78 * Revision 1.80 2005/08/24 15:31:36 kharlov
79 * Setting raw digits flag
81 * Revision 1.79 2005/07/25 15:53:53 kharlov
84 * Revision 1.78 2005/05/28 14:19:04 schutz
85 * Compilation warnings fixed by T.P.
89 //*-- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko (SUBATECH & Kurchatov Institute)
90 //////////////////////////////////////////////////////////////////////////////
91 // Clusterization class. Performs clusterization (collects neighbouring active cells) and
92 // unfolds the clusters having several local maxima.
93 // Results are stored in TreeR#, branches PHOSEmcRP (EMC recPoints),
94 // PHOSCpvRP (CPV RecPoints) and AliPHOSClusterizer (Clusterizer with all
95 // parameters including input digits branch title, thresholds etc.)
96 // This TTask is normally called from Reconstructioner, but can as well be used in
99 // root [0] AliPHOSClusterizerv1 * cl = new AliPHOSClusterizerv1("galice.root", "recpointsname", "digitsname")
100 // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
101 // // reads gAlice from header file "galice.root", uses digits stored in the branch names "digitsname" (default = "Default")
102 // // and saves recpoints in branch named "recpointsname" (default = "digitsname")
103 // root [1] cl->ExecuteTask()
104 // //finds RecPoints in all events stored in galice.root
105 // root [2] cl->SetDigitsBranch("digits2")
106 // //sets another title for Digitis (input) branch
107 // root [3] cl->SetRecPointsBranch("recp2")
108 // //sets another title four output branches
109 // root [4] cl->SetEmcLocalMaxCut(0.03)
110 // //set clusterization parameters
111 // root [5] cl->ExecuteTask("deb all time")
112 // //once more finds RecPoints options are
113 // // deb - print number of found rec points
114 // // deb all - print number of found RecPoints and some their characteristics
115 // // time - print benchmarking results
117 // --- ROOT system ---
122 #include "TBenchmark.h"
124 // --- Standard library ---
126 // --- AliRoot header files ---
128 #include "AliRunLoader.h"
129 #include "AliGenerator.h"
130 #include "AliPHOSGetter.h"
131 #include "AliPHOSGeometry.h"
132 #include "AliPHOSClusterizerv1.h"
133 #include "AliPHOSEmcRecPoint.h"
134 #include "AliPHOSCpvRecPoint.h"
135 #include "AliPHOSDigit.h"
136 #include "AliPHOSDigitizer.h"
137 #include "AliPHOSCalibrationDB.h"
138 #include "AliCDBManager.h"
139 #include "AliCDBStorage.h"
140 #include "AliCDBEntry.h"
142 ClassImp(AliPHOSClusterizerv1)
144 //____________________________________________________________________________
145 AliPHOSClusterizerv1::AliPHOSClusterizerv1() :
146 AliPHOSClusterizer(),
147 fDefaultInit(0), fEmcCrystals(0), fToUnfold(0),
148 fWrite(0), fNumberOfEmcClusters(0), fNumberOfCpvClusters(0),
149 fCalibData(0), fADCchanelEmc(0), fADCpedestalEmc(0),
150 fADCchanelCpv(0), fADCpedestalCpv(0), fEmcClusteringThreshold(0),
151 fCpvClusteringThreshold(0), fEmcMinE(0), fCpvMinE(0),
152 fEmcLocMaxCut(0), fW0(0), fCpvLocMaxCut(0),
153 fW0CPV(0), fRecPointsInRun(0), fEmcTimeGate(0),
156 // default ctor (to be used mainly by Streamer)
159 fDefaultInit = kTRUE ;
162 //____________________________________________________________________________
163 AliPHOSClusterizerv1::AliPHOSClusterizerv1(const TString alirunFileName, const TString eventFolderName) :
164 AliPHOSClusterizer(alirunFileName, eventFolderName),
165 fDefaultInit(0), fEmcCrystals(0), fToUnfold(0),
166 fWrite(0), fNumberOfEmcClusters(0), fNumberOfCpvClusters(0),
167 fCalibData(0), fADCchanelEmc(0), fADCpedestalEmc(0),
168 fADCchanelCpv(0), fADCpedestalCpv(0), fEmcClusteringThreshold(0),
169 fCpvClusteringThreshold(0), fEmcMinE(0), fCpvMinE(0),
170 fEmcLocMaxCut(0), fW0(0), fCpvLocMaxCut(0),
171 fW0CPV(0), fRecPointsInRun(0), fEmcTimeGate(0),
174 // ctor with the indication of the file where header Tree and digits Tree are stored
178 fDefaultInit = kFALSE ;
181 //____________________________________________________________________________
182 AliPHOSClusterizerv1::AliPHOSClusterizerv1(const AliPHOSClusterizerv1 & obj) :
183 AliPHOSClusterizer(obj),
184 fDefaultInit(0), fEmcCrystals(0), fToUnfold(0),
185 fWrite(0), fNumberOfEmcClusters(0), fNumberOfCpvClusters(0),
186 fCalibData(0), fADCchanelEmc(0), fADCpedestalEmc(0),
187 fADCchanelCpv(0), fADCpedestalCpv(0), fEmcClusteringThreshold(0),
188 fCpvClusteringThreshold(0), fEmcMinE(0), fCpvMinE(0),
189 fEmcLocMaxCut(0), fW0(0), fCpvLocMaxCut(0),
190 fW0CPV(0), fRecPointsInRun(0), fEmcTimeGate(0),
195 //____________________________________________________________________________
196 AliPHOSClusterizerv1::~AliPHOSClusterizerv1()
201 //____________________________________________________________________________
202 const TString AliPHOSClusterizerv1::BranchName() const
207 //____________________________________________________________________________
208 Float_t AliPHOSClusterizerv1::CalibrateEMC(Float_t amp, Int_t absId)
210 // Convert EMC measured amplitude into real energy.
211 // Calibration parameters are taken from calibration data base for raw data,
212 // or from digitizer parameters for simulated data.
216 AliPHOSGetter *gime = AliPHOSGetter::Instance();
217 gime->PHOSGeometry()->AbsToRelNumbering(absId,relId) ;
218 Int_t module = relId[0];
219 Int_t column = relId[3];
220 Int_t row = relId[2];
221 if(absId <= fEmcCrystals) { // this is EMC
222 fADCchanelEmc = fCalibData->GetADCchannelEmc (module,column,row);
223 return amp*fADCchanelEmc ;
227 if(absId <= fEmcCrystals) // this is EMC
228 return fADCpedestalEmc + amp*fADCchanelEmc ;
233 //____________________________________________________________________________
234 Float_t AliPHOSClusterizerv1::CalibrateCPV(Int_t amp, Int_t absId)
236 // Convert digitized CPV amplitude into charge.
237 // Calibration parameters are taken from calibration data base for raw data,
238 // or from digitizer parameters for simulated data.
242 AliPHOSGetter *gime = AliPHOSGetter::Instance();
243 gime->PHOSGeometry()->AbsToRelNumbering(absId,relId) ;
244 Int_t module = relId[0];
245 Int_t column = relId[3];
246 Int_t row = relId[2];
247 if(absId > fEmcCrystals) { // this is CPV
248 fADCchanelCpv = fCalibData->GetADCchannelCpv (module,column,row);
249 fADCpedestalCpv = fCalibData->GetADCpedestalCpv(module,column,row);
250 return fADCpedestalCpv + amp*fADCchanelCpv ;
254 if(absId > fEmcCrystals) // this is CPV
255 return fADCpedestalCpv+ amp*fADCchanelCpv ;
260 //____________________________________________________________________________
261 void AliPHOSClusterizerv1::Exec(Option_t *option)
263 // Steering method to perform clusterization for events
264 // in the range from fFirstEvent to fLastEvent.
265 // This range is optionally set by SetEventRange().
266 // if fLastEvent=-1 (by default), then process events until the end.
268 if(strstr(option,"tim"))
269 gBenchmark->Start("PHOSClusterizer");
271 if(strstr(option,"print")) {
276 GetCalibrationParameters() ;
278 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
280 gime->SetRawDigits(kFALSE);
282 gime->SetRawDigits(kTRUE);
284 if (fLastEvent == -1)
285 fLastEvent = gime->MaxEvent() - 1 ;
287 fLastEvent = TMath::Min(fFirstEvent, gime->MaxEvent()); // one event at the time
288 Int_t nEvents = fLastEvent - fFirstEvent + 1;
291 for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
293 gime->Event(ievent ,"D"); // Read digits from simulated data
295 AliRunLoader * rl = AliRunLoader::GetRunLoader(gime->PhosLoader()->GetTitle());
296 rl->GetEvent(ievent);
297 gime->Event(fRawReader,"W",fIsOldRCUFormat); // Read digits from raw data
299 fNumberOfEmcClusters = fNumberOfCpvClusters = 0 ;
303 AliDebug(2,Form(" ---- Printing clusters (%d) of event %d.\n",
304 gime->EmcRecPoints()->GetEntries(),ievent));
305 if(AliLog::GetGlobalDebugLevel()>1)
306 gime->EmcRecPoints()->Print();
313 if(strstr(option,"deb"))
314 PrintRecPoints(option) ;
316 //increment the total number of recpoints per run
317 fRecPointsInRun += gime->EmcRecPoints()->GetEntriesFast() ;
318 fRecPointsInRun += gime->CpvRecPoints()->GetEntriesFast() ;
321 if(fWrite) //do not unload in "on flight" mode
324 if(strstr(option,"tim")){
325 gBenchmark->Stop("PHOSClusterizer");
326 AliInfo(Form("took %f seconds for Clusterizing %f seconds per event \n",
327 gBenchmark->GetCpuTime("PHOSClusterizer"),
328 gBenchmark->GetCpuTime("PHOSClusterizer")/nEvents )) ;
332 //____________________________________________________________________________
333 Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, AliPHOSDigit ** maxAt, Float_t * maxAtEnergy,
334 Int_t nPar, Float_t * fitparameters) const
336 // Calls TMinuit to fit the energy distribution of a cluster with several maxima
337 // The initial values for fitting procedure are set equal to the positions of local maxima.
338 // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
341 AliPHOSGetter * gime = AliPHOSGetter::Instance();
342 TClonesArray * digits = gime->Digits();
344 gMinuit->mncler(); // Reset Minuit's list of paramters
345 gMinuit->SetPrintLevel(-1) ; // No Printout
346 gMinuit->SetFCN(AliPHOSClusterizerv1::UnfoldingChiSquare) ;
347 // To set the address of the minimization function
349 TList * toMinuit = new TList();
350 toMinuit->AddAt(emcRP,0) ;
351 toMinuit->AddAt(digits,1) ;
353 gMinuit->SetObjectFit(toMinuit) ; // To tranfer pointer to UnfoldingChiSquare
355 // filling initial values for fit parameters
356 AliPHOSDigit * digit ;
360 Int_t nDigits = (Int_t) nPar / 3 ;
364 const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
366 for(iDigit = 0; iDigit < nDigits; iDigit++){
367 digit = maxAt[iDigit];
372 geom->AbsToRelNumbering(digit->GetId(), relid) ;
373 geom->RelPosInModule(relid, x, z) ;
375 Float_t energy = maxAtEnergy[iDigit] ;
377 gMinuit->mnparm(index, "x", x, 0.1, 0, 0, ierflg) ;
380 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : x = %f\n", x ) ;
383 gMinuit->mnparm(index, "z", z, 0.1, 0, 0, ierflg) ;
386 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : z =%f\n", z ) ;
389 gMinuit->mnparm(index, "Energy", energy , 0.05*energy, 0., 4.*energy, ierflg) ;
392 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : energy = %f\n", energy ) ;
397 Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
402 gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ; // force TMinuit to reduce function calls
403 gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient
404 gMinuit->SetMaxIterations(5);
405 gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ; // No Warnings
407 gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize
409 if(ierflg == 4){ // Minimum not found
410 Warning("FindFit", "PHOS Unfolding fit not converged, cluster abandoned\n" );
413 for(index = 0; index < nPar; index++){
416 gMinuit->GetParameter(index, val, err) ; // Returns value and error of parameter index
417 fitparameters[index] = val ;
425 //____________________________________________________________________________
426 void AliPHOSClusterizerv1::GetCalibrationParameters()
428 // Set calibration parameters:
429 // if calibration database exists, they are read from database,
430 // otherwise, they are taken from digitizer.
432 // It is a user responsilibity to open CDB before reconstruction, for example:
433 // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
435 AliPHOSGetter * gime = AliPHOSGetter::Instance();
436 // fCalibData = new AliPHOSCalibData(gAlice->GetRunNumber()); //original
438 fCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number
442 if ( !gime->Digitizer() )
443 gime->LoadDigitizer();
444 AliPHOSDigitizer * dig = gime->Digitizer();
445 fADCchanelEmc = dig->GetEMCchannel() ;
446 fADCpedestalEmc = dig->GetEMCpedestal();
448 fADCchanelCpv = dig->GetCPVchannel() ;
449 fADCpedestalCpv = dig->GetCPVpedestal() ;
453 //____________________________________________________________________________
454 void AliPHOSClusterizerv1::Init()
456 // Make all memory allocations which can not be done in default constructor.
457 // Attach the Clusterizer task to the list of PHOS tasks
459 AliPHOSGetter* gime = AliPHOSGetter::Instance() ;
461 gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName.Data());
463 AliPHOSGeometry * geom = gime->PHOSGeometry();
465 fEmcCrystals = geom->GetNModules() * geom->GetNCristalsInModule() ;
468 gMinuit = new TMinuit(100);
470 if ( !gime->Clusterizer() ) {
471 gime->PostClusterizer(this);
475 //____________________________________________________________________________
476 void AliPHOSClusterizerv1::InitParameters()
479 fNumberOfCpvClusters = 0 ;
480 fNumberOfEmcClusters = 0 ;
482 fCpvClusteringThreshold = 0.0;
483 fEmcClusteringThreshold = 0.2;
485 fEmcLocMaxCut = 0.03 ;
486 fCpvLocMaxCut = 0.03 ;
494 fEmcTimeGate = 1.e-8 ;
498 fRecPointsInRun = 0 ;
504 SetEventRange(0,-1) ;
506 fIsOldRCUFormat = kFALSE;
509 //____________________________________________________________________________
510 Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)const
512 // Gives the neighbourness of two digits = 0 are not neighbour but continue searching
514 // = 2 are not neighbour but do not continue searching
515 // neighbours are defined as digits having at least a common vertex
516 // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
517 // which is compared to a digit (d2) not yet in a cluster
519 AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ;
524 geom->AbsToRelNumbering(d1->GetId(), relid1) ;
527 geom->AbsToRelNumbering(d2->GetId(), relid2) ;
529 if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module
530 Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
531 Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
533 if (( coldiff <= 1 ) && ( rowdiff <= 1 )){
534 if((relid1[1] != 0) || (TMath::Abs(d1->GetTime() - d2->GetTime() ) < fEmcTimeGate))
538 if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1))
539 rv = 2; // Difference in row numbers is too large to look further
545 if( (relid1[0] < relid2[0]) || (relid1[1] != relid2[1]) )
552 //____________________________________________________________________________
553 void AliPHOSClusterizerv1::CleanDigits(TClonesArray * digits)
555 // Remove digits with amplitudes below threshold
557 for(Int_t i=0; i<digits->GetEntriesFast(); i++){
558 AliPHOSDigit * digit = static_cast<AliPHOSDigit*>(digits->At(i)) ;
559 if ( (IsInEmc(digit) && CalibrateEMC(digit->GetEnergy(),digit->GetId()) < fEmcMinE) ||
560 (IsInCpv(digit) && CalibrateCPV(digit->GetAmp() ,digit->GetId()) < fCpvMinE) )
561 digits->RemoveAt(i) ;
564 for (Int_t i = 0 ; i < digits->GetEntriesFast() ; i++) {
565 AliPHOSDigit *digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
566 digit->SetIndexInList(i) ;
569 //Overwrite digits tree
570 AliPHOSGetter* gime = AliPHOSGetter::Instance();
571 TTree * treeD = gime->TreeD();
572 treeD->Branch("PHOS", &digits);
574 gime->WriteDigits("OVERWRITE");
575 gime->PhosLoader()->UnloadDigits() ;
577 //____________________________________________________________________________
578 Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit) const
580 // Tells if (true) or not (false) the digit is in a PHOS-EMC module
583 AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ;
585 Int_t nEMC = geom->GetNModules()*geom->GetNPhi()*geom->GetNZ();
587 if(digit->GetId() <= nEMC ) rv = kTRUE;
592 //____________________________________________________________________________
593 Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit) const
595 // Tells if (true) or not (false) the digit is in a PHOS-CPV module
599 AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ;
601 Int_t nEMC = geom->GetNModules()*geom->GetNPhi()*geom->GetNZ();
603 if(digit->GetId() > nEMC ) rv = kTRUE;
608 //____________________________________________________________________________
609 void AliPHOSClusterizerv1::Unload()
611 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
612 gime->PhosLoader()->UnloadDigits() ;
613 gime->PhosLoader()->UnloadRecPoints() ;
616 //____________________________________________________________________________
617 void AliPHOSClusterizerv1::WriteRecPoints()
620 // Creates new branches with given title
621 // fills and writes into TreeR.
623 AliPHOSGetter * gime = AliPHOSGetter::Instance();
625 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
626 TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
627 TClonesArray * digits = gime->Digits() ;
630 //Evaluate position, dispersion and other RecPoint properties..
631 Int_t nEmc = emcRecPoints->GetEntriesFast();
632 for(index = 0; index < nEmc; index++){
633 AliPHOSEmcRecPoint * rp = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) );
634 rp->Purify(fEmcMinE) ;
635 if(rp->GetMultiplicity()==0){
636 emcRecPoints->RemoveAt(index) ;
640 // No vertex is available now, calculate only vertex-independent quantities and
641 // postpone cluster properties calculation to TrackSegmentMaker
642 rp->EvalAll(fW0,digits) ;
644 emcRecPoints->Compress() ;
645 // emcRecPoints->Sort() ; //Can not sort until position is calculated!
646 // emcRecPoints->Expand(emcRecPoints->GetEntriesFast()) ;
647 for(index = 0; index < emcRecPoints->GetEntries(); index++){
648 dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) )->SetIndexInList(index) ;
651 //Now the same for CPV
652 for(index = 0; index < cpvRecPoints->GetEntries(); index++){
653 AliPHOSCpvRecPoint * rp = dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(index) );
654 rp->EvalAll(fW0CPV,digits) ;
656 // cpvRecPoints->Sort() ;
658 for(index = 0; index < cpvRecPoints->GetEntries(); index++)
659 dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(index) )->SetIndexInList(index) ;
661 cpvRecPoints->Expand(cpvRecPoints->GetEntriesFast()) ;
663 if(fWrite){ //We write TreeR
664 TTree * treeR = gime->TreeR();
666 Int_t bufferSize = 32000 ;
667 Int_t splitlevel = 0 ;
670 TBranch * emcBranch = treeR->Branch("PHOSEmcRP","TObjArray",&emcRecPoints,bufferSize,splitlevel);
671 emcBranch->SetTitle(BranchName());
674 TBranch * cpvBranch = treeR->Branch("PHOSCpvRP","TObjArray",&cpvRecPoints,bufferSize,splitlevel);
675 cpvBranch->SetTitle(BranchName());
680 gime->WriteRecPoints("OVERWRITE");
681 gime->WriteClusterizer("OVERWRITE");
685 //____________________________________________________________________________
686 void AliPHOSClusterizerv1::MakeClusters()
688 // Steering method to construct the clusters stored in a list of Reconstructed Points
689 // A cluster is defined as a list of neighbour digits
692 AliPHOSGetter * gime = AliPHOSGetter::Instance();
694 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
695 TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
697 emcRecPoints->Delete() ;
698 cpvRecPoints->Delete() ;
700 TClonesArray * digits = gime->Digits() ;
702 //Remove digits below threshold
703 CleanDigits(digits) ;
705 TClonesArray * digitsC = static_cast<TClonesArray*>( digits->Clone() ) ;
707 // Clusterization starts
709 TIter nextdigit(digitsC) ;
710 AliPHOSDigit * digit ;
711 Bool_t notremoved = kTRUE ;
713 while ( (digit = dynamic_cast<AliPHOSDigit *>( nextdigit()) ) ) { // scan over the list of digitsC
716 AliPHOSRecPoint * clu = 0 ;
718 TArrayI clusterdigitslist(1500) ;
721 if (( IsInEmc (digit) &&
722 CalibrateEMC(digit->GetEnergy(),digit->GetId()) > fEmcClusteringThreshold ) ||
724 CalibrateCPV(digit->GetAmp() ,digit->GetId()) > fCpvClusteringThreshold ) ) {
725 Int_t iDigitInCluster = 0 ;
727 if ( IsInEmc(digit) ) {
728 // start a new EMC RecPoint
729 if(fNumberOfEmcClusters >= emcRecPoints->GetSize())
730 emcRecPoints->Expand(2*fNumberOfEmcClusters+1) ;
732 emcRecPoints->AddAt(new AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ;
733 clu = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(fNumberOfEmcClusters) ) ;
734 fNumberOfEmcClusters++ ;
735 clu->AddDigit(*digit, CalibrateEMC(digit->GetEnergy(),digit->GetId())) ;
736 clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
738 digitsC->Remove(digit) ;
742 // start a new CPV cluster
743 if(fNumberOfCpvClusters >= cpvRecPoints->GetSize())
744 cpvRecPoints->Expand(2*fNumberOfCpvClusters+1);
746 cpvRecPoints->AddAt(new AliPHOSCpvRecPoint(""), fNumberOfCpvClusters) ;
748 clu = dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(fNumberOfCpvClusters) ) ;
749 fNumberOfCpvClusters++ ;
750 clu->AddDigit(*digit, CalibrateCPV(digit->GetAmp(),digit->GetId()) ) ;
751 clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
753 digitsC->Remove(digit) ;
756 // Here we remove remaining EMC digits, which cannot make a cluster
759 while( ( digit = dynamic_cast<AliPHOSDigit *>( nextdigit() ) ) ) {
761 digitsC->Remove(digit) ;
765 notremoved = kFALSE ;
772 AliPHOSDigit * digitN ;
774 while (index < iDigitInCluster){ // scan over digits already in cluster
775 digit = dynamic_cast<AliPHOSDigit*>( digits->At(clusterdigitslist[index]) ) ;
777 while ( (digitN = dynamic_cast<AliPHOSDigit *>( nextdigit() ) ) ) { // scan over the reduced list of digits
778 Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!!
780 case 0 : // not a neighbour
782 case 1 : // are neighbours
783 if (IsInEmc (digitN))
784 clu->AddDigit(*digitN, CalibrateEMC( digitN->GetEnergy(), digitN->GetId() ) );
786 clu->AddDigit(*digitN, CalibrateCPV( digitN->GetAmp() , digitN->GetId() ) );
788 clusterdigitslist[iDigitInCluster] = digitN->GetIndexInList() ;
790 digitsC->Remove(digitN) ;
792 case 2 : // too far from each other
801 } // loop over cluster
812 //____________________________________________________________________________
813 void AliPHOSClusterizerv1::MakeUnfolding()
815 // Unfolds clusters using the shape of an ElectroMagnetic shower
816 // Performs unfolding of all EMC/CPV clusters
818 AliPHOSGetter * gime = AliPHOSGetter::Instance();
820 const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
822 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
823 TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
824 TClonesArray * digits = gime->Digits() ;
826 // Unfold first EMC clusters
827 if(fNumberOfEmcClusters > 0){
829 Int_t nModulesToUnfold = geom->GetNModules() ;
831 Int_t numberofNotUnfolded = fNumberOfEmcClusters ;
833 for(index = 0 ; index < numberofNotUnfolded ; index++){
835 AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) ) ;
836 if(emcRecPoint->GetPHOSMod()> nModulesToUnfold)
839 Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
840 AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
841 Float_t * maxAtEnergy = new Float_t[nMultipl] ;
842 Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fEmcLocMaxCut,digits) ;
844 if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
845 UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
846 emcRecPoints->Remove(emcRecPoint);
847 emcRecPoints->Compress() ;
849 fNumberOfEmcClusters -- ;
850 numberofNotUnfolded-- ;
853 emcRecPoint->SetNExMax(1) ; //Only one local maximum
857 delete[] maxAtEnergy ;
860 // Unfolding of EMC clusters finished
863 // Unfold now CPV clusters
864 if(fNumberOfCpvClusters > 0){
866 Int_t nModulesToUnfold = geom->GetNModules() ;
868 Int_t numberofCpvNotUnfolded = fNumberOfCpvClusters ;
870 for(index = 0 ; index < numberofCpvNotUnfolded ; index++){
872 AliPHOSRecPoint * recPoint = dynamic_cast<AliPHOSRecPoint *>( cpvRecPoints->At(index) ) ;
874 if(recPoint->GetPHOSMod()> nModulesToUnfold)
877 AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint*>(recPoint) ;
879 Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
880 AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
881 Float_t * maxAtEnergy = new Float_t[nMultipl] ;
882 Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fCpvLocMaxCut,digits) ;
884 if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
885 UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
886 cpvRecPoints->Remove(emcRecPoint);
887 cpvRecPoints->Compress() ;
889 numberofCpvNotUnfolded-- ;
890 fNumberOfCpvClusters-- ;
894 delete[] maxAtEnergy ;
897 //Unfolding of Cpv clusters finished
901 //____________________________________________________________________________
902 Double_t AliPHOSClusterizerv1::ShowerShape(Double_t x, Double_t z)
904 // Shape of the shower (see PHOS TDR)
905 // If you change this function, change also the gradient evaluation in ChiSquare()
907 //for the moment we neglect dependence on the incident angle.
909 Double_t r2 = x*x + z*z ;
910 Double_t r4 = r2*r2 ;
911 Double_t r295 = TMath::Power(r2, 2.95/2.) ;
912 Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
916 //____________________________________________________________________________
917 void AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,
919 AliPHOSDigit ** maxAt,
920 Float_t * maxAtEnergy)
922 // Performs the unfolding of a cluster with nMax overlapping showers
924 AliPHOSGetter * gime = AliPHOSGetter::Instance();
926 const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
928 const TClonesArray * digits = gime->Digits() ;
929 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
930 TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
932 Int_t nPar = 3 * nMax ;
933 Float_t * fitparameters = new Float_t[nPar] ;
935 Bool_t rv = FindFit(iniEmc, maxAt, maxAtEnergy, nPar, fitparameters) ;
937 // Fit failed, return and remove cluster
938 iniEmc->SetNExMax(-1) ;
939 delete[] fitparameters ;
943 // create ufolded rec points and fill them with new energy lists
944 // First calculate energy deposited in each sell in accordance with fit (without fluctuations): efit[]
945 // and later correct this number in acordance with actual energy deposition
947 Int_t nDigits = iniEmc->GetMultiplicity() ;
948 Float_t * efit = new Float_t[nDigits] ;
949 Float_t xDigit=0.,zDigit=0. ;
950 Float_t xpar=0.,zpar=0.,epar=0. ;
952 AliPHOSDigit * digit = 0 ;
953 Int_t * emcDigits = iniEmc->GetDigitsList() ;
959 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
960 digit = dynamic_cast<AliPHOSDigit*>( digits->At(emcDigits[iDigit] ) ) ;
961 geom->AbsToRelNumbering(digit->GetId(), relid) ;
962 geom->RelPosInModule(relid, xDigit, zDigit) ;
966 while(iparam < nPar ){
967 xpar = fitparameters[iparam] ;
968 zpar = fitparameters[iparam+1] ;
969 epar = fitparameters[iparam+2] ;
971 // geom->GetIncidentVector(fVtx,relid[0],xpar,zpar,vIncid) ;
972 // efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) ;
973 efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
978 // Now create new RecPoints and fill energy lists with efit corrected to fluctuations
979 // so that energy deposited in each cell is distributed betwin new clusters proportionally
980 // to its contribution to efit
982 Float_t * emcEnergies = iniEmc->GetEnergiesList() ;
986 while(iparam < nPar ){
987 xpar = fitparameters[iparam] ;
988 zpar = fitparameters[iparam+1] ;
989 epar = fitparameters[iparam+2] ;
991 // geom->GetIncidentVector(fVtx,iniEmc->GetPHOSMod(),xpar,zpar,vIncid) ;
993 AliPHOSEmcRecPoint * emcRP = 0 ;
995 if(iniEmc->IsEmc()){ //create new entries in fEmcRecPoints...
997 if(fNumberOfEmcClusters >= emcRecPoints->GetSize())
998 emcRecPoints->Expand(2*fNumberOfEmcClusters) ;
1000 (*emcRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint("") ;
1001 emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(fNumberOfEmcClusters) ) ;
1002 fNumberOfEmcClusters++ ;
1003 emcRP->SetNExMax((Int_t)nPar/3) ;
1005 else{//create new entries in fCpvRecPoints
1006 if(fNumberOfCpvClusters >= cpvRecPoints->GetSize())
1007 cpvRecPoints->Expand(2*fNumberOfCpvClusters) ;
1009 (*cpvRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint("") ;
1010 emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( cpvRecPoints->At(fNumberOfCpvClusters) ) ;
1011 fNumberOfCpvClusters++ ;
1015 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
1016 digit = dynamic_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) ) ;
1017 geom->AbsToRelNumbering(digit->GetId(), relid) ;
1018 geom->RelPosInModule(relid, xDigit, zDigit) ;
1019 // ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) / efit[iDigit] ;
1020 ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar) / efit[iDigit] ;
1021 eDigit = emcEnergies[iDigit] * ratio ;
1022 emcRP->AddDigit( *digit, eDigit ) ;
1026 delete[] fitparameters ;
1031 //_____________________________________________________________________________
1032 void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
1034 // Calculates the Chi square for the cluster unfolding minimization
1035 // Number of parameters, Gradient, Chi squared, parameters, what to do
1037 TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;
1039 AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint*>( toMinuit->At(0) ) ;
1040 TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) ) ;
1041 // TVector3 * vtx = dynamic_cast<TVector3*>(toMinuit->At(2)) ; //Vertex position
1043 // AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( gMinuit->GetObjectFit() ) ; // EmcRecPoint to fit
1045 Int_t * emcDigits = emcRP->GetDigitsList() ;
1047 Int_t nOdigits = emcRP->GetDigitsMultiplicity() ;
1049 Float_t * emcEnergies = emcRP->GetEnergiesList() ;
1051 const AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ;
1057 for(iparam = 0 ; iparam < nPar ; iparam++)
1058 Grad[iparam] = 0 ; // Will evaluate gradient
1062 AliPHOSDigit * digit ;
1065 for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
1067 digit = dynamic_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) );
1073 geom->AbsToRelNumbering(digit->GetId(), relid) ;
1075 geom->RelPosInModule(relid, xDigit, zDigit) ;
1077 if(iflag == 2){ // calculate gradient
1080 while(iParam < nPar ){
1081 Double_t dx = (xDigit - x[iParam]) ;
1083 Double_t dz = (zDigit - x[iParam]) ;
1085 // geom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),x[iParam-2],x[iParam-1],vInc) ;
1086 // efit += x[iParam] * ShowerShape(dx,dz,vInc) ;
1087 efit += x[iParam] * ShowerShape(dx,dz) ;
1090 Double_t sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E)
1092 while(iParam < nPar ){
1093 Double_t xpar = x[iParam] ;
1094 Double_t zpar = x[iParam+1] ;
1095 Double_t epar = x[iParam+2] ;
1096 // geom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
1097 Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
1098 // Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
1099 Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar) ;
1100 //DP: No incident angle dependence in gradient yet!!!!!!
1101 Double_t r4 = dr*dr*dr*dr ;
1102 Double_t r295 = TMath::Power(dr,2.95) ;
1103 Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) +
1104 0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ;
1106 Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ; // Derivative over x
1108 Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ; // Derivative over z
1110 Grad[iParam] += shape ; // Derivative over energy
1117 while(iparam < nPar ){
1118 Double_t xpar = x[iparam] ;
1119 Double_t zpar = x[iparam+1] ;
1120 Double_t epar = x[iparam+2] ;
1122 // geom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
1123 // efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
1124 efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
1127 fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/emcEnergies[iDigit] ;
1128 // Here we assume, that sigma = sqrt(E)
1133 //____________________________________________________________________________
1134 void AliPHOSClusterizerv1::Print(const Option_t *)const
1136 // Print clusterizer parameters
1139 TString taskName(GetName()) ;
1140 taskName.ReplaceAll(Version(), "") ;
1142 if( strcmp(GetName(), "") !=0 ) {
1144 message = "\n--------------- %s %s -----------\n" ;
1145 message += "Clusterizing digits from the file: %s\n" ;
1146 message += " Branch: %s\n" ;
1147 message += " EMC Clustering threshold = %f\n" ;
1148 message += " EMC Local Maximum cut = %f\n" ;
1149 message += " EMC Logarothmic weight = %f\n" ;
1150 message += " CPV Clustering threshold = %f\n" ;
1151 message += " CPV Local Maximum cut = %f\n" ;
1152 message += " CPV Logarothmic weight = %f\n" ;
1154 message += " Unfolding on\n" ;
1156 message += " Unfolding off\n" ;
1158 message += "------------------------------------------------------------------" ;
1161 message = " AliPHOSClusterizerv1 not initialized " ;
1163 AliInfo(Form("%s, %s %s %s %s %s %s %s %s %s %s", message.Data(),
1168 fEmcClusteringThreshold,
1171 fCpvClusteringThreshold,
1175 //____________________________________________________________________________
1176 //void AliPHOSClusterizerv1::GetVertex(void)
1177 //{ //Extracts vertex posisition
1185 // if(gAlice && gAlice->GetMCApp() && gAlice->Generator()){
1187 // gAlice->Generator()->GetOrigin(x,y,z) ;
1188 // fVtx.SetXYZ(x,y,z) ;
1193 // fVtx[0]=fVtx[1]=fVtx[2]=0. ;
1196 //____________________________________________________________________________
1197 void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option)
1199 // Prints list of RecPoints produced at the current pass of AliPHOSClusterizer
1201 AliPHOSGetter * gime = AliPHOSGetter::Instance();
1203 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
1204 TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
1206 AliInfo(Form("\nevent %d \n Found %d EMC RecPoints and %d CPV RecPoints",
1207 gAlice->GetEvNumber(),
1208 emcRecPoints->GetEntriesFast(),
1209 cpvRecPoints->GetEntriesFast() )) ;
1211 fRecPointsInRun += emcRecPoints->GetEntriesFast() ;
1212 fRecPointsInRun += cpvRecPoints->GetEntriesFast() ;
1215 if(strstr(option,"all")) {
1216 printf("\n EMC clusters \n") ;
1217 printf("Index Ene(MeV) Multi Module X Y Z Lambdas_1 Lambda_2 # of prim Primaries list\n") ;
1219 for (index = 0 ; index < emcRecPoints->GetEntries() ; index++) {
1220 AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )emcRecPoints->At(index) ;
1222 rp->GetLocalPosition(locpos);
1224 rp->GetElipsAxis(lambda);
1227 primaries = rp->GetPrimaries(nprimaries);
1228 printf("\n%6d %8.2f %3d %2d %4.1f %4.1f %4.1f %4f %4f %2d : ",
1229 rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetPHOSMod(),
1230 locpos.X(), locpos.Y(), locpos.Z(), lambda[0], lambda[1], nprimaries) ;
1232 for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
1233 printf("%d ", primaries[iprimary] ) ;
1238 //Now plot CPV recPoints
1239 printf("\n CPV clusters \n") ;
1240 printf("Index Ene(MeV) Module X Y Z \n") ;
1241 for (index = 0 ; index < cpvRecPoints->GetEntries() ; index++) {
1242 AliPHOSCpvRecPoint * rp = (AliPHOSCpvRecPoint * )cpvRecPoints->At(index) ;
1245 rp->GetLocalPosition(locpos);
1247 printf("\n%6d %8.2f %2d %4.1f %4.1f %4.1f \n",
1248 rp->GetIndexInList(), rp->GetEnergy(), rp->GetPHOSMod(),
1249 locpos.X(), locpos.Y(), locpos.Z()) ;