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
  74 ClassImp(AliEMCALClusterizerv1)
  75
  76 AliEMCALGeometry * geom = 0;
  77 //____________________________________________________________________________
  78   AliEMCALClusterizerv1::AliEMCALClusterizerv1() : AliEMCALClusterizer()
  79 {
  80   // default ctor (to be used mainly by Streamer)
  81
  82   InitParameters() ;
  83   fDefaultInit = kTRUE ;
  84   geom = AliEMCALGeometry::GetInstance();
  85   geom->GetTransformationForSM(); // Global <-> Local
  86 }
  87
  88 //____________________________________________________________________________
  89 AliEMCALClusterizerv1::AliEMCALClusterizerv1(const TString alirunFileName, const TString eventFolderName)
  90 :AliEMCALClusterizer(alirunFileName, eventFolderName)
  91 {
  92   // ctor with the indication of the file where header Tree and digits Tree are stored
  93
  94   InitParameters() ;
  95   Init() ;
  96   fDefaultInit = kFALSE ;
  97 }
  98
  99 //____________________________________________________________________________
 100   AliEMCALClusterizerv1::~AliEMCALClusterizerv1()
 101 {
 102   // dtor
 103 }
 104
 105 //____________________________________________________________________________
 106 const TString AliEMCALClusterizerv1::BranchName() const
 107 {
 108    return GetName();
 109 }
 110
 111 //____________________________________________________________________________
 112 Float_t  AliEMCALClusterizerv1::Calibrate(Int_t amp) const
 113 {
 114   //To be replased later by the method, reading individual parameters from the database
 115   return -fADCpedestalECA + amp * fADCchannelECA ;
 116 }
 117
 118 //____________________________________________________________________________
 119 void AliEMCALClusterizerv1::Exec(Option_t * option)
 120 {
 121   // Steering method to perform clusterization for events
 122   // in the range from fFirstEvent to fLastEvent.
 123   // This range is optionally set by SetEventRange().
 124   // if fLastEvent=-1 (by default), then process events until the end.
 125
 126   if(strstr(option,"tim"))
 127     gBenchmark->Start("EMCALClusterizer");
 128
 129   if(strstr(option,"print"))
 130     Print("") ;
 131
 132   //AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
 133   AliRunLoader *rl = AliRunLoader::GetRunLoader();
 134   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
 135
 136   if (fLastEvent == -1)
 137     fLastEvent = rl->GetNumberOfEvents() - 1;
 138   Int_t nEvents   = fLastEvent - fFirstEvent + 1;
 139
 140   Int_t ievent ;
 141   rl->LoadDigits("EMCAL");
 142   for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
 143     rl->GetEvent(ievent);
 144     GetCalibrationParameters() ;
 145
 146     fNumberOfECAClusters = 0;
 147
 148     MakeClusters() ;
 149
 150     if(fToUnfold)
 151       MakeUnfolding() ;
 152
 153     WriteRecPoints() ;
 154
 155     if(strstr(option,"deb"))
 156       PrintRecPoints(option) ;
 157
 158     //increment the total number of recpoints per run
 159     fRecPointsInRun += emcalLoader->RecPoints()->GetEntriesFast() ;
 160   }
 161
 162   Unload();
 163
 164   if(strstr(option,"tim")){
 165     gBenchmark->Stop("EMCALClusterizer");
 166     printf("Exec took %f seconds for Clusterizing %f seconds per event",
 167          gBenchmark->GetCpuTime("EMCALClusterizer"), gBenchmark->GetCpuTime("EMCALClusterizer")/nEvents );
 168   }
 169
 170 }
 171
 172 //____________________________________________________________________________
 173 Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALRecPoint * emcRP, AliEMCALDigit ** maxAt, Float_t * maxAtEnergy,
 174                                     Int_t nPar, Float_t * fitparameters) const
 175 {
 176   // Calls TMinuit to fit the energy distribution of a cluster with several maxima
 177   // The initial values for fitting procedure are set equal to the positions of local maxima.
 178   // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
 179
 180   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
 181   TClonesArray *digits = emcalLoader->Digits();
 182   /*
 183   AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
 184   TClonesArray * digits = gime->Digits() ;
 185   */
 186
 187   gMinuit->mncler();                     // Reset Minuit's list of paramters
 188   gMinuit->SetPrintLevel(-1) ;           // No Printout
 189   gMinuit->SetFCN(AliEMCALClusterizerv1::UnfoldingChiSquare) ;
 190                                          // To set the address of the minimization function
 191   TList * toMinuit = new TList();
 192   toMinuit->AddAt(emcRP,0) ;
 193   toMinuit->AddAt(digits,1) ;
 194
 195   gMinuit->SetObjectFit(toMinuit) ;         // To tranfer pointer to UnfoldingChiSquare
 196
 197   // filling initial values for fit parameters
 198   AliEMCALDigit * digit ;
 199
 200   Int_t ierflg  = 0;
 201   Int_t index   = 0 ;
 202   Int_t nDigits = (Int_t) nPar / 3 ;
 203
 204   Int_t iDigit ;
 205
 206   for(iDigit = 0; iDigit < nDigits; iDigit++){
 207     digit = maxAt[iDigit];
 208
 209     Int_t relid[2] ;
 210     Float_t x = 0.;
 211     Float_t z = 0.;
 212     geom->AbsToRelNumbering(digit->GetId(), relid) ;
 213     geom->PosInAlice(relid, x, z) ;
 214
 215     Float_t energy = maxAtEnergy[iDigit] ;
 216
 217     gMinuit->mnparm(index, "x",  x, 0.1, 0, 0, ierflg) ;
 218     index++ ;
 219     if(ierflg != 0){
 220       Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : x = %f",  x ) ;
 221       return kFALSE;
 222     }
 223     gMinuit->mnparm(index, "z",  z, 0.1, 0, 0, ierflg) ;
 224     index++ ;
 225     if(ierflg != 0){
 226        Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : z = %f", z) ;
 227       return kFALSE;
 228     }
 229     gMinuit->mnparm(index, "Energy",  energy , 0.05*energy, 0., 4.*energy, ierflg) ;
 230     index++ ;
 231     if(ierflg != 0){
 232      Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : energy = %f", energy) ;
 233       return kFALSE;
 234     }
 235   }
 236
 237   Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
 238                       //  depends on it.
 239   Double_t p1 = 1.0 ;
 240   Double_t p2 = 0.0 ;
 241
 242   gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ;   // force TMinuit to reduce function calls
 243   gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ;   // force TMinuit to use my gradient
 244   gMinuit->SetMaxIterations(5);
 245   gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ;  // No Warnings
 246
 247   gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ;    // minimize
 248
 249   if(ierflg == 4){  // Minimum not found
 250     Error("FindFit", "EMCAL Unfolding  Fit not converged, cluster abandoned " ) ;
 251     return kFALSE ;
 252   }
 253   for(index = 0; index < nPar; index++){
 254     Double_t err ;
 255     Double_t val ;
 256     gMinuit->GetParameter(index, val, err) ;    // Returns value and error of parameter index
 257     fitparameters[index] = val ;
 258    }
 259
 260   delete toMinuit ;
 261   return kTRUE;
 262
 263 }
 264
 265 //____________________________________________________________________________
 266 void AliEMCALClusterizerv1::GetCalibrationParameters()
 267 {
 268   // Gets the parameters for the calibration from the digitizer
 269   //AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
 270   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
 271
 272   if ( !emcalLoader->Digitizer() )
 273     emcalLoader->LoadDigitizer();
 274   AliEMCALDigitizer * dig = dynamic_cast<AliEMCALDigitizer*>(emcalLoader->Digitizer());
 275
 276   fADCchannelECA   = dig->GetECAchannel() ;
 277   fADCpedestalECA  = dig->GetECApedestal();
 278 //PH  cout<<"ChannelECA, peds "<<fADCchannelECA<<" "<<fADCpedestalECA<<endl;
 279 }
 280
 281 //____________________________________________________________________________
 282 void AliEMCALClusterizerv1::Init()
 283 {
 284   // Make all memory allocations which can not be done in default constructor.
 285   // Attach the Clusterizer task to the list of EMCAL tasks
 286
 287   AliRunLoader *rl = AliRunLoader::GetRunLoader();
 288   geom = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
 289   geom->GetTransformationForSM(); // Global <-> Local
 290   AliInfo(Form("geom 0x%x",geom));
 291
 292 //Sub  fNTowers = geom->GetNZ() *  geom->GetNPhi() ;
 293   fNTowers =400; // ?
 294   if(!gMinuit)
 295     gMinuit = new TMinuit(100) ;
 296  //Sub if ( !gime->Clusterizer() )
 297  //Sub   gime->PostClusterizer(this);
 298  fHists = BookHists();
 299 }
 300
 301 //____________________________________________________________________________
 302 void AliEMCALClusterizerv1::InitParameters()
 303 {
 304   // Initializes the parameters for the Clusterizer
 305   fNumberOfECAClusters = 0;
 306
 307   fECAClusteringThreshold   = 0.5;  // value obtained from Alexei
 308   fECALocMaxCut = 0.03; // ??
 309
 310   fECAW0     = 4.5 ;
 311   fTimeGate = 1.e-8 ;
 312   fToUnfold = kFALSE ;
 313   fRecPointsInRun  = 0 ;
 314   fMinECut = 0.01; // have to be tune
 315 }
 316
 317 //____________________________________________________________________________
 318 Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2) const
 319 {
 320   // Gives the neighbourness of two digits = 0 are not neighbour but continue searching
 321   //                                       = 1 are neighbour
 322   //                                       = 2 is in different SM, quit from searching
 323   // neighbours are defined as digits having at least a common vertex
 324   // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
 325   //                                      which is compared to a digit (d2)  not yet in a cluster
 326
 327   static Int_t rv;
 328   static Int_t nSupMod1=0, nTower1=0, nIphi1=0, nIeta1=0, iphi1=0, ieta1=0;
 329   static Int_t nSupMod2=0, nTower2=0, nIphi2=0, nIeta2=0, iphi2=0, ieta2=0;
 330   static Int_t rowdiff, coldiff;
 331   rv = 0 ;
 332
 333   geom->GetCellIndex(d1->GetId(), nSupMod1,nTower1,nIphi1,nIeta1);
 334   geom->GetCellIndex(d2->GetId(), nSupMod2,nTower2,nIphi2,nIeta2);
 335   if(nSupMod1 != nSupMod2) return 2; // different SM
 336
 337   geom->GetCellPhiEtaIndexInSModule(nSupMod1,nTower1,nIphi1,nIeta1, iphi1,ieta1);
 338   geom->GetCellPhiEtaIndexInSModule(nSupMod2,nTower2,nIphi2,nIeta2, iphi2,ieta2);
 339
 340   rowdiff = TMath::Abs(iphi1 - iphi2);
 341   coldiff = TMath::Abs(ieta1 - ieta2) ;
 342
 343   if (( coldiff <= 1 )  && ( rowdiff <= 1 )) rv = 1;  // neighbours with at least commom vertex
 344
 345   if (gDebug == 2 && rv==1)
 346   printf("AreNeighbours: neighbours=%d, id1=%d, relid1=%d,%d \n id2=%d, relid2=%d,%d \n",
 347          rv, d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2);
 348
 349   return rv ;
 350 }
 351
 352 //____________________________________________________________________________
 353 void AliEMCALClusterizerv1::Unload()
 354 {
 355   // Unloads the Digits and RecPoints
 356   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
 357
 358   emcalLoader->UnloadDigits() ;
 359   emcalLoader->UnloadRecPoints() ;
 360 }
 361
 362 //____________________________________________________________________________
 363 void AliEMCALClusterizerv1::WriteRecPoints()
 364 {
 365
 366   // Creates new branches with given title
 367   // fills and writes into TreeR.
 368
 369   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
 370
 371   TObjArray * aECARecPoints = emcalLoader->RecPoints() ;
 372
 373   TClonesArray * digits = emcalLoader->Digits() ;
 374   TTree * treeR = emcalLoader->TreeR();
 375   if ( treeR==0 ) {
 376     emcalLoader->MakeRecPointsContainer();
 377     treeR = emcalLoader->TreeR();
 378   }
 379   Int_t index ;
 380
 381   //Evaluate position, dispersion and other RecPoint properties for EC section
 382   for(index = 0; index < aECARecPoints->GetEntries(); index++)
 383     (dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(index)))->EvalAll(fECAW0,digits) ;
 384
 385   aECARecPoints->Sort() ;
 386
 387   for(index = 0; index < aECARecPoints->GetEntries(); index++) {
 388     (dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(index)))->SetIndexInList(index) ;
 389     (dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(index)))->Print();
 390   }
 391
 392   Int_t bufferSize = 32000 ;
 393   Int_t splitlevel = 0 ;
 394
 395   //EC section branch
 396
 397   TBranch * branchECA = 0;
 398   if ((branchECA = treeR->GetBranch("EMCALECARP")))
 399     branchECA->SetAddress(&aECARecPoints);
 400   else
 401     treeR->Branch("EMCALECARP","TObjArray",&aECARecPoints,bufferSize,splitlevel);
 402   treeR->Fill() ;
 403
 404   emcalLoader->WriteRecPoints("OVERWRITE");
 405   //emcalLoader->WriteClusterizer("OVERWRITE");
 406   //emcalLoader->WriteReconstructioner("OVERWRITE");
 407 }
 408
 409 //____________________________________________________________________________
 410 void AliEMCALClusterizerv1::MakeClusters()
 411 {
 412   // Steering method to construct the clusters stored in a list of Reconstructed Points
 413   // A cluster is defined as a list of neighbour digits
 414
 415   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
 416
 417   TObjArray * aECARecPoints = emcalLoader->RecPoints() ;
 418
 419   if (geom==0)
 420      AliFatal("Did not get geometry from EMCALLoader");
 421
 422   aECARecPoints->Clear();
 423
 424   TClonesArray * digits  = emcalLoader->Digits() ;
 425   TClonesArray * digitsC =  dynamic_cast<TClonesArray*>(digits->Clone());
 426
 427   TIter nextdigit(digitsC) ;
 428   AliEMCALDigit * digit;
 429   double e=0.0, ehs = 0.0;
 430   while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigit())) ) { // clean up digits
 431     e = Calibrate(digit->GetAmp());
 432     AliEMCALHistoUtilities::FillH1(fHists, 10, digit->GetAmp());
 433     AliEMCALHistoUtilities::FillH1(fHists, 11, e);
 434     if(e < fMinECut ) digitsC->Remove(digit);
 435     else              ehs += e;
 436   }
 437   cout << " Number of digits " << digits->GetEntries() << " -> (e>" <<fMinECut <<")";
 438   cout << digitsC->GetEntries()<< " ehs "<<ehs<<endl;
 439
 440   // Clusterization starts
 441   //  cout << "Outer Loop" << endl;
 442   int ineb=0;
 443   nextdigit.Reset();
 444   AliEMCALRecPoint * recPoint = 0 ;
 445   while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigit())) ) { // scan over the list of digitsC
 446     recPoint = 0 ;
 447     TArrayI clusterECAdigitslist(1000); // what is this
 448
 449     if(geom->CheckAbsCellId(digit->GetId()) && (Calibrate(digit->GetAmp()) > fECAClusteringThreshold  ) ){
 450       Int_t iDigitInECACluster = 0; // start a new Tower RecPoint
 451       if(fNumberOfECAClusters >= aECARecPoints->GetSize()) aECARecPoints->Expand(2*fNumberOfECAClusters+1) ;
 452       AliEMCALRecPoint * rp = new  AliEMCALRecPoint("") ;
 453       aECARecPoints->AddAt(rp, fNumberOfECAClusters) ;
 454       recPoint = dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(fNumberOfECAClusters)) ;
 455       fNumberOfECAClusters++ ;
 456       recPoint->AddDigit(*digit, Calibrate(digit->GetAmp())) ;
 457       clusterECAdigitslist[iDigitInECACluster] = digit->GetIndexInList() ;
 458       iDigitInECACluster++ ;
 459       digitsC->Remove(digit) ;
 460       AliDebug(1,Form("MakeClusters: OK id = %d, ene = %f , thre = %f \n", digit->GetId(),
 461       Calibrate(digit->GetAmp()), fECAClusteringThreshold));
 462       nextdigit.Reset(); // will start from beggining
 463
 464       AliEMCALDigit * digitN = 0; // digi neighbor
 465       Int_t index = 0 ;
 466
 467       // Find the neighbours
 468       while (index < iDigitInECACluster){ // scan over digits already in cluster
 469         digit =  (AliEMCALDigit*)digits->At(clusterECAdigitslist[index]);
 470         index++ ;
 471         while ( (digitN = (AliEMCALDigit *)nextdigit())) { // scan over the reduced list of digits
 472           //          if( Calibrate(digitN->GetAmp()) < fMinECut  )  digitsC->Remove(digitN);
 473           ineb = AreNeighbours(digit, digitN);       // call (digit,digitN) in THAT oder !!!!!
 474           switch (ineb ) {
 475             case 0 :   // not a neighbour
 476               break ;
 477             case 1 :   // are neighbours
 478               recPoint->AddDigit(*digitN, Calibrate( digitN->GetAmp()) ) ;
 479               clusterECAdigitslist[iDigitInECACluster] = digitN->GetIndexInList() ;
 480               iDigitInECACluster++ ;
 481               digitsC->Remove(digitN) ;
 482                break ;
 483              case 2 :   // different SM
 484                break ;
 485           } // switch
 486         } // scan over the reduced list of digits
 487       } // scan over digits already in cluster
 488       nextdigit.Reset() ;  // will start from beggining
 489     }
 490   } // while digit
 491   if(recPoint) cout << "cl.e " << recPoint->GetEnergy() << endl;
 492   delete digitsC ;
 493   AliDebug(1,Form("total no of clusters %d from %d digits",fNumberOfECAClusters,digits->GetEntriesFast()));
 494 }
 495
 496 //____________________________________________________________________________
 497 void AliEMCALClusterizerv1::MakeUnfolding() const
 498 {
 499   Fatal("AliEMCALClusterizerv1::MakeUnfolding", "--> Unfolding not implemented") ;
 500
 501 }
 502
 503 //____________________________________________________________________________
 504 Double_t  AliEMCALClusterizerv1::ShowerShape(Double_t r)
 505 {
 506   // Shape of the shower (see EMCAL TDR)
 507   // If you change this function, change also the gradient evaluation in ChiSquare()
 508
 509   Double_t r4    = r*r*r*r ;
 510   Double_t r295  = TMath::Power(r, 2.95) ;
 511   Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
 512   return shape ;
 513 }
 514
 515 //____________________________________________________________________________
 516 void  AliEMCALClusterizerv1::UnfoldCluster(AliEMCALRecPoint * /*iniTower*/,
 517                                            Int_t /*nMax*/,
 518                                            AliEMCALDigit ** /*maxAt*/,
 519                                            Float_t * /*maxAtEnergy*/) const
 520 {
 521   // Performs the unfolding of a cluster with nMax overlapping showers
 522
 523   Fatal("UnfoldCluster", "--> Unfolding not implemented") ;
 524
 525 }
 526
 527 //_____________________________________________________________________________
 528 void AliEMCALClusterizerv1::UnfoldingChiSquare(Int_t & /*nPar*/, Double_t * /*Grad*/,
 529                                                Double_t & /*fret*/,
 530                                                Double_t * /*x*/, Int_t /*iflag*/)
 531 {
 532   // Calculates the Chi square for the cluster unfolding minimization
 533   // Number of parameters, Gradient, Chi squared, parameters, what to do
 534
 535   ::Fatal("UnfoldingChiSquare","Unfolding not implemented") ;
 536 }
 537 //____________________________________________________________________________
 538 void AliEMCALClusterizerv1::Print(Option_t * /*option*/)const
 539 {
 540   // Print clusterizer parameters
 541
 542   TString message("\n") ;
 543
 544   if( strcmp(GetName(), "") !=0 ){
 545
 546     // Print parameters
 547
 548     TString taskName(GetName()) ;
 549     taskName.ReplaceAll(Version(), "") ;
 550
 551     printf("--------------- ");
 552     printf(taskName.Data()) ;
 553     printf(" ");
 554     printf(GetTitle()) ;
 555     printf("-----------\n");
 556     printf("Clusterizing digits from the file: ");
 557     printf(taskName.Data());
 558     printf("\n                           Branch: ");
 559     printf(GetName());
 560     printf("\n                       ECA Local Maximum cut    = %f", fECALocMaxCut);
 561     printf("\n                       ECA Logarithmic weight   = %f", fECAW0);
 562     if(fToUnfold)
 563       printf("\nUnfolding on\n");
 564     else
 565       printf("\nUnfolding off\n");
 566
 567     printf("------------------------------------------------------------------");
 568   }
 569   else
 570     printf("AliEMCALClusterizerv1 not initialized ") ;
 571 }
 572
 573 //____________________________________________________________________________
 574 void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option)
 575 {
 576   // Prints list of RecPoints produced at the current pass of AliEMCALClusterizer
 577   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
 578   TObjArray * aECARecPoints = emcalLoader->RecPoints() ;
 579
 580   printf("PrintRecPoints: Clusterization result:") ;
 581
 582   printf("event # %d\n", emcalLoader->GetRunLoader()->GetEventNumber() ) ;
 583   printf("           Found %d ECA Rec Points\n ",
 584          aECARecPoints->GetEntriesFast()) ;
 585
 586   fRecPointsInRun +=  aECARecPoints->GetEntriesFast() ;
 587
 588   if(strstr(option,"all")) {
 589     Int_t index =0;
 590     printf("\n-----------------------------------------------------------------------\n") ;
 591     printf("Clusters in ECAL section\n") ;
 592     printf("Index    Ene(GeV) Multi Module     phi     r   theta    X    Y      Z   Dispersion Lambda 1   Lambda 2  # of prim  Primaries list\n") ;
 593    Float_t maxE=0;
 594    Float_t maxL1=0;
 595    Float_t maxL2=0;
 596    Float_t maxDis=0;
 597
 598     for (index = 0 ; index < aECARecPoints->GetEntries() ; index++) {
 599       AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint * >(aECARecPoints->At(index)) ;
 600       TVector3  globalpos;
 601       //rp->GetGlobalPosition(globalpos);
 602       TVector3  localpos;
 603       rp->GetLocalPosition(localpos);
 604       Float_t lambda[2];
 605       rp->GetElipsAxis(lambda);
 606       Int_t * primaries;
 607       Int_t nprimaries;
 608       primaries = rp->GetPrimaries(nprimaries);
 609       printf("\n%6d  %8.4f  %3d     %4.1f    %4.1f %4.1f  %4.1f %4.1f %4.1f    %4.1f   %4f  %4f    %2d     : ",
 610              rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(),
 611              globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(),
 612              rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ;
 613   /////////////
 614       if(rp->GetEnergy()>maxE){
 615               maxE=rp->GetEnergy();
 616               maxL1=lambda[0];
 617               maxL2=lambda[1];
 618               maxDis=rp->GetDispersion();
 619       }
 620       pointE->Fill(rp->GetEnergy());
 621       pointL1->Fill(lambda[0]);
 622       pointL2->Fill(lambda[1]);
 623       pointDis->Fill(rp->GetDispersion());
 624       pointMult->Fill(rp->GetMultiplicity());
 625       /////////////
 626       for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
 627         printf("%d ", primaries[iprimary] ) ;
 628       }
 629     }
 630
 631       MaxE->Fill(maxE);
 632       MaxL1->Fill(maxL1);
 633       MaxL2->Fill(maxL2);
 634       MaxDis->Fill(maxDis);
 635
 636
 637     printf("\n-----------------------------------------------------------------------\n");
 638   }
 639 }
 640 TList* AliEMCALClusterizerv1::BookHists()
 641 {
 642   gROOT->cd();
 643
 644         pointE = new TH1F("pointE","point energy", 2000, 0.0, 150.);
 645         pointL1 = new TH1F("pointL1","point L1", 1000, 0.0, 3.);
 646         pointL2 = new TH1F("pointL2","point L2", 1000, 0.0, 3.);
 647         pointDis = new TH1F("pointDis","point Dis", 1000, 0.0, 3.);
 648         pointMult = new TH1F("pointMult","point Mult", 100, 0.0, 100.);
 649         digitAmp = new TH1F("digitAmp","Digit Amplitude", 2000, 0.0, 5000.);
 650         MaxE = new TH1F("maxE","Max point energy", 2000, 0.0, 150.);
 651         MaxL1 = new TH1F("maxL1","Max point L1", 1000, 0.0, 3.);
 652         MaxL2 = new TH1F("maxL2","Max point L2", 1000, 0.0, 3.);
 653         MaxDis = new TH1F("maxDis","Max point Dis", 1000, 0.0, 3.); // 9
 654         //
 655         new TH1F("adcOfDigits","adc of digits(threshold control)", 1001, -0.5, 1000.5);   // 10
 656         new TH1F("energyOfDigits","energy of digits(threshold control)", 1000, 0.0, 1.);  // 11
 657
 658   return AliEMCALHistoUtilities::MoveHistsToList("EmcalClusterizerv1ControlHists", kFALSE);
 659 }
 660
 661 void AliEMCALClusterizerv1::SaveHists(const char *fn)
 662 {
 663   AliEMCALHistoUtilities::SaveListOfHists(fHists, fn, kTRUE);
 664 }
 665
 666 void AliEMCALClusterizerv1::Browse(TBrowser* b)
 667 {
 668   if(fHists) b->Add(fHists);
 669   TTask::Browse(b);
 670 }