]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGGA/CaloTrackCorrelations/AliAnaElectron.cxx
Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / PWGGA / CaloTrackCorrelations / AliAnaElectron.cxx
index 8604f45758e16f101a20b7284a7435c7f4e5d3e4..e7008ac2fdab6cbcfc9777a3998ad2d7448cf379 100755 (executable)
@@ -50,7 +50,7 @@ ClassImp(AliAnaElectron)
   
 //________________________________
 AliAnaElectron::AliAnaElectron() : 
-    AliAnaCaloTrackCorrBaseClass(),       fCalorimeter(""), 
+    AliAnaCaloTrackCorrBaseClass(),
     fMinDist(0.),                         fMinDist2(0.),                         fMinDist3(0.), 
     fTimeCutMin(-1),                      fTimeCutMax(999999),         
     fNCellsCut(0),                        fNLMCutMin(-1),                        fNLMCutMax(10),
@@ -59,6 +59,7 @@ AliAnaElectron::AliAnaElectron() :
     fdEdxMin(0.),                         fdEdxMax (200.), 
     fEOverPMin(0),                        fEOverPMax (2),
     fAODParticle(0),
+    fMomentum(),                          fMomentumMC(),                         fProdVertex(),
     // Histograms
     fhdEdxvsE(0),                         fhdEdxvsP(0),                 
     fhEOverPvsE(0),                       fhEOverPvsP(0),
@@ -158,49 +159,50 @@ AliAnaElectron::AliAnaElectron() :
 }
 
 //____________________________________________________________________________
-Bool_t  AliAnaElectron::ClusterSelected(AliVCluster* calo, TLorentzVector mom, Int_t nMaxima)
+Bool_t  AliAnaElectron::ClusterSelected(AliVCluster* calo, Int_t nMaxima)
 {
-  //Select clusters if they pass different cuts
-  if(GetDebug() > 2) 
-    printf("AliAnaElectron::ClusterSelected() Current Event %d; Before selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f\n",
-           GetReader()->GetEventNumber(),
-           mom.E(), mom.Pt(),calo->E(),mom.Phi()*TMath::RadToDeg(),mom.Eta());
+  // Select clusters if they pass different cuts
+  
+  AliDebug(1,Form("Current Event %d; Before selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f",
+           GetReader()->GetEventNumber(),fMomentum.E(),fMomentum.Pt(),calo->E(),fMomentum.Phi()*TMath::RadToDeg(),fMomentum.Eta()));
   
   //.......................................
   //If too small or big energy, skip it
-  if(mom.E() < GetMinEnergy() || mom.E() > GetMaxEnergy() ) return kFALSE ; 
-  if(GetDebug() > 2) printf("\t Cluster %d Pass E Cut \n",calo->GetID());
+  if(fMomentum.E() < GetMinEnergy() || fMomentum.E() > GetMaxEnergy() ) return kFALSE ; 
+  AliDebug(2,Form("\t Cluster %d Pass E Cut",calo->GetID()));
   
   //.......................................
   // TOF cut, BE CAREFUL WITH THIS CUT
   Double_t tof = calo->GetTOF()*1e9;
   if(tof < fTimeCutMin || tof > fTimeCutMax) return kFALSE;
-  if(GetDebug() > 2)  printf("\t Cluster %d Pass Time Cut \n",calo->GetID());
+  AliDebug(2,Form("\t Cluster %d Pass Time Cut",calo->GetID()));
   
   //.......................................
   if(calo->GetNCells() <= fNCellsCut && GetReader()->GetDataType() != AliCaloTrackReader::kMC) return kFALSE;
-  if(GetDebug() > 2) printf("\t Cluster %d Pass NCell Cut \n",calo->GetID());
+  AliDebug(2,Form("\t Cluster %d Pass NCell Cut",calo->GetID()));
   
   //.......................................
   //Check acceptance selection
-  if(IsFiducialCutOn()){
-    Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
+  if(IsFiducialCutOn())
+  {
+    Bool_t in = GetFiducialCut()->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),GetCalorimeter()) ;
     if(! in ) return kFALSE ;
   }
-  if(GetDebug() > 2) printf("Fiducial cut passed \n");
+  AliDebug(2,"\t Fiducial cut passed");
   
   //.......................................
   //Skip not matched clusters with tracks
-  if(!IsTrackMatched(calo, GetReader()->GetInputEvent())) {
-      if(GetDebug() > 2) printf("\t Reject non track-matched clusters\n");
+  if(!IsTrackMatched(calo, GetReader()->GetInputEvent()))
+  {
+      AliDebug(1,"\t Reject non track-matched clusters");
       return kFALSE ;
   }
-  else if(GetDebug() > 2)  printf(" Track-matching cut passed \n");
+  else AliDebug(2,"\t Track-matching cut passed");
   
   //...........................................
   // skip clusters with too many maxima
   if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax) return kFALSE ;
-  if(GetDebug() > 2) printf(" \t Cluster %d pass NLM %d of out of range \n",calo->GetID(), nMaxima);
+  AliDebug(2,Form("\t Cluster %d pass NLM %d of out of range",calo->GetID(), nMaxima));
   
   //.......................................
   //Check Distance to Bad channel, set bit.
@@ -209,24 +211,22 @@ Bool_t  AliAnaElectron::ClusterSelected(AliVCluster* calo, TLorentzVector mom, I
   if(distBad < fMinDist) {//In bad channel (PHOS cristal size 2.2x2.2 cm), EMCAL ( cell units )
     return kFALSE ;
   }
-  else if(GetDebug() > 2) printf("\t Bad channel cut passed %4.2f > %2.2f \n",distBad, fMinDist);
+  else AliDebug(2,Form("\t Bad channel cut passed %4.2f > %2.2f",distBad, fMinDist));
   //printf("Cluster %d Pass Bad Dist Cut \n",icalo);
 
-  if(GetDebug() > 0) 
-    printf("AliAnaElectron::ClusterSelected() Current Event %d; After  selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f\n",
+ AliDebug(1,Form("Current Event %d; After  selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f",
            GetReader()->GetEventNumber(), 
-           mom.E(), mom.Pt(),calo->E(),mom.Phi()*TMath::RadToDeg(),mom.Eta());
+           fMomentum.E(), fMomentum.Pt(),calo->E(),fMomentum.Phi()*TMath::RadToDeg(),fMomentum.Eta()));
   
   //All checks passed, cluster selected
   return kTRUE;
     
 }
 
-//__________________________________________________________________________________________________________
-void  AliAnaElectron::FillShowerShapeHistograms(AliVCluster* cluster, const Int_t mcTag, const Int_t pidTag)
+//______________________________________________________________________________________________
+void  AliAnaElectron::FillShowerShapeHistograms(AliVCluster* cluster, Int_t mcTag, Int_t pidTag)
 {
-  
-  //Fill cluster Shower Shape histograms
+  // Fill cluster Shower Shape histograms
   
   if(!fFillSSHistograms || GetMixedEvent()) return;
   
@@ -243,26 +243,18 @@ void  AliAnaElectron::FillShowerShapeHistograms(AliVCluster* cluster, const Int_
   
   Float_t l0   = 0., l1   = 0.;
   Float_t dispp= 0., dEta = 0., dPhi    = 0.; 
-  Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;    
+  Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
   
-  TLorentzVector mom;
-  if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
-  {
-    cluster->GetMomentum(mom,GetVertex(0)) ;}//Assume that come from vertex in straight line
-  else{
-    Double_t vertex[]={0,0,0};
-    cluster->GetMomentum(mom,vertex) ;
-  }
-  
-  Float_t eta = mom.Eta();
-  Float_t phi = mom.Phi();
+  Float_t eta = fMomentum.Eta();
+  Float_t phi = fMomentum.Phi();
   if(phi < 0) phi+=TMath::TwoPi();
   
   fhLam0E[pidIndex] ->Fill(energy,lambda0);
   fhLam1E[pidIndex] ->Fill(energy,lambda1);
   fhDispE[pidIndex] ->Fill(energy,disp);
   
-  if(fCalorimeter == "EMCAL" && GetModuleNumber(cluster) > 5)
+  if(GetCalorimeter() == kEMCAL &&  GetFirstSMCoveredByTRD() >= 0 &&
+     GetModuleNumber(cluster) >= GetFirstSMCoveredByTRD() )
   {
     fhLam0ETRD[pidIndex]->Fill(energy,lambda0);
     fhLam1ETRD[pidIndex]->Fill(energy,lambda1);
@@ -284,7 +276,7 @@ void  AliAnaElectron::FillShowerShapeHistograms(AliVCluster* cluster, const Int_
       fhPhiLam0HighE[pidIndex]   ->Fill(phi,   lambda0);
     }
     
-    if(fCalorimeter == "EMCAL")
+    if(GetCalorimeter() == kEMCAL)
     {
       GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), GetReader()->GetInputEvent()->GetEMCALCells(), cluster,
                                                                                    l0, l1, dispp, dEta, dPhi, sEta, sPhi, sEtaPhi);
@@ -308,8 +300,8 @@ void  AliAnaElectron::FillShowerShapeHistograms(AliVCluster* cluster, const Int_
   if(IsDataMC())
   {
     AliVCaloCells* cells = 0;
-    if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
-    else                        cells = GetPHOSCells();
+    if(GetCalorimeter() == kEMCAL) cells = GetEMCALCells();
+    else                           cells = GetPHOSCells();
     
     //Fill histograms to check shape of embedded clusters
     Float_t fraction = 0;
@@ -326,20 +318,16 @@ void  AliAnaElectron::FillShowerShapeHistograms(AliVCluster* cluster, const Int_
       //Fraction of total energy due to the embedded signal
       fraction/=clusterE;
       
-      if(GetDebug() > 1 ) printf("AliAnaElectron::FillShowerShapeHistogram() - Energy fraction of embedded signal %2.3f, Energy %2.3f\n",fraction, clusterE);
+      AliDebug(1,Form("Energy fraction of embedded signal %2.3f, Energy %2.3f",fraction, clusterE));
       
       fhEmbeddedSignalFractionEnergy->Fill(clusterE,fraction);
       
     }  // embedded fraction    
     
-    // Get the fraction of the cluster energy that carries the cell with highest energy
-    Int_t absID             =-1 ;
-    Float_t maxCellFraction = 0.;
-    Int_t index             = 0 ;
-    absID = GetCaloUtils()->GetMaxEnergyCell(cells, cluster,maxCellFraction);
-    
     // Check the origin and fill histograms
-    if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton) && 
+    Int_t index = -1;
+
+    if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton) &&
        !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) &&
        !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0) &&
        !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta))
@@ -358,11 +346,12 @@ void  AliAnaElectron::FillShowerShapeHistograms(AliVCluster* cluster, const Int_
         
         //Compare the primary depositing more energy with the rest, if no photon/electron as comon ancestor (conversions), count as other particle
         Int_t ancPDG = 0, ancStatus = -1;
-        TLorentzVector momentum; TVector3 prodVertex;
         Int_t ancLabel = 0;
         Int_t noverlaps = 1;      
-        for (UInt_t ilab = 0; ilab < cluster->GetNLabels(); ilab++ ) {
-          ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(cluster->GetLabels()[0],cluster->GetLabels()[ilab], GetReader(),ancPDG,ancStatus,momentum,prodVertex);
+        for (UInt_t ilab = 0; ilab < cluster->GetNLabels(); ilab++ )
+        {
+          ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(cluster->GetLabels()[0],cluster->GetLabels()[ilab], GetReader(),
+                                                               ancPDG,ancStatus,fMomentumMC,fProdVertex);
           if(ancPDG!=22 && TMath::Abs(ancPDG)!=11) noverlaps++;
         }
         
@@ -375,8 +364,9 @@ void  AliAnaElectron::FillShowerShapeHistograms(AliVCluster* cluster, const Int_
         else if(noverlaps > 2){          
           fhMCElectronELambda0NOverlap   ->Fill(energy, lambda0);
         }
-        else {
-          printf("AliAnaElectron::FillShowerShapeHistogram() - n overlaps = %d!!", noverlaps);
+        else
+        {
+          AliWarning(Form("N overlaps = %d for ancestor %d!!", noverlaps, ancLabel));
         }
       }//No embedding
       
@@ -421,17 +411,15 @@ void  AliAnaElectron::FillShowerShapeHistograms(AliVCluster* cluster, const Int_
     
     fhMCELambda0[pidIndex][index]    ->Fill(energy, lambda0);
     
-    if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
+    if(GetCalorimeter() == kEMCAL && !fFillOnlySimpleSSHisto)
     {
       fhMCEDispEta        [pidIndex][index]-> Fill(energy,dEta);
       fhMCEDispPhi        [pidIndex][index]-> Fill(energy,dPhi);
       fhMCESumEtaPhi      [pidIndex][index]-> Fill(energy,sEtaPhi);
       fhMCEDispEtaPhiDiff [pidIndex][index]-> Fill(energy,dPhi-dEta);
-      if(dEta+dPhi>0)fhMCESphericity     [pidIndex][index]-> Fill(energy,(dPhi-dEta)/(dEta+dPhi));  
-      
+      if(dEta+dPhi>0) fhMCESphericity[pidIndex][index]-> Fill(energy,(dPhi-dEta)/(dEta+dPhi));
     }
     
-    
   }//MC data
   
 }
@@ -444,19 +432,19 @@ TObjString *  AliAnaElectron::GetAnalysisCuts()
   const Int_t buffersize = 255;
   char onePar[buffersize] ;
   
-  snprintf(onePar,buffersize,"--- AliAnaElectron ---\n") ;
+  snprintf(onePar,buffersize,"--- AliAnaElectron ---") ;
   parList+=onePar ;    
-  snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
+  snprintf(onePar,buffersize,"Calorimeter: %s;",GetCalorimeterString().Data()) ;
   parList+=onePar ;
-  snprintf(onePar,buffersize," %2.2f < dEdx < %2.2f  \n",fdEdxMin,fdEdxMax) ;
+  snprintf(onePar,buffersize," %2.2f < dEdx < %2.2f;",fdEdxMin,fdEdxMax) ;
   parList+=onePar ;  
-  snprintf(onePar,buffersize," %2.2f <  E/P < %2.2f  \n",fEOverPMin, fEOverPMax) ;
+  snprintf(onePar,buffersize," %2.2f <  E/P < %2.2f;",fEOverPMin, fEOverPMax) ;
   parList+=onePar ;  
-  snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
+  snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster);",fMinDist) ;
   parList+=onePar ;
-  snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
+  snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation);",fMinDist2) ;
   parList+=onePar ;
-  snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
+  snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study);",fMinDist3) ;
   parList+=onePar ;  
   
   //Get parameters set in base class.
@@ -669,7 +657,7 @@ TList *  AliAnaElectron::GetCreateOutputObjects()
       fhDispE[pidIndex]->SetXTitle("E (GeV) ");
       outputContainer->Add(fhDispE[pidIndex]);
       
-      if(fCalorimeter == "EMCAL")
+      if(GetCalorimeter() == kEMCAL &&  GetFirstSMCoveredByTRD() >=0 )
       {
         fhLam0ETRD[pidIndex]  = new TH2F (Form("h%sLam0ETRD",pidParticle[pidIndex].Data()),
                                           Form("%s: #lambda_{0}^{2} vs E, EMCAL SM covered by TRD",pidParticle[pidIndex].Data()), 
@@ -739,7 +727,7 @@ TList *  AliAnaElectron::GetCreateOutputObjects()
         fhPhiLam0HighE[pidIndex]->SetXTitle("#phi");
         outputContainer->Add(fhPhiLam0HighE[pidIndex]);  
         
-        if(fCalorimeter == "EMCAL")
+        if(GetCalorimeter() == kEMCAL)
         {
           fhDispEtaE[pidIndex]  = new TH2F (Form("h%sDispEtaE",pidParticle[pidIndex].Data()),
                                             Form("%s: #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E",pidParticle[pidIndex].Data()),  
@@ -817,7 +805,7 @@ TList *  AliAnaElectron::GetCreateOutputObjects()
           fhMCELambda0[pidIndex][i]->SetXTitle("E (GeV)");
           outputContainer->Add(fhMCELambda0[pidIndex][i]) ; 
           
-          if(fCalorimeter=="EMCAL" && !fFillOnlySimpleSSHisto)
+          if(GetCalorimeter()==kEMCAL && !fFillOnlySimpleSSHisto)
           {
             fhMCEDispEta[pidIndex][i]  = new TH2F (Form("h%sEDispEtaE_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
                                                    Form("cluster from %s : %s like, #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E",ptypess[i].Data(),pidParticle[pidIndex].Data()),
@@ -1066,17 +1054,14 @@ TList *  AliAnaElectron::GetCreateOutputObjects()
 //_________________________
 void AliAnaElectron::Init()
 {
+  // Init
   
-  //Init
-  //Do some checks
-  if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
-    printf("AliAnaElectron::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
-    abort();
-  }
-  else  if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
-    printf("AliAnaElectron::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
-    abort();
-  }
+  // Do some checks
+  if       ( GetCalorimeter() == kPHOS  && !GetReader()->IsPHOSSwitchedOn()  && NewOutputAOD() )
+    AliFatal("STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!");
+  else  if ( GetCalorimeter() == kEMCAL && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD() )
+    AliFatal("STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!");
   
 }
 
@@ -1087,7 +1072,6 @@ void AliAnaElectron::InitParameters()
   //Initialize the parameters of the analysis.
   AddToHistogramsName("AnaElectron_");
   
-  fCalorimeter = "EMCAL" ;
   fMinDist     = 2.;
   fMinDist2    = 4.;
   fMinDist3    = 5.;
@@ -1115,88 +1099,69 @@ void  AliAnaElectron::MakeAnalysisFillAOD()
   
   //Select the Calorimeter of the photon
   TObjArray * pl = 0x0; 
-  if(fCalorimeter == "PHOS")
-    pl = GetPHOSClusters();
-  else if (fCalorimeter == "EMCAL")
-    pl = GetEMCALClusters();
+  if      (GetCalorimeter() == kPHOS ) pl = GetPHOSClusters ();
+  else if (GetCalorimeter() == kEMCAL) pl = GetEMCALClusters();
   
-  if(!pl) {
-    Info("MakeAnalysisFillAOD","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
+  if(!pl)
+  {
+    AliWarning(Form("TObjArray with %s clusters is NULL!",GetCalorimeterString().Data()));
     return;
   }
   
   //Init arrays, variables, get number of clusters
-  TLorentzVector mom, mom2 ;
   Int_t nCaloClusters = pl->GetEntriesFast();
   //List to be used in conversion analysis, to tag the cluster as candidate for conversion
   
-  if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillAOD() - input %s cluster entries %d\n", fCalorimeter.Data(), nCaloClusters);
+  AliDebug(1,Form("Input %s cluster entries %d", GetCalorimeterString().Data(), nCaloClusters));
   
   //----------------------------------------------------
   // Fill AOD with PHOS/EMCAL AliAODPWG4Particle objects
   //----------------------------------------------------
   // Loop on clusters
-  for(Int_t icalo = 0; icalo < nCaloClusters; icalo++){    
-         
+  for(Int_t icalo = 0; icalo < nCaloClusters; icalo++)
+  {
          AliVCluster * calo =  (AliVCluster*) (pl->At(icalo)); 
     //printf("calo %d, %f\n",icalo,calo->E());
     
     //Get the index where the cluster comes, to retrieve the corresponding vertex
     Int_t evtIndex = 0 ; 
-    if (GetMixedEvent()) {
+    if (GetMixedEvent())
+    {
       evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ; 
       //Get the vertex and check it is not too large in z
       if(TMath::Abs(GetVertex(evtIndex)[2])> GetZvertexCut()) continue;
     }
     
     //Cluster selection, not charged, with photon id and in fiducial cut         
-    if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
-      calo->GetMomentum(mom,GetVertex(evtIndex)) ;}//Assume that come from vertex in straight line
-    else{
+    if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
+    {
+      calo->GetMomentum(fMomentum,GetVertex(evtIndex)) ;
+    }//Assume that come from vertex in straight line
+    else
+    {
       Double_t vertex[]={0,0,0};
-      calo->GetMomentum(mom,vertex) ;
+      calo->GetMomentum(fMomentum,vertex) ;
     }
     
     //--------------------------------------
     // Cluster selection
     //--------------------------------------
     AliVCaloCells* cells    = 0;
-    if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
-    else                        cells = GetPHOSCells();
+    if(GetCalorimeter() == kEMCAL) cells = GetEMCALCells();
+    else                           cells = GetPHOSCells();
     
     Int_t nMaxima = GetCaloUtils()->GetNumberOfLocalMaxima(calo, cells); // NLM
-    if(!ClusterSelected(calo,mom,nMaxima)) continue;
-    
-    //----------------------------
-    //Create AOD for analysis
-    //----------------------------
-    AliAODPWG4Particle aodpart = AliAODPWG4Particle(mom);
-    
-    //...............................................
-    //Set the indeces of the original caloclusters (MC, ID), and calorimeter  
-    Int_t label = calo->GetLabel();
-    aodpart.SetLabel(label);
-    aodpart.SetCaloLabel(calo->GetID(),-1);
-    aodpart.SetDetector(fCalorimeter);
-    //printf("Index %d, Id %d, iaod %d\n",icalo, calo->GetID(),GetOutputAODBranch()->GetEntriesFast());
-    
-    //...............................................
-    //Set bad channel distance bit
-    Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
-    if     (distBad > fMinDist3) aodpart.SetDistToBad(2) ;
-    else if(distBad > fMinDist2) aodpart.SetDistToBad(1) ; 
-    else                         aodpart.SetDistToBad(0) ;
-    //printf("DistBad %f Bit %d\n",distBad, aodpart.DistToBad());
+    if(!ClusterSelected(calo,nMaxima)) continue;
     
     //-------------------------------------
-    //PID selection via dEdx
+    // PID selection via dE/dx
     //-------------------------------------
     
     AliVTrack *track = GetCaloUtils()->GetMatchedTrack(calo, GetReader()->GetInputEvent());
 
     if(!track)
     {
-      printf("AliAnaElectron::MakeAnalysisFillAOD() - Null track");
+      AliWarning("Null track");
       continue;
     }
     
@@ -1213,7 +1178,7 @@ void  AliAnaElectron::MakeAnalysisFillAOD()
       fhdEdxvsPCutEOverP  ->Fill(track->P(),dEdx);
     }
     
-    //Apply a mild cut on the cluster SS and check the value of dEdX and EOverP
+    // Apply a mild cut on the cluster SS and check the value of dEdX and EOverP
     Float_t m02 = calo->GetM02();
     if(m02 > 0.1 && m02 < 0.4)
     {
@@ -1239,27 +1204,24 @@ void  AliAnaElectron::MakeAnalysisFillAOD()
       if( eOverp < fEOverPMax && eOverp > fEOverPMin)
       {
         pid  = AliCaloPID::kElectron;
-      } //E/p
+      } // E/p
       
-    }// dEdx
-    
-    aodpart.SetIdentifiedParticleType(pid);    
+    }// dE/dx
     
     Int_t pidIndex = 0;// Electron
     if(pid == AliCaloPID::kChargedHadron) pidIndex = 1;
   
     //--------------------------------------------------------------------------------------
-    //Play with the MC stack if available
+    // Play with the MC stack if available
     //--------------------------------------------------------------------------------------
     
     //Check origin of the candidates
+    Int_t tag = -1 ;
     if(IsDataMC())
     {
-      Int_t tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader());
-      aodpart.SetTag(tag);
+      tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader(),GetCalorimeter());
       
-      if(GetDebug() > 0)
-        printf("AliAnaElectron::MakeAnalysisFillAOD() - Origin of candidate, bit map %d\n",aodpart.GetTag());
+      AliDebug(1,Form("Origin of candidate, bit map %d",tag));
          
       if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) && fhMCE[pidIndex][kmcPhoton])
       {
@@ -1340,13 +1302,13 @@ void  AliAnaElectron::MakeAnalysisFillAOD()
     //Fill some shower shape histograms
     //---------------------------------
 
-    FillShowerShapeHistograms(calo,aodpart.GetTag(),pid);
+    FillShowerShapeHistograms(calo,tag,pid);
   
     if(pid == AliCaloPID::kElectron)
       WeightHistograms(calo);
     
     //-----------------------------------------
-    //PID Shower Shape selection or bit setting
+    // PID Shower Shape selection or bit setting
     //-----------------------------------------
     
     // Data, PID check on
@@ -1357,34 +1319,72 @@ void  AliAnaElectron::MakeAnalysisFillAOD()
     
       if(GetCaloPID()->GetIdentifiedParticleType(calo)!=AliCaloPID::kPhoton)
       {
-        continue;
+        if(fAODParticle == AliCaloPID::kElectron)
+          continue;
+        
+        if(fAODParticle == 0 )
+          pid = AliCaloPID::kChargedHadron ;
       }
-      if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodpart.GetIdentifiedParticleType());
       
+      AliDebug(1,Form("PDG of identified particle %d",pid));
     }
         
-    if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() - Photon selection cuts passed: pT %3.2f, pdg %d\n",
-                              aodpart.Pt(), aodpart.GetIdentifiedParticleType());
+    AliDebug(1,Form("Photon selection cuts passed: pT %3.2f, pdg %d",fMomentum.Pt(),pid));
     
-    //FIXME, this to MakeAnalysisFillHistograms ...
-    Int_t absID             = 0; 
     Float_t maxCellFraction = 0;
+    Int_t absID = GetCaloUtils()->GetMaxEnergyCell(cells, calo,maxCellFraction);
+    if ( absID >= 0 )fhMaxCellDiffClusterE[pidIndex]->Fill(fMomentum.E(),maxCellFraction);
     
-    absID = GetCaloUtils()->GetMaxEnergyCell(cells, calo,maxCellFraction);
-    fhMaxCellDiffClusterE[pidIndex]->Fill(aodpart.E(),maxCellFraction);
-    fhNCellsE[pidIndex]            ->Fill(aodpart.E(),calo->GetNCells());
-    fhNLME[pidIndex]               ->Fill(aodpart.E(),nMaxima);
-    fhTimeE[pidIndex]              ->Fill(aodpart.E(),calo->GetTOF()*1.e9);
+    fhNCellsE[pidIndex] ->Fill(fMomentum.E(),calo->GetNCells());
+    fhNLME   [pidIndex] ->Fill(fMomentum.E(),nMaxima);
+    fhTimeE  [pidIndex] ->Fill(fMomentum.E(),calo->GetTOF()*1.e9);
+    
+    //----------------------------
+    // Create AOD for analysis
+    //----------------------------
 
     //Add AOD with electron/hadron object to aod branch
     if ( pid == fAODParticle || fAODParticle == 0 ) 
     {
+      AliAODPWG4Particle aodpart = AliAODPWG4Particle(fMomentum);
+      
+      //...............................................
+      //Set the indeces of the original caloclusters (MC, ID), and calorimeter
+      Int_t label = calo->GetLabel();
+      aodpart.SetLabel(label);
+      aodpart.SetCaloLabel (calo ->GetID(),-1);
+      aodpart.SetTrackLabel(track->GetID(),-1);
+
+      aodpart.SetDetectorTag(GetCalorimeter());
+      //printf("Index %d, Id %d, iaod %d\n",icalo, calo->GetID(),GetOutputAODBranch()->GetEntriesFast());
+      
+      aodpart.SetM02(calo->GetM02());
+      aodpart.SetNLM(nMaxima);
+      aodpart.SetTime(calo->GetTOF()*1e9);
+      aodpart.SetNCells(calo->GetNCells());
+      Int_t nSM = GetModuleNumber(calo);
+      aodpart.SetSModNumber(nSM);
+      
+      //...............................................
+      //Set bad channel distance bit
+      Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
+      if     (distBad > fMinDist3) aodpart.SetDistToBad(2) ;
+      else if(distBad > fMinDist2) aodpart.SetDistToBad(1) ;
+      else                         aodpart.SetDistToBad(0) ;
+      //printf("DistBad %f Bit %d\n",distBad, aodpart.DistToBad());
+
+      // MC tag
+      aodpart.SetTag(tag);
+      
+      // PID tag
+      aodpart.SetIdentifiedParticleType(pid);
+      
       AddAODParticle(aodpart);
     }
 
   }//loop
   
-  if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD()  End fill AODs, with %d entries \n",GetOutputAODBranch()->GetEntriesFast());  
+  AliDebug(1,Form("End fill AODs, with %d entries",GetOutputAODBranch()->GetEntriesFast()));
   
 }
 
@@ -1400,23 +1400,18 @@ void  AliAnaElectron::MakeAnalysisFillHistograms()
   TClonesArray     * mcparticles = 0x0;
   AliAODMCParticle * aodprimary  = 0x0; 
   
-  if(IsDataMC()){
-    
-    if(GetReader()->ReadStack()){
+  if(IsDataMC())
+  {
+    if(GetReader()->ReadStack())
+    {
       stack =  GetMCStack() ;
-      if(!stack) {
-        printf("AliAnaElectron::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
-        abort();
-      }
-      
+      if ( !stack )       AliFatal("Stack not available, is the MC handler called? STOP");
     }
-    else if(GetReader()->ReadAODMCParticles()){
-      
+    else if(GetReader()->ReadAODMCParticles())
+    {
       //Get the list of MC particles
       mcparticles = GetReader()->GetAODMCParticles();
-      if(!mcparticles && GetDebug() > 0)       {
-        printf("AliAnaElectron::MakeAnalysisFillHistograms() -  Standard MCParticles not available!\n");
-      }        
+      if ( !mcparticles ) AliFatal("Standard MCParticles not available! STOP");
     }
   }// is data and MC
   
@@ -1430,7 +1425,7 @@ void  AliAnaElectron::MakeAnalysisFillHistograms()
   //----------------------------------
   //Loop on stored AOD photons
   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
-  if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
+  AliDebug(1,Form("AOD branch entries %d", naod));
   
   for(Int_t iaod = 0; iaod < naod ; iaod++)
   {
@@ -1442,10 +1437,9 @@ void  AliAnaElectron::MakeAnalysisFillHistograms()
     else if(pdg == AliCaloPID::kChargedHadron) pidIndex = 1;
     else                                       continue    ;
           
-    if(ph->GetDetector() != fCalorimeter) continue;
+    if(((Int_t) ph->GetDetectorTag()) != GetCalorimeter()) continue;
     
-    if(GetDebug() > 2) 
-      printf("AliAnaElectron::MakeAnalysisFillHistograms() - ID Electron: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
+    AliDebug(1,Form("ID Electron: pt %f, phi %f, eta %f", ph->Pt(),ph->Phi(),ph->Eta())) ;
     
     //................................
     //Fill photon histograms 
@@ -1458,47 +1452,53 @@ void  AliAnaElectron::MakeAnalysisFillHistograms()
     fhPt[pidIndex]  ->Fill(ptcluster);
     fhPhi[pidIndex] ->Fill(ptcluster,phicluster);
     fhEta[pidIndex] ->Fill(ptcluster,etacluster);    
-    if     (ecluster > 0.5)   fhEtaPhi[pidIndex]  ->Fill(etacluster, phicluster);
+    if     (ecluster   > 0.5) fhEtaPhi[pidIndex]  ->Fill(etacluster, phicluster);
     else if(GetMinPt() < 0.5) fhEtaPhi05[pidIndex]->Fill(etacluster, phicluster);
   
     //.......................................
     //Play with the MC data if available
-    if(IsDataMC()){
-      
+    if(IsDataMC())
+    {
       //....................................................................
       // Access MC information in stack if requested, check that it exists.
       Int_t label =ph->GetLabel();
-      if(label < 0) {
-        if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillHistograms() *** bad label ***:  label %d \n", label);
+      if(label < 0)
+      {
+        AliDebug(1,Form("*** bad label ***:  label %d", label));
         continue;
       }
       
       Float_t eprim   = 0;
-      Float_t ptprim  = 0;
-      if(GetReader()->ReadStack()){
-        
-        if(label >=  stack->GetNtrack()) {
-          if(GetDebug() > 2)  printf("AliAnaElectron::MakeAnalysisFillHistograms() *** large label ***:  label %d, n tracks %d \n", label, stack->GetNtrack());
+      //Float_t ptprim  = 0;
+      if(GetReader()->ReadStack())
+      {
+        if(label >=  stack->GetNtrack())
+        {
+          AliDebug(1,Form("*** large label ***:  label %d, n tracks %d", label, stack->GetNtrack()));
           continue ;
         }
         
         primary = stack->Particle(label);
-        if(!primary){
-          printf("AliAnaElectron::MakeAnalysisFillHistograms() *** no primary ***:  label %d \n", label);
-          continue;
+        if(!primary)
+        {
+          AliWarning(Form("*** no primary ***:  label %d", label));
+          continue ;
         }
         
         eprim   = primary->Energy();
-        ptprim  = primary->Pt();               
+        //ptprim  = primary->Pt();
         
       }
-      else if(GetReader()->ReadAODMCParticles()){
+      else if(GetReader()->ReadAODMCParticles())
+      {
         //Check which is the input
-        if(ph->GetInputFileIndex() == 0){
+        if(ph->GetInputFileIndex() == 0)
+        {
           if(!mcparticles) continue;
-          if(label >=  mcparticles->GetEntriesFast()) {
-            if(GetDebug() > 2)  printf("AliAnaElectron::MakeAnalysisFillHistograms() *** large label ***:  label %d, n tracks %d \n", 
-                                       label, mcparticles->GetEntriesFast());
+          
+          if(label >=  mcparticles->GetEntriesFast())
+          {
+            AliDebug(1,Form("*** large label ***:  label %d, n tracks %d",label, mcparticles->GetEntriesFast()));
             continue ;
           }
           //Get the particle
@@ -1506,14 +1506,14 @@ void  AliAnaElectron::MakeAnalysisFillHistograms()
           
         }
         
-        if(!aodprimary){
-          printf("AliAnaElectron::MakeAnalysisFillHistograms() *** no primary ***:  label %d \n", label);
+        if(!aodprimary)
+        {
+          AliWarning(Form("*** no primary ***:  label %d", label));
           continue;
         }
         
         eprim   = aodprimary->E();
-        ptprim  = aodprimary->Pt();
-        
+        //ptprim  = aodprimary->Pt();
       }
       
       Int_t tag =ph->GetTag();
@@ -1617,7 +1617,8 @@ void  AliAnaElectron::MakeAnalysisFillHistograms()
         fhMCDeltaE[pidIndex][kmcElectron] ->Fill(ecluster,eprim-ecluster);
         
       }     
-      else if( fhMCE[pidIndex][kmcOther]){
+      else if( fhMCE[pidIndex][kmcOther])
+      {
         fhMCE  [pidIndex][kmcOther] ->Fill(ecluster);
         fhMCPt [pidIndex][kmcOther] ->Fill(ptcluster);
         fhMCPhi[pidIndex][kmcOther] ->Fill(ecluster,phicluster);
@@ -1645,7 +1646,7 @@ void AliAnaElectron::Print(const Option_t * opt) const
   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
   AliAnaCaloTrackCorrBaseClass::Print(" ");
 
-  printf("Calorimeter            =     %s\n", fCalorimeter.Data()) ;
+  printf("Calorimeter            =     %s\n", GetCalorimeterString().Data()) ;
   printf(" %2.2f < dEdx < %2.2f  \n",fdEdxMin,fdEdxMax) ;
   printf(" %2.2f <  E/P < %2.2f  \n",fEOverPMin,fEOverPMax) ;
   printf("Min Distance to Bad Channel   = %2.1f\n",fMinDist);
@@ -1665,8 +1666,8 @@ void AliAnaElectron::WeightHistograms(AliVCluster *clus)
   if(!fFillWeightHistograms || GetMixedEvent()) return;
   
   AliVCaloCells* cells = 0;
-  if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
-  else                        cells = GetPHOSCells();
+  if(GetCalorimeter() == kEMCAL) cells = GetEMCALCells();
+  else                           cells = GetPHOSCells();
   
   // First recalculate energy in case non linearity was applied
   Float_t  energy = 0;
@@ -1677,7 +1678,7 @@ void AliAnaElectron::WeightHistograms(AliVCluster *clus)
     
     //Recalibrate cell energy if needed
     Float_t amp = cells->GetCellAmplitude(id);
-    GetCaloUtils()->RecalibrateCellAmplitude(amp,fCalorimeter, id);
+    GetCaloUtils()->RecalibrateCellAmplitude(amp,GetCalorimeter(), id);
     
     energy    += amp;
     
@@ -1686,8 +1687,9 @@ void AliAnaElectron::WeightHistograms(AliVCluster *clus)
     
   } // energy loop       
   
-  if(energy <=0 ) {
-    printf("AliAnaCalorimeterQA::WeightHistograms()- Wrong calculated energy %f\n",energy);
+  if ( energy <= 0 )
+  {
+    AliWarning(Form("Wrong calculated energy %f",energy));
     return;
   }
   
@@ -1701,7 +1703,7 @@ void AliAnaElectron::WeightHistograms(AliVCluster *clus)
     
     //Recalibrate cell energy if needed
     Float_t amp = cells->GetCellAmplitude(id);
-    GetCaloUtils()->RecalibrateCellAmplitude(amp, fCalorimeter, id);
+    GetCaloUtils()->RecalibrateCellAmplitude(amp, GetCalorimeter(), id);
 
     //printf("energy %f, amp %f, rat %f, lograt %f\n",energy,amp,amp/energy,TMath::Log(amp/energy));
     fhECellClusterRatio   ->Fill(energy,amp/energy);
@@ -1709,8 +1711,8 @@ void AliAnaElectron::WeightHistograms(AliVCluster *clus)
   }        
   
   //Recalculate shower shape for different W0
-  if(fCalorimeter=="EMCAL"){
-    
+  if(GetCalorimeter()==kEMCAL)
+  {
     Float_t l0org = clus->GetM02();
     Float_t l1org = clus->GetM20();
     Float_t dorg  = clus->GetDispersion();