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