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