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