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