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