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