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.1;  // value obtained from Aleksei
 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.01; // have to be tune
 411
 412   fCalibData               = 0 ;
 413 }
 414
 415 //____________________________________________________________________________
 416 Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2) const
 417 {
 418   // Gives the neighbourness of two digits = 0 are not neighbour ; continue searching
 419   //                                       = 1 are neighbour
 420   //                                       = 2 is in different SM; continue searching
 421   // neighbours are defined as digits having at least a common vertex
 422   // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
 423   //                                      which is compared to a digit (d2)  not yet in a cluster
 424
 425   static Int_t rv;
 426   static Int_t nSupMod1=0, nModule1=0, nIphi1=0, nIeta1=0, iphi1=0, ieta1=0;
 427   static Int_t nSupMod2=0, nModule2=0, nIphi2=0, nIeta2=0, iphi2=0, ieta2=0;
 428   static Int_t rowdiff, coldiff;
 429   rv = 0 ;
 430
 431   fGeom->GetCellIndex(d1->GetId(), nSupMod1,nModule1,nIphi1,nIeta1);
 432   fGeom->GetCellIndex(d2->GetId(), nSupMod2,nModule2,nIphi2,nIeta2);
 433   if(nSupMod1 != nSupMod2) return 2; // different SM
 434
 435   fGeom->GetCellPhiEtaIndexInSModule(nSupMod1,nModule1,nIphi1,nIeta1, iphi1,ieta1);
 436   fGeom->GetCellPhiEtaIndexInSModule(nSupMod2,nModule2,nIphi2,nIeta2, iphi2,ieta2);
 437
 438   rowdiff = TMath::Abs(iphi1 - iphi2);
 439   coldiff = TMath::Abs(ieta1 - ieta2) ;
 440
 441   // neighbours with at least commom side; May 11, 2007
 442   if ((coldiff==0 && abs(rowdiff)==1) || (rowdiff==0 && abs(coldiff)==1)) rv = 1;
 443
 444   if (gDebug == 2 && rv==1)
 445   printf("AreNeighbours: neighbours=%d, id1=%d, relid1=%d,%d \n id2=%d, relid2=%d,%d \n",
 446          rv, d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2);
 447
 448   return rv ;
 449 }
 450
 451 //____________________________________________________________________________
 452 Int_t AliEMCALClusterizerv1::AreInGroup(AliEMCALDigit * d1, AliEMCALDigit * d2) const
 453 {
 454   // Tells whether two digits fall within the same supermodule and
 455   // tower grouping.  The number of towers in a group is controlled by
 456   // the parameter nTowersInGroup
 457   //    = 0 are not in same group but continue searching
 458   //    = 1 same group
 459   //    = 2 is in different SM, quit from searching
 460   //    = 3 different tower group, quit from searching
 461   //
 462   // The order of d1 and d2 is important: first (d1) should be a digit
 463   // already in a cluster which is compared to a digit (d2)  not yet in a cluster
 464
 465   //JLK Question: does the quit from searching assume that the digits
 466   //are ordered, so that once you are in a different SM, you'll not
 467   //find another in the list that will match?  How about my TowerGroup search?
 468
 469   static Int_t rv;
 470   static Int_t nSupMod1=0, nModule1=0, nIphi1=0, nIeta1=0, iphi1=0, ieta1=0;
 471   static Int_t nSupMod2=0, nModule2=0, nIphi2=0, nIeta2=0, iphi2=0, ieta2=0;
 472   static Int_t towerGroup1 = -1, towerGroup2 = -1;
 473   rv = 0 ;
 474
 475   fGeom->GetCellIndex(d1->GetId(), nSupMod1,nModule1,nIphi1,nIeta1);
 476   fGeom->GetCellIndex(d2->GetId(), nSupMod2,nModule2,nIphi2,nIeta2);
 477   if(nSupMod1 != nSupMod2) return 2; // different SM
 478
 479   static Int_t nTowerInSM = fGeom->GetNCellsInSupMod()/fGeom->GetNCellsInModule();
 480
 481   //figure out which tower grouping each digit belongs to
 482   for(int it = 0; it < nTowerInSM/fNTowerInGroup; it++) {
 483     if(nModule1 <= nTowerInSM - it*fNTowerInGroup) towerGroup1 = it;
 484     if(nModule2 <= nTowerInSM - it*fNTowerInGroup) towerGroup2 = it;
 485   }
 486   if(towerGroup1 != towerGroup2) return 3; //different Towergroup
 487
 488   //same SM, same towergroup, we're happy
 489   if(towerGroup1 == towerGroup2 && towerGroup2 >= 0)
 490     rv = 1;
 491
 492   if (gDebug == 2 && rv==1)
 493   printf("AreInGroup: neighbours=%d, id1=%d, relid1=%d,%d \n id2=%d, relid2=%d,%d \n",
 494          rv, d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2);
 495
 496   return rv ;
 497 }
 498
 499 //____________________________________________________________________________
 500 void AliEMCALClusterizerv1::WriteRecPoints()
 501 {
 502
 503   // Creates new branches with given title
 504   // fills and writes into TreeR.
 505
 506   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
 507
 508   TObjArray * aECARecPoints = emcalLoader->RecPoints() ;
 509
 510   TClonesArray * digits = emcalLoader->Digits() ;
 511   TTree * treeR = emcalLoader->TreeR();
 512   if ( treeR==0 ) {
 513     emcalLoader->MakeRecPointsContainer();
 514     treeR = emcalLoader->TreeR();
 515   }
 516   else if (treeR->GetEntries() > 0) {
 517     Warning("WriteRecPoints","RecPoints already exist in output file. New Recpoitns will not be visible.");
 518   }
 519   Int_t index ;
 520
 521   //Evaluate position, dispersion and other RecPoint properties for EC section
 522   for(index = 0; index < aECARecPoints->GetEntries(); index++) {
 523     if (dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(index))->GetClusterType() != AliESDCaloCluster::kPseudoCluster)
 524        dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(index))->EvalAll(fECAW0,digits) ;
 525   }
 526
 527   aECARecPoints->Sort() ;
 528
 529   for(index = 0; index < aECARecPoints->GetEntries(); index++) {
 530     (dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(index)))->SetIndexInList(index) ;
 531     (dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(index)))->Print();
 532   }
 533
 534   Int_t bufferSize = 32000 ;
 535   Int_t splitlevel = 0 ;
 536
 537   //EC section branch
 538
 539   TBranch * branchECA = 0;
 540   if ((branchECA = treeR->GetBranch("EMCALECARP")))
 541     branchECA->SetAddress(&aECARecPoints);
 542   else
 543     treeR->Branch("EMCALECARP","TObjArray",&aECARecPoints,bufferSize,splitlevel);
 544   treeR->Fill() ;
 545
 546   emcalLoader->WriteRecPoints("OVERWRITE");
 547 }
 548
 549 //____________________________________________________________________________
 550 void AliEMCALClusterizerv1::MakeClusters(char* option)
 551 {
 552   // Steering method to construct the clusters stored in a list of Reconstructed Points
 553   // A cluster is defined as a list of neighbour digits
 554   // Mar 03, 2007 by PAI
 555
 556   if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader");
 557
 558   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
 559   TObjArray * aECARecPoints = emcalLoader->RecPoints() ;
 560   aECARecPoints->Clear();
 561
 562   TClonesArray *digits = emcalLoader->Digits();
 563
 564   // Set up TObjArray with pointers to digits to work on
 565   TObjArray *digitsC = new TObjArray();
 566   TIter nextdigit(digits);
 567   AliEMCALDigit *digit;
 568   while ( (digit = dynamic_cast<AliEMCALDigit*>(nextdigit())) ) {
 569     digitsC->AddLast(digit);
 570   }
 571
 572   //Start with pseudoclusters, if option
 573   if(strstr(option,"pseudo")) {
 574     //
 575     // New algorithm : will be created one pseudo cluster per module
 576     //
 577     AliDebug(1,Form("Pseudo clustering #digits : %i ",digits->GetEntries()));
 578
 579     AliEMCALRecPoint *recPoints[12]; // max size is 12 : see fGeom->GetNumberOfSuperModules();
 580     for(int i=0; i<12; i++) recPoints[i] = 0;
 581     TIter nextdigitC(digitsC) ;
 582
 583     // PseudoClusterization starts
 584     int nSM = 0; // # of SM
 585     while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigitC())) ) { // scan over the list of digitsC
 586       if(fGeom->CheckAbsCellId(digit->GetId()) ) { //Is this an EMCAL digit? Just maing sure...
 587         nSM = fGeom->GetSuperModuleNumber(digit->GetId());
 588         if(recPoints[nSM] == 0) {
 589           recPoints[nSM] = new AliEMCALRecPoint(Form("PC%2.2i", nSM));
 590           recPoints[nSM]->SetClusterType(AliESDCaloCluster::kPseudoCluster);
 591         }
 592         recPoints[nSM]->AddDigit(*digit, Calibrate(digit->GetAmp(), digit->GetId()));
 593       }
 594     }
 595     fNumberOfECAClusters = 0;
 596     for(int i=0; i<fGeom->GetNumberOfSuperModules(); i++) { // put non empty rec.points to container
 597       if(recPoints[i]) aECARecPoints->AddAt(recPoints[i], fNumberOfECAClusters++);
 598     }
 599     AliDebug(1,Form(" Number of PC %d ", fNumberOfECAClusters));
 600   }
 601
 602   //
 603   // Now do real clusters
 604   //
 605
 606   double e = 0.0, ehs = 0.0;
 607   TIter nextdigitC(digitsC);
 608
 609   while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigitC())) ) { // clean up digits
 610     e = Calibrate(digit->GetAmp(), digit->GetId());
 611     AliEMCALHistoUtilities::FillH1(fHists, 10, digit->GetAmp());
 612     AliEMCALHistoUtilities::FillH1(fHists, 11, e);
 613     if ( e < fMinECut || digit->GetTimeR() > fTimeCut )
 614       digitsC->Remove(digit);
 615     else
 616       ehs += e;
 617   }
 618   AliDebug(1,Form("MakeClusters: Number of digits %d  -> (e %f), ehs %d\n",
 619                   digits->GetEntries(),fMinECut,ehs));
 620
 621   nextdigitC.Reset();
 622
 623   while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigitC())) ) { // scan over the list of digitsC
 624     TArrayI clusterECAdigitslist(digits->GetEntries());
 625
 626     if(fGeom->CheckAbsCellId(digit->GetId()) && (Calibrate(digit->GetAmp(), digit->GetId()) > fECAClusteringThreshold  ) ){
 627       // start a new Tower RecPoint
 628       if(fNumberOfECAClusters >= aECARecPoints->GetSize()) aECARecPoints->Expand(2*fNumberOfECAClusters+1) ;
 629       AliEMCALRecPoint *recPoint = new  AliEMCALRecPoint("") ;
 630       aECARecPoints->AddAt(recPoint, fNumberOfECAClusters) ;
 631       recPoint = dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(fNumberOfECAClusters)) ;
 632       fNumberOfECAClusters++ ;
 633
 634       recPoint->SetClusterType(AliESDCaloCluster::kClusterv1);
 635
 636       recPoint->AddDigit(*digit, Calibrate(digit->GetAmp(), digit->GetId())) ;
 637       TObjArray clusterDigits;
 638       clusterDigits.AddLast(digit);
 639       digitsC->Remove(digit) ;
 640
 641       AliDebug(1,Form("MakeClusters: OK id = %d, ene = %f , cell.th. = %f \n", digit->GetId(),
 642       Calibrate(digit->GetAmp(),digit->GetId()), fECAClusteringThreshold));
 643
 644       // Grow cluster by finding neighbours
 645       TIter nextClusterDigit(&clusterDigits);
 646       while ( (digit = dynamic_cast<AliEMCALDigit*>(nextClusterDigit())) ) { // scan over digits in cluster
 647         TIter nextdigitN(digitsC);
 648         AliEMCALDigit *digitN = 0; // digi neighbor
 649         while ( (digitN = (AliEMCALDigit *)nextdigitN()) ) { // scan over all digits to look for neighbours
 650           if (AreNeighbours(digit, digitN)==1) {      // call (digit,digitN) in THAT oder !!!!!
 651             recPoint->AddDigit(*digitN, Calibrate(digitN->GetAmp(),digitN->GetId()) ) ;
 652             clusterDigits.AddLast(digitN) ;
 653             digitsC->Remove(digitN) ;
 654           } // if(ineb==1)
 655         } // scan over digits
 656       } // scan over digits already in cluster
 657       if(recPoint)
 658         AliDebug(2,Form("MakeClusters: %d digitd, energy %f \n", clusterDigits.GetEntries(), recPoint->GetEnergy()));
 659     } // If seed found
 660   } // while digit
 661
 662   delete digitsC ;
 663
 664   AliDebug(1,Form("total no of clusters %d from %d digits",fNumberOfECAClusters,digits->GetEntriesFast()));
 665 }
 666
 667 void AliEMCALClusterizerv1::MakeUnfolding() const
 668 {
 669   Fatal("AliEMCALClusterizerv1::MakeUnfolding", "--> Unfolding not implemented") ;
 670 }
 671
 672 //____________________________________________________________________________
 673 Double_t  AliEMCALClusterizerv1::ShowerShape(Double_t r)
 674 {
 675   // Shape of the shower (see EMCAL TDR)
 676   // If you change this function, change also the gradient evaluation in ChiSquare()
 677
 678   Double_t r4    = r*r*r*r ;
 679   Double_t r295  = TMath::Power(r, 2.95) ;
 680   Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
 681   return shape ;
 682 }
 683
 684 //____________________________________________________________________________
 685 void  AliEMCALClusterizerv1::UnfoldCluster(AliEMCALRecPoint * /*iniTower*/,
 686                                            Int_t /*nMax*/,
 687                                            AliEMCALDigit ** /*maxAt*/,
 688                                            Float_t * /*maxAtEnergy*/) const
 689 {
 690   // Performs the unfolding of a cluster with nMax overlapping showers
 691
 692   Fatal("UnfoldCluster", "--> Unfolding not implemented") ;
 693
 694 }
 695
 696 //_____________________________________________________________________________
 697 void AliEMCALClusterizerv1::UnfoldingChiSquare(Int_t & /*nPar*/, Double_t * /*Grad*/,
 698                                                Double_t & /*fret*/,
 699                                                Double_t * /*x*/, Int_t /*iflag*/)
 700 {
 701   // Calculates the Chi square for the cluster unfolding minimization
 702   // Number of parameters, Gradient, Chi squared, parameters, what to do
 703
 704   ::Fatal("UnfoldingChiSquare","Unfolding not implemented") ;
 705 }
 706 //____________________________________________________________________________
 707 void AliEMCALClusterizerv1::Print(Option_t * /*option*/)const
 708 {
 709   // Print clusterizer parameters
 710
 711   TString message("\n") ;
 712
 713   if( strcmp(GetName(), "") !=0 ){
 714
 715     // Print parameters
 716
 717     TString taskName(GetName()) ;
 718     taskName.ReplaceAll(Version(), "") ;
 719
 720     printf("--------------- ");
 721     printf(taskName.Data()) ;
 722     printf(" ");
 723     printf(GetTitle()) ;
 724     printf("-----------\n");
 725     printf("Clusterizing digits from the file: ");
 726     printf(taskName.Data());
 727     printf("\n                           Branch: ");
 728     printf(GetName());
 729     printf("\n                       ECA Local Maximum cut    = %f", fECALocMaxCut);
 730     printf("\n                       ECA Logarithmic weight   = %f", fECAW0);
 731     if(fToUnfold)
 732       printf("\nUnfolding on\n");
 733     else
 734       printf("\nUnfolding off\n");
 735
 736     printf("------------------------------------------------------------------");
 737   }
 738   else
 739     printf("AliEMCALClusterizerv1 not initialized ") ;
 740 }
 741
 742 //____________________________________________________________________________
 743 void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option)
 744 {
 745   // Prints list of RecPoints produced at the current pass of AliEMCALClusterizer
 746   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
 747   TObjArray * aECARecPoints = emcalLoader->RecPoints() ;
 748
 749   if(strstr(option,"deb")) {
 750     printf("PrintRecPoints: Clusterization result:") ;
 751
 752     printf("event # %d\n", emcalLoader->GetRunLoader()->GetEventNumber() ) ;
 753     printf("           Found %d ECA Rec Points\n ",
 754          aECARecPoints->GetEntriesFast()) ;
 755   }
 756
 757   fRecPointsInRun +=  aECARecPoints->GetEntriesFast() ;
 758
 759   if(strstr(option,"all")) {
 760     if(strstr(option,"deb")) {
 761       printf("\n-----------------------------------------------------------------------\n") ;
 762       printf("Clusters in ECAL section\n") ;
 763       printf("Index    Ene(GeV) Multi Module     GX    GY   GZ  lX    lY   lZ   Dispersion Lambda 1   Lambda 2  # of prim  Primaries list\n") ;
 764     }
 765    Int_t index =0;
 766    Float_t maxE=0;
 767    Float_t maxL1=0;
 768    Float_t maxL2=0;
 769    Float_t maxDis=0;
 770
 771     AliEMCALHistoUtilities::FillH1(fHists, 12, double(aECARecPoints->GetEntries()));
 772
 773     for (index = 0 ; index < aECARecPoints->GetEntries() ; index++) {
 774       AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint * >(aECARecPoints->At(index)) ;
 775       TVector3  globalpos;
 776       //rp->GetGlobalPosition(globalpos);
 777       TVector3  localpos;
 778       rp->GetLocalPosition(localpos);
 779       Float_t lambda[2];
 780       rp->GetElipsAxis(lambda);
 781       Int_t * primaries;
 782       Int_t nprimaries;
 783       primaries = rp->GetPrimaries(nprimaries);
 784       if(strstr(option,"deb"))
 785       printf("\n%6d  %8.4f  %3d     %4.1f    %4.1f %4.1f  %4.1f %4.1f %4.1f    %4.1f   %4f  %4f    %2d     : ",
 786              rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(),
 787              globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(),
 788              rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ;
 789   /////////////
 790       if(rp->GetEnergy()>maxE){
 791               maxE=rp->GetEnergy();
 792               maxL1=lambda[0];
 793               maxL2=lambda[1];
 794               maxDis=rp->GetDispersion();
 795       }
 796       fPointE->Fill(rp->GetEnergy());
 797       fPointL1->Fill(lambda[0]);
 798       fPointL2->Fill(lambda[1]);
 799       fPointDis->Fill(rp->GetDispersion());
 800       fPointMult->Fill(rp->GetMultiplicity());
 801       /////////////
 802       if(strstr(option,"deb")){
 803         for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
 804           printf("%d ", primaries[iprimary] ) ;
 805         }
 806       }
 807     }
 808
 809       fMaxE->Fill(maxE);
 810       fMaxL1->Fill(maxL1);
 811       fMaxL2->Fill(maxL2);
 812       fMaxDis->Fill(maxDis);
 813
 814     if(strstr(option,"deb"))
 815     printf("\n-----------------------------------------------------------------------\n");
 816   }
 817 }
 818 TList* AliEMCALClusterizerv1::BookHists()
 819 {
 820   //set up histograms for monitoring clusterizer performance
 821
 822   gROOT->cd();
 823
 824         fPointE = new TH1F("00_pointE","point energy", 2000, 0.0, 150.);
 825         fPointL1 = new TH1F("01_pointL1","point L1", 1000, 0.0, 3.);
 826         fPointL2 = new TH1F("02_pointL2","point L2", 1000, 0.0, 3.);
 827         fPointDis = new TH1F("03_pointDisp","point dispersion", 1000, 0.0, 10.);
 828         fPointMult = new TH1F("04_pointMult","#cell in point(cluster)", 101, -0.5, 100.5);
 829         fDigitAmp = new TH1F("05_digitAmp","Digit Amplitude", 2000, 0.0, 5000.);
 830         fMaxE = new TH1F("06_maxE","Max point energy", 2000, 0.0, 150.);
 831         fMaxL1 = new TH1F("07_maxL1","Largest (first) of eigenvalue of covariance matrix", 1000, 0.0, 3.);
 832         fMaxL2 = new TH1F("08_maxL2","Smalest (second) of eigenvalue of covariace matrix", 1000, 0.0, 3.);
 833         fMaxDis = new TH1F("09_maxDis","Point dispersion", 1000, 0.0, 10.); // 9
 834         //
 835         new TH1F("10_adcOfDigits","adc of digits(threshold control)", 1001, -0.5, 1000.5);   // 10
 836         new TH1F("11_energyOfDigits","energy of digits(threshold control)", 1000, 0.0, 1.);  // 11
 837         new TH1F("12_numberOfPoints","number of points(clusters)", 101, -0.5, 100.5);        // 12
 838
 839   return AliEMCALHistoUtilities::MoveHistsToList("EmcalClusterizerv1ControlHists", kFALSE);
 840 }
 841
 842 void AliEMCALClusterizerv1::SaveHists(const char *fn)
 843 {
 844   AliEMCALHistoUtilities::SaveListOfHists(fHists, fn, kTRUE);
 845 }
 846
 847 void  AliEMCALClusterizerv1::PrintRecoInfo()
 848 {
 849   printf(" AliEMCALClusterizerv1::PrintRecoInfo() : version %s \n", Version() );
 850   TH1F *h = (TH1F*)fHists->At(12);
 851   if(h) {
 852     printf(" ## Multiplicity of RecPoints ## \n");
 853     for(int i=1; i<=h->GetNbinsX(); i++) {
 854       int nbin = int((*h)[i]);
 855       int mult = int(h->GetBinCenter(i));
 856       if(nbin > 0) printf(" %i : %5.5i %6.3f %% \n", mult, nbin, 100.*nbin/h->GetEntries());
 857     }
 858   }
 859 }
 860
 861 void AliEMCALClusterizerv1::DrawLambdasHists()
 862 {
 863   if(fMaxL1) {
 864     fMaxL1->Draw();
 865     if(fMaxL2) fMaxL2->Draw("same");
 866     if(fMaxDis) {
 867       fMaxDis->Draw("same");
 868     }
 869   }
 870 }
 871
 872 void AliEMCALClusterizerv1::Browse(TBrowser* b)
 873 {
 874   if(fHists) b->Add(fHists);
 875   if(fGeom)  b->Add(fGeom);
 876   TTask::Browse(b);
 877 }