Coding conventions (Gustavo)
[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
18 //*-- Author: Yves Schutz (SUBATECH)  & Dmitri Peressounko (SUBATECH & Kurchatov Institute)
19 //  August 2002 Yves Schutz: clone PHOS as closely as possible and intoduction
20 //                           of new  IO (à la PHOS)
21 //////////////////////////////////////////////////////////////////////////////
22 //  Clusterization class. Performs clusterization (collects neighbouring active cells) and 
23 //  unfolds the clusters having several local maxima.  
24 //  Results are stored in TreeR#, branches EMCALTowerRP (EMC recPoints),
25 //  EMCALPreShoRP (CPV RecPoints) and AliEMCALClusterizer (Clusterizer with all 
26 //  parameters including input digits branch title, thresholds etc.)
27 //  This TTask is normally called from Reconstructioner, but can as well be used in 
28 //  standalone mode.
29 // Use Case:
30 //  root [0] AliEMCALClusterizerv1 * cl = new AliEMCALClusterizerv1("galice.root")  
31 //  Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
32 //               //reads gAlice from header file "..."                      
33 //  root [1] cl->ExecuteTask()  
34 //               //finds RecPoints in all events stored in galice.root
35 //  root [2] cl->SetDigitsBranch("digits2") 
36 //               //sets another title for Digitis (input) branch
37 //  root [3] cl->SetRecPointsBranch("recp2")  
38 //               //sets another title four output branches
39 //  root [4] cl->SetTowerLocalMaxCut(0.03)  
40 //               //set clusterization parameters
41 //  root [5] cl->ExecuteTask("deb all time")  
42 //               //once more finds RecPoints options are 
43 //               // deb - print number of found rec points
44 //               // deb all - print number of found RecPoints and some their characteristics 
45 //               // time - print benchmarking results
46
47 // --- ROOT system ---
48
49 class TROOT;
50 #include <TH1.h>
51 #include <TFile.h> 
52 class TFolder;
53 #include <TMath.h> 
54 #include <TMinuit.h>
55 #include <TTree.h> 
56 class TSystem; 
57 #include <TBenchmark.h>
58 #include <TBrowser.h>
59
60 // --- Standard library ---
61
62
63 // --- AliRoot header files ---
64 #include "AliRunLoader.h"
65 #include "AliRun.h"
66 #include "AliESD.h"
67 #include "AliEMCALLoader.h"
68 #include "AliEMCALClusterizerv1.h"
69 #include "AliEMCALRecPoint.h"
70 #include "AliEMCALDigit.h"
71 #include "AliEMCALDigitizer.h"
72 #include "AliEMCAL.h"
73 #include "AliEMCALGeometry.h"
74 #include "AliEMCALHistoUtilities.h"
75 #include "AliCDBManager.h"
76
77 class AliCDBStorage;
78 #include "AliCDBEntry.h"
79
80 ClassImp(AliEMCALClusterizerv1)
81
82 //____________________________________________________________________________
83   AliEMCALClusterizerv1::AliEMCALClusterizerv1() : AliEMCALClusterizer()
84 {
85   // default ctor (to be used mainly by Streamer)
86   
87   InitParameters() ; 
88   fDefaultInit = kTRUE ;
89   fGeom = AliEMCALGeometry::GetInstance();
90   fGeom->GetTransformationForSM(); // Global <-> Local
91 }
92
93 //____________________________________________________________________________
94 AliEMCALClusterizerv1::AliEMCALClusterizerv1(const TString alirunFileName, const TString eventFolderName)
95 :AliEMCALClusterizer(alirunFileName, eventFolderName)
96 {
97   // ctor with the indication of the file where header Tree and digits Tree are stored
98   
99   InitParameters() ; 
100   Init() ;
101   fDefaultInit = kFALSE ; 
102 }
103
104 //____________________________________________________________________________
105   AliEMCALClusterizerv1::~AliEMCALClusterizerv1()
106 {
107   // dtor
108 }
109
110 //____________________________________________________________________________
111 const TString AliEMCALClusterizerv1::BranchName() const 
112
113    return GetName();
114 }
115
116 //____________________________________________________________________________
117 Float_t  AliEMCALClusterizerv1::Calibrate(Int_t amp, Int_t AbsId) 
118 {
119  
120   // Convert digitized amplitude into energy.
121   // Calibration parameters are taken from calibration data base for raw data,
122   // or from digitizer parameters for simulated data.
123
124   if(fCalibData){
125
126     //JLK 13-Mar-2006
127     //We now get geometry at a higher level
128     //
129     // Loader
130     //    AliRunLoader *rl = AliRunLoader::GetRunLoader();
131     
132     // Load EMCAL Geomtry
133     //    rl->LoadgAlice(); 
134     //AliRun   * gAlice = rl->GetAliRun(); 
135     //AliEMCAL * emcal  = (AliEMCAL*)gAlice->GetDetector("EMCAL");
136     //AliEMCALGeometry * geom = emcal->GetGeometry();
137     
138     if (fGeom==0)
139       AliFatal("Did not get geometry from EMCALLoader") ;
140     
141     Int_t iSupMod = -1;
142     Int_t nTower  = -1;
143     Int_t nIphi   = -1;
144     Int_t nIeta   = -1;
145     Int_t iphi    = -1;
146     Int_t ieta    = -1;
147     
148     Bool_t bCell = fGeom->GetCellIndex(AbsId, iSupMod, nTower, nIphi, nIeta) ;
149     if(!bCell)
150       Error("DigitizeEnergy","Wrong cell id number") ;
151     fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nTower,nIphi, nIeta,iphi,ieta);
152
153     fADCchannelECA  = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
154     fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
155     return -fADCpedestalECA + amp * fADCchannelECA ;        
156  
157   }
158   else //Return energy with default parameters if calibration is not available
159     return -fADCpedestalECA + amp * fADCchannelECA ; 
160   
161 }
162
163 //____________________________________________________________________________
164 void AliEMCALClusterizerv1::Exec(Option_t * option)
165 {
166   // Steering method to perform clusterization for events
167   // in the range from fFirstEvent to fLastEvent.
168   // This range is optionally set by SetEventRange().
169   // if fLastEvent=-1 (by default), then process events until the end.
170
171   if(strstr(option,"tim"))
172     gBenchmark->Start("EMCALClusterizer"); 
173   
174   if(strstr(option,"print"))
175     Print("") ; 
176
177   AliRunLoader *rl = AliRunLoader::GetRunLoader();
178   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
179
180   //Get calibration parameters from file or digitizer default values.
181   GetCalibrationParameters() ;
182
183   if (fLastEvent == -1) 
184     fLastEvent = rl->GetNumberOfEvents() - 1;
185   Int_t nEvents   = fLastEvent - fFirstEvent + 1;
186
187   Int_t ievent ;
188   rl->LoadDigits("EMCAL");
189   for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
190     rl->GetEvent(ievent);
191
192     fNumberOfECAClusters = 0;
193
194     if(strstr(option,"pseudo"))
195       MakeClusters("pseudo") ;  //both types
196     else
197       MakeClusters("") ;  //only the real clusters
198
199     if(fToUnfold)
200       MakeUnfolding() ;
201
202     WriteRecPoints() ;
203
204     if(strstr(option,"deb"))  
205       PrintRecPoints(option) ;
206
207     //increment the total number of recpoints per run   
208     fRecPointsInRun += emcalLoader->RecPoints()->GetEntriesFast() ;  
209   }
210   
211   Unload();
212
213   if(strstr(option,"tim")){
214     gBenchmark->Stop("EMCALClusterizer");
215     printf("Exec took %f seconds for Clusterizing %f seconds per event", 
216          gBenchmark->GetCpuTime("EMCALClusterizer"), gBenchmark->GetCpuTime("EMCALClusterizer")/nEvents );
217   }
218
219 }
220
221 //____________________________________________________________________________
222 Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALRecPoint * emcRP, AliEMCALDigit ** maxAt, Float_t * maxAtEnergy,
223                                     Int_t nPar, Float_t * fitparameters) const
224
225   // Calls TMinuit to fit the energy distribution of a cluster with several maxima 
226   // The initial values for fitting procedure are set equal to the positions of local maxima.
227   // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
228
229   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
230   TClonesArray *digits = emcalLoader->Digits();
231
232   gMinuit->mncler();                     // Reset Minuit's list of paramters
233   gMinuit->SetPrintLevel(-1) ;           // No Printout
234   gMinuit->SetFCN(AliEMCALClusterizerv1::UnfoldingChiSquare) ;  
235                                          // To set the address of the minimization function 
236   TList * toMinuit = new TList();
237   toMinuit->AddAt(emcRP,0) ;
238   toMinuit->AddAt(digits,1) ;
239   
240   gMinuit->SetObjectFit(toMinuit) ;         // To tranfer pointer to UnfoldingChiSquare
241
242   // filling initial values for fit parameters
243   AliEMCALDigit * digit ;
244
245   Int_t ierflg  = 0; 
246   Int_t index   = 0 ;
247   Int_t nDigits = (Int_t) nPar / 3 ;
248
249   Int_t iDigit ;
250
251   for(iDigit = 0; iDigit < nDigits; iDigit++){
252     digit = maxAt[iDigit]; 
253
254     Int_t relid[2] ;
255     Float_t x = 0.;
256     Float_t z = 0.;
257     fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
258     fGeom->PosInAlice(relid, x, z) ;
259
260     Float_t energy = maxAtEnergy[iDigit] ;
261
262     gMinuit->mnparm(index, "x",  x, 0.1, 0, 0, ierflg) ;
263     index++ ;   
264     if(ierflg != 0){ 
265       Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : x = %f",  x ) ;
266       return kFALSE;
267     }
268     gMinuit->mnparm(index, "z",  z, 0.1, 0, 0, ierflg) ;
269     index++ ;   
270     if(ierflg != 0){
271        Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : z = %f", z) ;
272       return kFALSE;
273     }
274     gMinuit->mnparm(index, "Energy",  energy , 0.05*energy, 0., 4.*energy, ierflg) ;
275     index++ ;   
276     if(ierflg != 0){
277      Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : energy = %f", energy) ;      
278       return kFALSE;
279     }
280   }
281
282   Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
283                       //  depends on it. 
284   Double_t p1 = 1.0 ;
285   Double_t p2 = 0.0 ;
286
287   gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ;   // force TMinuit to reduce function calls  
288   gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ;   // force TMinuit to use my gradient  
289   gMinuit->SetMaxIterations(5);
290   gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ;  // No Warnings
291
292   gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ;    // minimize 
293
294   if(ierflg == 4){  // Minimum not found   
295     Error("FindFit", "EMCAL Unfolding  Fit not converged, cluster abandoned " ) ;      
296     return kFALSE ;
297   }            
298   for(index = 0; index < nPar; index++){
299     Double_t err ;
300     Double_t val ;
301     gMinuit->GetParameter(index, val, err) ;    // Returns value and error of parameter index
302     fitparameters[index] = val ;
303    }
304
305   delete toMinuit ;
306   return kTRUE;
307
308 }
309
310 //____________________________________________________________________________
311 void AliEMCALClusterizerv1::GetCalibrationParameters() 
312 {
313   // Set calibration parameters:
314   // if calibration database exists, they are read from database,
315   // otherwise, they are taken from digitizer.
316   //
317   // It is a user responsilibity to open CDB before reconstruction, 
318   // for example: 
319   // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
320
321   //Check if calibration is stored in data base
322    if(AliCDBManager::Instance()->IsDefaultStorageSet()){
323      AliCDBEntry *entry = (AliCDBEntry*) 
324      AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
325      if (entry) fCalibData = (AliEMCALCalibData*) entry->GetObject();
326    }
327    if(!fCalibData)
328      {
329        //If calibration is not available use default parameters
330        //Loader
331        AliEMCALLoader *emcalLoader = 
332          dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
333        if ( !emcalLoader->Digitizer() ) 
334          emcalLoader->LoadDigitizer();
335        AliEMCALDigitizer * dig = dynamic_cast<AliEMCALDigitizer*>(emcalLoader->Digitizer());
336        
337        fADCchannelECA   = dig->GetECAchannel() ;
338        fADCpedestalECA  = dig->GetECApedestal();
339      }
340 }
341
342 //____________________________________________________________________________
343 void AliEMCALClusterizerv1::Init()
344 {
345   // Make all memory allocations which can not be done in default constructor.
346   // Attach the Clusterizer task to the list of EMCAL tasks
347   
348   AliRunLoader *rl = AliRunLoader::GetRunLoader();
349   fGeom = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
350   fGeom->GetTransformationForSM(); // Global <-> Local
351   AliInfo(Form("geom 0x%x",fGeom));
352
353   if(!gMinuit) 
354     gMinuit = new TMinuit(100) ;
355
356   fHists = BookHists();
357 }
358
359 //____________________________________________________________________________
360 void AliEMCALClusterizerv1::InitParameters()
361
362   // Initializes the parameters for the Clusterizer
363   fNumberOfECAClusters = 0;
364
365   fNTowerInGroup = 36;  //Produces maximum of 80 pseudoclusters per event
366
367   fECAClusteringThreshold   = 0.5;  // value obtained from Alexei
368   fECALocMaxCut = 0.03; // ??
369
370   fECAW0     = 4.5 ;
371   fTimeGate = 1.e-8 ; 
372   fToUnfold = kFALSE ;
373   fRecPointsInRun  = 0 ;
374   fMinECut = 0.01; // have to be tune
375
376   fCalibData               = 0 ;
377 }
378
379 //____________________________________________________________________________
380 Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2) const
381 {
382   // Gives the neighbourness of two digits = 0 are not neighbour but continue searching 
383   //                                       = 1 are neighbour
384   //                                       = 2 is in different SM, quit from searching
385   // neighbours are defined as digits having at least a common vertex 
386   // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster 
387   //                                      which is compared to a digit (d2)  not yet in a cluster  
388
389   static Int_t rv; 
390   static Int_t nSupMod1=0, nTower1=0, nIphi1=0, nIeta1=0, iphi1=0, ieta1=0;
391   static Int_t nSupMod2=0, nTower2=0, nIphi2=0, nIeta2=0, iphi2=0, ieta2=0;
392   static Int_t rowdiff, coldiff;
393   rv = 0 ; 
394
395   fGeom->GetCellIndex(d1->GetId(), nSupMod1,nTower1,nIphi1,nIeta1);
396   fGeom->GetCellIndex(d2->GetId(), nSupMod2,nTower2,nIphi2,nIeta2);
397   if(nSupMod1 != nSupMod2) return 2; // different SM
398
399   fGeom->GetCellPhiEtaIndexInSModule(nSupMod1,nTower1,nIphi1,nIeta1, iphi1,ieta1);
400   fGeom->GetCellPhiEtaIndexInSModule(nSupMod2,nTower2,nIphi2,nIeta2, iphi2,ieta2);
401
402   rowdiff = TMath::Abs(iphi1 - iphi2);  
403   coldiff = TMath::Abs(ieta1 - ieta2) ;  
404   
405   if (( coldiff <= 1 )  && ( rowdiff <= 1 )) rv = 1;  // neighbours with at least commom vertex
406  
407   if (gDebug == 2 && rv==1) 
408   printf("AreNeighbours: neighbours=%d, id1=%d, relid1=%d,%d \n id2=%d, relid2=%d,%d \n", 
409          rv, d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2);   
410   
411   return rv ; 
412 }
413
414 //____________________________________________________________________________
415 Int_t AliEMCALClusterizerv1::AreInGroup(AliEMCALDigit * d1, AliEMCALDigit * d2) const
416 {
417   // Tells whether two digits fall within the same supermodule and
418   // tower grouping.  The number of towers in a group is controlled by
419   // the parameter nTowersInGroup
420   //    = 0 are not in same group but continue searching 
421   //    = 1 same group
422   //    = 2 is in different SM, quit from searching
423   //    = 3 different tower group, quit from searching
424   //
425   // The order of d1 and d2 is important: first (d1) should be a digit 
426   // already in a cluster which is compared to a digit (d2)  not yet in a cluster  
427
428   //JLK Question: does the quit from searching assume that the digits
429   //are ordered, so that once you are in a different SM, you'll not
430   //find another in the list that will match?  How about my TowerGroup search?
431
432   static Int_t rv;
433   static Int_t nSupMod1=0, nTower1=0, nIphi1=0, nIeta1=0, iphi1=0, ieta1=0;
434   static Int_t nSupMod2=0, nTower2=0, nIphi2=0, nIeta2=0, iphi2=0, ieta2=0;
435   static Int_t towerGroup1 = -1, towerGroup2 = -1;
436   rv = 0 ;
437
438   fGeom->GetCellIndex(d1->GetId(), nSupMod1,nTower1,nIphi1,nIeta1);
439   fGeom->GetCellIndex(d2->GetId(), nSupMod2,nTower2,nIphi2,nIeta2);
440   if(nSupMod1 != nSupMod2) return 2; // different SM
441
442   static Int_t nTowerInSM = fGeom->GetNCellsInSupMod()/fGeom->GetNCellsInTower();
443
444   //figure out which tower grouping each digit belongs to
445   for(int it = 0; it < nTowerInSM/fNTowerInGroup; it++) {
446     if(nTower1 <= nTowerInSM - it*fNTowerInGroup) towerGroup1 = it;
447     if(nTower2 <= nTowerInSM - it*fNTowerInGroup) towerGroup2 = it;
448   }
449   if(towerGroup1 != towerGroup2) return 3; //different Towergroup
450
451   //same SM, same towergroup, we're happy
452   if(towerGroup1 == towerGroup2 && towerGroup2 >= 0)
453     rv = 1;
454
455   if (gDebug == 2 && rv==1)
456   printf("AreInGroup: neighbours=%d, id1=%d, relid1=%d,%d \n id2=%d, relid2=%d,%d \n",
457          rv, d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2);
458
459   return rv ;
460 }
461
462 //____________________________________________________________________________
463 void AliEMCALClusterizerv1::Unload() 
464 {
465   // Unloads the Digits and RecPoints
466   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
467     
468   emcalLoader->UnloadDigits() ; 
469   emcalLoader->UnloadRecPoints() ; 
470 }
471  
472 //____________________________________________________________________________
473 void AliEMCALClusterizerv1::WriteRecPoints()
474 {
475
476   // Creates new branches with given title
477   // fills and writes into TreeR.
478
479   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
480
481   TObjArray * aECARecPoints = emcalLoader->RecPoints() ; 
482
483   TClonesArray * digits = emcalLoader->Digits() ; 
484   TTree * treeR = emcalLoader->TreeR();  
485   if ( treeR==0 ) {
486     emcalLoader->MakeRecPointsContainer();
487     treeR = emcalLoader->TreeR();  
488   }
489   Int_t index ;
490
491   //Evaluate position, dispersion and other RecPoint properties for EC section
492   for(index = 0; index < aECARecPoints->GetEntries(); index++)
493     (dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(index)))->EvalAll(fECAW0,digits) ;
494   
495   aECARecPoints->Sort() ;
496
497   for(index = 0; index < aECARecPoints->GetEntries(); index++) {
498     (dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(index)))->SetIndexInList(index) ;
499     (dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(index)))->Print();
500   }
501
502   Int_t bufferSize = 32000 ;    
503   Int_t splitlevel = 0 ; 
504
505   //EC section branch
506   
507   TBranch * branchECA = 0;
508   if ((branchECA = treeR->GetBranch("EMCALECARP")))
509     branchECA->SetAddress(&aECARecPoints);
510   else
511     treeR->Branch("EMCALECARP","TObjArray",&aECARecPoints,bufferSize,splitlevel);
512   treeR->Fill() ;
513
514   emcalLoader->WriteRecPoints("OVERWRITE");
515
516 }
517
518 //____________________________________________________________________________
519 void AliEMCALClusterizerv1::MakeClusters(char* option)
520 {
521   // Steering method to construct the clusters stored in a list of Reconstructed Points
522   // A cluster is defined as a list of neighbour digits
523
524   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
525
526   TObjArray * aECARecPoints = emcalLoader->RecPoints() ; 
527     
528   if (fGeom==0) 
529      AliFatal("Did not get geometry from EMCALLoader");
530
531   aECARecPoints->Clear();
532
533   //Start with pseudoclusters, if option
534   if(strstr(option,"pseudo")) {
535     TClonesArray * digits  = emcalLoader->Digits() ; 
536     TClonesArray * digitsC =  dynamic_cast<TClonesArray*>(digits->Clone());
537
538     TIter nextdigit(digitsC) ;
539     AliEMCALDigit * digit;
540
541     AliEMCALRecPoint * recPoint = 0 ; 
542     int ineb=0;
543     nextdigit.Reset();
544
545     // PseudoClusterization starts    
546     while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigit())) ) { // scan over the list of digitsC
547       recPoint = 0 ; 
548       TArrayI clusterECAdigitslist(1000); // what is this   
549
550       if(fGeom->CheckAbsCellId(digit->GetId()) ) { //Is this an EMCAL digit? Just maing sure...
551         Int_t iDigitInECACluster = 0; // start a new Tower RecPoint
552         if(fNumberOfECAClusters >= aECARecPoints->GetSize()) aECARecPoints->Expand(2*fNumberOfECAClusters+1) ;
553         AliEMCALRecPoint * rp = new  AliEMCALRecPoint("") ; 
554         aECARecPoints->AddAt(rp, fNumberOfECAClusters) ;
555         recPoint = dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(fNumberOfECAClusters)) ; 
556         fNumberOfECAClusters++ ; 
557
558         recPoint->SetClusterType(AliESDCaloCluster::kPseudoCluster);
559
560         recPoint->AddDigit(*digit, digit->GetAmp()) ; 
561         clusterECAdigitslist[iDigitInECACluster] = digit->GetIndexInList() ;    
562         iDigitInECACluster++ ; 
563         digitsC->Remove(digit) ; 
564         AliDebug(1,Form("MakePseudoClusters: OK id = %d, adc = %f \n", digit->GetId(), digit->GetAmp()));
565         nextdigit.Reset(); // will start from beggining
566       
567         AliEMCALDigit * digitN = 0; // digi neighbor
568         Int_t index = 0 ;
569
570         // Find the neighbours
571         while (index < iDigitInECACluster){ // scan over digits already in cluster 
572           digit =  (AliEMCALDigit*)digits->At(clusterECAdigitslist[index]);
573           index++ ; 
574           while ( (digitN = (AliEMCALDigit *)nextdigit())) { // scan over the reduced list of digits 
575             ineb = AreInGroup(digit, digitN);       // call (digit,digitN) in THAT oder !!!!! 
576             switch (ineb ) {
577               case 0 :   // not a neighbour
578                break ;
579               case 1 :   // are neighbours 
580                recPoint->AddDigit(*digitN, digitN->GetAmp() ) ;
581                clusterECAdigitslist[iDigitInECACluster] = digitN->GetIndexInList() ; 
582                iDigitInECACluster++ ; 
583                digitsC->Remove(digitN) ;
584                break ;
585               case 2 :   // different SM
586                break ;
587               case 3 :   // different Tower group
588                break ;
589             } // switch
590           } // scan over the reduced list of digits
591         } // scan over digits already in cluster
592         nextdigit.Reset() ;  // will start from beggining
593       }
594     }
595     if(recPoint)
596       AliDebug(1,Form("MakeClusters: cl.e %f \n", recPoint->GetEnergy())); 
597     //if(recPoint) cout << "cl.e " << recPoint->GetEnergy() << endl; 
598     delete digitsC ;
599   }
600
601   //Now do real clusters
602   TClonesArray * digits  = emcalLoader->Digits() ; 
603   TClonesArray * digitsC =  dynamic_cast<TClonesArray*>(digits->Clone());
604
605   TIter nextdigit(digitsC) ;
606   AliEMCALDigit * digit;
607
608   AliEMCALRecPoint * recPoint = 0 ; 
609   int ineb=0;
610   nextdigit.Reset();
611
612   double e=0.0, ehs = 0.0;
613   while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigit())) ) { // clean up digits
614     e = Calibrate(digit->GetAmp(), digit->GetId());
615     AliEMCALHistoUtilities::FillH1(fHists, 10, digit->GetAmp());
616     AliEMCALHistoUtilities::FillH1(fHists, 11, e);
617     if(e < fMinECut ) digitsC->Remove(digit);
618     else              ehs += e;
619   } 
620   AliDebug(1,Form("MakeClusters: Number of digits %d  -> (e %f), ehs %d\n",
621                   digits->GetEntries(),fMinECut,ehs));
622   //cout << " Number of digits " << digits->GetEntries() << " -> (e>" <<fMinECut <<")";
623   //cout << digitsC->GetEntries()<< " ehs "<<ehs<<endl; 
624
625   // Clusterization starts    
626   //  cout << "Outer Loop" << endl;
627   ineb=0;
628   nextdigit.Reset();
629   recPoint = 0 ; 
630   while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigit())) ) { // scan over the list of digitsC
631     recPoint = 0 ; 
632     TArrayI clusterECAdigitslist(1000); // what is this   
633
634     if(fGeom->CheckAbsCellId(digit->GetId()) && (Calibrate(digit->GetAmp(), digit->GetId()) > fECAClusteringThreshold  ) ){
635       Int_t iDigitInECACluster = 0; // start a new Tower RecPoint
636       if(fNumberOfECAClusters >= aECARecPoints->GetSize()) aECARecPoints->Expand(2*fNumberOfECAClusters+1) ;
637       AliEMCALRecPoint * rp = new  AliEMCALRecPoint("") ; 
638       aECARecPoints->AddAt(rp, fNumberOfECAClusters) ;
639       recPoint = dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(fNumberOfECAClusters)) ; 
640       fNumberOfECAClusters++ ; 
641
642       recPoint->SetClusterType(AliESDCaloCluster::kClusterv1);
643
644       recPoint->AddDigit(*digit, Calibrate(digit->GetAmp(), digit->GetId())) ; 
645       clusterECAdigitslist[iDigitInECACluster] = digit->GetIndexInList() ;      
646       iDigitInECACluster++ ; 
647       digitsC->Remove(digit) ; 
648       AliDebug(1,Form("MakeClusters: OK id = %d, ene = %f , thre = %f \n", digit->GetId(),
649       Calibrate(digit->GetAmp(),digit->GetId()), fECAClusteringThreshold));  
650       nextdigit.Reset(); // will start from beggining
651       
652       AliEMCALDigit * digitN = 0; // digi neighbor
653       Int_t index = 0 ;
654
655       // Find the neighbours
656       while (index < iDigitInECACluster){ // scan over digits already in cluster 
657         digit =  (AliEMCALDigit*)digits->At(clusterECAdigitslist[index]);
658         index++ ; 
659         while ( (digitN = (AliEMCALDigit *)nextdigit())) { // scan over the reduced list of digits 
660           ineb = AreNeighbours(digit, digitN);       // call (digit,digitN) in THAT oder !!!!! 
661           switch (ineb ) {
662             case 0 :   // not a neighbour
663               break ;
664             case 1 :   // are neighbours 
665               recPoint->AddDigit(*digitN, Calibrate(digitN->GetAmp(),digitN->GetId()) ) ;
666               clusterECAdigitslist[iDigitInECACluster] = digitN->GetIndexInList() ; 
667               iDigitInECACluster++ ; 
668               digitsC->Remove(digitN) ;
669                break ;
670              case 2 :   // different SM
671                break ;
672           } // switch
673         } // scan over the reduced list of digits
674       } // scan over digits already in cluster
675       nextdigit.Reset() ;  // will start from beggining
676     }
677   } // while digit 
678   if(recPoint)
679     AliDebug(1,Form("MakeClusters: cl.e %f \n", recPoint->GetEnergy())); 
680   //if(recPoint) cout << "cl.e " << recPoint->GetEnergy() << endl; 
681   delete digitsC ;
682
683   AliDebug(1,Form("total no of clusters %d from %d digits",fNumberOfECAClusters,digits->GetEntriesFast())); 
684 }
685
686 //____________________________________________________________________________
687 void AliEMCALClusterizerv1::MakeUnfolding() const
688 {
689   Fatal("AliEMCALClusterizerv1::MakeUnfolding", "--> Unfolding not implemented") ;
690  
691 }
692
693 //____________________________________________________________________________
694 Double_t  AliEMCALClusterizerv1::ShowerShape(Double_t r)
695
696   // Shape of the shower (see EMCAL TDR)
697   // If you change this function, change also the gradient evaluation in ChiSquare()
698
699   Double_t r4    = r*r*r*r ;
700   Double_t r295  = TMath::Power(r, 2.95) ;
701   Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
702   return shape ;
703 }
704
705 //____________________________________________________________________________
706 void  AliEMCALClusterizerv1::UnfoldCluster(AliEMCALRecPoint * /*iniTower*/, 
707                                            Int_t /*nMax*/, 
708                                            AliEMCALDigit ** /*maxAt*/, 
709                                            Float_t * /*maxAtEnergy*/) const
710 {
711   // Performs the unfolding of a cluster with nMax overlapping showers 
712   
713   Fatal("UnfoldCluster", "--> Unfolding not implemented") ;
714
715 }
716
717 //_____________________________________________________________________________
718 void AliEMCALClusterizerv1::UnfoldingChiSquare(Int_t & /*nPar*/, Double_t * /*Grad*/,
719                                                Double_t & /*fret*/,
720                                                Double_t * /*x*/, Int_t /*iflag*/)
721 {
722   // Calculates the Chi square for the cluster unfolding minimization
723   // Number of parameters, Gradient, Chi squared, parameters, what to do
724   
725   ::Fatal("UnfoldingChiSquare","Unfolding not implemented") ;
726 }
727 //____________________________________________________________________________
728 void AliEMCALClusterizerv1::Print(Option_t * /*option*/)const
729 {
730   // Print clusterizer parameters
731
732   TString message("\n") ; 
733   
734   if( strcmp(GetName(), "") !=0 ){
735     
736     // Print parameters
737  
738     TString taskName(GetName()) ;
739     taskName.ReplaceAll(Version(), "") ;
740     
741     printf("--------------- "); 
742     printf(taskName.Data()) ; 
743     printf(" "); 
744     printf(GetTitle()) ; 
745     printf("-----------\n");  
746     printf("Clusterizing digits from the file: "); 
747     printf(taskName.Data());  
748     printf("\n                           Branch: "); 
749     printf(GetName()); 
750     printf("\n                       ECA Local Maximum cut    = %f", fECALocMaxCut); 
751     printf("\n                       ECA Logarithmic weight   = %f", fECAW0); 
752     if(fToUnfold)
753       printf("\nUnfolding on\n");
754     else
755       printf("\nUnfolding off\n");
756     
757     printf("------------------------------------------------------------------"); 
758   }
759   else
760     printf("AliEMCALClusterizerv1 not initialized ") ;
761 }
762
763 //____________________________________________________________________________
764 void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option)
765 {
766   // Prints list of RecPoints produced at the current pass of AliEMCALClusterizer
767   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
768   TObjArray * aECARecPoints = emcalLoader->RecPoints() ; 
769     
770   printf("PrintRecPoints: Clusterization result:") ; 
771   
772   printf("event # %d\n", emcalLoader->GetRunLoader()->GetEventNumber() ) ;
773   printf("           Found %d ECA Rec Points\n ", 
774          aECARecPoints->GetEntriesFast()) ; 
775
776   fRecPointsInRun +=  aECARecPoints->GetEntriesFast() ; 
777   
778   if(strstr(option,"all")) {
779     Int_t index =0;
780     printf("\n-----------------------------------------------------------------------\n") ;
781     printf("Clusters in ECAL section\n") ;
782     printf("Index    Ene(GeV) Multi Module     phi     r   theta    X    Y      Z   Dispersion Lambda 1   Lambda 2  # of prim  Primaries list\n") ;      
783    Float_t maxE=0; 
784    Float_t maxL1=0; 
785    Float_t maxL2=0; 
786    Float_t maxDis=0; 
787
788     for (index = 0 ; index < aECARecPoints->GetEntries() ; index++) {
789       AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint * >(aECARecPoints->At(index)) ; 
790       TVector3  globalpos;  
791       //rp->GetGlobalPosition(globalpos);
792       TVector3  localpos;  
793       rp->GetLocalPosition(localpos);
794       Float_t lambda[2]; 
795       rp->GetElipsAxis(lambda);
796       Int_t * primaries; 
797       Int_t nprimaries;
798       primaries = rp->GetPrimaries(nprimaries);
799       printf("\n%6d  %8.4f  %3d     %4.1f    %4.1f %4.1f  %4.1f %4.1f %4.1f    %4.1f   %4f  %4f    %2d     : ", 
800              rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(),
801              globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(), 
802              rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ; 
803   /////////////
804       if(rp->GetEnergy()>maxE){
805               maxE=rp->GetEnergy();
806               maxL1=lambda[0];
807               maxL2=lambda[1];
808               maxDis=rp->GetDispersion();
809       }
810       fPointE->Fill(rp->GetEnergy());
811       fPointL1->Fill(lambda[0]);
812       fPointL2->Fill(lambda[1]);
813       fPointDis->Fill(rp->GetDispersion());
814       fPointMult->Fill(rp->GetMultiplicity());
815       ///////////// 
816       for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
817         printf("%d ", primaries[iprimary] ) ; 
818       } 
819     }
820
821       fMaxE->Fill(maxE);
822       fMaxL1->Fill(maxL1);
823       fMaxL2->Fill(maxL2);
824       fMaxDis->Fill(maxDis);
825
826
827     printf("\n-----------------------------------------------------------------------\n");
828   }
829 }
830 TList* AliEMCALClusterizerv1::BookHists()
831 {
832   //set up histograms for monitoring clusterizer performance
833
834   gROOT->cd();
835
836         fPointE = new TH1F("pointE","point energy", 2000, 0.0, 150.);
837         fPointL1 = new TH1F("pointL1","point L1", 1000, 0.0, 3.);
838         fPointL2 = new TH1F("pointL2","point L2", 1000, 0.0, 3.);
839         fPointDis = new TH1F("pointDis","point Dis", 1000, 0.0, 3.);
840         fPointMult = new TH1F("pointMult","point Mult", 100, 0.0, 100.);
841         fDigitAmp = new TH1F("digitAmp","Digit Amplitude", 2000, 0.0, 5000.);
842         fMaxE = new TH1F("maxE","Max point energy", 2000, 0.0, 150.);
843         fMaxL1 = new TH1F("maxL1","Max point L1", 1000, 0.0, 3.);
844         fMaxL2 = new TH1F("maxL2","Max point L2", 1000, 0.0, 3.);
845         fMaxDis = new TH1F("maxDis","Max point Dis", 1000, 0.0, 3.); // 9
846         //
847         new TH1F("adcOfDigits","adc of digits(threshold control)", 1001, -0.5, 1000.5);   // 10
848         new TH1F("energyOfDigits","energy of digits(threshold control)", 1000, 0.0, 1.);  // 11
849
850   return AliEMCALHistoUtilities::MoveHistsToList("EmcalClusterizerv1ControlHists", kFALSE);
851 }
852
853 void AliEMCALClusterizerv1::SaveHists(const char *fn)
854 {
855   AliEMCALHistoUtilities::SaveListOfHists(fHists, fn, kTRUE);
856 }
857
858 void AliEMCALClusterizerv1::Browse(TBrowser* b)
859 {
860   if(fHists) b->Add(fHists);
861   TTask::Browse(b);
862 }