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