removed iostream
[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
72 #include "AliEMCALClusterizerv1.h"
73 #include "AliEMCALDigit.h"
74 #include "AliEMCALDigitizer.h"
75 #include "AliEMCALTowerRecPoint.h"
76 #include "AliEMCAL.h"
77 #include "AliEMCALGetter.h"
78 #include "AliEMCALGeometry.h"
79 #include "AliRun.h"
80
81 ClassImp(AliEMCALClusterizerv1)
82   
83 //____________________________________________________________________________
84   AliEMCALClusterizerv1::AliEMCALClusterizerv1() : AliEMCALClusterizer()
85 {
86   // default ctor (to be used mainly by Streamer)
87   
88   InitParameters() ; 
89   fDefaultInit = kTRUE ; 
90 }
91
92 //____________________________________________________________________________
93 AliEMCALClusterizerv1::AliEMCALClusterizerv1(const char* headerFile, const char* name, const Bool_t toSplit)
94 :AliEMCALClusterizer(headerFile, name, toSplit)
95 {
96   // ctor with the indication of the file where header Tree and digits Tree are stored
97   
98   InitParameters() ; 
99   Init() ;
100   fDefaultInit = kFALSE ; 
101
102 }
103
104 //____________________________________________________________________________
105   AliEMCALClusterizerv1::~AliEMCALClusterizerv1()
106 {
107   // dtor
108   fSplitFile = 0 ; 
109   
110 }
111
112 //____________________________________________________________________________
113 const TString AliEMCALClusterizerv1::BranchName() const 
114
115   TString branchName(GetName() ) ;
116   branchName.Remove(branchName.Index(Version())-1) ;
117   return branchName ;
118 }
119
120 //____________________________________________________________________________
121 Float_t  AliEMCALClusterizerv1::Calibrate(Int_t amp, Bool_t inpresho) const
122 {
123   //To be replased later by the method, reading individual parameters from the database
124   if ( inpresho ) // calibrate as pre shower
125     return -fADCpedestalPreSho + amp * fADCchannelPreSho ; 
126   else //calibrate as tower 
127     return -fADCpedestalTower + amp * fADCchannelTower ;                
128 }
129
130 //____________________________________________________________________________
131 void AliEMCALClusterizerv1::Exec(Option_t * option)
132 {
133   // Steering method
134
135   if( strcmp(GetName(), "")== 0 ) 
136     Init() ;
137
138   if(strstr(option,"tim"))
139     gBenchmark->Start("EMCALClusterizer"); 
140   
141   if(strstr(option,"print"))
142     Print("") ; 
143
144   AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
145   if(gime->BranchExists("RecPoints"))
146     return ;
147   Int_t nevents = gime->MaxEvent() ;
148   Int_t ievent ;
149
150   for(ievent = 0; ievent < nevents; ievent++){
151
152     gime->Event(ievent,"D") ;
153
154     if(ievent == 0)
155       GetCalibrationParameters() ;
156
157     fNumberOfTowerClusters = fNumberOfPreShoClusters = 0 ;
158            
159     MakeClusters() ;
160     
161     if(fToUnfold)
162       MakeUnfolding() ;
163
164     WriteRecPoints(ievent) ;
165
166     if(strstr(option,"deb"))  
167       PrintRecPoints(option) ;
168
169     //increment the total number of digits per run 
170     fRecPointsInRun += gime->TowerRecPoints()->GetEntriesFast() ;  
171     fRecPointsInRun += gime->PreShowerRecPoints()->GetEntriesFast() ;  
172  }
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 //____________________________________________________________________________
183 Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALTowerRecPoint * emcRP, AliEMCALDigit ** maxAt, Float_t * maxAtEnergy,
184                                     Int_t nPar, Float_t * fitparameters) const
185
186   // Calls TMinuit to fit the energy distribution of a cluster with several maxima 
187   // The initial values for fitting procedure are set equal to the positions of local maxima.
188   // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
189
190   AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; 
191   TClonesArray * digits = gime->Digits() ; 
192   
193
194   gMinuit->mncler();                     // Reset Minuit's list of paramters
195   gMinuit->SetPrintLevel(-1) ;           // No Printout
196   gMinuit->SetFCN(AliEMCALClusterizerv1::UnfoldingChiSquare) ;  
197                                          // To set the address of the minimization function 
198   TList * toMinuit = new TList();
199   toMinuit->AddAt(emcRP,0) ;
200   toMinuit->AddAt(digits,1) ;
201   
202   gMinuit->SetObjectFit(toMinuit) ;         // To tranfer pointer to UnfoldingChiSquare
203
204   // filling initial values for fit parameters
205   AliEMCALDigit * digit ;
206
207   Int_t ierflg  = 0; 
208   Int_t index   = 0 ;
209   Int_t nDigits = (Int_t) nPar / 3 ;
210
211   Int_t iDigit ;
212
213   AliEMCALGeometry * geom = gime->EMCALGeometry() ; 
214
215   for(iDigit = 0; iDigit < nDigits; iDigit++){
216     digit = maxAt[iDigit]; 
217
218     Int_t relid[4] ;
219     Float_t x = 0.;
220     Float_t z = 0.;
221     geom->AbsToRelNumbering(digit->GetId(), relid) ;
222     geom->PosInAlice(relid, x, z) ;
223
224     Float_t energy = maxAtEnergy[iDigit] ;
225
226     gMinuit->mnparm(index, "x",  x, 0.1, 0, 0, ierflg) ;
227     index++ ;   
228     if(ierflg != 0){ 
229       Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : x = %f",  x ) ;
230       return kFALSE;
231     }
232     gMinuit->mnparm(index, "z",  z, 0.1, 0, 0, ierflg) ;
233     index++ ;   
234     if(ierflg != 0){
235        Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : z = %f", z) ;
236       return kFALSE;
237     }
238     gMinuit->mnparm(index, "Energy",  energy , 0.05*energy, 0., 4.*energy, ierflg) ;
239     index++ ;   
240     if(ierflg != 0){
241      Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : energy = %f", energy) ;      
242       return kFALSE;
243     }
244   }
245
246   Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
247                       //  depends on it. 
248   Double_t p1 = 1.0 ;
249   Double_t p2 = 0.0 ;
250
251   gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ;   // force TMinuit to reduce function calls  
252   gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ;   // force TMinuit to use my gradient  
253   gMinuit->SetMaxIterations(5);
254   gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ;  // No Warnings
255
256   gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ;    // minimize 
257
258   if(ierflg == 4){  // Minimum not found   
259     Error("FindFit", "EMCAL Unfolding  Fit not converged, cluster abandoned " ) ;      
260     return kFALSE ;
261   }            
262   for(index = 0; index < nPar; index++){
263     Double_t err ;
264     Double_t val ;
265     gMinuit->GetParameter(index, val, err) ;    // Returns value and error of parameter index
266     fitparameters[index] = val ;
267    }
268
269   delete toMinuit ;
270   return kTRUE;
271
272 }
273
274 //____________________________________________________________________________
275 void AliEMCALClusterizerv1::GetCalibrationParameters() 
276 {
277   AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
278   const AliEMCALDigitizer * dig = gime->Digitizer(BranchName()) ;
279
280   fADCchannelTower   = dig->GetTowerchannel() ;
281   fADCpedestalTower  = dig->GetTowerpedestal();
282
283   fADCchannelPreSho  = dig->GetPreShochannel() ;
284   fADCpedestalPreSho = dig->GetPreShopedestal() ; 
285 }
286
287 //____________________________________________________________________________
288 void AliEMCALClusterizerv1::Init()
289 {
290   // Make all memory allocations which can not be done in default constructor.
291   // Attach the Clusterizer task to the list of EMCAL tasks
292   
293   if ( strcmp(GetTitle(), "") == 0 )
294     SetTitle("galice.root") ;
295
296   TString branchname = GetName() ;
297   branchname.Remove(branchname.Index(Version())-1) ;
298
299   AliEMCALGetter * gime = AliEMCALGetter::GetInstance(GetTitle(), branchname.Data(), fToSplit ) ; 
300   if ( gime == 0 ) {
301     Error("Init", "Could not obtain the Getter object !" ) ; 
302     return ;
303   } 
304
305   fSplitFile = 0 ;
306   if(fToSplit){
307     // construct the name of the file as /path/EMCAL.SDigits.root
308     //First - extract full path if necessary
309     TString fileName(GetTitle()) ;
310     Ssiz_t islash = fileName.Last('/') ;
311     if(islash<fileName.Length())
312       fileName.Remove(islash+1,fileName.Length()) ;
313     else
314       fileName="" ;
315     // Next - append the file name 
316     fileName+="EMCAL.RecData." ;
317     if((strcmp(branchname.Data(),"Default")!=0)&&(strcmp(branchname.Data(),"")!=0)){
318       fileName+=branchname ;
319       fileName+="." ;
320     }
321     fileName+="root" ;
322     // Finally - check if the file already opened or open the file
323     fSplitFile = static_cast<TFile*>(gROOT->GetFile(fileName.Data()));   
324     if(!fSplitFile)
325       fSplitFile =  TFile::Open(fileName.Data(),"update") ;
326   }
327
328   const AliEMCALGeometry * geom = gime->EMCALGeometry() ;
329   fNTowers = geom->GetNZ() *  geom->GetNPhi() ;
330   if(!gMinuit) 
331     gMinuit = new TMinuit(100) ;
332
333   gime->PostClusterizer(this) ;
334   gime->PostRecPoints(branchname ) ;
335  
336 }
337
338 //____________________________________________________________________________
339 void AliEMCALClusterizerv1::InitParameters()
340 {
341   fNumberOfPreShoClusters = fNumberOfTowerClusters = 0 ;   
342   fPreShoClusteringThreshold  = 0.0001;
343   fTowerClusteringThreshold   = 0.2;    
344   fTowerLocMaxCut  = 0.03 ;
345   fPreShoLocMaxCut = 0.03 ;
346   
347   fW0     = 4.5 ;
348   fW0CPV  = 4.0 ;
349
350   fTimeGate = 1.e-8 ; 
351   
352   fToUnfold = kFALSE ;
353   
354   TString clusterizerName( GetName()) ; 
355   if (clusterizerName.IsNull() ) 
356     clusterizerName = "Default" ; 
357   clusterizerName.Append(":") ; 
358   clusterizerName.Append(Version()) ; 
359   SetName(clusterizerName) ;
360   fRecPointsInRun          = 0 ; 
361
362 }
363
364 //____________________________________________________________________________
365 Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2)const
366 {
367   // Gives the neighbourness of two digits = 0 are not neighbour but continue searching 
368   //                                       = 1 are neighbour
369   //                                       = 2 are not neighbour but do not continue searching
370   // neighbours are defined as digits having at least a common vertex 
371   // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster 
372   //                                      which is compared to a digit (d2)  not yet in a cluster  
373
374    AliEMCALGeometry * geom = AliEMCALGetter::GetInstance()->EMCALGeometry() ;
375
376   Int_t rv = 0 ; 
377
378   Int_t relid1[4] ; 
379   geom->AbsToRelNumbering(d1->GetId(), relid1) ; 
380
381   Int_t relid2[4] ; 
382   geom->AbsToRelNumbering(d2->GetId(), relid2) ; 
383  
384   if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same EMCAL Arm 
385     Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;  
386     Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;  
387     
388     if (( coldiff <= 1 )  && ( rowdiff <= 1 )){
389       if((relid1[1] != 0) || (TMath::Abs(d1->GetTime() - d2->GetTime() ) < fTimeGate))
390       rv = 1 ; 
391     }
392     else {
393       if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1)) 
394         rv = 2; //  Difference in row numbers is too large to look further 
395     }
396
397   } 
398   else {
399     
400     if( (relid1[0] < relid2[0]) || (relid1[1] != relid2[1]) )  
401       rv=2 ;
402   }
403   return rv ; 
404 }
405
406 //____________________________________________________________________________
407 Bool_t AliEMCALClusterizerv1::IsInTower(AliEMCALDigit * digit) const
408 {
409   // Tells if (true) or not (false) the digit is in a EMCAL-Tower 
410  
411   Bool_t rv = kFALSE ; 
412   if (!digit->IsInPreShower()) 
413     rv = kTRUE; 
414   return rv ; 
415 }
416
417 //____________________________________________________________________________
418 Bool_t AliEMCALClusterizerv1::IsInPreShower(AliEMCALDigit * digit) const
419 {
420   // Tells if (true) or not (false) the digit is in a EMCAL-PreShower
421  
422   Bool_t rv = kFALSE ; 
423   if (digit->IsInPreShower()) 
424     rv = kTRUE; 
425   return rv ; 
426 }
427
428 //____________________________________________________________________________
429 void AliEMCALClusterizerv1::WriteRecPoints(Int_t event)
430 {
431
432   // Creates new branches with given title
433   // fills and writes into TreeR.
434
435   AliEMCALGetter *gime = AliEMCALGetter::GetInstance() ; 
436   TObjArray * towerRecPoints = gime->TowerRecPoints() ; 
437   TObjArray * preshoRecPoints = gime->PreShowerRecPoints() ; 
438   TClonesArray * digits = gime->Digits() ; 
439   TTree * treeR ; 
440   
441   if(fToSplit){
442     if(!fSplitFile)
443       return ;
444     fSplitFile->cd() ;
445     TString name("TreeR") ;
446     name += event ; 
447     treeR = dynamic_cast<TTree*>(fSplitFile->Get(name)); 
448   }
449   else{
450     treeR = gAlice->TreeR();
451   }
452
453   if(!treeR){
454     gAlice->MakeTree("R", fSplitFile);
455     treeR = gAlice->TreeR() ;
456   }
457  
458   Int_t index ;
459   //Evaluate position, dispersion and other RecPoint properties...
460   for(index = 0; index < towerRecPoints->GetEntries(); index++)
461     (dynamic_cast<AliEMCALTowerRecPoint *>(towerRecPoints->At(index)))->EvalAll(fW0,digits) ;
462   
463   towerRecPoints->Sort() ;
464
465   for(index = 0; index < towerRecPoints->GetEntries(); index++)
466     (dynamic_cast<AliEMCALTowerRecPoint *>(towerRecPoints->At(index)))->SetIndexInList(index) ;
467
468   towerRecPoints->Expand(towerRecPoints->GetEntriesFast()) ; 
469   //Now the same for pre shower
470   for(index = 0; index < preshoRecPoints->GetEntries(); index++)
471     (dynamic_cast<AliEMCALRecPoint *>(preshoRecPoints->At(index)))->EvalAll(fW0CPV,digits)  ;
472   preshoRecPoints->Sort() ;
473
474   for(index = 0; index < preshoRecPoints->GetEntries(); index++)
475     (dynamic_cast<AliEMCALRecPoint *>(preshoRecPoints->At(index)))->SetIndexInList(index) ;
476
477   preshoRecPoints->Expand(preshoRecPoints->GetEntriesFast()) ;
478   
479   Int_t bufferSize = 32000 ;    
480   Int_t splitlevel = 0 ; 
481
482   //First Tower branch
483   TBranch * towerBranch = treeR->Branch("EMCALTowerRP","TObjArray",&towerRecPoints,bufferSize,splitlevel);
484   towerBranch->SetTitle(BranchName());
485   
486   //Now Pre Shower branch 
487   TBranch * preshoBranch = treeR->Branch("EMCALPreShoRP","TObjArray",&preshoRecPoints,bufferSize,splitlevel);
488   preshoBranch->SetTitle(BranchName());
489     
490   //And Finally  clusterizer branch
491   AliEMCALClusterizerv1 * cl = (AliEMCALClusterizerv1*)gime->Clusterizer(BranchName()) ;
492   TBranch * clusterizerBranch = treeR->Branch("AliEMCALClusterizer","AliEMCALClusterizerv1",
493                                               &cl,bufferSize,splitlevel);
494   clusterizerBranch->SetTitle(BranchName());
495
496   towerBranch        ->Fill() ;
497   preshoBranch        ->Fill() ;
498   clusterizerBranch->Fill() ;
499
500   treeR->AutoSave() ; //Write(0,kOverwrite) ;  
501   if(gAlice->TreeR()!=treeR)
502     treeR->Delete(); 
503 }
504
505 //____________________________________________________________________________
506 void AliEMCALClusterizerv1::MakeClusters()
507 {
508   // Steering method to construct the clusters stored in a list of Reconstructed Points
509   // A cluster is defined as a list of neighbour digits
510     
511   AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; 
512   
513   TObjArray * towerRecPoints  = gime->TowerRecPoints(BranchName()) ; 
514   TObjArray * preshoRecPoints = gime->PreShowerRecPoints(BranchName()) ; 
515   towerRecPoints->Delete() ;
516   preshoRecPoints->Delete() ;
517   
518   TClonesArray * digits = gime->Digits() ; 
519   if ( !digits ) {
520     Fatal("MakeClusters -> Digits with name %s not found", GetName() ) ; 
521   } 
522   TClonesArray * digitsC =  dynamic_cast<TClonesArray*>(digits->Clone()) ;
523   
524   
525   // Clusterization starts  
526   
527   TIter nextdigit(digitsC) ; 
528   AliEMCALDigit * digit ; 
529   Bool_t notremoved = kTRUE ;
530   
531   while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigit())) ) { // scan over the list of digitsC
532     AliEMCALRecPoint * clu = 0 ; 
533     
534     TArrayI clusterdigitslist(1500) ;   
535     Int_t index ;
536     if (( IsInTower (digit)  && Calibrate(digit->GetAmp(),digit->IsInPreShower()) > fTowerClusteringThreshold  ) || 
537         ( IsInPreShower (digit) && Calibrate(digit->GetAmp(),digit->IsInPreShower()) > fPreShoClusteringThreshold  ) ) {
538       
539       Int_t iDigitInCluster = 0 ; 
540       
541       if  ( IsInTower(digit) ) {   
542         // start a new Tower RecPoint
543         if(fNumberOfTowerClusters >= towerRecPoints->GetSize()) 
544           towerRecPoints->Expand(2*fNumberOfTowerClusters+1) ;
545         
546         towerRecPoints->AddAt(new  AliEMCALTowerRecPoint(""), fNumberOfTowerClusters) ;
547         clu = dynamic_cast<AliEMCALTowerRecPoint *>(towerRecPoints->At(fNumberOfTowerClusters)) ; 
548         fNumberOfTowerClusters++ ; 
549         clu->AddDigit(*digit, Calibrate(digit->GetAmp(),digit->IsInPreShower())) ; 
550         clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;  
551         iDigitInCluster++ ; 
552         digitsC->Remove(digit) ; 
553         
554       } else { 
555         
556         // start a new Pre Shower cluster
557         if(fNumberOfPreShoClusters >= preshoRecPoints->GetSize()) 
558           preshoRecPoints->Expand(2*fNumberOfPreShoClusters+1);
559         
560         preshoRecPoints->AddAt(new AliEMCALTowerRecPoint(""), fNumberOfPreShoClusters) ;
561         
562         clu =  dynamic_cast<AliEMCALTowerRecPoint *>(preshoRecPoints->At(fNumberOfPreShoClusters))  ;  
563         fNumberOfPreShoClusters++ ; 
564         clu->AddDigit(*digit, Calibrate(digit->GetAmp(),digit->IsInPreShower() ) );     
565         clusterdigitslist[iDigitInCluster] = digit->GetIndexInList()  ; 
566         iDigitInCluster++ ; 
567         digitsC->Remove(digit) ; 
568         nextdigit.Reset() ;
569         
570         // Here we remove remaining Tower digits, which cannot make a cluster
571         
572         if( notremoved ) { 
573           while( ( digit = dynamic_cast<AliEMCALDigit *>(nextdigit()) ) ) {
574             if( IsInTower(digit) )
575               digitsC->Remove(digit) ;
576             else 
577               break ; 
578           }
579           notremoved = kFALSE ;
580         }
581         
582       } // else        
583       
584       nextdigit.Reset() ;
585       
586       AliEMCALDigit * digitN ; 
587       index = 0 ;
588       while (index < iDigitInCluster){ // scan over digits already in cluster 
589         digit =  (AliEMCALDigit*)digits->At(clusterdigitslist[index])  ;      
590         index++ ; 
591         while ( (digitN = (AliEMCALDigit *)nextdigit()) ) { // scan over the reduced list of digits 
592           Int_t ineb = AreNeighbours(digit, digitN);       // call (digit,digitN) in THAT oder !!!!!
593          switch (ineb ) {
594           case 0 :   // not a neighbour
595             break ;
596           case 1 :   // are neighbours 
597             clu->AddDigit(*digitN, Calibrate( digitN->GetAmp(), digitN->IsInPreShower() ) ) ;
598             clusterdigitslist[iDigitInCluster] = digitN->GetIndexInList() ; 
599             iDigitInCluster++ ; 
600             digitsC->Remove(digitN) ;
601             break ;
602           case 2 :   // too far from each other
603             goto endofloop;   
604           } // switch
605           
606         } // while digitN
607         
608       endofloop: ;
609         nextdigit.Reset() ; 
610       } // loop over cluster     
611     } // energy theshold     
612   } // while digit  
613   delete digitsC ;
614 }
615
616 //____________________________________________________________________________
617 void AliEMCALClusterizerv1::MakeUnfolding()
618 {
619   Fatal("AliEMCALClusterizerv1::MakeUnfolding", "--> Unfolding not implemented") ;
620  
621 }
622
623 //____________________________________________________________________________
624 Double_t  AliEMCALClusterizerv1::ShowerShape(Double_t r)
625
626   // Shape of the shower (see EMCAL TDR)
627   // If you change this function, change also the gradient evaluation in ChiSquare()
628
629   Double_t r4    = r*r*r*r ;
630   Double_t r295  = TMath::Power(r, 2.95) ;
631   Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
632   return shape ;
633 }
634
635 //____________________________________________________________________________
636 void  AliEMCALClusterizerv1::UnfoldCluster(AliEMCALTowerRecPoint * iniTower, 
637                                                  Int_t nMax, 
638                                                  AliEMCALDigit ** maxAt, 
639                                                  Float_t * maxAtEnergy)
640 {
641   // Performs the unfolding of a cluster with nMax overlapping showers 
642   
643   Fatal("UnfoldCluster", "--> Unfolding not implemented") ;
644
645 }
646
647 //_____________________________________________________________________________
648 void AliEMCALClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
649 {
650   // Calculates the Chi square for the cluster unfolding minimization
651   // Number of parameters, Gradient, Chi squared, parameters, what to do
652   
653   ::Fatal("UnfoldingChiSquare","Unfolding not implemented") ;
654 }
655 //____________________________________________________________________________
656 void AliEMCALClusterizerv1::Print(Option_t * option)const
657 {
658   // Print clusterizer parameters
659
660   TString message("\n") ; 
661   
662   if( strcmp(GetName(), "") !=0 ){
663     
664     // Print parameters
665  
666     TString taskName(GetName()) ; 
667     taskName.ReplaceAll(Version(), "") ;
668     
669     message += "--------------- " ; 
670     message += taskName.Data() ; 
671     message += " " ; 
672     message += GetTitle() ; 
673     message += "-----------\n" ;  
674     message += "Clusterizing digits from the file: " ; 
675     message += taskName.Data() ;  
676     message += "\n                           Branch: " ; 
677     message += GetName() ;  
678     message += "\n                       EMC Clustering threshold = " ; 
679     message += fTowerClusteringThreshold ; 
680     message += "\n                       EMC Local Maximum cut    = " ;
681     message += fTowerLocMaxCut ; 
682     message += "\n                       EMC Logarothmic weight   = " ;
683     message += fW0 ;
684     message += "\n                       CPV Clustering threshold = " ; 
685     message += fPreShoClusteringThreshold ;
686     message += "\n                       CPV Local Maximum cut    = " ;
687     message += fPreShoLocMaxCut ;
688     message += "\n                       CPV Logarothmic weight   = " ; 
689     message += fW0CPV ;
690     if(fToUnfold)
691       message +="\nUnfolding on\n" ;
692     else
693       message += "\nUnfolding off\n";
694     
695     message += "------------------------------------------------------------------" ; 
696   }
697   else
698     message += "AliEMCALClusterizerv1 not initialized " ;
699   
700   Info("Print", message.Data() ) ; 
701 }
702
703 //____________________________________________________________________________
704 void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option)
705 {
706   // Prints list of RecPoints produced at the current pass of AliEMCALClusterizer
707
708   TObjArray * towerRecPoints = AliEMCALGetter::GetInstance()->TowerRecPoints() ; 
709   TObjArray * preshoRecPoints = AliEMCALGetter::GetInstance()->PreShowerRecPoints() ; 
710
711   TString message("\n")  ;
712
713   message += "event " ; 
714   message += gAlice->GetEvNumber() ;
715   message += "\n       Found " ; 
716   message += towerRecPoints->GetEntriesFast() ; 
717   message += " TOWER Rec Points and " ;
718   message += preshoRecPoints->GetEntriesFast() ; 
719   message += " PRE SHOWER RecPoints\n" ;
720
721   fRecPointsInRun +=  towerRecPoints->GetEntriesFast() ; 
722   fRecPointsInRun +=  preshoRecPoints->GetEntriesFast() ; 
723
724   if(strstr(option,"all")) {
725
726     message += "Tower clusters\n" ;
727     message += "Index  Ene(MeV)   Multi  Module     phi     r  theta    Lambda 1   Lambda 2  # of prim  Primaries list\n" ;      
728     
729     Int_t index ;
730     for (index = 0 ; index < towerRecPoints->GetEntries() ; index++) {
731       AliEMCALTowerRecPoint * rp = dynamic_cast<AliEMCALTowerRecPoint * >(towerRecPoints->At(index)) ; 
732       TVector3  globalpos;  
733       rp->GetGlobalPosition(globalpos);
734       Float_t lambda[2]; 
735       rp->GetElipsAxis(lambda);
736       Int_t * primaries; 
737       Int_t nprimaries;
738       primaries = rp->GetPrimaries(nprimaries);
739
740       message += rp->GetIndexInList() ; 
741       message += "   " ; 
742       message += rp->GetEnergy() ;
743       message += "      " ; 
744       message += rp->GetMultiplicity() ;
745       message += "      " ; 
746       message += rp->GetEMCALArm() ;
747       message += "     " ;
748       message += globalpos.X() ; 
749       message += "  " ;
750       message += globalpos.Y() ;  
751       message += globalpos.Z() ; 
752       message += "     " ;
753       message += lambda[0] ; 
754       message += "  " ;
755       message += lambda[1] ; 
756       message += "  " ;
757       message += nprimaries ; 
758       message += " :  " ;
759      
760       for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
761         message += primaries[iprimary] ; 
762         message += " "  ;
763       } 
764     }
765
766     //Now plot Pre shower recPoints
767
768     message += "\n-----------------------------------------------------------------------\n" ;
769
770     message += "PreShower clusters\n" ;
771     message += " Index  Ene(MeV)   Multi  Module     phi     r  theta    Lambda 1   Lambda 2  # of prim  Primaries list\n" ;      
772     
773     for (index = 0 ; index < preshoRecPoints->GetEntries() ; index++) {
774       AliEMCALTowerRecPoint * rp = dynamic_cast<AliEMCALTowerRecPoint *>(preshoRecPoints->At(index)) ; 
775       TVector3  globalpos;  
776       rp->GetGlobalPosition(globalpos);
777       Float_t lambda[2]; 
778       rp->GetElipsAxis(lambda);
779       Int_t * primaries;
780       Int_t nprimaries;
781       primaries = rp->GetPrimaries(nprimaries);
782
783       message += rp->GetIndexInList() ; 
784       message += "   " ; 
785       message += rp->GetEnergy() ; 
786       message += "   " ; 
787       message += rp->GetMultiplicity() ;
788       message += "   " ; 
789       message += rp->GetEMCALArm() ;
790       message += "   " ; 
791       message += globalpos.X() ;
792       message += "   " ; 
793       message += globalpos.Y() ;
794       message += "   " ; 
795       message += globalpos.Z() ;
796       message += "   " ; 
797       message += lambda[0] ;
798       message += "   " ; 
799       message += lambda[1] ;
800       message += "   " ; 
801       message += nprimaries ;
802       message += " : " ; 
803     
804       for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
805         message += primaries[iprimary] ; 
806         message += "  " ;
807       }          
808     }
809
810     message += "\n-----------------------------------------------------------------------" ;
811   }
812
813   Info("PrintRecPoints", message.Data() ) ; 
814   
815 }