PHOS/AliPHOSClusterizerv1.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 /* $Log:
  19    1 October 2000. Yuri Kharlov:
  20      AreNeighbours()
  21      PPSD upper layer is considered if number of layers>1
  22
  23    18 October 2000. Yuri Kharlov:
  24      AliPHOSClusterizerv1()
  25      CPV clusterizing parameters added
  26
  27      MakeClusters()
  28      After first PPSD digit remove EMC digits only once
  29 */
  30 //*-- Author: Yves Schutz (SUBATECH)  & Dmitri Peressounko (SUBATECH & Kurchatov Institute)
  31 //////////////////////////////////////////////////////////////////////////////
  32 //  Clusterization class. Performs clusterization (collects neighbouring active cells) and
  33 //  unfolds the clusters having several local maxima.
  34 //  Results are stored in TreeR#, branches PHOSEmcRP (EMC recPoints),
  35 //  PHOSCpvRP (CPV RecPoints) and AliPHOSClusterizer (Clusterizer with all
  36 //  parameters including input digits branch title, thresholds etc.)
  37 //  This TTask is normally called from Reconstructioner, but can as well be used in
  38 //  standalone mode.
  39 // Use Case:
  40 //  root [0] AliPHOSClusterizerv1 * cl = new AliPHOSClusterizerv1("galice.root", "recpointsname", "digitsname")
  41 //  Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
  42 //               // reads gAlice from header file "galice.root", uses digits stored in the branch names "digitsname" (default = "Default")
  43 //               // and saves recpoints in branch named "recpointsname" (default = "digitsname")
  44 //  root [1] cl->ExecuteTask()
  45 //               //finds RecPoints in all events stored in galice.root
  46 //  root [2] cl->SetDigitsBranch("digits2")
  47 //               //sets another title for Digitis (input) branch
  48 //  root [3] cl->SetRecPointsBranch("recp2")
  49 //               //sets another title four output branches
  50 //  root [4] cl->SetEmcLocalMaxCut(0.03)
  51 //               //set clusterization parameters
  52 //  root [5] cl->ExecuteTask("deb all time")
  53 //               //once more finds RecPoints options are
  54 //               // deb - print number of found rec points
  55 //               // deb all - print number of found RecPoints and some their characteristics
  56 //               // time - print benchmarking results
  57
  58 // --- ROOT system ---
  59
  60 #include "TROOT.h"
  61 #include "TFile.h"
  62 #include "TFolder.h"
  63 #include "TMath.h"
  64 #include "TMinuit.h"
  65 #include "TTree.h"
  66 #include "TSystem.h"
  67 #include "TBenchmark.h"
  68
  69 // --- Standard library ---
  70
  71 // --- AliRoot header files ---
  72 #include "AliPHOSClusterizerv1.h"
  73 #include "AliPHOSCpvRecPoint.h"
  74 #include "AliPHOSDigit.h"
  75 #include "AliPHOSDigitizer.h"
  76 #include "AliPHOSEmcRecPoint.h"
  77 #include "AliPHOS.h"
  78 #include "AliPHOSGetter.h"
  79 #include "AliRun.h"
  80 #include "AliPHOSCalibrManager.h"
  81 #include "AliPHOSCalibrationData.h"
  82
  83 ClassImp(AliPHOSClusterizerv1)
  84
  85 //____________________________________________________________________________
  86   AliPHOSClusterizerv1::AliPHOSClusterizerv1() : AliPHOSClusterizer()
  87 {
  88   // default ctor (to be used mainly by Streamer)
  89
  90   InitParameters() ;
  91   fDefaultInit = kTRUE ;
  92 }
  93
  94 //____________________________________________________________________________
  95 AliPHOSClusterizerv1::AliPHOSClusterizerv1(const char* headerFile,const char* name, const Bool_t toSplit)
  96 :AliPHOSClusterizer(headerFile, name, toSplit)
  97 {
  98   // ctor with the indication of the file where header Tree and digits Tree are stored
  99
 100   InitParameters() ;
 101   Init() ;
 102   fDefaultInit = kFALSE ;
 103
 104 }
 105
 106 //____________________________________________________________________________
 107   AliPHOSClusterizerv1::~AliPHOSClusterizerv1()
 108 {
 109   delete fPedestals ;
 110   delete fGains ;
 111
 112   fSplitFile = 0 ;
 113
 114 }
 115
 116 //____________________________________________________________________________
 117 const TString AliPHOSClusterizerv1::BranchName() const
 118 {
 119   TString branchName(GetName() ) ;
 120   branchName.Remove(branchName.Index(Version())-1) ;
 121   return branchName ;
 122 }
 123
 124 //____________________________________________________________________________
 125 Float_t  AliPHOSClusterizerv1::Calibrate(Int_t amp, Int_t absId) const
 126 {
 127   if(fPedestals || fGains ){  //use calibration data
 128     if(!fPedestals || !fGains ){
 129       Error("Calibrate","Either Pedestals of Gains not set!") ;
 130       return 0 ;
 131     }
 132     Float_t en=(amp - fPedestals->Data(absId))*fGains->Data(absId) ;
 133     if(en>0)
 134       return en ;
 135     else
 136       return 0. ;
 137   }
 138   else{ //simulation
 139     if(absId <= fEmcCrystals) //calibrate as EMC
 140       return fADCpedestalEmc + amp*fADCchanelEmc ;
 141     else //calibrate as CPV
 142       return fADCpedestalCpv+ amp*fADCchanelCpv ;
 143   }
 144 }
 145
 146 //____________________________________________________________________________
 147 void AliPHOSClusterizerv1::Exec(Option_t * option)
 148 {
 149   // Steering method
 150   if( strcmp(GetName(), "")== 0 )
 151     Init() ;
 152
 153   if(strstr(option,"tim"))
 154     gBenchmark->Start("PHOSClusterizer");
 155
 156   if(strstr(option,"print"))
 157     Print("") ;
 158
 159   AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
 160   if(gime->BranchExists("RecPoints"))
 161     return ;
 162
 163   //test, if this is real data:
 164   TBranch * eh = gAlice->TreeE()->GetBranch("AliPHOSBeamTestEvent") ;
 165   if(eh){ //Read CalibrationData
 166     AliPHOSCalibrManager * cmngr=AliPHOSCalibrManager::GetInstance() ;
 167     if(!cmngr){ //not defined yet
 168       cmngr=AliPHOSCalibrManager::GetInstance("PHOSBTCalibration.root") ;
 169     }
 170     if(fPedestals)
 171       delete fPedestals ;
 172     fPedestals = new AliPHOSCalibrationData("Pedestals",fCalibrVersion) ;
 173     cmngr->ReadFromRoot(*fPedestals,fCalibrRun) ;
 174     if(fGains)
 175       delete fGains ;
 176     fGains = new AliPHOSCalibrationData("Gains",fCalibrVersion) ;
 177     cmngr->ReadFromRoot(*fGains,fCalibrRun) ;
 178   }
 179
 180   Int_t nevents = gime->MaxEvent() ;
 181   Int_t ievent ;
 182
 183   for(ievent = 0; ievent < nevents; ievent++){
 184     gime->Event(ievent,"D") ;
 185     Int_t pattern = gime->EventPattern() ;
 186     if(pattern){
 187             printf("pattern %d \n",pattern) ;
 188       if(!fPatterns){
 189         Error("Exec","Patterns for reconstruction of events are not set, use SetPatterns() ") ;
 190         return ;
 191       }
 192       Bool_t skip = kTRUE ;
 193       Int_t npat = fPatterns->GetSize() ;
 194       for(Int_t i = 0; i<npat;i++)
 195         if(pattern==fPatterns->At(i)){
 196           skip = kFALSE ;
 197           break ;
 198         }
 199       if(skip)
 200         continue ;
 201     }
 202
 203     if(ievent == 0)
 204       GetCalibrationParameters() ;
 205
 206     fNumberOfEmcClusters  = fNumberOfCpvClusters  = 0 ;
 207
 208     MakeClusters() ;
 209
 210     if(fToUnfold)
 211       MakeUnfolding() ;
 212
 213     WriteRecPoints(ievent) ;
 214
 215     if(strstr(option,"deb"))
 216       PrintRecPoints(option) ;
 217
 218     //increment the total number of digits per run
 219     fRecPointsInRun += gime->EmcRecPoints()->GetEntriesFast() ;
 220     fRecPointsInRun += gime->CpvRecPoints()->GetEntriesFast() ;
 221  }
 222
 223   if(strstr(option,"tim")){
 224     gBenchmark->Stop("PHOSClusterizer");
 225     Info("Exec", "  took %f seconds for Clusterizing %f seconds per event \n",
 226          gBenchmark->GetCpuTime("PHOSClusterizer"),
 227          gBenchmark->GetCpuTime("PHOSClusterizer")/nevents ) ;
 228   }
 229 }
 230
 231 //____________________________________________________________________________
 232 Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, AliPHOSDigit ** maxAt, Float_t * maxAtEnergy,
 233                                     Int_t nPar, Float_t * fitparameters) const
 234 {
 235   // Calls TMinuit to fit the energy distribution of a cluster with several maxima
 236   // The initial values for fitting procedure are set equal to the positions of local maxima.
 237   // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
 238
 239   AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
 240   TClonesArray * digits = gime->Digits() ;
 241
 242
 243   gMinuit->mncler();                     // Reset Minuit's list of paramters
 244   gMinuit->SetPrintLevel(-1) ;           // No Printout
 245   gMinuit->SetFCN(AliPHOSClusterizerv1::UnfoldingChiSquare) ;
 246                                          // To set the address of the minimization function
 247
 248   TList * toMinuit = new TList();
 249   toMinuit->AddAt(emcRP,0) ;
 250   toMinuit->AddAt(digits,1) ;
 251
 252   gMinuit->SetObjectFit(toMinuit) ;         // To tranfer pointer to UnfoldingChiSquare
 253
 254   // filling initial values for fit parameters
 255   AliPHOSDigit * digit ;
 256
 257   Int_t ierflg  = 0;
 258   Int_t index   = 0 ;
 259   Int_t nDigits = (Int_t) nPar / 3 ;
 260
 261   Int_t iDigit ;
 262
 263   const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
 264
 265   for(iDigit = 0; iDigit < nDigits; iDigit++){
 266     digit = maxAt[iDigit];
 267
 268     Int_t relid[4] ;
 269     Float_t x = 0.;
 270     Float_t z = 0.;
 271     geom->AbsToRelNumbering(digit->GetId(), relid) ;
 272     geom->RelPosInModule(relid, x, z) ;
 273
 274     Float_t energy = maxAtEnergy[iDigit] ;
 275
 276     gMinuit->mnparm(index, "x",  x, 0.1, 0, 0, ierflg) ;
 277     index++ ;
 278     if(ierflg != 0){
 279       Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : x = %f\n", x ) ;
 280       return kFALSE;
 281     }
 282     gMinuit->mnparm(index, "z",  z, 0.1, 0, 0, ierflg) ;
 283     index++ ;
 284     if(ierflg != 0){
 285        Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : z =%f\n", z ) ;
 286       return kFALSE;
 287     }
 288     gMinuit->mnparm(index, "Energy",  energy , 0.05*energy, 0., 4.*energy, ierflg) ;
 289     index++ ;
 290     if(ierflg != 0){
 291       Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : energy = %f\n", energy ) ;
 292       return kFALSE;
 293     }
 294   }
 295
 296   Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
 297                       //  depends on it.
 298   Double_t p1 = 1.0 ;
 299   Double_t p2 = 0.0 ;
 300
 301   gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ;   // force TMinuit to reduce function calls
 302   gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ;   // force TMinuit to use my gradient
 303   gMinuit->SetMaxIterations(5);
 304   gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ;  // No Warnings
 305
 306   gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ;    // minimize
 307
 308   if(ierflg == 4){  // Minimum not found
 309     Warning("FindFit", "PHOS Unfolding fit not converged, cluster abandoned\n" );
 310     return kFALSE ;
 311   }
 312   for(index = 0; index < nPar; index++){
 313     Double_t err ;
 314     Double_t val ;
 315     gMinuit->GetParameter(index, val, err) ;    // Returns value and error of parameter index
 316     fitparameters[index] = val ;
 317    }
 318
 319   delete toMinuit ;
 320   return kTRUE;
 321
 322 }
 323
 324 //____________________________________________________________________________
 325 void AliPHOSClusterizerv1::GetCalibrationParameters()
 326 {
 327   AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
 328
 329   const TTask * task = gime->Digitizer(BranchName()) ;
 330   if(strcmp(task->IsA()->GetName(),"AliPHOSDigitizer")==0){
 331     const AliPHOSDigitizer * dig = static_cast<const AliPHOSDigitizer *>(task) ;
 332
 333     fADCchanelEmc   = dig->GetEMCchannel() ;
 334     fADCpedestalEmc = dig->GetEMCpedestal();
 335
 336     fADCchanelCpv   = dig->GetCPVchannel() ;
 337     fADCpedestalCpv = dig->GetCPVpedestal() ;
 338   }
 339   else{
 340  ;
 341 //    fCalibrationDB = gime->CalibrationDB();
 342   }
 343 }
 344
 345 //____________________________________________________________________________
 346 void AliPHOSClusterizerv1::Init()
 347 {
 348   // Make all memory allocations which can not be done in default constructor.
 349   // Attach the Clusterizer task to the list of PHOS tasks
 350
 351   if ( strcmp(GetTitle(), "") == 0 )
 352     SetTitle("galice.root") ;
 353
 354   TString branchname = GetName() ;
 355   branchname.Remove(branchname.Index(Version())-1) ;
 356
 357   AliPHOSGetter * gime = AliPHOSGetter::GetInstance(GetTitle(),branchname.Data(), fToSplit ) ;
 358   if ( gime == 0 ) {
 359     Error("Init", "Could not obtain the Getter object !") ;
 360     return ;
 361   }
 362
 363   fSplitFile = 0 ;
 364   if(fToSplit){
 365     // construct the name of the file as /path/EMCAL.SDigits.root
 366     //First - extract full path if necessary
 367     TString fileName(GetTitle()) ;
 368     Ssiz_t islash = fileName.Last('/') ;
 369     if(islash<fileName.Length())
 370       fileName.Remove(islash+1,fileName.Length()) ;
 371     else
 372       fileName="" ;
 373     // Next - append the file name
 374     fileName+="PHOS.RecData." ;
 375     if((strcmp(branchname.Data(),"Default")!=0)&&(strcmp(branchname.Data(),"")!=0)){
 376       fileName+=branchname ;
 377       fileName+="." ;
 378     }
 379     fileName+="root" ;
 380     // Finally - check if the file already opened or open the file
 381     fSplitFile = static_cast<TFile*>(gROOT->GetFile(fileName.Data()));
 382     if(!fSplitFile)
 383       fSplitFile =  TFile::Open(fileName.Data(),"update") ;
 384   }
 385
 386
 387
 388   const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
 389   fEmcCrystals = geom->GetNModules() *  geom->GetNCristalsInModule() ;
 390
 391   if(!gMinuit)
 392     gMinuit = new TMinuit(100) ;
 393
 394   gime->PostClusterizer(this) ;
 395   gime->PostRecPoints(branchname) ;
 396
 397 }
 398
 399 //____________________________________________________________________________
 400 void AliPHOSClusterizerv1::InitParameters()
 401 {
 402
 403   fNumberOfCpvClusters     = 0 ;
 404   fNumberOfEmcClusters     = 0 ;
 405
 406   fCpvClusteringThreshold  = 0.0;
 407   fEmcClusteringThreshold  = 0.2;
 408
 409   fEmcLocMaxCut            = 0.03 ;
 410   fCpvLocMaxCut            = 0.03 ;
 411
 412   fW0                      = 4.5 ;
 413   fW0CPV                   = 4.0 ;
 414
 415   fEmcTimeGate             = 1.e-8 ;
 416
 417   fToUnfold                = kTRUE ;
 418
 419   fPurifyThreshold         = 0.012 ;
 420
 421   fPedestals = 0 ;
 422   fGains     = 0 ;
 423   fCalibrVersion= "v1" ;
 424   fCalibrRun = 1 ;
 425   fPatterns  = 0 ;
 426
 427   TString clusterizerName( GetName()) ;
 428   if (clusterizerName.IsNull() )
 429     clusterizerName = "Default" ;
 430   clusterizerName.Append(":") ;
 431   clusterizerName.Append(Version()) ;
 432   SetName(clusterizerName) ;
 433   fRecPointsInRun          = 0 ;
 434
 435
 436 }
 437
 438 //____________________________________________________________________________
 439 Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)const
 440 {
 441   // Gives the neighbourness of two digits = 0 are not neighbour but continue searching
 442   //                                       = 1 are neighbour
 443   //                                       = 2 are not neighbour but do not continue searching
 444   // neighbours are defined as digits having at least a common vertex
 445   // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
 446   //                                      which is compared to a digit (d2)  not yet in a cluster
 447
 448   const AliPHOSGeometry * geom = AliPHOSGetter::GetInstance()->PHOSGeometry() ;
 449
 450   Int_t rv = 0 ;
 451
 452   Int_t relid1[4] ;
 453   geom->AbsToRelNumbering(d1->GetId(), relid1) ;
 454
 455   Int_t relid2[4] ;
 456   geom->AbsToRelNumbering(d2->GetId(), relid2) ;
 457
 458   if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module
 459     Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
 460     Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
 461
 462     if (( coldiff <= 1 )  && ( rowdiff <= 1 )){
 463       if((relid1[1] != 0) || (TMath::Abs(d1->GetTime() - d2->GetTime() ) < fEmcTimeGate))
 464       rv = 1 ;
 465     }
 466     else {
 467       if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1))
 468         rv = 2; //  Difference in row numbers is too large to look further
 469     }
 470
 471   }
 472   else {
 473
 474     if( (relid1[0] < relid2[0]) || (relid1[1] != relid2[1]) )
 475       rv=2 ;
 476
 477   }
 478
 479   return rv ;
 480 }
 481
 482
 483 //____________________________________________________________________________
 484 Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit) const
 485 {
 486   // Tells if (true) or not (false) the digit is in a PHOS-EMC module
 487
 488   Bool_t rv = kFALSE ;
 489   const AliPHOSGeometry * geom = AliPHOSGetter::GetInstance()->PHOSGeometry() ;
 490
 491   Int_t nEMC = geom->GetNModules()*geom->GetNPhi()*geom->GetNZ();
 492
 493   if(digit->GetId() <= nEMC )   rv = kTRUE;
 494
 495   return rv ;
 496 }
 497
 498 //____________________________________________________________________________
 499 Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit) const
 500 {
 501   // Tells if (true) or not (false) the digit is in a PHOS-CPV module
 502
 503   Bool_t rv = kFALSE ;
 504   const AliPHOSGeometry * geom = AliPHOSGetter::GetInstance()->PHOSGeometry() ;
 505
 506   Int_t nEMC = geom->GetNModules()*geom->GetNPhi()*geom->GetNZ();
 507
 508   if(digit->GetId() > nEMC )   rv = kTRUE;
 509
 510   return rv ;
 511 }
 512
 513 //____________________________________________________________________________
 514 void AliPHOSClusterizerv1::WriteRecPoints(Int_t event)
 515 {
 516
 517   // Creates new branches with given title
 518   // fills and writes into TreeR.
 519
 520
 521   AliPHOSGetter *gime = AliPHOSGetter::GetInstance() ;
 522   TObjArray * emcRecPoints = gime->EmcRecPoints() ;
 523   TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
 524   TClonesArray * digits = gime->Digits() ;
 525   TTree * treeR ;
 526
 527   if(fToSplit){
 528     if(!fSplitFile)
 529       return ;
 530     fSplitFile->cd() ;
 531     TString name("TreeR") ;
 532     name += event ;
 533     treeR = dynamic_cast<TTree*>(fSplitFile->Get(name));
 534   }
 535   else{
 536     treeR = gAlice->TreeR();
 537   }
 538
 539   if(!treeR){
 540     if(fSplitFile)
 541       gAlice->MakeTree("R", fSplitFile);
 542     else
 543       gAlice->MakeTree("R", gROOT->GetFile(GetTitle()));
 544     treeR = gAlice->TreeR() ;
 545   }
 546
 547   Int_t index ;
 548   //Evaluate position, dispersion and other RecPoint properties...
 549   for(index = 0; index < emcRecPoints->GetEntries(); index++){
 550     AliPHOSEmcRecPoint * emcrp = static_cast<AliPHOSEmcRecPoint *>(emcRecPoints->At(index)) ;
 551     emcrp->Purify(fPurifyThreshold) ;
 552     if(emcrp->GetMultiplicity())
 553       emcrp->EvalAll(fW0,digits) ;
 554     else
 555       emcRecPoints->Remove(emcrp) ;
 556  }
 557
 558   emcRecPoints->Compress() ;
 559   emcRecPoints->Sort() ;
 560   for(index = 0; index < emcRecPoints->GetEntries(); index++)
 561     ((AliPHOSEmcRecPoint *)emcRecPoints->At(index))->SetIndexInList(index) ;
 562
 563   emcRecPoints->Expand(emcRecPoints->GetEntriesFast()) ;
 564
 565   //Now the same for CPV
 566   for(index = 0; index < cpvRecPoints->GetEntries(); index++)
 567     ((AliPHOSRecPoint *)cpvRecPoints->At(index))->EvalAll(fW0CPV,digits)  ;
 568
 569   cpvRecPoints->Sort() ;
 570
 571   for(index = 0; index < cpvRecPoints->GetEntries(); index++)
 572     ((AliPHOSRecPoint *)cpvRecPoints->At(index))->SetIndexInList(index) ;
 573
 574   cpvRecPoints->Expand(cpvRecPoints->GetEntriesFast()) ;
 575
 576   Int_t bufferSize = 32000 ;
 577   Int_t splitlevel = 0 ;
 578
 579   //First EMC
 580   TBranch * emcBranch = treeR->Branch("PHOSEmcRP","TObjArray",&emcRecPoints,bufferSize,splitlevel);
 581   emcBranch->SetTitle(BranchName());
 582
 583   //Now CPV branch
 584   TBranch * cpvBranch = treeR->Branch("PHOSCpvRP","TObjArray",&cpvRecPoints,bufferSize,splitlevel);
 585   cpvBranch->SetTitle(BranchName());
 586
 587   //And Finally  clusterizer branch
 588   AliPHOSClusterizerv1 * cl = this ;
 589   TBranch * clusterizerBranch = treeR->Branch("AliPHOSClusterizer","AliPHOSClusterizerv1",
 590                                               &cl,bufferSize,splitlevel);
 591   clusterizerBranch->SetTitle(BranchName());
 592
 593   emcBranch        ->Fill() ;
 594   cpvBranch        ->Fill() ;
 595   clusterizerBranch->Fill() ;
 596
 597   treeR->AutoSave() ;
 598   if(gAlice->TreeR()!=treeR)
 599     treeR->Delete();
 600 }
 601
 602 //____________________________________________________________________________
 603 void AliPHOSClusterizerv1::MakeClusters()
 604 {
 605   // Steering method to construct the clusters stored in a list of Reconstructed Points
 606   // A cluster is defined as a list of neighbour digits
 607
 608   AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
 609
 610   TObjArray * emcRecPoints = gime->EmcRecPoints(BranchName()) ;
 611   TObjArray * cpvRecPoints = gime->CpvRecPoints(BranchName()) ;
 612   emcRecPoints->Delete() ;
 613   cpvRecPoints->Delete() ;
 614
 615   TClonesArray * digits = gime->Digits() ;
 616    if ( !digits ) {
 617    Fatal("MakeClusters", "Digits with name %s not found !", BranchName().Data() ) ;
 618   }
 619   TClonesArray * digitsC =  (TClonesArray*)digits->Clone() ;
 620
 621
 622   // Clusterization starts
 623
 624   TIter nextdigit(digitsC) ;
 625   AliPHOSDigit * digit ;
 626   Bool_t notremoved = kTRUE ;
 627
 628   while ( (digit = (AliPHOSDigit *)nextdigit()) ) { // scan over the list of digitsC
 629     AliPHOSRecPoint * clu = 0 ;
 630
 631     TArrayI clusterdigitslist(1500) ;
 632     Int_t index ;
 633    // cout << "digit " << digit << " " << Calibrate(digit->GetAmp(),digit->GetId()) << endl ;
 634
 635     if (( IsInEmc (digit) && Calibrate(digit->GetAmp(),digit->GetId()) > fEmcClusteringThreshold  ) ||
 636         ( IsInCpv (digit) && Calibrate(digit->GetAmp(),digit->GetId()) > fCpvClusteringThreshold  ) ) {
 637       Int_t iDigitInCluster = 0 ;
 638
 639       if  ( IsInEmc(digit) ) {
 640         // start a new EMC RecPoint
 641         if(fNumberOfEmcClusters >= emcRecPoints->GetSize())
 642           emcRecPoints->Expand(2*fNumberOfEmcClusters+1) ;
 643
 644         emcRecPoints->AddAt(new  AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ;
 645         clu = (AliPHOSEmcRecPoint *) emcRecPoints->At(fNumberOfEmcClusters) ;
 646         fNumberOfEmcClusters++ ;
 647         clu->AddDigit(*digit, Calibrate(digit->GetAmp(),digit->GetId())) ;
 648         clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
 649         iDigitInCluster++ ;
 650         digitsC->Remove(digit) ;
 651
 652       } else {
 653
 654         // start a new CPV cluster
 655         if(fNumberOfCpvClusters >= cpvRecPoints->GetSize())
 656           cpvRecPoints->Expand(2*fNumberOfCpvClusters+1);
 657
 658         cpvRecPoints->AddAt(new AliPHOSCpvRecPoint(""), fNumberOfCpvClusters) ;
 659
 660         clu =  (AliPHOSCpvRecPoint *) cpvRecPoints->At(fNumberOfCpvClusters)  ;
 661         fNumberOfCpvClusters++ ;
 662         clu->AddDigit(*digit, Calibrate(digit->GetAmp(),digit->GetId()) ) ;
 663         clusterdigitslist[iDigitInCluster] = digit->GetIndexInList()  ;
 664         iDigitInCluster++ ;
 665         digitsC->Remove(digit) ;
 666         nextdigit.Reset() ;
 667
 668         // Here we remove remaining EMC digits, which cannot make a cluster
 669
 670         if( notremoved ) {
 671           while( ( digit = (AliPHOSDigit *)nextdigit() ) ) {
 672             if( IsInEmc(digit) )
 673               digitsC->Remove(digit) ;
 674             else
 675               break ;
 676           }
 677           notremoved = kFALSE ;
 678         }
 679
 680       } // else
 681
 682       nextdigit.Reset() ;
 683
 684       AliPHOSDigit * digitN ;
 685       index = 0 ;
 686       while (index < iDigitInCluster){ // scan over digits already in cluster
 687         digit =  (AliPHOSDigit*)digits->At(clusterdigitslist[index])  ;
 688         index++ ;
 689         while ( (digitN = (AliPHOSDigit *)nextdigit()) ) { // scan over the reduced list of digits
 690           Int_t ineb = AreNeighbours(digit, digitN);       // call (digit,digitN) in THAT oder !!!!!
 691           switch (ineb ) {
 692           case 0 :   // not a neighbour
 693             break ;
 694           case 1 :   // are neighbours
 695             clu->AddDigit(*digitN, Calibrate( digitN->GetAmp(), digitN->GetId() ) ) ;
 696             clusterdigitslist[iDigitInCluster] = digitN->GetIndexInList() ;
 697             iDigitInCluster++ ;
 698             digitsC->Remove(digitN) ;
 699             break ;
 700           case 2 :   // too far from each other
 701             goto endofloop;
 702           } // switch
 703
 704         } // while digitN
 705
 706       endofloop: ;
 707         nextdigit.Reset() ;
 708
 709       } // loop over cluster
 710
 711     } // energy theshold
 712
 713
 714   } // while digit
 715
 716   delete digitsC ;
 717
 718 }
 719
 720 //____________________________________________________________________________
 721 void AliPHOSClusterizerv1::MakeUnfolding()
 722 {
 723   // Unfolds clusters using the shape of an ElectroMagnetic shower
 724   // Performs unfolding of all EMC/CPV clusters
 725
 726   AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
 727
 728
 729   const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
 730   TObjArray * emcRecPoints = gime->EmcRecPoints(BranchName()) ;
 731   TObjArray * cpvRecPoints = gime->CpvRecPoints(BranchName()) ;
 732   TClonesArray * digits = gime->Digits() ;
 733
 734   // Unfold first EMC clusters
 735   if(fNumberOfEmcClusters > 0){
 736
 737     Int_t nModulesToUnfold = geom->GetNModules() ;
 738
 739     Int_t numberofNotUnfolded = fNumberOfEmcClusters ;
 740     Int_t index ;
 741     for(index = 0 ; index < numberofNotUnfolded ; index++){
 742
 743       AliPHOSEmcRecPoint * emcRecPoint = (AliPHOSEmcRecPoint *) emcRecPoints->At(index) ;
 744       if(emcRecPoint->GetPHOSMod()> nModulesToUnfold)
 745         break ;
 746
 747       Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
 748       AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
 749       Float_t * maxAtEnergy = new Float_t[nMultipl] ;
 750       Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fEmcLocMaxCut,digits) ;
 751
 752       if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0
 753         UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
 754         if(emcRecPoint->GetNExMax()!=-1){ //Do not remove clusters which failed to unfold
 755           emcRecPoints->Remove(emcRecPoint);
 756           emcRecPoints->Compress() ;
 757           index-- ;
 758           fNumberOfEmcClusters -- ;
 759           numberofNotUnfolded-- ;
 760         }
 761       }
 762       else
 763         emcRecPoint->SetNExMax(1) ; //Only one local maximum
 764
 765       delete[] maxAt ;
 766       delete[] maxAtEnergy ;
 767     }
 768   }
 769   // Unfolding of EMC clusters finished
 770
 771
 772   // Unfold now CPV clusters
 773   if(fNumberOfCpvClusters > 0){
 774
 775     Int_t nModulesToUnfold = geom->GetNModules() ;
 776
 777     Int_t numberofCpvNotUnfolded = fNumberOfCpvClusters ;
 778     Int_t index ;
 779     for(index = 0 ; index < numberofCpvNotUnfolded ; index++){
 780
 781       AliPHOSRecPoint * recPoint = (AliPHOSRecPoint *) cpvRecPoints->At(index) ;
 782
 783       if(recPoint->GetPHOSMod()> nModulesToUnfold)
 784         break ;
 785
 786       AliPHOSEmcRecPoint * emcRecPoint = (AliPHOSEmcRecPoint*) recPoint ;
 787
 788       Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
 789       AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
 790       Float_t * maxAtEnergy = new Float_t[nMultipl] ;
 791       Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fCpvLocMaxCut,digits) ;
 792
 793       if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0
 794         UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
 795         cpvRecPoints->Remove(emcRecPoint);
 796         cpvRecPoints->Compress() ;
 797         index-- ;
 798         numberofCpvNotUnfolded-- ;
 799         fNumberOfCpvClusters-- ;
 800       }
 801
 802       delete[] maxAt ;
 803       delete[] maxAtEnergy ;
 804     }
 805   }
 806   //Unfolding of Cpv clusters finished
 807
 808 }
 809
 810 //____________________________________________________________________________
 811 Double_t  AliPHOSClusterizerv1::ShowerShape(Double_t r)
 812 {
 813   // Shape of the shower (see PHOS TDR)
 814   // If you change this function, change also the gradient evaluation in ChiSquare()
 815
 816   Double_t r4    = r*r*r*r ;
 817   Double_t r295  = TMath::Power(r, 2.95) ;
 818   Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
 819   return shape ;
 820 }
 821
 822 //____________________________________________________________________________
 823 void  AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,
 824                                                  Int_t nMax,
 825                                                  AliPHOSDigit ** maxAt,
 826                                                  Float_t * maxAtEnergy)
 827 {
 828   // Performs the unfolding of a cluster with nMax overlapping showers
 829
 830   AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
 831
 832   const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
 833   const TClonesArray * digits = gime->Digits() ;
 834   TObjArray * emcRecPoints = gime->EmcRecPoints() ;
 835   TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
 836
 837   Int_t nPar = 3 * nMax ;
 838   Float_t * fitparameters = new Float_t[nPar] ;
 839
 840   Bool_t rv = FindFit(iniEmc, maxAt, maxAtEnergy, nPar, fitparameters) ;
 841   if( !rv ) {
 842     // Fit failed, return
 843     iniEmc->SetNExMax(-1) ;
 844     delete[] fitparameters ;
 845     return ;
 846   }
 847
 848   // create ufolded rec points and fill them with new energy lists
 849   // First calculate energy deposited in each sell in accordance with fit (without fluctuations): efit[]
 850   // and later correct this number in acordance with actual energy deposition
 851
 852   Int_t nDigits = iniEmc->GetMultiplicity() ;
 853   Float_t * efit = new Float_t[nDigits] ;
 854   Float_t xDigit=0.,zDigit=0.,distance=0. ;
 855   Float_t xpar=0.,zpar=0.,epar=0.  ;
 856   Int_t relid[4] ;
 857   AliPHOSDigit * digit = 0 ;
 858   Int_t * emcDigits = iniEmc->GetDigitsList() ;
 859
 860   Int_t iparam ;
 861   Int_t iDigit ;
 862   for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
 863     digit = (AliPHOSDigit*) digits->At(emcDigits[iDigit] ) ;
 864     geom->AbsToRelNumbering(digit->GetId(), relid) ;
 865     geom->RelPosInModule(relid, xDigit, zDigit) ;
 866     efit[iDigit] = 0;
 867
 868     iparam = 0 ;
 869     while(iparam < nPar ){
 870       xpar = fitparameters[iparam] ;
 871       zpar = fitparameters[iparam+1] ;
 872       epar = fitparameters[iparam+2] ;
 873       iparam += 3 ;
 874       distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar)  ;
 875       distance =  TMath::Sqrt(distance) ;
 876       efit[iDigit] += epar * ShowerShape(distance) ;
 877     }
 878   }
 879
 880
 881   // Now create new RecPoints and fill energy lists with efit corrected to fluctuations
 882   // so that energy deposited in each cell is distributed betwin new clusters proportionally
 883   // to its contribution to efit
 884
 885   Float_t * emcEnergies = iniEmc->GetEnergiesList() ;
 886   Float_t ratio ;
 887
 888   iparam = 0 ;
 889   while(iparam < nPar ){
 890     xpar = fitparameters[iparam] ;
 891     zpar = fitparameters[iparam+1] ;
 892     epar = fitparameters[iparam+2] ;
 893     iparam += 3 ;
 894
 895     AliPHOSEmcRecPoint * emcRP = 0 ;
 896
 897     if(iniEmc->IsEmc()){ //create new entries in fEmcRecPoints...
 898
 899       if(fNumberOfEmcClusters >= emcRecPoints->GetSize())
 900         emcRecPoints->Expand(2*fNumberOfEmcClusters) ;
 901
 902       (*emcRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint("") ;
 903       emcRP = (AliPHOSEmcRecPoint *) emcRecPoints->At(fNumberOfEmcClusters);
 904       fNumberOfEmcClusters++ ;
 905       emcRP->SetNExMax((Int_t)nPar/3) ;
 906     }
 907     else{//create new entries in fCpvRecPoints
 908       if(fNumberOfCpvClusters >= cpvRecPoints->GetSize())
 909         cpvRecPoints->Expand(2*fNumberOfCpvClusters) ;
 910
 911       (*cpvRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint("") ;
 912       emcRP = (AliPHOSEmcRecPoint *) cpvRecPoints->At(fNumberOfCpvClusters);
 913       fNumberOfCpvClusters++ ;
 914     }
 915
 916     Float_t eDigit ;
 917     for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
 918       digit = (AliPHOSDigit*) digits->At( emcDigits[iDigit] ) ;
 919       geom->AbsToRelNumbering(digit->GetId(), relid) ;
 920       geom->RelPosInModule(relid, xDigit, zDigit) ;
 921       distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar)  ;
 922       distance =  TMath::Sqrt(distance) ;
 923       ratio = epar * ShowerShape(distance) / efit[iDigit] ;
 924       eDigit = emcEnergies[iDigit] * ratio ;
 925       emcRP->AddDigit( *digit, eDigit ) ;
 926     }
 927   }
 928
 929   delete[] fitparameters ;
 930   delete[] efit ;
 931
 932 }
 933
 934 //_____________________________________________________________________________
 935 void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
 936 {
 937   // Calculates the Chi square for the cluster unfolding minimization
 938   // Number of parameters, Gradient, Chi squared, parameters, what to do
 939
 940   TList * toMinuit = (TList*) gMinuit->GetObjectFit() ;
 941
 942   AliPHOSEmcRecPoint * emcRP = (AliPHOSEmcRecPoint*) toMinuit->At(0)  ;
 943   TClonesArray * digits = (TClonesArray*)toMinuit->At(1)  ;
 944
 945
 946
 947   //  AliPHOSEmcRecPoint * emcRP = (AliPHOSEmcRecPoint *) gMinuit->GetObjectFit() ; // EmcRecPoint to fit
 948
 949   Int_t * emcDigits     = emcRP->GetDigitsList() ;
 950
 951   Int_t nOdigits = emcRP->GetDigitsMultiplicity() ;
 952
 953   Float_t * emcEnergies = emcRP->GetEnergiesList() ;
 954
 955   const AliPHOSGeometry * geom = AliPHOSGetter::GetInstance()->PHOSGeometry() ;
 956   fret = 0. ;
 957   Int_t iparam ;
 958
 959   if(iflag == 2)
 960     for(iparam = 0 ; iparam < nPar ; iparam++)
 961       Grad[iparam] = 0 ; // Will evaluate gradient
 962
 963   Double_t efit ;
 964
 965   AliPHOSDigit * digit ;
 966   Int_t iDigit ;
 967
 968   for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
 969
 970     digit = (AliPHOSDigit*) digits->At( emcDigits[iDigit] ) ;
 971
 972     Int_t relid[4] ;
 973     Float_t xDigit ;
 974     Float_t zDigit ;
 975
 976     geom->AbsToRelNumbering(digit->GetId(), relid) ;
 977
 978     geom->RelPosInModule(relid, xDigit, zDigit) ;
 979
 980      if(iflag == 2){  // calculate gradient
 981        Int_t iParam = 0 ;
 982        efit = 0 ;
 983        while(iParam < nPar ){
 984          Double_t distance = (xDigit - x[iParam]) * (xDigit - x[iParam]) ;
 985          iParam++ ;
 986          distance += (zDigit - x[iParam]) * (zDigit - x[iParam]) ;
 987          distance = TMath::Sqrt( distance ) ;
 988          iParam++ ;
 989          efit += x[iParam] * ShowerShape(distance) ;
 990          iParam++ ;
 991        }
 992        Double_t sum=0 ;
 993        if(emcEnergies[iDigit] > 0)
 994         sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ;
 995        // Here we assume, that sigma = sqrt(E)
 996        iParam = 0 ;
 997        while(iParam < nPar ){
 998          Double_t xpar = x[iParam] ;
 999          Double_t zpar = x[iParam+1] ;
1000          Double_t epar = x[iParam+2] ;
1001          Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
1002          Double_t shape = sum * ShowerShape(dr) ;
1003          Double_t r4 = dr*dr*dr*dr ;
1004          Double_t r295 = TMath::Power(dr,2.95) ;
1005          Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) +
1006                                          0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ;
1007
1008          Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ;  // Derivative over x
1009          iParam++ ;
1010          Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ;  // Derivative over z
1011          iParam++ ;
1012          Grad[iParam] += shape ;                                  // Derivative over energy
1013          iParam++ ;
1014        }
1015      }
1016      efit = 0;
1017      iparam = 0 ;
1018
1019      while(iparam < nPar ){
1020        Double_t xpar = x[iparam] ;
1021        Double_t zpar = x[iparam+1] ;
1022        Double_t epar = x[iparam+2] ;
1023        iparam += 3 ;
1024        Double_t distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar)  ;
1025        distance =  TMath::Sqrt(distance) ;
1026        efit += epar * ShowerShape(distance) ;
1027      }
1028
1029      if(emcEnergies[iDigit] > 0) //for the case if rounding error
1030        fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/
1031                emcEnergies[iDigit] ;
1032      // Here we assume, that sigma = sqrt(E)
1033   }
1034
1035 }
1036
1037 //____________________________________________________________________________
1038 void AliPHOSClusterizerv1::Print(Option_t * option)const
1039 {
1040   // Print clusterizer parameters
1041
1042   TString message ;
1043   TString taskName(GetName()) ;
1044   taskName.ReplaceAll(Version(), "") ;
1045
1046   if( strcmp(GetName(), "") !=0 ) {
1047     // Print parameters
1048     message  = "\n--------------- %s %s -----------\n" ;
1049     message += "Clusterizing digits from the file: %s\n" ;
1050     message += "                           Branch: %s\n" ;
1051     message += "                       EMC Clustering threshold = %f\n" ;
1052     message += "                       EMC Local Maximum cut    = %f\n" ;
1053     message += "                       EMC Logarothmic weight   = %f\n" ;
1054     message += "                       CPV Clustering threshold = %f\n" ;
1055     message += "                       CPV Local Maximum cut    = %f\n" ;
1056     message += "                       CPV Logarothmic weight   = %f\n" ;
1057     if(fToUnfold)
1058       message += " Unfolding on\n" ;
1059     else
1060       message += " Unfolding off\n" ;
1061
1062     message += "------------------------------------------------------------------" ;
1063   }
1064   else
1065     message = " AliPHOSClusterizerv1 not initialized " ;
1066
1067   Info("Print", message.Data(),
1068        taskName.Data(),
1069        GetTitle(),
1070        taskName.Data(),
1071        GetName(),
1072        fEmcClusteringThreshold,
1073        fEmcLocMaxCut,
1074        fW0,
1075        fCpvClusteringThreshold,
1076        fCpvLocMaxCut,
1077        fW0CPV ) ;
1078 }
1079
1080
1081 //____________________________________________________________________________
1082 void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option)
1083 {
1084   // Prints list of RecPoints produced at the current pass of AliPHOSClusterizer
1085
1086   TObjArray * emcRecPoints = AliPHOSGetter::GetInstance()->EmcRecPoints(BranchName()) ;
1087   TObjArray * cpvRecPoints = AliPHOSGetter::GetInstance()->CpvRecPoints(BranchName()) ;
1088
1089   if(strstr(option,"deb")){
1090     printf("event %d, found %d EmcRecPoints and %d CPV RecPoints \n",
1091            gAlice->GetEvNumber(),
1092            emcRecPoints->GetEntriesFast(),
1093            cpvRecPoints->GetEntriesFast()) ;
1094   }
1095
1096
1097   if(strstr(option,"all")) {
1098     printf("\nEMC clusters\n" );
1099     printf("Index    Ene(MeV) Multi Module     X     Y   Z  Lambda 1 Lambda 2  # of prim  Primaries list\n" );
1100     Int_t index ;
1101     for (index = 0 ; index < emcRecPoints->GetEntries() ; index++) {
1102       AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )emcRecPoints->At(index) ;
1103       TVector3  locpos;
1104       rp->GetLocalPosition(locpos);
1105       Float_t lambda[2];
1106       rp->GetElipsAxis(lambda);
1107       Int_t * primaries;
1108       Int_t nprimaries;
1109       primaries = rp->GetPrimaries(nprimaries);
1110       printf("\n%6d  %8.2f  %3d     %2d     %4.1f    %4.1f %4.1f %4f  %4f    %2d     : ",
1111               rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetPHOSMod(),
1112               locpos.X(), locpos.Y(), locpos.Z(), lambda[0], lambda[1], nprimaries) ;
1113
1114       for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
1115         printf("%d ", primaries[iprimary] ) ;
1116       }
1117     }
1118     //Now plot CPV recPoints
1119     printf("\n CPV RecPoints \n") ;
1120     printf("Index    Ene(MeV) Module     X     Y   Z  # of prim  Primaries list\n") ;
1121     for (index = 0 ; index < cpvRecPoints->GetEntries() ; index++) {
1122         AliPHOSRecPoint * rp = (AliPHOSRecPoint * )cpvRecPoints->At(index) ;
1123
1124         TVector3  locpos;
1125         rp->GetLocalPosition(locpos);
1126
1127         Int_t * primaries;
1128         Int_t nprimaries ;
1129         primaries = rp->GetPrimaries(nprimaries);
1130         printf("\n%6d  %8.2f  %2d     %4.1f    %4.1f %4.1f %2d     : ",
1131                 rp->GetIndexInList(), rp->GetEnergy(), rp->GetPHOSMod(),
1132                 locpos.X(), locpos.Y(), locpos.Z(), nprimaries) ;
1133         primaries = rp->GetPrimaries(nprimaries);
1134         for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
1135           printf("%d ", primaries[iprimary] ) ;
1136         }
1137       }
1138   }
1139 }
1140