EMCAL/AliEMCALClusterizerv1.cxx

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