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