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