Transition to NewIO
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALClusterizerv1.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17 /* $Log:
18    1 October 2000. Yuri Kharlov:
19      AreNeighbours()
20      PPSD upper layer is considered if number of layers>1
21    18 October 2000. Yuri Kharlov:
22      AliEMCALClusterizerv1()
23      CPV clusterizing parameters added
24      MakeClusters()
25      After first PPSD digit remove EMC digits only once
26 */
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 
37 //  standalone mode.
38 // Use Case:
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
55
56 // --- ROOT system ---
57
58 #include "TROOT.h" 
59 #include "TFile.h" 
60 #include "TFolder.h" 
61 #include "TMath.h" 
62 #include "TMinuit.h"
63 #include "TTree.h" 
64 #include "TSystem.h" 
65 #include "TBenchmark.h"
66
67 // --- Standard library ---
68
69
70 // --- AliRoot header files ---
71 #include "AliEMCALGetter.h"
72 #include "AliEMCALClusterizerv1.h"
73 #include "AliEMCALTowerRecPoint.h"
74 #include "AliEMCALDigit.h"
75 #include "AliEMCALDigitizer.h"
76 #include "AliEMCAL.h"
77 #include "AliEMCALGeometry.h"
78
79 ClassImp(AliEMCALClusterizerv1)
80   
81 //____________________________________________________________________________
82   AliEMCALClusterizerv1::AliEMCALClusterizerv1() : AliEMCALClusterizer()
83 {
84   // default ctor (to be used mainly by Streamer)
85   
86   InitParameters() ; 
87   fDefaultInit = kTRUE ; 
88 }
89
90 //____________________________________________________________________________
91 AliEMCALClusterizerv1::AliEMCALClusterizerv1(const TString alirunFileName, const TString eventFolderName)
92 :AliEMCALClusterizer(alirunFileName, eventFolderName)
93 {
94   // ctor with the indication of the file where header Tree and digits Tree are stored
95   
96   InitParameters() ; 
97   Init() ;
98   fDefaultInit = kFALSE ; 
99
100 }
101
102 //____________________________________________________________________________
103   AliEMCALClusterizerv1::~AliEMCALClusterizerv1()
104 {
105   // dtor
106   
107 }
108
109 //____________________________________________________________________________
110 const TString AliEMCALClusterizerv1::BranchName() const 
111
112    return GetName();
113
114 }
115
116 //____________________________________________________________________________
117 Float_t  AliEMCALClusterizerv1::Calibrate(Int_t amp, Int_t where) const
118 {
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 ;
127   else 
128     Fatal("Calibrate", "Something went wrong!") ;
129   return -9999999. ; 
130 }
131
132 //____________________________________________________________________________
133 void AliEMCALClusterizerv1::Exec(Option_t * option)
134 {
135   // Steering method
136
137   if(strstr(option,"tim"))
138     gBenchmark->Start("EMCALClusterizer"); 
139   
140   if(strstr(option,"print"))
141     Print("") ; 
142
143   AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
144
145   Int_t nevents = gime->MaxEvent() ;
146   Int_t ievent ;
147
148   for(ievent = 0; ievent < nevents; ievent++){
149
150     gime->Event(ievent,"D") ;
151
152     GetCalibrationParameters() ;
153
154     fNumberOfPREClusters = fNumberOfECAClusters = fNumberOfHCAClusters = 0 ;
155            
156     MakeClusters() ;
157
158     if(fToUnfold)
159       MakeUnfolding() ;
160
161     WriteRecPoints(ievent) ;
162
163     if(strstr(option,"deb"))  
164       PrintRecPoints(option) ;
165
166     //increment the total number of recpoints per run 
167     fRecPointsInRun += gime->PRERecPoints()->GetEntriesFast() ;  
168     fRecPointsInRun += gime->ECARecPoints()->GetEntriesFast() ;  
169     fRecPointsInRun += gime->HCARecPoints()->GetEntriesFast() ;  
170   }
171   
172   Unload();
173
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 ) ;
178   }
179 }
180
181 //____________________________________________________________________________
182 Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALTowerRecPoint * emcRP, AliEMCALDigit ** maxAt, Float_t * maxAtEnergy,
183                                     Int_t nPar, Float_t * fitparameters) const
184
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
188
189   AliEMCALGetter * gime = AliEMCALGetter::Instance() ; 
190   TClonesArray * digits = gime->Digits() ; 
191   
192
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) ;
200   
201   gMinuit->SetObjectFit(toMinuit) ;         // To tranfer pointer to UnfoldingChiSquare
202
203   // filling initial values for fit parameters
204   AliEMCALDigit * digit ;
205
206   Int_t ierflg  = 0; 
207   Int_t index   = 0 ;
208   Int_t nDigits = (Int_t) nPar / 3 ;
209
210   Int_t iDigit ;
211
212   AliEMCALGeometry * geom = gime->EMCALGeometry() ; 
213
214   for(iDigit = 0; iDigit < nDigits; iDigit++){
215     digit = maxAt[iDigit]; 
216
217     Int_t relid[4] ;
218     Float_t x = 0.;
219     Float_t z = 0.;
220     geom->AbsToRelNumbering(digit->GetId(), relid) ;
221     geom->PosInAlice(relid, x, z) ;
222
223     Float_t energy = maxAtEnergy[iDigit] ;
224
225     gMinuit->mnparm(index, "x",  x, 0.1, 0, 0, ierflg) ;
226     index++ ;   
227     if(ierflg != 0){ 
228       Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : x = %f",  x ) ;
229       return kFALSE;
230     }
231     gMinuit->mnparm(index, "z",  z, 0.1, 0, 0, ierflg) ;
232     index++ ;   
233     if(ierflg != 0){
234        Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : z = %f", z) ;
235       return kFALSE;
236     }
237     gMinuit->mnparm(index, "Energy",  energy , 0.05*energy, 0., 4.*energy, ierflg) ;
238     index++ ;   
239     if(ierflg != 0){
240      Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : energy = %f", energy) ;      
241       return kFALSE;
242     }
243   }
244
245   Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
246                       //  depends on it. 
247   Double_t p1 = 1.0 ;
248   Double_t p2 = 0.0 ;
249
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
254
255   gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ;    // minimize 
256
257   if(ierflg == 4){  // Minimum not found   
258     Error("FindFit", "EMCAL Unfolding  Fit not converged, cluster abandoned " ) ;      
259     return kFALSE ;
260   }            
261   for(index = 0; index < nPar; index++){
262     Double_t err ;
263     Double_t val ;
264     gMinuit->GetParameter(index, val, err) ;    // Returns value and error of parameter index
265     fitparameters[index] = val ;
266    }
267
268   delete toMinuit ;
269   return kTRUE;
270
271 }
272
273 //____________________________________________________________________________
274 void AliEMCALClusterizerv1::GetCalibrationParameters() 
275 {
276   AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
277
278   if ( !gime->Digitizer() ) 
279     gime->LoadDigitizer();
280   AliEMCALDigitizer * dig = gime->Digitizer(); 
281    
282   fADCchannelPRE  = dig->GetPREchannel() ;
283   fADCpedestalPRE = dig->GetPREpedestal() ; 
284
285   fADCchannelECA   = dig->GetECAchannel() ;
286   fADCpedestalECA  = dig->GetECApedestal();
287
288   fADCchannelHCA   = dig->GetHCAchannel() ;
289   fADCpedestalHCA  = dig->GetHCApedestal();
290 }
291
292 //____________________________________________________________________________
293 void AliEMCALClusterizerv1::Init()
294 {
295   // Make all memory allocations which can not be done in default constructor.
296   // Attach the Clusterizer task to the list of EMCAL tasks
297   
298   AliEMCALGetter * gime = AliEMCALGetter::Instance(GetTitle(), fEventFolderName.Data());
299
300   AliEMCALGeometry * geom = gime->EMCALGeometry() ;
301
302   fNTowers = geom->GetNZ() *  geom->GetNPhi() ;
303   if(!gMinuit) 
304     gMinuit = new TMinuit(100) ;
305
306  if ( !gime->Clusterizer() ) 
307     gime->PostClusterizer(this); 
308 }
309
310 //____________________________________________________________________________
311 void AliEMCALClusterizerv1::InitParameters()
312 {
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 ;
320   
321   fPREW0    = 4.0 ;
322   fECAW0     = 4.5 ;
323   fHCAW0     = 4.5 ;
324
325   fTimeGate = 1.e-8 ; 
326   
327   fToUnfold = kFALSE ;
328    
329   fRecPointsInRun          = 0 ;
330  
331 }
332
333 //____________________________________________________________________________
334 Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2)const
335 {
336   // Gives the neighbourness of two digits = 0 are not neighbour but continue searching 
337   //                                       = 1 are neighbour
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  
342
343    AliEMCALGeometry * geom = AliEMCALGetter::Instance()->EMCALGeometry() ;
344
345   Int_t rv = 0 ; 
346
347   Int_t relid1[4] ; 
348   geom->AbsToRelNumbering(d1->GetId(), relid1) ; 
349
350   Int_t relid2[4] ; 
351   geom->AbsToRelNumbering(d2->GetId(), relid2) ; 
352   
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] ) ;  
357     
358     if (( coldiff <= 1 )  && ( rowdiff <= 1 )){
359       if((relid1[1] != 0) || (TMath::Abs(d1->GetTime() - d2->GetTime() ) < fTimeGate))
360       rv = 1 ; 
361     }
362     else {
363       if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1)) 
364         rv = 2; //  Difference in row numbers is too large to look further 
365     }
366
367   } 
368   else {
369     
370     if( (relid1[0] < relid2[0]) || (relid1[1] != relid2[1]) )  
371       rv=0 ;
372   }
373
374   if (gDebug == 2 ) 
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]) ;   
378   
379   return rv ; 
380 }
381
382 //____________________________________________________________________________
383 void AliEMCALClusterizerv1::Unload() 
384 {
385   AliEMCALGetter * gime = AliEMCALGetter::Instance() ; 
386   gime->EmcalLoader()->UnloadDigits() ; 
387   gime->EmcalLoader()->UnloadRecPoints() ; 
388 }
389  
390 //____________________________________________________________________________
391 void AliEMCALClusterizerv1::WriteRecPoints(Int_t event)
392 {
393
394   // Creates new branches with given title
395   // fills and writes into TreeR.
396
397   AliEMCALGetter *gime = AliEMCALGetter::Instance() ; 
398
399   TObjArray * aPRERecPoints = gime->PRERecPoints() ; 
400   TObjArray * aECARecPoints = gime->ECARecPoints() ; 
401   TObjArray * aHCARecPoints = gime->HCARecPoints() ; 
402
403   TClonesArray * digits = gime->Digits() ; 
404   TTree * treeR = gime->TreeR(); ; 
405   
406   Int_t index ;
407
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() ;
412   
413   for(index = 0; index < aPRERecPoints->GetEntries(); index++)
414     (dynamic_cast<AliEMCALRecPoint *>(aPRERecPoints->At(index)))->SetIndexInList(index) ;
415   
416   aPRERecPoints->Expand(aPRERecPoints->GetEntriesFast()) ;
417   
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) ;
421   
422   aECARecPoints->Sort() ;
423
424   for(index = 0; index < aECARecPoints->GetEntries(); index++)
425     (dynamic_cast<AliEMCALTowerRecPoint *>(aECARecPoints->At(index)))->SetIndexInList(index) ;
426
427   aECARecPoints->Expand(aECARecPoints->GetEntriesFast()) ; 
428   
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) ;
432   
433   aHCARecPoints->Sort() ;
434
435   for(index = 0; index < aHCARecPoints->GetEntries(); index++)
436     (dynamic_cast<AliEMCALTowerRecPoint *>(aHCARecPoints->At(index)))->SetIndexInList(index) ;
437
438   aHCARecPoints->Expand(aHCARecPoints->GetEntriesFast()) ; 
439  
440   Int_t bufferSize = 32000 ;    
441   Int_t splitlevel = 0 ; 
442   
443   //PRE section branch 
444   TBranch * branchPRE = treeR->Branch("EMCALPRERP","TObjArray",&aPRERecPoints,bufferSize,splitlevel);
445   branchPRE->SetTitle(BranchName());
446
447   //EC section branch
448   TBranch * branchECA = treeR->Branch("EMCALECARP","TObjArray",&aECARecPoints,bufferSize,splitlevel);
449   branchECA->SetTitle(BranchName());
450
451   //HCA section branch
452   TBranch * branchHCA = treeR->Branch("EMCALHCARP","TObjArray",&aHCARecPoints,bufferSize,splitlevel);
453   branchHCA->SetTitle(BranchName());
454
455   branchPRE->Fill() ;
456   branchECA->Fill() ;
457   branchHCA->Fill() ;
458
459   gime->WriteRecPoints("OVERWRITE");
460   gime->WriteClusterizer("OVERWRITE");
461 }
462
463 //____________________________________________________________________________
464 void AliEMCALClusterizerv1::MakeClusters()
465 {
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
468     
469   AliEMCALGetter * gime = AliEMCALGetter::Instance() ; 
470
471   AliEMCALGeometry * geom = gime->EMCALGeometry() ; 
472
473
474   TObjArray * aPRERecPoints = gime->PRERecPoints() ; 
475   TObjArray * aECARecPoints  = gime->ECARecPoints() ; 
476   TObjArray * aHCARecPoints  = gime->HCARecPoints() ; 
477
478   aPRERecPoints->Delete() ;
479   aECARecPoints->Delete() ;
480   aHCARecPoints->Delete() ;
481
482   TClonesArray * digits = gime->Digits() ; 
483
484   TIter next(digits) ; 
485   AliEMCALDigit * digit ; 
486   Int_t ndigECA=0, ndigPRE=0, ndigHCA=0 ; 
487
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())) 
491       ndigECA++ ; 
492     else if (geom->IsInPRE(digit->GetId()))
493       ndigPRE++; 
494     else if (geom->IsInHCA(digit->GetId()))
495       ndigHCA++;
496     else {
497       Error("MakeClusters", "id = %d is a wrong ID!", digit->GetId()) ; 
498       abort() ;
499     }
500   }
501
502   // add amplitude of PRE and ECA sections
503   Int_t digECA ; 
504   for (digECA = 0 ; digECA < ndigECA ; digECA++) {
505     digit = dynamic_cast<AliEMCALDigit *>(digits->At(digECA)) ;
506     Int_t digPRE ;
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)) ; 
512         break ; 
513       }
514     }
515     if (gDebug) 
516       Info("MakeClusters","id = %d amp = %d", digit->GetId(), digit->GetAmp()) ; 
517   }  
518
519   TClonesArray * digitsC =  dynamic_cast<TClonesArray*>(digits->Clone()) ;
520   
521   
522   // Clusterization starts  
523   
524   TIter nextdigit(digitsC) ; 
525   Bool_t notremovedECA = kTRUE, notremovedPRE = kTRUE ;
526   
527   while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigit())) ) { // scan over the list of digitsC
528     AliEMCALRecPoint * clu = 0 ; 
529     
530     TArrayI clusterPREdigitslist(50), clusterECAdigitslist(50), clusterHCAdigitslist(50);   
531  
532     Bool_t inPRE = kFALSE, inECA = kFALSE, inHCA = kFALSE ;
533     if( geom->IsInPRE(digit->GetId()) ) {
534       inPRE = kTRUE ; 
535     }
536     else if( geom->IsInECA(digit->GetId()) ) {
537       inECA = kTRUE ;
538     }
539     else if( geom->IsInHCA(digit->GetId()) ) {
540       inHCA = kTRUE ;
541     }
542     
543     if (gDebug == 2) { 
544       if (inPRE)
545         Info("MakeClusters","id = %d, ene = %f , thre = %f ", 
546              digit->GetId(),Calibrate(digit->GetAmp(), 0), fPREClusteringThreshold) ;  
547       if (inECA)
548         Info("MakeClusters","id = %d, ene = %f , thre = %f", 
549              digit->GetId(),Calibrate(digit->GetAmp(), 1), fECAClusteringThreshold) ;  
550       if (inHCA)
551         Info("MakeClusters","id = %d, ene = %f , thre = %f", 
552              digit->GetId(),Calibrate(digit->GetAmp(), 2), fHCAClusteringThreshold ) ;  
553     }
554     
555     if ( (inPRE  && (Calibrate(digit->GetAmp(), 0) > fPREClusteringThreshold  )) || 
556          (inECA && (Calibrate(digit->GetAmp(), 1) > fECAClusteringThreshold  ))  || 
557          (inHCA && (Calibrate(digit->GetAmp(), 2) > fHCAClusteringThreshold  )) ) {
558       
559       Int_t  iDigitInPRECluster = 0, iDigitInECACluster = 0, iDigitInHCACluster = 0; 
560       Int_t where ; // PRE = 0, ECAl = 1, HCAL = 2
561
562       // Find the seed in each of the section ECAL/PRE/HCAL
563
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("") ; 
570         rp->SetECA() ; 
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) ; 
578         if (gDebug == 2 ) 
579           Info("MakeClusters","OK id = %d, ene = %f , thre = %f ", digit->GetId(),Calibrate(digit->GetAmp(), 1), fECAClusteringThreshold) ;  
580         
581       } 
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("") ;    
588         rp->SetPRE() ; 
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) ;
596         if (gDebug == 2 ) 
597           Info("MakeClusters","OK id = %d, ene = %f , thre = %f ", digit->GetId(),Calibrate(digit->GetAmp(), 0), fPREClusteringThreshold) ;  
598         
599         nextdigit.Reset() ;
600
601         // Here we remove remaining ECA digits, which cannot make a cluster
602         
603         if( notremovedECA ) { 
604           while( ( digit = dynamic_cast<AliEMCALDigit *>(nextdigit()) ) ) {
605             if( geom->IsInECA(digit->GetId()) )
606               digitsC->Remove(digit) ;
607             else 
608               break ; 
609           }
610           notremovedECA = kFALSE ;
611         }
612
613       } 
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("") ;    
620         rp->SetHCA() ; 
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) ;
628         if (gDebug == 2 ) 
629           Info("MakeClusters","OK id = %d, ene = %f , thre = %f ", digit->GetId(),Calibrate(digit->GetAmp(), 2), fHCAClusteringThreshold) ;  
630  
631         nextdigit.Reset() ;
632    
633         // Here we remove remaining PRE digits, which cannot make a cluster
634         
635         if( notremovedPRE ) { 
636           while( ( digit = dynamic_cast<AliEMCALDigit *>(nextdigit()) ) ) {
637             if( geom->IsInPRE(digit->GetId()) )
638               digitsC->Remove(digit) ;
639             else 
640               break ; 
641           }
642           notremovedPRE = kFALSE ;
643         }       
644       }        
645       
646       nextdigit.Reset() ;
647       
648       AliEMCALDigit * digitN ; 
649       Int_t index = 0 ;
650
651       // Do the Clustering in each of the three section ECAL/PRE/HCAL
652
653       while (index < iDigitInECACluster){ // scan over digits already in cluster 
654         digit =  (AliEMCALDigit*)digits->At(clusterECAdigitslist[index])  ;      
655         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) ;  
659          switch (ineb ) {
660           case 0 :   // not a neighbour
661             break ;
662           case 1 :   // are neighbours 
663             clu->AddDigit(*digitN, Calibrate( digitN->GetAmp(), 1) ) ;
664             clusterECAdigitslist[iDigitInECACluster] = digitN->GetIndexInList() ; 
665             iDigitInECACluster++ ; 
666             digitsC->Remove(digitN) ;
667             break ;
668           case 2 :   // too far from each other
669             goto endofloop1;   
670           } // switch
671           
672         } // while digitN
673         
674       endofloop1: ;
675         nextdigit.Reset() ; 
676       } // loop over ECA cluster
677       
678       index = 0 ; 
679       while (index < iDigitInPRECluster){ // scan over digits already in cluster 
680         digit =  (AliEMCALDigit*)digits->At(clusterPREdigitslist[index])  ;      
681         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) ;  
685           switch (ineb ) {
686           case 0 :   // not a neighbour
687             break ;
688           case 1 :   // are neighbours 
689             clu->AddDigit(*digitN, Calibrate( digitN->GetAmp(), 0) ) ;
690             clusterPREdigitslist[iDigitInPRECluster] = digitN->GetIndexInList() ; 
691             iDigitInPRECluster++ ; 
692             digitsC->Remove(digitN) ;
693             break ;
694           case 2 :   // too far from each other
695             goto endofloop2;   
696           } // switch
697           
698         } // while digitN
699         
700       endofloop2: ;
701         nextdigit.Reset() ; 
702       } // loop over PRE cluster
703     
704       index = 0 ; 
705       while (index < iDigitInHCACluster){ // scan over digits already in cluster 
706         digit =  (AliEMCALDigit*)digits->At(clusterHCAdigitslist[index])  ;      
707         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) ;  
711           switch (ineb ) {
712           case 0 :   // not a neighbour
713             break ;
714           case 1 :   // are neighbours 
715             clu->AddDigit(*digitN, Calibrate( digitN->GetAmp(), 2) ) ;
716             clusterHCAdigitslist[iDigitInHCACluster] = digitN->GetIndexInList() ; 
717             iDigitInHCACluster++ ; 
718             digitsC->Remove(digitN) ;
719             break ;
720           case 2 :   // too far from each other
721             goto endofloop3;   
722           } // switch  
723         } // while digitN
724         
725       endofloop3: ;
726         nextdigit.Reset() ; 
727       } // loop over HCA cluster
728
729     } // energy theshold     
730   } // while digit  
731   delete digitsC ;
732 }
733
734 //____________________________________________________________________________
735 void AliEMCALClusterizerv1::MakeUnfolding()
736 {
737   Fatal("AliEMCALClusterizerv1::MakeUnfolding", "--> Unfolding not implemented") ;
738  
739 }
740
741 //____________________________________________________________________________
742 Double_t  AliEMCALClusterizerv1::ShowerShape(Double_t r)
743
744   // Shape of the shower (see EMCAL TDR)
745   // If you change this function, change also the gradient evaluation in ChiSquare()
746
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) ) ) ;
750   return shape ;
751 }
752
753 //____________________________________________________________________________
754 void  AliEMCALClusterizerv1::UnfoldCluster(AliEMCALTowerRecPoint * iniTower, 
755                                                  Int_t nMax, 
756                                                  AliEMCALDigit ** maxAt, 
757                                                  Float_t * maxAtEnergy)
758 {
759   // Performs the unfolding of a cluster with nMax overlapping showers 
760   
761   Fatal("UnfoldCluster", "--> Unfolding not implemented") ;
762
763 }
764
765 //_____________________________________________________________________________
766 void AliEMCALClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
767 {
768   // Calculates the Chi square for the cluster unfolding minimization
769   // Number of parameters, Gradient, Chi squared, parameters, what to do
770   
771   ::Fatal("UnfoldingChiSquare","Unfolding not implemented") ;
772 }
773 //____________________________________________________________________________
774 void AliEMCALClusterizerv1::Print(Option_t * option)const
775 {
776   // Print clusterizer parameters
777
778   TString message("\n") ; 
779   
780   if( strcmp(GetName(), "") !=0 ){
781     
782     // Print parameters
783  
784     TString taskName(GetName()) ; 
785     taskName.ReplaceAll(Version(), "") ;
786     
787     message += "--------------- " ; 
788     message += taskName.Data() ; 
789     message += " " ; 
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   = " ; 
801     message += fPREW0 ;
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   = " ;
807     message += fECAW0 ;
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   = " ;
813     message += fHCAW0 ;
814     if(fToUnfold)
815       message +="\nUnfolding on\n" ;
816     else
817       message += "\nUnfolding off\n";
818     
819     message += "------------------------------------------------------------------" ; 
820   }
821   else
822     message += "AliEMCALClusterizerv1 not initialized " ;
823   
824   Info("Print", message.Data() ) ; 
825 }
826
827 //____________________________________________________________________________
828 void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option)
829 {
830   // Prints list of RecPoints produced at the current pass of AliEMCALClusterizer
831
832   TObjArray * aPRERecPoints = AliEMCALGetter::Instance()->PRERecPoints() ; 
833   TObjArray * aECARecPoints = AliEMCALGetter::Instance()->ECARecPoints() ; 
834   TObjArray * aHCARecPoints = AliEMCALGetter::Instance()->HCARecPoints() ; 
835
836   Info("PrintRecPoints", "Clusterization result:") ; 
837   
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() ) ; 
841
842   fRecPointsInRun +=  aPRERecPoints->GetEntriesFast() ; 
843   fRecPointsInRun +=  aECARecPoints->GetEntriesFast() ; 
844   fRecPointsInRun +=  aHCARecPoints->GetEntriesFast() ; 
845   
846   if(strstr(option,"all")) {
847
848     //Pre shower recPoints
849
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") ;      
853
854     Int_t index ;
855     
856     for (index = 0 ; index < aPRERecPoints->GetEntries() ; index++) {
857       AliEMCALTowerRecPoint * rp = dynamic_cast<AliEMCALTowerRecPoint *>(aPRERecPoints->At(index)) ; 
858       TVector3  globalpos;  
859       rp->GetGlobalPosition(globalpos);
860       TVector3  localpos;  
861       rp->GetLocalPosition(localpos);
862       Float_t lambda[2]; 
863       rp->GetElipsAxis(lambda);
864       Int_t * primaries;
865       Int_t nprimaries;
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] ) ; 
873       }          
874     }
875     
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") ;      
879     
880     for (index = 0 ; index < aECARecPoints->GetEntries() ; index++) {
881       AliEMCALTowerRecPoint * rp = dynamic_cast<AliEMCALTowerRecPoint * >(aECARecPoints->At(index)) ; 
882       TVector3  globalpos;  
883       rp->GetGlobalPosition(globalpos);
884       TVector3  localpos;  
885       rp->GetLocalPosition(localpos);
886       Float_t lambda[2]; 
887       rp->GetElipsAxis(lambda);
888       Int_t * primaries; 
889       Int_t nprimaries;
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] ) ; 
897       } 
898     }
899
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") ;      
903     
904     for (index = 0 ; index < aHCARecPoints->GetEntries() ; index++) {
905       AliEMCALTowerRecPoint * rp = dynamic_cast<AliEMCALTowerRecPoint * >(aHCARecPoints->At(index)) ; 
906       TVector3  globalpos;  
907       rp->GetGlobalPosition(globalpos);
908       TVector3  localpos;  
909       rp->GetLocalPosition(localpos);
910       Float_t lambda[2]; 
911       rp->GetElipsAxis(lambda);
912       Int_t * primaries; 
913       Int_t nprimaries;
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] ) ; 
921       } 
922     }
923
924     printf("\n-----------------------------------------------------------------------\n");
925   }
926 }