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 ---
71 #include "AliEMCALGetter.h"
72 #include "AliEMCALClusterizerv1.h"
73 #include "AliEMCALTowerRecPoint.h"
74 #include "AliEMCALDigit.h"
75 #include "AliEMCALDigitizer.h"
77 #include "AliEMCALGeometry.h"
79 ClassImp(AliEMCALClusterizerv1)
81 //____________________________________________________________________________
82 AliEMCALClusterizerv1::AliEMCALClusterizerv1() : AliEMCALClusterizer()
84 // default ctor (to be used mainly by Streamer)
87 fDefaultInit = kTRUE ;
90 //____________________________________________________________________________
91 AliEMCALClusterizerv1::AliEMCALClusterizerv1(const TString alirunFileName, const TString eventFolderName)
92 :AliEMCALClusterizer(alirunFileName, eventFolderName)
94 // ctor with the indication of the file where header Tree and digits Tree are stored
98 fDefaultInit = kFALSE ;
102 //____________________________________________________________________________
103 AliEMCALClusterizerv1::~AliEMCALClusterizerv1()
109 //____________________________________________________________________________
110 const TString AliEMCALClusterizerv1::BranchName() const
116 //____________________________________________________________________________
117 Float_t AliEMCALClusterizerv1::Calibrate(Int_t amp, Int_t where) const
119 //To be replased later by the method, reading individual parameters from the database
120 // where = 0 == PRE ; where = 1 == ECAL ; where = 2 == HCAL
121 if ( where == 0 ) // calibrate as PRE section
122 return -fADCpedestalPRE + amp * fADCchannelPRE ;
123 else if (where == 1) //calibrate as ECA section
124 return -fADCpedestalECA + amp * fADCchannelECA ;
125 else if (where == 2) //calibrate as HCA section
126 return -fADCpedestalHCA + amp * fADCchannelHCA ;
128 Fatal("Calibrate", "Something went wrong!") ;
132 //____________________________________________________________________________
133 void AliEMCALClusterizerv1::Exec(Option_t * option)
137 if(strstr(option,"tim"))
138 gBenchmark->Start("EMCALClusterizer");
140 if(strstr(option,"print"))
143 AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
145 Int_t nevents = gime->MaxEvent() ;
148 for(ievent = 0; ievent < nevents; ievent++){
150 gime->Event(ievent,"D") ;
152 GetCalibrationParameters() ;
154 fNumberOfPREClusters = fNumberOfECAClusters = fNumberOfHCAClusters = 0 ;
161 WriteRecPoints(ievent) ;
163 if(strstr(option,"deb"))
164 PrintRecPoints(option) ;
166 //increment the total number of recpoints per run
167 fRecPointsInRun += gime->PRERecPoints()->GetEntriesFast() ;
168 fRecPointsInRun += gime->ECARecPoints()->GetEntriesFast() ;
169 fRecPointsInRun += gime->HCARecPoints()->GetEntriesFast() ;
174 if(strstr(option,"tim")){
175 gBenchmark->Stop("EMCALClusterizer");
176 Info("Exec", "took %f seconds for Clusterizing %f seconds per event",
177 gBenchmark->GetCpuTime("EMCALClusterizer"), gBenchmark->GetCpuTime("EMCALClusterizer")/nevents ) ;
181 //____________________________________________________________________________
182 Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALTowerRecPoint * emcRP, AliEMCALDigit ** maxAt, Float_t * maxAtEnergy,
183 Int_t nPar, Float_t * fitparameters) const
185 // Calls TMinuit to fit the energy distribution of a cluster with several maxima
186 // The initial values for fitting procedure are set equal to the positions of local maxima.
187 // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
189 AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
190 TClonesArray * digits = gime->Digits() ;
193 gMinuit->mncler(); // Reset Minuit's list of paramters
194 gMinuit->SetPrintLevel(-1) ; // No Printout
195 gMinuit->SetFCN(AliEMCALClusterizerv1::UnfoldingChiSquare) ;
196 // To set the address of the minimization function
197 TList * toMinuit = new TList();
198 toMinuit->AddAt(emcRP,0) ;
199 toMinuit->AddAt(digits,1) ;
201 gMinuit->SetObjectFit(toMinuit) ; // To tranfer pointer to UnfoldingChiSquare
203 // filling initial values for fit parameters
204 AliEMCALDigit * digit ;
208 Int_t nDigits = (Int_t) nPar / 3 ;
212 AliEMCALGeometry * geom = gime->EMCALGeometry() ;
214 for(iDigit = 0; iDigit < nDigits; iDigit++){
215 digit = maxAt[iDigit];
220 geom->AbsToRelNumbering(digit->GetId(), relid) ;
221 geom->PosInAlice(relid, x, z) ;
223 Float_t energy = maxAtEnergy[iDigit] ;
225 gMinuit->mnparm(index, "x", x, 0.1, 0, 0, ierflg) ;
228 Error("FindFit", "EMCAL Unfolding Unable to set initial value for fit procedure : x = %f", x ) ;
231 gMinuit->mnparm(index, "z", z, 0.1, 0, 0, ierflg) ;
234 Error("FindFit", "EMCAL Unfolding Unable to set initial value for fit procedure : z = %f", z) ;
237 gMinuit->mnparm(index, "Energy", energy , 0.05*energy, 0., 4.*energy, ierflg) ;
240 Error("FindFit", "EMCAL Unfolding Unable to set initial value for fit procedure : energy = %f", energy) ;
245 Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
250 gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ; // force TMinuit to reduce function calls
251 gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient
252 gMinuit->SetMaxIterations(5);
253 gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ; // No Warnings
255 gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize
257 if(ierflg == 4){ // Minimum not found
258 Error("FindFit", "EMCAL Unfolding Fit not converged, cluster abandoned " ) ;
261 for(index = 0; index < nPar; index++){
264 gMinuit->GetParameter(index, val, err) ; // Returns value and error of parameter index
265 fitparameters[index] = val ;
273 //____________________________________________________________________________
274 void AliEMCALClusterizerv1::GetCalibrationParameters()
276 AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
278 if ( !gime->Digitizer() )
279 gime->LoadDigitizer();
280 AliEMCALDigitizer * dig = gime->Digitizer();
282 fADCchannelPRE = dig->GetPREchannel() ;
283 fADCpedestalPRE = dig->GetPREpedestal() ;
285 fADCchannelECA = dig->GetECAchannel() ;
286 fADCpedestalECA = dig->GetECApedestal();
288 fADCchannelHCA = dig->GetHCAchannel() ;
289 fADCpedestalHCA = dig->GetHCApedestal();
292 //____________________________________________________________________________
293 void AliEMCALClusterizerv1::Init()
295 // Make all memory allocations which can not be done in default constructor.
296 // Attach the Clusterizer task to the list of EMCAL tasks
298 AliEMCALGetter * gime = AliEMCALGetter::Instance(GetTitle(), fEventFolderName.Data());
300 AliEMCALGeometry * geom = gime->EMCALGeometry() ;
302 fNTowers = geom->GetNZ() * geom->GetNPhi() ;
304 gMinuit = new TMinuit(100) ;
306 if ( !gime->Clusterizer() )
307 gime->PostClusterizer(this);
310 //____________________________________________________________________________
311 void AliEMCALClusterizerv1::InitParameters()
313 fNumberOfPREClusters = fNumberOfECAClusters = fNumberOfHCAClusters = 0 ;
314 fPREClusteringThreshold = 0.0001; // must be adjusted according to the noise leve set by digitizer
315 fECAClusteringThreshold = 0.0045; // must be adjusted according to the noise leve set by digitizer
316 fHCAClusteringThreshold = 0.001; // must be adjusted according to the noise leve set by digitizer
317 fPRELocMaxCut = 0.03 ;
318 fECALocMaxCut = 0.03 ;
319 fHCALocMaxCut = 0.03 ;
329 fRecPointsInRun = 0 ;
333 //____________________________________________________________________________
334 Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2)const
336 // Gives the neighbourness of two digits = 0 are not neighbour but continue searching
338 // = 2 are not neighbour but do not continue searching
339 // neighbours are defined as digits having at least a common vertex
340 // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
341 // which is compared to a digit (d2) not yet in a cluster
343 AliEMCALGeometry * geom = AliEMCALGetter::Instance()->EMCALGeometry() ;
348 geom->AbsToRelNumbering(d1->GetId(), relid1) ;
351 geom->AbsToRelNumbering(d2->GetId(), relid2) ;
353 if ( (relid1[0] == relid2[0]) && // inside the same EMCAL Arm
354 (relid1[1]==relid2[1]) ) { // and same tower section
355 Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
356 Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
358 if (( coldiff <= 1 ) && ( rowdiff <= 1 )){
359 if((relid1[1] != 0) || (TMath::Abs(d1->GetTime() - d2->GetTime() ) < fTimeGate))
363 if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1))
364 rv = 2; // Difference in row numbers is too large to look further
370 if( (relid1[0] < relid2[0]) || (relid1[1] != relid2[1]) )
375 Info("AreNeighbours", "neighbours=%d, id1=%d, relid1=%d,%d,%d,%d \n id2=%d, relid2=%d,%d,%d,%d",
376 rv, d1->GetId(), relid1[0], relid1[1], relid1[2], relid1[3],
377 d2->GetId(), relid2[0], relid2[1], relid2[2], relid2[3]) ;
382 //____________________________________________________________________________
383 void AliEMCALClusterizerv1::Unload()
385 AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
386 gime->EmcalLoader()->UnloadDigits() ;
387 gime->EmcalLoader()->UnloadRecPoints() ;
390 //____________________________________________________________________________
391 void AliEMCALClusterizerv1::WriteRecPoints(Int_t event)
394 // Creates new branches with given title
395 // fills and writes into TreeR.
397 AliEMCALGetter *gime = AliEMCALGetter::Instance() ;
399 TObjArray * aPRERecPoints = gime->PRERecPoints() ;
400 TObjArray * aECARecPoints = gime->ECARecPoints() ;
401 TObjArray * aHCARecPoints = gime->HCARecPoints() ;
403 TClonesArray * digits = gime->Digits() ;
404 TTree * treeR = gime->TreeR(); ;
408 //Evaluate position, dispersion and other RecPoint properties for PRE section
409 for(index = 0; index < aPRERecPoints->GetEntries(); index++)
410 (dynamic_cast<AliEMCALRecPoint *>(aPRERecPoints->At(index)))->EvalAll(fPREW0,digits) ;
411 aPRERecPoints->Sort() ;
413 for(index = 0; index < aPRERecPoints->GetEntries(); index++)
414 (dynamic_cast<AliEMCALRecPoint *>(aPRERecPoints->At(index)))->SetIndexInList(index) ;
416 aPRERecPoints->Expand(aPRERecPoints->GetEntriesFast()) ;
418 //Evaluate position, dispersion and other RecPoint properties for EC section
419 for(index = 0; index < aECARecPoints->GetEntries(); index++)
420 (dynamic_cast<AliEMCALTowerRecPoint *>(aECARecPoints->At(index)))->EvalAll(fECAW0,digits) ;
422 aECARecPoints->Sort() ;
424 for(index = 0; index < aECARecPoints->GetEntries(); index++)
425 (dynamic_cast<AliEMCALTowerRecPoint *>(aECARecPoints->At(index)))->SetIndexInList(index) ;
427 aECARecPoints->Expand(aECARecPoints->GetEntriesFast()) ;
429 //Evaluate position, dispersion and other RecPoint properties for HCA section
430 for(index = 0; index < aHCARecPoints->GetEntries(); index++)
431 (dynamic_cast<AliEMCALTowerRecPoint *>(aHCARecPoints->At(index)))->EvalAll(fHCAW0,digits) ;
433 aHCARecPoints->Sort() ;
435 for(index = 0; index < aHCARecPoints->GetEntries(); index++)
436 (dynamic_cast<AliEMCALTowerRecPoint *>(aHCARecPoints->At(index)))->SetIndexInList(index) ;
438 aHCARecPoints->Expand(aHCARecPoints->GetEntriesFast()) ;
440 Int_t bufferSize = 32000 ;
441 Int_t splitlevel = 0 ;
444 TBranch * branchPRE = treeR->Branch("EMCALPRERP","TObjArray",&aPRERecPoints,bufferSize,splitlevel);
445 branchPRE->SetTitle(BranchName());
448 TBranch * branchECA = treeR->Branch("EMCALECARP","TObjArray",&aECARecPoints,bufferSize,splitlevel);
449 branchECA->SetTitle(BranchName());
452 TBranch * branchHCA = treeR->Branch("EMCALHCARP","TObjArray",&aHCARecPoints,bufferSize,splitlevel);
453 branchHCA->SetTitle(BranchName());
459 gime->WriteRecPoints("OVERWRITE");
460 gime->WriteClusterizer("OVERWRITE");
463 //____________________________________________________________________________
464 void AliEMCALClusterizerv1::MakeClusters()
466 // Steering method to construct the clusters stored in a list of Reconstructed Points
467 // A cluster is defined as a list of neighbour digits
469 AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
471 AliEMCALGeometry * geom = gime->EMCALGeometry() ;
474 TObjArray * aPRERecPoints = gime->PRERecPoints() ;
475 TObjArray * aECARecPoints = gime->ECARecPoints() ;
476 TObjArray * aHCARecPoints = gime->HCARecPoints() ;
478 aPRERecPoints->Delete() ;
479 aECARecPoints->Delete() ;
480 aHCARecPoints->Delete() ;
482 TClonesArray * digits = gime->Digits() ;
485 AliEMCALDigit * digit ;
486 Int_t ndigECA=0, ndigPRE=0, ndigHCA=0 ;
488 // count the number of digits in ECA section
489 while ( (digit = dynamic_cast<AliEMCALDigit *>(next())) ) { // scan over the list of digits
490 if (geom->IsInECA(digit->GetId()))
492 else if (geom->IsInPRE(digit->GetId()))
494 else if (geom->IsInHCA(digit->GetId()))
497 Error("MakeClusters", "id = %d is a wrong ID!", digit->GetId()) ;
502 // add amplitude of PRE and ECA sections
504 for (digECA = 0 ; digECA < ndigECA ; digECA++) {
505 digit = dynamic_cast<AliEMCALDigit *>(digits->At(digECA)) ;
507 for (digPRE = ndigECA ; digPRE < ndigECA+ndigPRE ; digPRE++) {
508 AliEMCALDigit * digitPRE = dynamic_cast<AliEMCALDigit *>(digits->At(digPRE)) ;
509 if ( geom->AreInSameTower(digit->GetId(), digitPRE->GetId()) ){
510 Float_t amp = static_cast<Float_t>(digit->GetAmp()) + geom->GetSummationFraction() * static_cast<Float_t>(digitPRE->GetAmp()) + 0.5 ;
511 digit->SetAmp(static_cast<Int_t>(amp)) ;
516 Info("MakeClusters","id = %d amp = %d", digit->GetId(), digit->GetAmp()) ;
519 TClonesArray * digitsC = dynamic_cast<TClonesArray*>(digits->Clone()) ;
522 // Clusterization starts
524 TIter nextdigit(digitsC) ;
525 Bool_t notremovedECA = kTRUE, notremovedPRE = kTRUE ;
527 while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigit())) ) { // scan over the list of digitsC
528 AliEMCALRecPoint * clu = 0 ;
530 TArrayI clusterPREdigitslist(50), clusterECAdigitslist(50), clusterHCAdigitslist(50);
532 Bool_t inPRE = kFALSE, inECA = kFALSE, inHCA = kFALSE ;
533 if( geom->IsInPRE(digit->GetId()) ) {
536 else if( geom->IsInECA(digit->GetId()) ) {
539 else if( geom->IsInHCA(digit->GetId()) ) {
545 Info("MakeClusters","id = %d, ene = %f , thre = %f ",
546 digit->GetId(),Calibrate(digit->GetAmp(), 0), fPREClusteringThreshold) ;
548 Info("MakeClusters","id = %d, ene = %f , thre = %f",
549 digit->GetId(),Calibrate(digit->GetAmp(), 1), fECAClusteringThreshold) ;
551 Info("MakeClusters","id = %d, ene = %f , thre = %f",
552 digit->GetId(),Calibrate(digit->GetAmp(), 2), fHCAClusteringThreshold ) ;
555 if ( (inPRE && (Calibrate(digit->GetAmp(), 0) > fPREClusteringThreshold )) ||
556 (inECA && (Calibrate(digit->GetAmp(), 1) > fECAClusteringThreshold )) ||
557 (inHCA && (Calibrate(digit->GetAmp(), 2) > fHCAClusteringThreshold )) ) {
559 Int_t iDigitInPRECluster = 0, iDigitInECACluster = 0, iDigitInHCACluster = 0;
560 Int_t where ; // PRE = 0, ECAl = 1, HCAL = 2
562 // Find the seed in each of the section ECAL/PRE/HCAL
564 if( geom->IsInECA(digit->GetId()) ) {
565 where = 1 ; // to tell we are in ECAL
566 // start a new Tower RecPoint
567 if(fNumberOfECAClusters >= aECARecPoints->GetSize())
568 aECARecPoints->Expand(2*fNumberOfECAClusters+1) ;
569 AliEMCALTowerRecPoint * rp = new AliEMCALTowerRecPoint("") ;
571 aECARecPoints->AddAt(rp, fNumberOfECAClusters) ;
572 clu = dynamic_cast<AliEMCALTowerRecPoint *>(aECARecPoints->At(fNumberOfECAClusters)) ;
573 fNumberOfECAClusters++ ;
574 clu->AddDigit(*digit, Calibrate(digit->GetAmp(), where)) ;
575 clusterECAdigitslist[iDigitInECACluster] = digit->GetIndexInList() ;
576 iDigitInECACluster++ ;
577 digitsC->Remove(digit) ;
579 Info("MakeClusters","OK id = %d, ene = %f , thre = %f ", digit->GetId(),Calibrate(digit->GetAmp(), 1), fECAClusteringThreshold) ;
582 else if( geom->IsInPRE(digit->GetId()) ) {
583 where = 0 ; // to tell we are in PRE
584 // start a new Pre Shower cluster
585 if(fNumberOfPREClusters >= aPRERecPoints->GetSize())
586 aPRERecPoints->Expand(2*fNumberOfPREClusters+1);
587 AliEMCALTowerRecPoint * rp = new AliEMCALTowerRecPoint("") ;
589 aPRERecPoints->AddAt(rp, fNumberOfPREClusters) ;
590 clu = dynamic_cast<AliEMCALTowerRecPoint *>(aPRERecPoints->At(fNumberOfPREClusters)) ;
591 fNumberOfPREClusters++ ;
592 clu->AddDigit(*digit, Calibrate(digit->GetAmp(), where));
593 clusterPREdigitslist[iDigitInPRECluster] = digit->GetIndexInList() ;
594 iDigitInPRECluster++ ;
595 digitsC->Remove(digit) ;
597 Info("MakeClusters","OK id = %d, ene = %f , thre = %f ", digit->GetId(),Calibrate(digit->GetAmp(), 0), fPREClusteringThreshold) ;
601 // Here we remove remaining ECA digits, which cannot make a cluster
603 if( notremovedECA ) {
604 while( ( digit = dynamic_cast<AliEMCALDigit *>(nextdigit()) ) ) {
605 if( geom->IsInECA(digit->GetId()) )
606 digitsC->Remove(digit) ;
610 notremovedECA = kFALSE ;
614 else if( geom->IsInHCA(digit->GetId()) ) {
615 where = 2 ; // to tell we are in HCAL
616 // start a new HCAL cluster
617 if(fNumberOfHCAClusters >= aHCARecPoints->GetSize())
618 aHCARecPoints->Expand(2*fNumberOfHCAClusters+1);
619 AliEMCALTowerRecPoint * rp = new AliEMCALTowerRecPoint("") ;
621 aHCARecPoints->AddAt(rp, fNumberOfHCAClusters) ;
622 clu = dynamic_cast<AliEMCALTowerRecPoint *>(aHCARecPoints->At(fNumberOfHCAClusters)) ;
623 fNumberOfHCAClusters++ ;
624 clu->AddDigit(*digit, Calibrate(digit->GetAmp(), where));
625 clusterHCAdigitslist[iDigitInHCACluster] = digit->GetIndexInList() ;
626 iDigitInHCACluster++ ;
627 digitsC->Remove(digit) ;
629 Info("MakeClusters","OK id = %d, ene = %f , thre = %f ", digit->GetId(),Calibrate(digit->GetAmp(), 2), fHCAClusteringThreshold) ;
633 // Here we remove remaining PRE digits, which cannot make a cluster
635 if( notremovedPRE ) {
636 while( ( digit = dynamic_cast<AliEMCALDigit *>(nextdigit()) ) ) {
637 if( geom->IsInPRE(digit->GetId()) )
638 digitsC->Remove(digit) ;
642 notremovedPRE = kFALSE ;
648 AliEMCALDigit * digitN ;
651 // Do the Clustering in each of the three section ECAL/PRE/HCAL
653 while (index < iDigitInECACluster){ // scan over digits already in cluster
654 digit = (AliEMCALDigit*)digits->At(clusterECAdigitslist[index]) ;
656 while ( (digitN = (AliEMCALDigit *)nextdigit()) ) { // scan over the reduced list of digits
657 Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!!
658 // Info("MakeClusters","id1 = %d, id2 = %d , neighbours = %d", digit->GetId(), digitN->GetId(), ineb) ;
660 case 0 : // not a neighbour
662 case 1 : // are neighbours
663 clu->AddDigit(*digitN, Calibrate( digitN->GetAmp(), 1) ) ;
664 clusterECAdigitslist[iDigitInECACluster] = digitN->GetIndexInList() ;
665 iDigitInECACluster++ ;
666 digitsC->Remove(digitN) ;
668 case 2 : // too far from each other
676 } // loop over ECA cluster
679 while (index < iDigitInPRECluster){ // scan over digits already in cluster
680 digit = (AliEMCALDigit*)digits->At(clusterPREdigitslist[index]) ;
682 while ( (digitN = (AliEMCALDigit *)nextdigit()) ) { // scan over the reduced list of digits
683 Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!!
684 // Info("MakeClusters","id1 = %d, id2 = %d , neighbours = %d", digit->GetId(), digitN->GetId(), ineb) ;
686 case 0 : // not a neighbour
688 case 1 : // are neighbours
689 clu->AddDigit(*digitN, Calibrate( digitN->GetAmp(), 0) ) ;
690 clusterPREdigitslist[iDigitInPRECluster] = digitN->GetIndexInList() ;
691 iDigitInPRECluster++ ;
692 digitsC->Remove(digitN) ;
694 case 2 : // too far from each other
702 } // loop over PRE cluster
705 while (index < iDigitInHCACluster){ // scan over digits already in cluster
706 digit = (AliEMCALDigit*)digits->At(clusterHCAdigitslist[index]) ;
708 while ( (digitN = (AliEMCALDigit *)nextdigit()) ) { // scan over the reduced list of digits
709 Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!!
710 //Info("MakeClusters","id1 = %d, id2 = %d , neighbours = %d", digit->GetId(), digitN->GetId(), ineb) ;
712 case 0 : // not a neighbour
714 case 1 : // are neighbours
715 clu->AddDigit(*digitN, Calibrate( digitN->GetAmp(), 2) ) ;
716 clusterHCAdigitslist[iDigitInHCACluster] = digitN->GetIndexInList() ;
717 iDigitInHCACluster++ ;
718 digitsC->Remove(digitN) ;
720 case 2 : // too far from each other
727 } // loop over HCA cluster
734 //____________________________________________________________________________
735 void AliEMCALClusterizerv1::MakeUnfolding()
737 Fatal("AliEMCALClusterizerv1::MakeUnfolding", "--> Unfolding not implemented") ;
741 //____________________________________________________________________________
742 Double_t AliEMCALClusterizerv1::ShowerShape(Double_t r)
744 // Shape of the shower (see EMCAL TDR)
745 // If you change this function, change also the gradient evaluation in ChiSquare()
747 Double_t r4 = r*r*r*r ;
748 Double_t r295 = TMath::Power(r, 2.95) ;
749 Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
753 //____________________________________________________________________________
754 void AliEMCALClusterizerv1::UnfoldCluster(AliEMCALTowerRecPoint * iniTower,
756 AliEMCALDigit ** maxAt,
757 Float_t * maxAtEnergy)
759 // Performs the unfolding of a cluster with nMax overlapping showers
761 Fatal("UnfoldCluster", "--> Unfolding not implemented") ;
765 //_____________________________________________________________________________
766 void AliEMCALClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
768 // Calculates the Chi square for the cluster unfolding minimization
769 // Number of parameters, Gradient, Chi squared, parameters, what to do
771 ::Fatal("UnfoldingChiSquare","Unfolding not implemented") ;
773 //____________________________________________________________________________
774 void AliEMCALClusterizerv1::Print(Option_t * option)const
776 // Print clusterizer parameters
778 TString message("\n") ;
780 if( strcmp(GetName(), "") !=0 ){
784 TString taskName(GetName()) ;
785 taskName.ReplaceAll(Version(), "") ;
787 message += "--------------- " ;
788 message += taskName.Data() ;
790 message += GetTitle() ;
791 message += "-----------\n" ;
792 message += "Clusterizing digits from the file: " ;
793 message += taskName.Data() ;
794 message += "\n Branch: " ;
795 message += GetName() ;
796 message += "\n Pre Shower Clustering threshold = " ;
797 message += fPREClusteringThreshold ;
798 message += "\n Pre Shower Local Maximum cut = " ;
799 message += fPRELocMaxCut ;
800 message += "\n Pre Shower Logarothmic weight = " ;
802 message += "\n ECA Clustering threshold = " ;
803 message += fECAClusteringThreshold ;
804 message += "\n ECA Local Maximum cut = " ;
805 message += fECALocMaxCut ;
806 message += "\n ECA Logarothmic weight = " ;
808 message += "\n Pre Shower Clustering threshold = " ;
809 message += fHCAClusteringThreshold ;
810 message += "\n HCA Local Maximum cut = " ;
811 message += fHCALocMaxCut ;
812 message += "\n HCA Logarothmic weight = " ;
815 message +="\nUnfolding on\n" ;
817 message += "\nUnfolding off\n";
819 message += "------------------------------------------------------------------" ;
822 message += "AliEMCALClusterizerv1 not initialized " ;
824 Info("Print", message.Data() ) ;
827 //____________________________________________________________________________
828 void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option)
830 // Prints list of RecPoints produced at the current pass of AliEMCALClusterizer
832 TObjArray * aPRERecPoints = AliEMCALGetter::Instance()->PRERecPoints() ;
833 TObjArray * aECARecPoints = AliEMCALGetter::Instance()->ECARecPoints() ;
834 TObjArray * aHCARecPoints = AliEMCALGetter::Instance()->HCARecPoints() ;
836 Info("PrintRecPoints", "Clusterization result:") ;
838 printf("event # %d\n", gAlice->GetEvNumber() ) ;
839 printf(" Found %d PRE SHOWER RecPoints, %d ECA Rec Points and %d HCA Rec Points\n ",
840 aPRERecPoints->GetEntriesFast(), aECARecPoints->GetEntriesFast(), aHCARecPoints->GetEntriesFast() ) ;
842 fRecPointsInRun += aPRERecPoints->GetEntriesFast() ;
843 fRecPointsInRun += aECARecPoints->GetEntriesFast() ;
844 fRecPointsInRun += aHCARecPoints->GetEntriesFast() ;
846 if(strstr(option,"all")) {
848 //Pre shower recPoints
850 printf("-----------------------------------------------------------------------\n") ;
851 printf("Clusters in PRE section\n") ;
852 printf("Index Ene(GeV) Multi Module phi r theta X Y Z Dispersion Lambda 1 Lambda 2 # of prim Primaries list\n") ;
856 for (index = 0 ; index < aPRERecPoints->GetEntries() ; index++) {
857 AliEMCALTowerRecPoint * rp = dynamic_cast<AliEMCALTowerRecPoint *>(aPRERecPoints->At(index)) ;
859 rp->GetGlobalPosition(globalpos);
861 rp->GetLocalPosition(localpos);
863 rp->GetElipsAxis(lambda);
866 primaries = rp->GetPrimaries(nprimaries);
867 printf("\n%6d %8.4f %3d %2d %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4f %4f %2d : ",
868 rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetEMCALArm(),
869 globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(),
870 rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ;
871 for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
872 printf("%d ", primaries[iprimary] ) ;
876 printf("\n-----------------------------------------------------------------------\n") ;
877 printf("Clusters in ECAL section\n") ;
878 printf("Index Ene(GeV) Multi Module phi r theta X Y Z Dispersion Lambda 1 Lambda 2 # of prim Primaries list\n") ;
880 for (index = 0 ; index < aECARecPoints->GetEntries() ; index++) {
881 AliEMCALTowerRecPoint * rp = dynamic_cast<AliEMCALTowerRecPoint * >(aECARecPoints->At(index)) ;
883 rp->GetGlobalPosition(globalpos);
885 rp->GetLocalPosition(localpos);
887 rp->GetElipsAxis(lambda);
890 primaries = rp->GetPrimaries(nprimaries);
891 printf("\n%6d %8.4f %3d %2d %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4f %4f %2d : ",
892 rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetEMCALArm(),
893 globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(),
894 rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ;
895 for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
896 printf("%d ", primaries[iprimary] ) ;
900 printf("\n-----------------------------------------------------------------------\n") ;
901 printf("Clusters in HCAL section\n") ;
902 printf("Index Ene(GeV) Multi Module phi r theta X Y Z Dispersion Lambda 1 Lambda 2 # of prim Primaries list\n") ;
904 for (index = 0 ; index < aHCARecPoints->GetEntries() ; index++) {
905 AliEMCALTowerRecPoint * rp = dynamic_cast<AliEMCALTowerRecPoint * >(aHCARecPoints->At(index)) ;
907 rp->GetGlobalPosition(globalpos);
909 rp->GetLocalPosition(localpos);
911 rp->GetElipsAxis(lambda);
914 primaries = rp->GetPrimaries(nprimaries);
915 printf("\n%6d %8.4f %3d %2d %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4f %4f %2d : ",
916 rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetEMCALArm(),
917 globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(),
918 rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ;
919 for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
920 printf("%d ", primaries[iprimary] ) ;
924 printf("\n-----------------------------------------------------------------------\n");