PWG4/PartCorrDep/AliAnaInsideClusterInvariantMass.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 //_________________________________________________________________________
  17 //
  18 // Split clusters with some criteria and calculate invariant mass
  19 // to identify them as pi0 or conversion
  20 //
  21 //
  22 //-- Author: Gustavo Conesa (LPSC-Grenoble)
  23 //_________________________________________________________________________
  24
  25 //////////////////////////////////////////////////////////////////////////////
  26
  27
  28 // --- ROOT system ---
  29 #include <TList.h>
  30 #include <TClonesArray.h>
  31 #include <TObjString.h>
  32 #include <TH3F.h>
  33
  34 // --- Analysis system ---
  35 #include "AliAnaInsideClusterInvariantMass.h"
  36 #include "AliCaloTrackReader.h"
  37 #include "AliMCAnalysisUtils.h"
  38 #include "AliStack.h"
  39 #include "AliFiducialCut.h"
  40 #include "TParticle.h"
  41 #include "AliVCluster.h"
  42 #include "AliAODEvent.h"
  43 #include "AliAODMCParticle.h"
  44 #include "AliEMCALGeoParams.h"
  45
  46 // --- Detectors ---
  47 //#include "AliPHOSGeoUtils.h"
  48 #include "AliEMCALGeometry.h"
  49
  50 ClassImp(AliAnaInsideClusterInvariantMass)
  51
  52 //__________________________________________________________________
  53 AliAnaInsideClusterInvariantMass::AliAnaInsideClusterInvariantMass() :
  54   AliAnaPartCorrBaseClass(),
  55   fCalorimeter(""),
  56   fM02Cut(0),
  57   fMinNCells(0)
  58 {
  59   //default ctor
  60
  61   // Init array of histograms
  62   for(Int_t i = 0; i < 7; i++){
  63     fhMassNLocMax1[i]  = 0;
  64     fhMassNLocMax2[i]  = 0;
  65     fhMassNLocMaxN[i]  = 0;
  66     fhNLocMax[i]       = 0;
  67     fhNLocMaxM02Cut[i] = 0;
  68     fhM02NLocMax1[i]   = 0;
  69     fhM02NLocMax2[i]   = 0;
  70     fhM02NLocMaxN[i]   = 0;
  71     fhNCellNLocMax1[i] = 0;
  72     fhNCellNLocMax2[i] = 0;
  73     fhNCellNLocMaxN[i] = 0;
  74     fhM02Pi0[i]        = 0;
  75     fhM02Eta[i]        = 0;
  76     fhM02Con[i]        = 0;
  77     fhInvMassAllCells[i]=0;
  78   }
  79
  80   InitParameters();
  81
  82 }
  83
  84 //_______________________________________________________________
  85 TObjString *  AliAnaInsideClusterInvariantMass::GetAnalysisCuts()
  86 {
  87         //Save parameters used for analysis
  88   TString parList ; //this will be list of parameters used for this analysis.
  89   const Int_t buffersize = 255;
  90   char onePar[buffersize] ;
  91
  92   snprintf(onePar,buffersize,"--- AliAnaInsideClusterInvariantMass ---\n") ;
  93   parList+=onePar ;
  94
  95   snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
  96   parList+=onePar ;
  97   snprintf(onePar,buffersize,"fM02Cut =%f \n"   ,fM02Cut) ;
  98   parList+=onePar ;
  99   snprintf(onePar,buffersize,"fMinNCells =%f \n",fMinNCells) ;
 100   parList+=onePar ;
 101
 102   return new TObjString(parList) ;
 103
 104 }
 105
 106 //____________________________________________________________________________________________________
 107 Bool_t AliAnaInsideClusterInvariantMass::AreNeighbours( const Int_t absId1, const Int_t absId2 ) const
 108 {
 109   // Tells if (true) or not (false) two digits are neighbours
 110   // A neighbour is defined as being two digits which share a corner
 111
 112   Bool_t areNeighbours = kFALSE ;
 113   Int_t nSupMod =0, nModule =0, nIphi =0, nIeta =0;
 114   Int_t nSupMod1=0, nModule1=0, nIphi1=0, nIeta1=0;
 115   Int_t relid1[2],  relid2[2] ;
 116   Int_t rowdiff=0,  coldiff=0;
 117
 118   areNeighbours = kFALSE ;
 119
 120   GetEMCALGeometry()->GetCellIndex(absId1, nSupMod,nModule,nIphi,nIeta);
 121   GetEMCALGeometry()->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, relid1[0],relid1[1]);
 122
 123   GetEMCALGeometry()->GetCellIndex(absId2, nSupMod1,nModule1,nIphi1,nIeta1);
 124   GetEMCALGeometry()->GetCellPhiEtaIndexInSModule(nSupMod1,nModule1,nIphi1,nIeta1, relid2[0],relid2[1]);
 125
 126   // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2-1
 127   // C Side impair SM, nSupMod%2=1; A side pair SM nSupMod%2=0
 128   if(nSupMod1!=nSupMod){
 129     if(nSupMod1%2) relid1[1]+=AliEMCALGeoParams::fgkEMCALCols;
 130     else           relid2[1]+=AliEMCALGeoParams::fgkEMCALCols;
 131   }
 132
 133   rowdiff = TMath::Abs( relid1[0] - relid2[0] ) ;
 134   coldiff = TMath::Abs( relid1[1] - relid2[1] ) ;
 135
 136   if (( coldiff <= 1 )  && ( rowdiff <= 1 ) && (coldiff + rowdiff > 0))
 137     areNeighbours = kTRUE ;
 138
 139   return areNeighbours;
 140 }
 141
 142 //_____________________________________________________________________________________
 143 TLorentzVector AliAnaInsideClusterInvariantMass::GetCellMomentum(const Int_t absId,
 144                                                                  AliVCaloCells * cells)
 145 {
 146
 147   // Cell momentum calculation for invariant mass
 148
 149   Double_t cellpos[] = {0, 0, 0};
 150   GetEMCALGeometry()->GetGlobal(absId, cellpos);
 151
 152   if(GetVertex(0)){//calculate direction from vertex
 153     cellpos[0]-=GetVertex(0)[0];
 154     cellpos[1]-=GetVertex(0)[1];
 155     cellpos[2]-=GetVertex(0)[2];
 156   }
 157
 158   Double_t r = TMath::Sqrt(cellpos[0]*cellpos[0]+cellpos[1]*cellpos[1]+cellpos[2]*cellpos[2] ) ;
 159
 160   Float_t en = cells->GetCellAmplitude(absId);
 161   RecalibrateCellAmplitude(en,absId);
 162   TLorentzVector cellMom ;
 163   cellMom.SetPxPyPzE( en*cellpos[0]/r,  en*cellpos[1]/r, en*cellpos[2]/r,  en) ;
 164
 165   return cellMom;
 166
 167 }
 168
 169 //________________________________________________________________
 170 TList * AliAnaInsideClusterInvariantMass::GetCreateOutputObjects()
 171 {
 172   // Create histograms to be saved in output file and
 173   // store them in outputContainer
 174   TList * outputContainer = new TList() ;
 175   outputContainer->SetName("InsideClusterHistos") ;
 176
 177   Int_t nptbins  = GetHistoPtBins();           Float_t ptmax  = GetHistoPtMax();           Float_t ptmin  = GetHistoPtMin();
 178   Int_t ssbins   = GetHistoShowerShapeBins();  Float_t ssmax  = GetHistoShowerShapeMax();  Float_t ssmin  = GetHistoShowerShapeMin();
 179   Int_t mbins    = GetHistoMassBins();         Float_t mmax   = GetHistoMassMax();         Float_t mmin   = GetHistoMassMin();
 180   Int_t ncbins   = GetHistoNClusterCellBins(); Int_t   ncmax  = GetHistoNClusterCellMax(); Int_t   ncmin  = GetHistoNClusterCellMin();
 181
 182   TString ptype[] ={"","#gamma","#gamma->e^{#pm}","#pi^{0}","#eta","e^{#pm}", "hadron"};
 183   TString pname[] ={"","Photon","Conversion",     "Pi0",    "Eta", "Electron","Hadron"};
 184
 185   Int_t n = 1;
 186
 187   if(IsDataMC()) n = 7;
 188
 189   for(Int_t i = 0; i < n; i++){
 190
 191     fhInvMassAllCells[i]  = new TH2F(Form("hInvMassAllCells%s",pname[i].Data()),
 192                                   Form("Invariant mass of all cells %s",ptype[i].Data()),
 193                                   nptbins,ptmin,ptmax,mbins,mmin,mmax);
 194     fhInvMassAllCells[i]->SetYTitle("M (MeV/c^2)");
 195     fhInvMassAllCells[i]->SetXTitle("E (GeV)");
 196     outputContainer->Add(fhInvMassAllCells[i]) ;
 197
 198
 199     fhMassNLocMax1[i]  = new TH2F(Form("hMassNLocMax1%s",pname[i].Data()),
 200                                   Form("Invariant mass of 2 highest energy cells %s",ptype[i].Data()),
 201                                   nptbins,ptmin,ptmax,mbins,mmin,mmax);
 202     fhMassNLocMax1[i]->SetYTitle("M (MeV/c^2)");
 203     fhMassNLocMax1[i]->SetXTitle("E (GeV)");
 204     outputContainer->Add(fhMassNLocMax1[i]) ;
 205
 206     fhMassNLocMax2[i]  = new TH2F(Form("hMassNLocMax2%s",pname[i].Data()),
 207                                   Form("Invariant mass of 2 local maxima cells %s",ptype[i].Data()),
 208                                   nptbins,ptmin,ptmax,mbins,mmin,mmax);
 209     fhMassNLocMax2[i]->SetYTitle("M (MeV/c^2)");
 210     fhMassNLocMax2[i]->SetXTitle("E (GeV)");
 211     outputContainer->Add(fhMassNLocMax2[i]) ;
 212
 213     fhMassNLocMaxN[i]  = new TH2F(Form("hMassNLocMaxN%s",pname[i].Data()),
 214                                   Form("Invariant mass of 2 local maxima cells %s",ptype[i].Data()),
 215                                   nptbins,ptmin,ptmax,mbins,mmin,mmax);
 216     fhMassNLocMaxN[i]->SetYTitle("M (MeV/c^2)");
 217     fhMassNLocMax2[i]->SetXTitle("E (GeV)");
 218     outputContainer->Add(fhMassNLocMaxN[i]) ;
 219
 220
 221     fhNLocMax[i]     = new TH2F(Form("hNLocMax%s",pname[i].Data()),
 222                              Form("Number of local maxima in cluster %s",ptype[i].Data()),
 223                              nptbins,ptmin,ptmax,5,0,5);
 224     fhNLocMax[i]   ->SetYTitle("N maxima");
 225     fhNLocMax[i]   ->SetXTitle("E (GeV)");
 226     outputContainer->Add(fhNLocMax[i]) ;
 227
 228     fhNLocMaxM02Cut[i] = new TH2F(Form("hNLocMaxM02Cut%s",pname[i].Data()),
 229                               Form("Number of local maxima in cluster %s for M02 > %2.2f",ptype[i].Data(),fM02Cut),
 230                               nptbins,ptmin,ptmax,5,0,5);
 231     fhNLocMaxM02Cut[i]->SetYTitle("N maxima");
 232     fhNLocMaxM02Cut[i]->SetXTitle("E (GeV)");
 233     outputContainer->Add(fhNLocMaxM02Cut[i]) ;
 234
 235
 236     fhM02NLocMax1[i]     = new TH2F(Form("hM02NLocMax1%s",pname[i].Data()),
 237                                      Form("#lambda_{0}^{2} vs E for N max  = 1 %s",ptype[i].Data()),
 238                                      nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 239     fhM02NLocMax1[i]   ->SetYTitle("#lambda_{0}^{2}");
 240     fhM02NLocMax1[i]   ->SetXTitle("E (GeV)");
 241     outputContainer->Add(fhM02NLocMax1[i]) ;
 242
 243     fhM02NLocMax2[i]     = new TH2F(Form("hM02NLocMax2%s",pname[i].Data()),
 244                                      Form("#lambda_{0}^{2} vs E for N max  = 2 %s",ptype[i].Data()),
 245                                      nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 246     fhM02NLocMax2[i]   ->SetYTitle("#lambda_{0}^{2}");
 247     fhM02NLocMax2[i]   ->SetXTitle("E (GeV)");
 248     outputContainer->Add(fhM02NLocMax2[i]) ;
 249
 250
 251     fhM02NLocMaxN[i]    = new TH2F(Form("hM02NLocMaxN%s",pname[i].Data()),
 252                                    Form("#lambda_{0}^{2} vs E for N max  > 2 %s",ptype[i].Data()),
 253                                    nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 254     fhM02NLocMaxN[i]   ->SetYTitle("#lambda_{0}^{2}");
 255     fhM02NLocMaxN[i]   ->SetXTitle("E (GeV)");
 256     outputContainer->Add(fhM02NLocMaxN[i]) ;
 257
 258
 259     fhNCellNLocMax1[i]  = new TH2F(Form("hNCellNLocMax1%s",pname[i].Data()),
 260                                    Form("#lambda_{0}^{2} vs E for N max  = 1 %s",ptype[i].Data()),
 261                                    nptbins,ptmin,ptmax,ncbins,ncmin,ncmax);
 262     fhNCellNLocMax1[i] ->SetYTitle("N cells");
 263     fhNCellNLocMax1[i] ->SetXTitle("E (GeV)");
 264     outputContainer->Add(fhNCellNLocMax1[i]) ;
 265
 266     fhNCellNLocMax2[i]     = new TH2F(Form("hNCellNLocMax2%s",pname[i].Data()),
 267                                     Form("#lambda_{0}^{2} vs E for N max  = 2 %s",ptype[i].Data()),
 268                                     nptbins,ptmin,ptmax,ncbins,ncmin,ncmax);
 269     fhNCellNLocMax2[i]   ->SetYTitle("N cells");
 270     fhNCellNLocMax2[i]   ->SetXTitle("E (GeV)");
 271     outputContainer->Add(fhNCellNLocMax2[i]) ;
 272
 273
 274     fhNCellNLocMaxN[i]     = new TH2F(Form("hNCellNLocMaxN%s",pname[i].Data()),
 275                                     Form("#lambda_{0}^{2} vs E for N max  > 2 %s",ptype[i].Data()),
 276                                     nptbins,ptmin,ptmax,ncbins,ncmin,ncmax);
 277     fhNCellNLocMaxN[i]   ->SetYTitle("N cells");
 278     fhNCellNLocMaxN[i]   ->SetXTitle("E (GeV)");
 279     outputContainer->Add(fhNCellNLocMaxN[i]) ;
 280
 281
 282     fhM02Pi0[i]     = new TH2F(Form("hM02Pi0%s",pname[i].Data()),
 283                                     Form("#lambda_{0}^{2} vs E ffor mass [0.1-0.2] MeV/c2 %s",ptype[i].Data()),
 284                                     nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 285     fhM02Pi0[i]   ->SetYTitle("#lambda_{0}^{2}");
 286     fhM02Pi0[i]   ->SetXTitle("E (GeV)");
 287     outputContainer->Add(fhM02Pi0[i]) ;
 288
 289     fhM02Eta[i]     = new TH2F(Form("hM02Eta%s",pname[i].Data()),
 290                                     Form("#lambda_{0}^{2} vs E for mass [0.4-0.7] MeV/c2, %s",ptype[i].Data()),
 291                                     nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 292     fhM02Eta[i]   ->SetYTitle("#lambda_{0}^{2}");
 293     fhM02Eta[i]   ->SetXTitle("E (GeV)");
 294     outputContainer->Add(fhM02Eta[i]) ;
 295
 296
 297     fhM02Con[i]    = new TH2F(Form("hM02Con%s",pname[i].Data()),
 298                                    Form("#lambda_{0}^{2} vs E for mass < 0.5 MeV/c2, %s",ptype[i].Data()),
 299                                    nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 300     fhM02Con[i]   ->SetYTitle("#lambda_{0}^{2}");
 301     fhM02Con[i]   ->SetXTitle("E (GeV)");
 302     outputContainer->Add(fhM02Con[i]) ;
 303
 304   }
 305
 306   return outputContainer ;
 307
 308 }
 309
 310 //________________________________________________________________________________________________________
 311 Int_t AliAnaInsideClusterInvariantMass::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCaloCells* cells,
 312                                                                Int_t *absIdList,     Float_t *maxEList)
 313 {
 314   // Find local maxima in cluster
 315
 316   Float_t locMaxCut = 0; // not used for the moment
 317
 318   Int_t iDigitN = 0 ;
 319   Int_t iDigit  = 0 ;
 320   Int_t absId1 = -1 ;
 321   Int_t absId2 = -1 ;
 322   const Int_t nCells = cluster->GetNCells();
 323   //printf("cluster :");
 324   for(iDigit = 0; iDigit < nCells ; iDigit++){
 325     absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit]  ;
 326
 327     //Float_t en = cells->GetCellAmplitude(absIdList[iDigit]);
 328     //RecalibrateCellAmplitude(en,absIdList[iDigit]);
 329     //Int_t icol = -1, irow = -1, iRCU = -1;
 330     //Int_t sm = GetCaloUtils()->GetModuleNumberCellIndexes(absIdList[iDigit], fCalorimeter, icol, irow, iRCU) ;
 331
 332     //printf("\t cell %d, id %d, sm %d, col %d, row %d, e %f\n", iDigit, absIdList[iDigit], sm, icol, irow, en );
 333   }
 334
 335   for(iDigit = 0 ; iDigit < nCells; iDigit++) {
 336     if(absIdList[iDigit]>=0) {
 337
 338       absId1 = absIdList[iDigit] ;
 339       Float_t en1 = cells->GetCellAmplitude(absId1);
 340       RecalibrateCellAmplitude(en1,absId1);
 341
 342       for(iDigitN = 0; iDigitN < nCells; iDigitN++) {
 343         absId2 = absIdList[iDigitN] ;
 344
 345         Float_t en2 = cells->GetCellAmplitude(absId2);
 346         RecalibrateCellAmplitude(en2,absId2);
 347
 348         if ( AreNeighbours(absId1, absId2) ) {
 349
 350           if (en1 > en2 ) {
 351             absIdList[iDigitN] = -1 ;
 352             // but may be digit too is not local max ?
 353             if(en1 < en2 + locMaxCut)
 354               absIdList[iDigit] = -1 ;
 355           }
 356           else {
 357             absIdList[iDigit] = -1 ;
 358             // but may be digitN too is not local max ?
 359             if(en1 > en2 - locMaxCut)
 360               absIdList[iDigitN] = -1 ;
 361           }
 362         } // if Areneighbours
 363       } // while digitN
 364     } // slot not empty
 365   } // while digit
 366
 367   iDigitN = 0 ;
 368   for(iDigit = 0; iDigit < nCells; iDigit++) {
 369     if(absIdList[iDigit]>=0 ){
 370       absIdList[iDigitN] = absIdList[iDigit] ;
 371       Float_t en = cells->GetCellAmplitude(absIdList[iDigit]);
 372       RecalibrateCellAmplitude(en,absIdList[iDigit]);
 373       if(en < 0.1) continue; // Maxima only with seed energy at least
 374       maxEList[iDigitN] = en ;
 375       //printf("Local max %d, id %d, en %f\n", iDigit,absIdList[iDigitN],en);
 376       iDigitN++ ;
 377     }
 378   }
 379
 380   return iDigitN ;
 381
 382 }
 383
 384
 385 //___________________________________________
 386 void AliAnaInsideClusterInvariantMass::Init()
 387 {
 388   //Init
 389   //Do some checks
 390   if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
 391     printf("AliAnaInsideClusterInvariantMass::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
 392     abort();
 393   }
 394   else  if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
 395     printf("AliAnaInsideClusterInvariantMass::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
 396     abort();
 397   }
 398
 399   if( GetReader()->GetDataType() == AliCaloTrackReader::kMC ){
 400     printf("AliAnaInsideClusterInvariantMass::Init() - !!STOP: You want to use pure MC data!!\n");
 401     abort();
 402
 403   }
 404
 405 }
 406
 407 //_____________________________________________________
 408 void AliAnaInsideClusterInvariantMass::InitParameters()
 409 {
 410   //Initialize the parameters of the analysis.
 411   AddToHistogramsName("AnaPi0InsideClusterInvariantMass_");
 412
 413   fCalorimeter = "EMCAL" ;
 414   fM02Cut      = 0.26 ;
 415   fMinNCells   = 4 ;
 416 }
 417
 418
 419 //__________________________________________________________________
 420 void  AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms()
 421 {
 422   //Search for pi0 in fCalorimeter with shower shape analysis
 423
 424   TObjArray * pl       = 0x0;
 425   AliVCaloCells* cells = 0x0;
 426
 427   //Select the Calorimeter of the photon
 428   if(fCalorimeter == "PHOS"){
 429     pl    = GetPHOSClusters();
 430     cells = GetPHOSCells();
 431   }
 432   else if (fCalorimeter == "EMCAL"){
 433     pl    = GetEMCALClusters();
 434     cells = GetEMCALCells();
 435   }
 436
 437   if(!pl || !cells) {
 438     Info("MakeAnalysisFillHistograms","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
 439     return;
 440   }
 441
 442         if(fCalorimeter == "PHOS") return; // Not implemented for PHOS yet
 443
 444   for(Int_t icluster = 0; icluster < pl->GetEntriesFast(); icluster++){
 445     AliVCluster * cluster = (AliVCluster*) (pl->At(icluster));
 446
 447     // Study clusters with large shape parameter
 448     Float_t en = cluster->E();
 449     Float_t l0 = cluster->GetM02();
 450     Int_t   nc = cluster->GetNCells();
 451
 452     //If too small or big E or low number of cells, skip it
 453     if( ( en < GetMinEnergy() || en > GetMaxEnergy() ) && nc < fMinNCells) continue ;
 454
 455     Int_t    absId1    = -1; Int_t absId2 = -1;
 456     Int_t   *absIdList = new Int_t  [nc];
 457     Float_t *maxEList  = new Float_t[nc];
 458     Int_t    nMax      = GetNumberOfLocalMaxima(cluster, cells, absIdList, maxEList) ;
 459
 460     if (nMax <= 0) {
 461       printf("AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms() - No local maximum found!");
 462       delete [] absIdList ;
 463       delete [] maxEList  ;
 464       return;
 465     }
 466
 467     fhNLocMax[0]->Fill(en,nMax);
 468
 469     if     ( nMax == 1  ) { fhM02NLocMax1[0]->Fill(en,l0) ; fhNCellNLocMax1[0]->Fill(en,nc) ; }
 470     else if( nMax == 2  ) { fhM02NLocMax2[0]->Fill(en,l0) ; fhNCellNLocMax2[0]->Fill(en,nc) ; }
 471     else if( nMax >= 3  ) { fhM02NLocMaxN[0]->Fill(en,l0) ; fhNCellNLocMaxN[0]->Fill(en,nc) ; }
 472     else printf("N max smaller than 1 -> %d \n",nMax);
 473
 474     // Play with the MC stack if available
 475     // Check origin of the candidates
 476     Int_t mcindex = -1;
 477     if(IsDataMC()){
 478
 479       Int_t tag = GetMCAnalysisUtils()->CheckOrigin(cluster->GetLabels(),cluster->GetNLabels(), GetReader(), 0);
 480
 481       if      ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)  )      mcindex = kmcPi0;
 482       else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)  )      mcindex = kmcEta;
 483       else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
 484                !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) mcindex = kmcPhoton;
 485       else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
 486                 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) mcindex = kmcConversion;
 487       else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron))   mcindex = kmcElectron;
 488       else                                                                                mcindex = kmcHadron;
 489
 490 /*      printf("AliAnaInsideClusterInvariantMass::FillAnalysisMakeHistograms() - tag %d, photon %d, prompt %d, frag %d, isr %d, pi0 decay %d, eta decay %d, other decay %d \n conv %d, pi0 %d, hadron %d, electron %d, unk %d, muon %d,pion %d, proton %d, neutron %d, kaon %d, antiproton %d, antineutron %d, bad %d\n",tag,
 491              GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton),
 492              GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt),
 493              GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation),
 494              GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR),
 495              GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay),
 496              GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay),
 497              GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay),
 498
 499              GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion),
 500              GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0),
 501              GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOther),
 502              GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron),
 503              GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCUnknown),
 504              GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCMuon),
 505              GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPion),
 506              GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCProton),
 507              GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron),
 508              GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCKaon),
 509              GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton),
 510              GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron),
 511              GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCBadLabel)
 512 );
 513 */
 514
 515       fhNLocMax[mcindex]->Fill(en,nMax);
 516
 517       if     (nMax == 1 ) { fhM02NLocMax1[mcindex]->Fill(en,l0) ; fhNCellNLocMax1[mcindex]->Fill(en,nc) ; }
 518       else if(nMax == 2 ) { fhM02NLocMax2[mcindex]->Fill(en,l0) ; fhNCellNLocMax2[mcindex]->Fill(en,nc) ; }
 519       else if(nMax >= 3 ) { fhM02NLocMaxN[mcindex]->Fill(en,l0) ; fhNCellNLocMaxN[mcindex]->Fill(en,nc) ; }
 520
 521     }
 522
 523     if( l0 < fM02Cut) {
 524       delete [] absIdList ;
 525       delete [] maxEList  ;
 526       continue;
 527     }
 528
 529     // Get the 2 max indeces and do inv mass
 530
 531     absId1 = absIdList[0];
 532     TLorentzVector cellMomi = GetCellMomentum(absId1, cells);
 533
 534     if     ( nMax == 2 ) absId2 = absIdList[1];
 535     else if( nMax == 1 ){
 536       //Find second highest energy cell
 537       Float_t enmax = 0 ;
 538       for(Int_t iDigit = 0 ; iDigit < cluster->GetNCells() ; iDigit++){
 539         Int_t absId = cluster->GetCellsAbsId()[iDigit];
 540         if( absId == absId1 ) continue ;
 541         Float_t endig = cells->GetCellAmplitude(absId);
 542         RecalibrateCellAmplitude(endig,absId);
 543         if(endig > enmax) {
 544           enmax  = endig ;
 545           absId2 = absId ;
 546         }
 547
 548         //Get mass of all cell in cluster combinations
 549
 550
 551         Float_t enj = cells->GetCellAmplitude(absId);
 552         RecalibrateCellAmplitude(enj,absId);
 553
 554         if(enj < 0.3) continue;
 555
 556         TLorentzVector cellMomj = GetCellMomentum(absId, cells);
 557
 558         fhInvMassAllCells[0]->Fill(en,(cellMomj+cellMomi).M());
 559
 560         if(IsDataMC()) fhInvMassAllCells[mcindex]->Fill(en,(cellMomj+cellMomi).M());
 561
 562       }// cell loop
 563
 564     }
 565     else { // loop on maxima, find 2 highest
 566
 567       Int_t enmax = 0 ;
 568       for(Int_t iDigit = 0 ; iDigit < nMax ; iDigit++){
 569         Float_t endig = maxEList[iDigit];
 570         if(endig > enmax) {
 571           endig  = enmax ;
 572           absId2 = absIdList[iDigit];
 573         }
 574       }// maxima loop
 575
 576       // First max is not highest, check if there is other higher
 577       if(maxEList[0] < enmax){
 578
 579         for(Int_t iDigit = 0 ; iDigit < nMax ; iDigit++){
 580           if(absId2 == absIdList[iDigit]) continue;
 581           Float_t endig = maxEList[iDigit];
 582           if(endig > enmax) {
 583             endig  = enmax ;
 584             absId1 = absIdList[iDigit];
 585           }
 586         }// maxima loop
 587
 588       }
 589
 590     }
 591
 592     fhNLocMaxM02Cut[0]->Fill(en,nMax);
 593
 594     TLorentzVector cellMom1 = GetCellMomentum(absId1, cells);
 595     TLorentzVector cellMom2 = GetCellMomentum(absId2, cells);
 596
 597     Float_t mass = (cellMom1+cellMom2).M();
 598
 599     if     (nMax==1) fhMassNLocMax1[0]->Fill(en,mass);
 600     else if(nMax==2) fhMassNLocMax2[0]->Fill(en,mass);
 601     else if(nMax >2) fhMassNLocMaxN[0]->Fill(en,mass);
 602
 603     if     (mass < 0.1)               fhM02Con[0]->Fill(en,l0);
 604     else if(mass < 0.2)               fhM02Pi0[0]->Fill(en,l0);
 605     else if(mass < 0.6 && mass > 0.4) fhM02Eta[0]->Fill(en,l0);
 606
 607     if(IsDataMC()){
 608
 609       fhNLocMaxM02Cut[mcindex]->Fill(en,nMax);
 610
 611       if     (nMax==1) fhMassNLocMax1[mcindex]->Fill(en,mass);
 612       else if(nMax==2) fhMassNLocMax2[mcindex]->Fill(en,mass);
 613       else if(nMax >2) fhMassNLocMaxN[mcindex]->Fill(en,mass);
 614
 615       if     (mass < 0.1)               fhM02Con[mcindex]->Fill(en,l0);
 616       else if(mass < 0.2)               fhM02Pi0[mcindex]->Fill(en,l0);
 617       else if(mass < 0.6 && mass > 0.4) fhM02Eta[mcindex]->Fill(en,l0);
 618
 619
 620     }//Work with MC truth first
 621
 622     delete [] absIdList ;
 623     delete [] maxEList  ;
 624
 625   }//loop
 626
 627   if(GetDebug() > 1) printf("AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms() - END \n");
 628
 629 }
 630
 631 //______________________________________________________________________
 632 void AliAnaInsideClusterInvariantMass::Print(const Option_t * opt) const
 633 {
 634   //Print some relevant parameters set for the analysis
 635   if(! opt)
 636     return;
 637
 638   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
 639   AliAnaPartCorrBaseClass::Print("");
 640   printf("Calorimeter     =     %s\n", fCalorimeter.Data()) ;
 641   printf("lambda 0 sqared >  %2.1f\n", fM02Cut);
 642   printf("    \n") ;
 643
 644 }
 645
 646 //____________________________________________________________________________________________
 647 void AliAnaInsideClusterInvariantMass::RecalibrateCellAmplitude(Float_t & amp, const Int_t id)
 648 {
 649   //Recaculate cell energy if recalibration factor
 650
 651   Int_t icol     = -1; Int_t irow     = -1; Int_t iRCU     = -1;
 652   Int_t nModule  = GetModuleNumberCellIndexes(id,fCalorimeter, icol, irow, iRCU);
 653
 654   if (GetCaloUtils()->IsRecalibrationOn()) {
 655     if(fCalorimeter == "PHOS") {
 656       amp *= GetCaloUtils()->GetPHOSChannelRecalibrationFactor(nModule,icol,irow);
 657     }
 658     else                                   {
 659       amp *= GetCaloUtils()->GetEMCALChannelRecalibrationFactor(nModule,icol,irow);
 660     }
 661   }
 662 }
 663
 664 //________________________________________________________________________________________
 665 void AliAnaInsideClusterInvariantMass::SplitEnergy(const Int_t absId1, const Int_t absId2,
 666                                                    const AliVCaloCells* cells,
 667                                                    Float_t & e1, Float_t & e2 )
 668 {
 669
 670   // Split energy of cluster between the 2 local maxima.
 671
 672   Int_t icol1 = -1, irow1 = -1, iRCU1 = -1;
 673   Int_t sm1 = GetCaloUtils()->GetModuleNumberCellIndexes(absId1, fCalorimeter, icol1, irow1, iRCU1) ;
 674   Int_t icol2 = -1, irow2 = -1, iRCU2 = -1;
 675   Int_t sm2 = GetCaloUtils()->GetModuleNumberCellIndexes(absId2, fCalorimeter, icol2, irow2, iRCU2) ;
 676
 677   if(sm1!=sm2) {
 678     if(sm1%2) icol1+=AliEMCALGeoParams::fgkEMCALCols;
 679     else      icol2+=AliEMCALGeoParams::fgkEMCALCols;
 680   }
 681
 682   // just to avoid compilation warning
 683   Int_t nTotCells = cells->GetNumberOfCells();
 684   if(GetDebug() > 2)printf("N cells %d, e1 %f, e2 %f\n", nTotCells,e1, e2);
 685 /// continue here
 686
 687 }
 688
 689