EMCAL/AliEMCALClusterizerv1.cxx

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