EMCAL/AliEMCALClusterizerv1.cxx

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