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 1 October 2000. Yuri Kharlov:
20 PPSD upper layer is considered if number of layers>1
21 18 October 2000. Yuri Kharlov:
22 AliEMCALClusterizerv1()
23 CPV clusterizing parameters added
25 After first PPSD digit remove EMC digits only once
27 //*-- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko (SUBATECH & Kurchatov Institute)
28 // August 2002 Yves Schutz: clone PHOS as closely as possible and intoduction
29 // of new IO (à la PHOS)
30 //////////////////////////////////////////////////////////////////////////////
31 // Clusterization class. Performs clusterization (collects neighbouring active cells) and
32 // unfolds the clusters having several local maxima.
33 // Results are stored in TreeR#, branches EMCALTowerRP (EMC recPoints),
34 // EMCALPreShoRP (CPV RecPoints) and AliEMCALClusterizer (Clusterizer with all
35 // parameters including input digits branch title, thresholds etc.)
36 // This TTask is normally called from Reconstructioner, but can as well be used in
39 // root [0] AliEMCALClusterizerv1 * cl = new AliEMCALClusterizerv1("galice.root")
40 // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
41 // //reads gAlice from header file "..."
42 // root [1] cl->ExecuteTask()
43 // //finds RecPoints in all events stored in galice.root
44 // root [2] cl->SetDigitsBranch("digits2")
45 // //sets another title for Digitis (input) branch
46 // root [3] cl->SetRecPointsBranch("recp2")
47 // //sets another title four output branches
48 // root [4] cl->SetTowerLocalMaxCut(0.03)
49 // //set clusterization parameters
50 // root [5] cl->ExecuteTask("deb all time")
51 // //once more finds RecPoints options are
52 // // deb - print number of found rec points
53 // // deb all - print number of found RecPoints and some their characteristics
54 // // time - print benchmarking results
56 // --- ROOT system ---
65 #include "TBenchmark.h"
67 // --- Standard library ---
70 // --- AliRoot header files ---
72 #include "AliEMCALClusterizerv1.h"
73 #include "AliEMCALDigit.h"
74 #include "AliEMCALDigitizer.h"
75 #include "AliEMCALTowerRecPoint.h"
77 #include "AliEMCALGetter.h"
78 #include "AliEMCALGeometry.h"
81 ClassImp(AliEMCALClusterizerv1)
83 //____________________________________________________________________________
84 AliEMCALClusterizerv1::AliEMCALClusterizerv1() : AliEMCALClusterizer()
86 // default ctor (to be used mainly by Streamer)
89 fDefaultInit = kTRUE ;
92 //____________________________________________________________________________
93 AliEMCALClusterizerv1::AliEMCALClusterizerv1(const char* headerFile, const char* name, const Bool_t toSplit)
94 :AliEMCALClusterizer(headerFile, name, toSplit)
96 // ctor with the indication of the file where header Tree and digits Tree are stored
100 fDefaultInit = kFALSE ;
104 //____________________________________________________________________________
105 AliEMCALClusterizerv1::~AliEMCALClusterizerv1()
112 //____________________________________________________________________________
113 const TString AliEMCALClusterizerv1::BranchName() const
115 TString branchName(GetName() ) ;
116 branchName.Remove(branchName.Index(Version())-1) ;
120 //____________________________________________________________________________
121 Float_t AliEMCALClusterizerv1::Calibrate(Int_t amp, Int_t where) const
123 //To be replased later by the method, reading individual parameters from the database
124 // where = 0 == PRE ; where = 1 == ECAL ; where = 2 == HCAL
125 if ( where == 0 ) // calibrate as PRE section
126 return -fADCpedestalPRE + amp * fADCchannelPRE ;
127 else if (where == 1) //calibrate as EC section
128 return -fADCpedestalEC + amp * fADCchannelEC ;
129 else if (where == 2) //calibrate as HC section
130 return -fADCpedestalHC + amp * fADCchannelHC ;
132 Fatal("Calibrate", "Something went wrong!") ;
136 //____________________________________________________________________________
137 void AliEMCALClusterizerv1::Exec(Option_t * option)
141 if( strcmp(GetName(), "")== 0 )
144 if(strstr(option,"tim"))
145 gBenchmark->Start("EMCALClusterizer");
147 if(strstr(option,"print"))
150 AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
151 if(gime->BranchExists("RecPoints"))
153 Int_t nevents = gime->MaxEvent() ;
156 for(ievent = 0; ievent < nevents; ievent++){
158 gime->Event(ievent,"D") ;
161 GetCalibrationParameters() ;
163 fNumberOfPREClusters = fNumberOfECClusters = fNumberOfHCClusters = 0 ;
170 WriteRecPoints(ievent) ;
172 if(strstr(option,"deb"))
173 PrintRecPoints(option) ;
175 //increment the total number of digits per run
176 fRecPointsInRun += gime->PRERecPoints()->GetEntriesFast() ;
177 fRecPointsInRun += gime->ECALRecPoints()->GetEntriesFast() ;
178 fRecPointsInRun += gime->HCALRecPoints()->GetEntriesFast() ;
181 if(strstr(option,"tim")){
182 gBenchmark->Stop("EMCALClusterizer");
183 Info("Exec", "took %f seconds for Clusterizing %f seconds per event",
184 gBenchmark->GetCpuTime("EMCALClusterizer"), gBenchmark->GetCpuTime("EMCALClusterizer")/nevents ) ;
189 //____________________________________________________________________________
190 Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALTowerRecPoint * emcRP, AliEMCALDigit ** maxAt, Float_t * maxAtEnergy,
191 Int_t nPar, Float_t * fitparameters) const
193 // Calls TMinuit to fit the energy distribution of a cluster with several maxima
194 // The initial values for fitting procedure are set equal to the positions of local maxima.
195 // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
197 AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
198 TClonesArray * digits = gime->Digits() ;
201 gMinuit->mncler(); // Reset Minuit's list of paramters
202 gMinuit->SetPrintLevel(-1) ; // No Printout
203 gMinuit->SetFCN(AliEMCALClusterizerv1::UnfoldingChiSquare) ;
204 // To set the address of the minimization function
205 TList * toMinuit = new TList();
206 toMinuit->AddAt(emcRP,0) ;
207 toMinuit->AddAt(digits,1) ;
209 gMinuit->SetObjectFit(toMinuit) ; // To tranfer pointer to UnfoldingChiSquare
211 // filling initial values for fit parameters
212 AliEMCALDigit * digit ;
216 Int_t nDigits = (Int_t) nPar / 3 ;
220 AliEMCALGeometry * geom = gime->EMCALGeometry() ;
222 for(iDigit = 0; iDigit < nDigits; iDigit++){
223 digit = maxAt[iDigit];
228 geom->AbsToRelNumbering(digit->GetId(), relid) ;
229 geom->PosInAlice(relid, x, z) ;
231 Float_t energy = maxAtEnergy[iDigit] ;
233 gMinuit->mnparm(index, "x", x, 0.1, 0, 0, ierflg) ;
236 Error("FindFit", "EMCAL Unfolding Unable to set initial value for fit procedure : x = %f", x ) ;
239 gMinuit->mnparm(index, "z", z, 0.1, 0, 0, ierflg) ;
242 Error("FindFit", "EMCAL Unfolding Unable to set initial value for fit procedure : z = %f", z) ;
245 gMinuit->mnparm(index, "Energy", energy , 0.05*energy, 0., 4.*energy, ierflg) ;
248 Error("FindFit", "EMCAL Unfolding Unable to set initial value for fit procedure : energy = %f", energy) ;
253 Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
258 gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ; // force TMinuit to reduce function calls
259 gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient
260 gMinuit->SetMaxIterations(5);
261 gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ; // No Warnings
263 gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize
265 if(ierflg == 4){ // Minimum not found
266 Error("FindFit", "EMCAL Unfolding Fit not converged, cluster abandoned " ) ;
269 for(index = 0; index < nPar; index++){
272 gMinuit->GetParameter(index, val, err) ; // Returns value and error of parameter index
273 fitparameters[index] = val ;
281 //____________________________________________________________________________
282 void AliEMCALClusterizerv1::GetCalibrationParameters()
284 AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
285 const AliEMCALDigitizer * dig = gime->Digitizer(BranchName()) ;
287 fADCchannelPRE = dig->GetPREchannel() ;
288 fADCpedestalPRE = dig->GetPREpedestal() ;
290 fADCchannelEC = dig->GetECchannel() ;
291 fADCpedestalEC = dig->GetECpedestal();
293 fADCchannelHC = dig->GetHCchannel() ;
294 fADCpedestalHC = dig->GetHCpedestal();
297 //____________________________________________________________________________
298 void AliEMCALClusterizerv1::Init()
300 // Make all memory allocations which can not be done in default constructor.
301 // Attach the Clusterizer task to the list of EMCAL tasks
303 if ( strcmp(GetTitle(), "") == 0 )
304 SetTitle("galice.root") ;
306 TString branchname = GetName() ;
307 branchname.Remove(branchname.Index(Version())-1) ;
309 AliEMCALGetter * gime = AliEMCALGetter::GetInstance(GetTitle(), branchname.Data(), fToSplit ) ;
311 Error("Init", "Could not obtain the Getter object !" ) ;
317 // construct the name of the file as /path/EMCAL.SDigits.root
318 //First - extract full path if necessary
319 TString fileName(GetTitle()) ;
320 Ssiz_t islash = fileName.Last('/') ;
321 if(islash<fileName.Length())
322 fileName.Remove(islash+1,fileName.Length()) ;
325 // Next - append the file name
326 fileName+="EMCAL.RecData." ;
327 if((strcmp(branchname.Data(),"Default")!=0)&&(strcmp(branchname.Data(),"")!=0)){
328 fileName+=branchname ;
332 // Finally - check if the file already opened or open the file
333 fSplitFile = static_cast<TFile*>(gROOT->GetFile(fileName.Data()));
335 fSplitFile = TFile::Open(fileName.Data(),"update") ;
338 const AliEMCALGeometry * geom = gime->EMCALGeometry() ;
339 fNTowers = geom->GetNZ() * geom->GetNPhi() ;
341 gMinuit = new TMinuit(100) ;
343 gime->PostClusterizer(this) ;
344 gime->PostRecPoints(branchname ) ;
348 //____________________________________________________________________________
349 void AliEMCALClusterizerv1::InitParameters()
351 fNumberOfPREClusters = fNumberOfECClusters = fNumberOfHCClusters = 0 ;
352 fPREClusteringThreshold = 0.0001; // must be adjusted according to the noise leve set by digitizer
353 fECClusteringThreshold = 0.0045; // must be adjusted according to the noise leve set by digitizer
354 fHCClusteringThreshold = 0.001; // must be adjusted according to the noise leve set by digitizer
355 fPRELocMaxCut = 0.03 ;
356 fECLocMaxCut = 0.03 ;
357 fHCLocMaxCut = 0.03 ;
367 TString clusterizerName( GetName()) ;
368 if (clusterizerName.IsNull() )
369 clusterizerName = "Default" ;
370 clusterizerName.Append(":") ;
371 clusterizerName.Append(Version()) ;
372 SetName(clusterizerName) ;
373 fRecPointsInRun = 0 ;
377 //____________________________________________________________________________
378 Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2)const
380 // Gives the neighbourness of two digits = 0 are not neighbour but continue searching
382 // = 2 are not neighbour but do not continue searching
383 // neighbours are defined as digits having at least a common vertex
384 // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
385 // which is compared to a digit (d2) not yet in a cluster
387 AliEMCALGeometry * geom = AliEMCALGetter::GetInstance()->EMCALGeometry() ;
392 geom->AbsToRelNumbering(d1->GetId(), relid1) ;
395 geom->AbsToRelNumbering(d2->GetId(), relid2) ;
397 if ( (relid1[0] == relid2[0]) && // inside the same EMCAL Arm
398 (relid1[1]==relid2[1]) ) { // and same tower section
399 Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
400 Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
402 if (( coldiff <= 1 ) && ( rowdiff <= 1 )){
403 if((relid1[1] != 0) || (TMath::Abs(d1->GetTime() - d2->GetTime() ) < fTimeGate))
407 if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1))
408 rv = 2; // Difference in row numbers is too large to look further
414 if( (relid1[0] < relid2[0]) || (relid1[1] != relid2[1]) )
419 Info("AreNeighbours", "neighbours=%d, id1=%d, relid1=%d,%d,%d,%d \n id2=%d, relid2=%d,%d,%d,%d",
420 rv, d1->GetId(), relid1[0], relid1[1], relid1[2], relid1[3],
421 d2->GetId(), relid2[0], relid2[1], relid2[2], relid2[3]) ;
427 //____________________________________________________________________________
428 void AliEMCALClusterizerv1::WriteRecPoints(Int_t event)
431 // Creates new branches with given title
432 // fills and writes into TreeR.
434 AliEMCALGetter *gime = AliEMCALGetter::GetInstance() ;
436 TObjArray * aPRERecPoints = gime->PRERecPoints() ;
437 TObjArray * aECRecPoints = gime->ECALRecPoints() ;
438 TObjArray * aHCRecPoints = gime->HCALRecPoints() ;
440 TClonesArray * digits = gime->Digits() ;
447 TString name("TreeR") ;
449 treeR = dynamic_cast<TTree*>(fSplitFile->Get(name));
452 treeR = gAlice->TreeR();
456 gAlice->MakeTree("R", fSplitFile);
457 treeR = gAlice->TreeR() ;
462 //Evaluate position, dispersion and other RecPoint properties for PRE section
463 for(index = 0; index < aPRERecPoints->GetEntries(); index++)
464 (dynamic_cast<AliEMCALRecPoint *>(aPRERecPoints->At(index)))->EvalAll(fPREW0,digits) ;
465 aPRERecPoints->Sort() ;
467 for(index = 0; index < aPRERecPoints->GetEntries(); index++)
468 (dynamic_cast<AliEMCALRecPoint *>(aPRERecPoints->At(index)))->SetIndexInList(index) ;
470 aPRERecPoints->Expand(aPRERecPoints->GetEntriesFast()) ;
472 //Evaluate position, dispersion and other RecPoint properties for EC section
473 for(index = 0; index < aECRecPoints->GetEntries(); index++)
474 (dynamic_cast<AliEMCALTowerRecPoint *>(aECRecPoints->At(index)))->EvalAll(fECW0,digits) ;
476 aECRecPoints->Sort() ;
478 for(index = 0; index < aECRecPoints->GetEntries(); index++)
479 (dynamic_cast<AliEMCALTowerRecPoint *>(aECRecPoints->At(index)))->SetIndexInList(index) ;
481 aECRecPoints->Expand(aECRecPoints->GetEntriesFast()) ;
483 //Evaluate position, dispersion and other RecPoint properties for HC section
484 for(index = 0; index < aHCRecPoints->GetEntries(); index++)
485 (dynamic_cast<AliEMCALTowerRecPoint *>(aHCRecPoints->At(index)))->EvalAll(fHCW0,digits) ;
487 aHCRecPoints->Sort() ;
489 for(index = 0; index < aHCRecPoints->GetEntries(); index++)
490 (dynamic_cast<AliEMCALTowerRecPoint *>(aHCRecPoints->At(index)))->SetIndexInList(index) ;
492 aHCRecPoints->Expand(aHCRecPoints->GetEntriesFast()) ;
494 Int_t bufferSize = 32000 ;
495 Int_t splitlevel = 0 ;
498 TBranch * branchPRE = treeR->Branch("EMCALPRERP","TObjArray",&aPRERecPoints,bufferSize,splitlevel);
499 branchPRE->SetTitle(BranchName());
502 TBranch * branchEC = treeR->Branch("EMCALECRP","TObjArray",&aECRecPoints,bufferSize,splitlevel);
503 branchEC->SetTitle(BranchName());
506 TBranch * branchHC = treeR->Branch("EMCALHCRP","TObjArray",&aHCRecPoints,bufferSize,splitlevel);
507 branchHC->SetTitle(BranchName());
509 //And Finally clusterizer branch
510 AliEMCALClusterizerv1 * cl = (AliEMCALClusterizerv1*)gime->Clusterizer(BranchName()) ;
511 TBranch * clusterizerBranch = treeR->Branch("AliEMCALClusterizer","AliEMCALClusterizerv1",
512 &cl,bufferSize,splitlevel);
513 clusterizerBranch->SetTitle(BranchName());
518 clusterizerBranch->Fill() ;
520 treeR->AutoSave() ; //Write(0,kOverwrite) ;
521 if(gAlice->TreeR()!=treeR)
525 //____________________________________________________________________________
526 void AliEMCALClusterizerv1::MakeClusters()
528 // Steering method to construct the clusters stored in a list of Reconstructed Points
529 // A cluster is defined as a list of neighbour digits
531 AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
533 AliEMCALGeometry * geom = gime->EMCALGeometry() ;
536 TObjArray * aPRERecPoints = gime->PRERecPoints(BranchName()) ;
537 TObjArray * aECRecPoints = gime->ECALRecPoints(BranchName()) ;
538 TObjArray * aHCRecPoints = gime->HCALRecPoints(BranchName()) ;
540 aPRERecPoints->Delete() ;
541 aECRecPoints->Delete() ;
542 aHCRecPoints->Delete() ;
544 TClonesArray * digits = gime->Digits() ;
546 Fatal("MakeClusters -> Digits with name %s not found", GetName() ) ;
550 AliEMCALDigit * digit ;
551 Int_t ndigEC=0, ndigPRE=0, ndigHC=0 ;
553 // count the number of digits in EC section
554 while ( (digit = dynamic_cast<AliEMCALDigit *>(next())) ) { // scan over the list of digits
555 if (geom->IsInECAL(digit->GetId()))
557 else if (geom->IsInPRE(digit->GetId()))
559 else if (geom->IsInHCAL(digit->GetId()))
562 Error("MakeClusters", "id = %d is a wrong ID!", digit->GetId()) ;
567 // add amplitude of PRE and EC sections
569 for (digEC = 0 ; digEC < ndigEC ; digEC++) {
570 digit = dynamic_cast<AliEMCALDigit *>(digits->At(digEC)) ;
572 for (digPRE = ndigEC ; digPRE < ndigEC+ndigPRE ; digPRE++) {
573 AliEMCALDigit * digitPRE = dynamic_cast<AliEMCALDigit *>(digits->At(digPRE)) ;
574 if ( geom->AreInSameTower(digit->GetId(), digitPRE->GetId()) ){
575 Float_t amp = static_cast<Float_t>(digit->GetAmp()) + geom->GetSummationFraction() * static_cast<Float_t>(digitPRE->GetAmp()) + 0.5 ;
576 digit->SetAmp(static_cast<Int_t>(amp)) ;
581 Info("MakeClusters","id = %d amp = %d", digit->GetId(), digit->GetAmp()) ;
584 TClonesArray * digitsC = dynamic_cast<TClonesArray*>(digits->Clone()) ;
587 // Clusterization starts
589 TIter nextdigit(digitsC) ;
590 Bool_t notremovedEC = kTRUE, notremovedPRE = kTRUE ;
592 while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigit())) ) { // scan over the list of digitsC
593 AliEMCALRecPoint * clu = 0 ;
595 TArrayI clusterPREdigitslist(50), clusterECdigitslist(50), clusterHCdigitslist(50);
597 Bool_t inPRE = kFALSE, inECAL = kFALSE, inHCAL = kFALSE ;
598 if( geom->IsInPRE(digit->GetId()) ) {
601 else if( geom->IsInECAL(digit->GetId()) ) {
604 else if( geom->IsInHCAL(digit->GetId()) ) {
610 Info("MakeClusters","id = %d, ene = %f , thre = %f ",
611 digit->GetId(),Calibrate(digit->GetAmp(), 0), fPREClusteringThreshold) ;
613 Info("MakeClusters","id = %d, ene = %f , thre = %f",
614 digit->GetId(),Calibrate(digit->GetAmp(), 1), fECClusteringThreshold) ;
616 Info("MakeClusters","id = %d, ene = %f , thre = %f",
617 digit->GetId(),Calibrate(digit->GetAmp(), 2), fHCClusteringThreshold ) ;
620 if ( (inPRE && (Calibrate(digit->GetAmp(), 0) > fPREClusteringThreshold )) ||
621 (inECAL && (Calibrate(digit->GetAmp(), 1) > fECClusteringThreshold )) ||
622 (inHCAL && (Calibrate(digit->GetAmp(), 2) > fHCClusteringThreshold )) ) {
624 Int_t iDigitInPRECluster = 0, iDigitInECCluster = 0, iDigitInHCCluster = 0;
625 Int_t where ; // PRE = 0, ECAl = 1, HCAL = 2
627 // Find the seed in each of the section ECAL/PRE/HCAL
629 if( geom->IsInECAL(digit->GetId()) ) {
630 where = 1 ; // to tell we are in ECAL
631 // start a new Tower RecPoint
632 if(fNumberOfECClusters >= aECRecPoints->GetSize())
633 aECRecPoints->Expand(2*fNumberOfECClusters+1) ;
634 AliEMCALTowerRecPoint * rp = new AliEMCALTowerRecPoint("") ;
636 aECRecPoints->AddAt(rp, fNumberOfECClusters) ;
637 clu = dynamic_cast<AliEMCALTowerRecPoint *>(aECRecPoints->At(fNumberOfECClusters)) ;
638 fNumberOfECClusters++ ;
639 clu->AddDigit(*digit, Calibrate(digit->GetAmp(), where)) ;
640 clusterECdigitslist[iDigitInECCluster] = digit->GetIndexInList() ;
641 iDigitInECCluster++ ;
642 digitsC->Remove(digit) ;
644 Info("MakeClusters","OK id = %d, ene = %f , thre = %f ", digit->GetId(),Calibrate(digit->GetAmp(), 1), fECClusteringThreshold) ;
647 else if( geom->IsInPRE(digit->GetId()) ) {
648 where = 0 ; // to tell we are in PRE
649 // start a new Pre Shower cluster
650 if(fNumberOfPREClusters >= aPRERecPoints->GetSize())
651 aPRERecPoints->Expand(2*fNumberOfPREClusters+1);
652 AliEMCALTowerRecPoint * rp = new AliEMCALTowerRecPoint("") ;
654 aPRERecPoints->AddAt(rp, fNumberOfPREClusters) ;
655 clu = dynamic_cast<AliEMCALTowerRecPoint *>(aPRERecPoints->At(fNumberOfPREClusters)) ;
656 fNumberOfPREClusters++ ;
657 clu->AddDigit(*digit, Calibrate(digit->GetAmp(), where));
658 clusterPREdigitslist[iDigitInPRECluster] = digit->GetIndexInList() ;
659 iDigitInPRECluster++ ;
660 digitsC->Remove(digit) ;
662 Info("MakeClusters","OK id = %d, ene = %f , thre = %f ", digit->GetId(),Calibrate(digit->GetAmp(), 0), fPREClusteringThreshold) ;
666 // Here we remove remaining EC digits, which cannot make a cluster
669 while( ( digit = dynamic_cast<AliEMCALDigit *>(nextdigit()) ) ) {
670 if( geom->IsInECAL(digit->GetId()) )
671 digitsC->Remove(digit) ;
675 notremovedEC = kFALSE ;
679 else if( geom->IsInHCAL(digit->GetId()) ) {
680 where = 2 ; // to tell we are in HCAL
681 // start a new HCAL cluster
682 if(fNumberOfHCClusters >= aHCRecPoints->GetSize())
683 aHCRecPoints->Expand(2*fNumberOfHCClusters+1);
684 AliEMCALTowerRecPoint * rp = new AliEMCALTowerRecPoint("") ;
686 aHCRecPoints->AddAt(rp, fNumberOfHCClusters) ;
687 clu = dynamic_cast<AliEMCALTowerRecPoint *>(aHCRecPoints->At(fNumberOfHCClusters)) ;
688 fNumberOfHCClusters++ ;
689 clu->AddDigit(*digit, Calibrate(digit->GetAmp(), where));
690 clusterHCdigitslist[iDigitInHCCluster] = digit->GetIndexInList() ;
691 iDigitInHCCluster++ ;
692 digitsC->Remove(digit) ;
694 Info("MakeClusters","OK id = %d, ene = %f , thre = %f ", digit->GetId(),Calibrate(digit->GetAmp(), 2), fHCClusteringThreshold) ;
698 // Here we remove remaining PRE digits, which cannot make a cluster
700 if( notremovedPRE ) {
701 while( ( digit = dynamic_cast<AliEMCALDigit *>(nextdigit()) ) ) {
702 if( geom->IsInPRE(digit->GetId()) )
703 digitsC->Remove(digit) ;
707 notremovedPRE = kFALSE ;
713 AliEMCALDigit * digitN ;
716 // Do the Clustering in each of the three section ECAL/PRE/HCAL
718 while (index < iDigitInECCluster){ // scan over digits already in cluster
719 digit = (AliEMCALDigit*)digits->At(clusterECdigitslist[index]) ;
721 while ( (digitN = (AliEMCALDigit *)nextdigit()) ) { // scan over the reduced list of digits
722 Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!!
723 // Info("MakeClusters","id1 = %d, id2 = %d , neighbours = %d", digit->GetId(), digitN->GetId(), ineb) ;
725 case 0 : // not a neighbour
727 case 1 : // are neighbours
728 clu->AddDigit(*digitN, Calibrate( digitN->GetAmp(), 1) ) ;
729 clusterECdigitslist[iDigitInECCluster] = digitN->GetIndexInList() ;
730 iDigitInECCluster++ ;
731 digitsC->Remove(digitN) ;
733 case 2 : // too far from each other
741 } // loop over EC cluster
744 while (index < iDigitInPRECluster){ // scan over digits already in cluster
745 digit = (AliEMCALDigit*)digits->At(clusterPREdigitslist[index]) ;
747 while ( (digitN = (AliEMCALDigit *)nextdigit()) ) { // scan over the reduced list of digits
748 Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!!
749 // Info("MakeClusters","id1 = %d, id2 = %d , neighbours = %d", digit->GetId(), digitN->GetId(), ineb) ;
751 case 0 : // not a neighbour
753 case 1 : // are neighbours
754 clu->AddDigit(*digitN, Calibrate( digitN->GetAmp(), 0) ) ;
755 clusterPREdigitslist[iDigitInPRECluster] = digitN->GetIndexInList() ;
756 iDigitInPRECluster++ ;
757 digitsC->Remove(digitN) ;
759 case 2 : // too far from each other
767 } // loop over PRE cluster
770 while (index < iDigitInHCCluster){ // scan over digits already in cluster
771 digit = (AliEMCALDigit*)digits->At(clusterHCdigitslist[index]) ;
773 while ( (digitN = (AliEMCALDigit *)nextdigit()) ) { // scan over the reduced list of digits
774 Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!!
775 //Info("MakeClusters","id1 = %d, id2 = %d , neighbours = %d", digit->GetId(), digitN->GetId(), ineb) ;
777 case 0 : // not a neighbour
779 case 1 : // are neighbours
780 clu->AddDigit(*digitN, Calibrate( digitN->GetAmp(), 2) ) ;
781 clusterHCdigitslist[iDigitInHCCluster] = digitN->GetIndexInList() ;
782 iDigitInHCCluster++ ;
783 digitsC->Remove(digitN) ;
785 case 2 : // too far from each other
792 } // loop over HC cluster
799 //____________________________________________________________________________
800 void AliEMCALClusterizerv1::MakeUnfolding()
802 Fatal("AliEMCALClusterizerv1::MakeUnfolding", "--> Unfolding not implemented") ;
806 //____________________________________________________________________________
807 Double_t AliEMCALClusterizerv1::ShowerShape(Double_t r)
809 // Shape of the shower (see EMCAL TDR)
810 // If you change this function, change also the gradient evaluation in ChiSquare()
812 Double_t r4 = r*r*r*r ;
813 Double_t r295 = TMath::Power(r, 2.95) ;
814 Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
818 //____________________________________________________________________________
819 void AliEMCALClusterizerv1::UnfoldCluster(AliEMCALTowerRecPoint * iniTower,
821 AliEMCALDigit ** maxAt,
822 Float_t * maxAtEnergy)
824 // Performs the unfolding of a cluster with nMax overlapping showers
826 Fatal("UnfoldCluster", "--> Unfolding not implemented") ;
830 //_____________________________________________________________________________
831 void AliEMCALClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
833 // Calculates the Chi square for the cluster unfolding minimization
834 // Number of parameters, Gradient, Chi squared, parameters, what to do
836 ::Fatal("UnfoldingChiSquare","Unfolding not implemented") ;
838 //____________________________________________________________________________
839 void AliEMCALClusterizerv1::Print(Option_t * option)const
841 // Print clusterizer parameters
843 TString message("\n") ;
845 if( strcmp(GetName(), "") !=0 ){
849 TString taskName(GetName()) ;
850 taskName.ReplaceAll(Version(), "") ;
852 message += "--------------- " ;
853 message += taskName.Data() ;
855 message += GetTitle() ;
856 message += "-----------\n" ;
857 message += "Clusterizing digits from the file: " ;
858 message += taskName.Data() ;
859 message += "\n Branch: " ;
860 message += GetName() ;
861 message += "\n Pre Shower Clustering threshold = " ;
862 message += fPREClusteringThreshold ;
863 message += "\n Pre Shower Local Maximum cut = " ;
864 message += fPRELocMaxCut ;
865 message += "\n Pre Shower Logarothmic weight = " ;
867 message += "\n EC Clustering threshold = " ;
868 message += fECClusteringThreshold ;
869 message += "\n EC Local Maximum cut = " ;
870 message += fECLocMaxCut ;
871 message += "\n EC Logarothmic weight = " ;
873 message += "\n Pre Shower Clustering threshold = " ;
874 message += fHCClusteringThreshold ;
875 message += "\n HC Local Maximum cut = " ;
876 message += fHCLocMaxCut ;
877 message += "\n HC Logarothmic weight = " ;
880 message +="\nUnfolding on\n" ;
882 message += "\nUnfolding off\n";
884 message += "------------------------------------------------------------------" ;
887 message += "AliEMCALClusterizerv1 not initialized " ;
889 Info("Print", message.Data() ) ;
892 //____________________________________________________________________________
893 void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option)
895 // Prints list of RecPoints produced at the current pass of AliEMCALClusterizer
897 TObjArray * aPRERecPoints = AliEMCALGetter::GetInstance()->PRERecPoints() ;
898 TObjArray * aECRecPoints = AliEMCALGetter::GetInstance()->ECALRecPoints() ;
899 TObjArray * aHCRecPoints = AliEMCALGetter::GetInstance()->HCALRecPoints() ;
901 Info("PrintRecPoints", "Clusterization result:") ;
903 printf("event # %d\n", gAlice->GetEvNumber() ) ;
904 printf(" Found %d PRE SHOWER RecPoints, %d EC Rec Points and %d HC Rec Points\n ",
905 aPRERecPoints->GetEntriesFast(), aECRecPoints->GetEntriesFast(), aHCRecPoints->GetEntriesFast() ) ;
907 fRecPointsInRun += aPRERecPoints->GetEntriesFast() ;
908 fRecPointsInRun += aECRecPoints->GetEntriesFast() ;
909 fRecPointsInRun += aHCRecPoints->GetEntriesFast() ;
911 if(strstr(option,"all")) {
913 //Pre shower recPoints
915 printf("-----------------------------------------------------------------------\n") ;
916 printf("Clusters in PRE section\n") ;
917 printf("Index Ene(GeV) Multi Module phi r theta X Y Z Dispersion Lambda 1 Lambda 2 # of prim Primaries list\n") ;
921 for (index = 0 ; index < aPRERecPoints->GetEntries() ; index++) {
922 AliEMCALTowerRecPoint * rp = dynamic_cast<AliEMCALTowerRecPoint *>(aPRERecPoints->At(index)) ;
924 rp->GetGlobalPosition(globalpos);
926 rp->GetLocalPosition(localpos);
928 rp->GetElipsAxis(lambda);
931 primaries = rp->GetPrimaries(nprimaries);
932 printf("\n%6d %8.4f %3d %2d %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4f %4f %2d : ",
933 rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetEMCALArm(),
934 globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(),
935 rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ;
936 for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
937 printf("%d ", primaries[iprimary] ) ;
941 printf("\n-----------------------------------------------------------------------\n") ;
942 printf("Clusters in ECAL section\n") ;
943 printf("Index Ene(GeV) Multi Module phi r theta X Y Z Dispersion Lambda 1 Lambda 2 # of prim Primaries list\n") ;
945 for (index = 0 ; index < aECRecPoints->GetEntries() ; index++) {
946 AliEMCALTowerRecPoint * rp = dynamic_cast<AliEMCALTowerRecPoint * >(aECRecPoints->At(index)) ;
948 rp->GetGlobalPosition(globalpos);
950 rp->GetLocalPosition(localpos);
952 rp->GetElipsAxis(lambda);
955 primaries = rp->GetPrimaries(nprimaries);
956 printf("\n%6d %8.4f %3d %2d %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4f %4f %2d : ",
957 rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetEMCALArm(),
958 globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(),
959 rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ;
960 for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
961 printf("%d ", primaries[iprimary] ) ;
965 printf("\n-----------------------------------------------------------------------\n") ;
966 printf("Clusters in HCAL section\n") ;
967 printf("Index Ene(GeV) Multi Module phi r theta X Y Z Dispersion Lambda 1 Lambda 2 # of prim Primaries list\n") ;
969 for (index = 0 ; index < aHCRecPoints->GetEntries() ; index++) {
970 AliEMCALTowerRecPoint * rp = dynamic_cast<AliEMCALTowerRecPoint * >(aHCRecPoints->At(index)) ;
972 rp->GetGlobalPosition(globalpos);
974 rp->GetLocalPosition(localpos);
976 rp->GetElipsAxis(lambda);
979 primaries = rp->GetPrimaries(nprimaries);
980 printf("\n%6d %8.4f %3d %2d %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4f %4f %2d : ",
981 rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetEMCALArm(),
982 globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(),
983 rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ;
984 for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
985 printf("%d ", primaries[iprimary] ) ;
989 printf("\n-----------------------------------------------------------------------\n");