EMCAL/AliEMCALClusterizerv1.cxx

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