PWG4/GammaConv/AliAnalysisTaskCaloConv.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: Dmitri Peressounko (RRC KI)                                    *
   5  *                                                                        *
   6  * Permission to use, copy, modify and distribute this software and its   *
   7  * documentation strictly for non-commercial purposes is hereby granted   *
   8  * without fee, provided that the above copyright notice appears in all   *
   9  * copies and that both the copyright notice and this permission notice   *
  10  * appear in the supporting documentation. The authors make no claims     *
  11  * about the suitability of this software for any purpose. It is          *
  12  * provided "as is" without express or implied warranty.                  *
  13  **************************************************************************/
  14
  15 ////////////////////////////////////////////////
  16 //---------------------------------------------
  17 // Class used to do analysis on conversion + calorimeter pairs
  18 // A lot of code cut-and-pasted from AliV0Reader and GammaConversion classes
  19 //
  20 //---------------------------------------------
  21 ////////////////////////////////////////////////
  22
  23 // root
  24 #include "TChain.h"
  25 #include "TH3.h"
  26 #include "TH2.h"
  27
  28 // analysis
  29 #include "AliAnalysisManager.h"
  30 #include "AliAnalysisTaskCaloConv.h"
  31 #include "AliStack.h"
  32 #include "AliLog.h"
  33 #include "AliESDEvent.h"
  34 #include "AliESDpid.h"
  35 #include "AliESDtrackCuts.h"
  36 #include "AliCFContainer.h"   // for CF
  37 #include "AliESDCaloCluster.h"
  38 #include "AliPHOSGeoUtils.h"
  39 #include "AliEMCALGeoUtils.h"
  40 #include "AliCFContainer.h"
  41 #include "AliMCEventHandler.h"
  42 #include "AliMCEvent.h"
  43 #include "AliESDv0.h"
  44 #include "AliKFParticle.h"
  45 #include "AliKFVertex.h"
  46
  47 class AliESDInputHandler;
  48 class Riostream;
  49 class TFile;
  50
  51 ClassImp(AliAnalysisTaskCaloConv)
  52
  53
  54 AliAnalysisTaskCaloConv::AliAnalysisTaskCaloConv():
  55 AliAnalysisTaskSE(),
  56   fESDEvent(NULL),
  57   fESDpid(NULL),
  58   fStack(NULL),
  59   fOutputContainer(NULL),
  60   fCFOutputContainer(NULL),
  61   fConvCFCont(0x0),
  62   fPHOSCFCont(0x0),
  63   fEMCALCFCont(0x0),
  64   fPi0CFCont(0x0),
  65   fTriggerCINT1B(kFALSE),
  66   fMinOpeningAngleGhostCut(0.),
  67   fPHOSgeom(0x0),
  68   fEMCALgeom(0x0),
  69   fPi0Thresh1(0.5),
  70   fPi0Thresh2(1.),
  71   fConvEvent(NULL) ,
  72   fPHOSEvent(NULL),
  73   fEMCALEvent(NULL),
  74   fnSigmaAboveElectronLine(5.),
  75   fnSigmaBelowElectronLine(-3.),
  76   fnSigmaAbovePionLine(0.),
  77   fpnSigmaAbovePionLine(1.),
  78   fprobCut(0.),
  79   fmaxR(180.),
  80   fmaxZ(240.),
  81   fetaCut(0.9),
  82   fptCut(0.02),
  83   fchi2CutConversion(30.)
  84 {
  85   // Default constructor
  86   Int_t nBin=10 ;
  87   for(Int_t i=0;i<nBin;i++){
  88     fPHOSEvents[i]=0 ;
  89     fEMCALEvents[i]=0;
  90     fConvEvents[i]=0;
  91   }
  92 }
  93 AliAnalysisTaskCaloConv::AliAnalysisTaskCaloConv(const char* name):
  94   AliAnalysisTaskSE(name),
  95   fESDEvent(NULL),
  96   fESDpid(NULL),
  97   fStack(NULL),
  98   fOutputContainer(NULL),
  99   fCFOutputContainer(NULL),
 100   fConvCFCont(0x0),
 101   fPHOSCFCont(0x0),
 102   fEMCALCFCont(0x0),
 103   fPi0CFCont(0x0),
 104   fTriggerCINT1B(kFALSE),
 105   fMinOpeningAngleGhostCut(0.),
 106   fPHOSgeom(0x0),
 107   fEMCALgeom(0x0),
 108   fPi0Thresh1(0.5),
 109   fPi0Thresh2(1.),
 110   fConvEvent(NULL) ,
 111   fPHOSEvent(NULL),
 112   fEMCALEvent(NULL),
 113   fnSigmaAboveElectronLine(5.),
 114   fnSigmaBelowElectronLine(-3.),
 115   fnSigmaAbovePionLine(0.),
 116   fpnSigmaAbovePionLine(1.),
 117   fprobCut(0.),
 118   fmaxR(180.),
 119   fmaxZ(240.),
 120   fetaCut(0.9),
 121   fptCut(0.02),
 122   fchi2CutConversion(30.)
 123 {
 124   // Common I/O in slot 0
 125   DefineInput (0, TChain::Class());
 126   DefineOutput(0, TTree::Class());
 127
 128   // Your private output
 129   DefineOutput(1, TList::Class());
 130   DefineOutput(2, TList::Class());  // for CF
 131
 132   Int_t nBin=10 ;
 133   for(Int_t i=0;i<nBin;i++){
 134     fPHOSEvents[i]=0 ;
 135     fEMCALEvents[i]=0;
 136     fConvEvents[i]=0;
 137   }
 138   fESDpid = new AliESDpid;
 139 }
 140 //_____________________________________________________
 141 AliAnalysisTaskCaloConv::~AliAnalysisTaskCaloConv()
 142 {
 143   // Remove all pointers
 144
 145   if(fOutputContainer){
 146     fOutputContainer->Clear() ;
 147     delete fOutputContainer ;
 148   }
 149   if(fCFOutputContainer){
 150     fCFOutputContainer->Clear() ;
 151     delete fCFOutputContainer ;
 152   }
 153
 154   if(fPHOSgeom){
 155     delete fPHOSgeom ;
 156     fPHOSgeom=0x0 ;
 157   }
 158
 159   if(fEMCALgeom){
 160     delete fEMCALgeom ;
 161     fEMCALgeom=0x0;
 162   }
 163
 164   for(Int_t ivtx=0; ivtx<10; ivtx++){
 165     if(fPHOSEvents[ivtx]){
 166       delete fPHOSEvents[ivtx] ;
 167       fPHOSEvents[ivtx]=0x0 ;
 168     }
 169     if(fEMCALEvents[ivtx]){
 170       delete fEMCALEvents[ivtx] ;
 171       fEMCALEvents[ivtx]=0x0 ;
 172     }
 173     if(fConvEvents[ivtx]){
 174       delete fConvEvents[ivtx] ;
 175       fConvEvents[ivtx]=0x0 ;
 176     }
 177   }
 178
 179 }
 180 //_____________________________________________________
 181 void AliAnalysisTaskCaloConv::Init()
 182 {
 183   // Initialization
 184   // AliLog::SetGlobalLogLevel(AliLog::kError);
 185 }
 186 //_____________________________________________________
 187 void AliAnalysisTaskCaloConv::UserExec(Option_t */*option*/)
 188 {
 189   // Execute analysis for current event
 190   // First select conversion and calorimeter photons
 191   // then construct inv. mass distributions
 192
 193   //First try to find Stack information.
 194   if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()){
 195     if(static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent())
 196       fStack = static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent()->Stack();
 197   }
 198
 199
 200   fESDEvent=(AliESDEvent*)InputEvent();
 201 //  //Take Only events with proper trigger
 202 //  if(fTriggerCINT1B){
 203 //    if(!fESDEvent->IsTriggerClassFired("CINT1B-ABCE-NOPF-ALL")) return;
 204 //  }
 205
 206   //Init geometry if not done yet
 207   InitGeometry();
 208
 209   //Select conversion and calorimeter photons
 210   //Conversion photons should go first since they are used in calibration in SelectCALOPhotons()
 211   SelectConvPhotons() ;
 212   SelectPHOSPhotons() ;
 213   SelectEMCALPhotons() ;
 214   FillRealMixed() ;
 215   //Fill MC histograms if MC is present
 216   ProcessMC();
 217
 218   PostData(1, fOutputContainer);
 219 }
 220 //____________________________________________________________
 221 void AliAnalysisTaskCaloConv::ConnectInputData(Option_t *option){
 222   // see header file for documentation
 223
 224   AliAnalysisTaskSE::ConnectInputData(option);
 225
 226 }
 227 //____________________________________________________________
 228 void AliAnalysisTaskCaloConv::UserCreateOutputObjects()
 229 {
 230   // Create the output container
 231   if(fOutputContainer != NULL){
 232     delete fOutputContainer;
 233     fOutputContainer = NULL;
 234   }
 235   fOutputContainer = new TList();
 236   fOutputContainer->SetOwner(kTRUE);
 237
 238   if(fCFOutputContainer != NULL){
 239     delete fCFOutputContainer;
 240     fCFOutputContainer = NULL;
 241   }
 242   fCFOutputContainer = new TList();
 243   fCFOutputContainer->SetOwner(kTRUE);
 244
 245   //===========Correction Framework ======================
 246   //bins: pt,eta,mass
 247   Int_t iBin[3]={500,40,100};
 248   fConvCFCont = new AliCFContainer("ConvContainer","container for converted photons", 23,3,iBin);
 249   fConvCFCont->SetBinLimits(0,0.,50.);
 250   fConvCFCont->SetBinLimits(1,-2.,2.) ;
 251   fConvCFCont->SetBinLimits(2,0.,1.);
 252   fCFOutputContainer->Add(fConvCFCont) ;
 253
 254   fPHOSCFCont = new AliCFContainer("PHOSContainer","container for PHOS photons", 10,2,iBin);
 255   fPHOSCFCont->SetBinLimits(0,0.,50.);
 256   fPHOSCFCont->SetBinLimits(1,-2.,2.) ;
 257   fCFOutputContainer->Add(fPHOSCFCont) ;
 258
 259   fEMCALCFCont = new AliCFContainer("EMCALContainer","container for EMCAL photons", 10,2,iBin);
 260   fEMCALCFCont->SetBinLimits(0,0.,50.);
 261   fEMCALCFCont->SetBinLimits(1,-2.,2.) ;
 262   fCFOutputContainer->Add(fEMCALCFCont) ;
 263
 264   fPi0CFCont = new AliCFContainer("Pi0Container","container for EMCAL photons", 10,2,iBin);
 265   fPi0CFCont->SetBinLimits(0,0.,50.);
 266   fPi0CFCont->SetBinLimits(1,-2.,2.) ;
 267   fCFOutputContainer->Add(fPi0CFCont) ;
 268
 269   //========================================================
 270
 271
 272   //Adding the histograms to the output container
 273   Int_t firstRun= 114700 ;
 274   Int_t lastRun = 128000 ;
 275   Int_t nRuns =lastRun-firstRun+1 ;
 276
 277   //Run QA histigrams
 278   fOutputContainer->Add(new TH2F("hRunTrigger","Triggers fired",nRuns,float(firstRun),float(lastRun),2,0.,2.)) ;
 279   fOutputContainer->Add(new TH1F("hRunEvents","Events per run",nRuns,float(firstRun),float(lastRun))) ;
 280   fOutputContainer->Add(new TH1F("hRunConvs","Conversion photons per run",nRuns,float(firstRun),float(lastRun))) ;
 281   fOutputContainer->Add(new TH1F("hRunPHOS","PHOS photons per run",nRuns,float(firstRun),float(lastRun))) ;
 282   fOutputContainer->Add(new TH1F("hRunEMCAL","EMCAL photons per run",nRuns,float(firstRun),float(lastRun))) ;
 283   fOutputContainer->Add(new TH1F("hVtxBin","Vtx distribution",10,0.,10.)) ;
 284   fOutputContainer->Add(new TH1F("hEvents","Events processed",1,0.,1.)) ;
 285
 286   //Check of geometry
 287   fOutputContainer->Add(new TH3F("PHOS_beyond","PHOS clusters not in PHOS",200,-300.,300.,200,-500.,0.,200,-100.,100.)) ;
 288
 289   Int_t npt=200 ;
 290   Double_t ptmax=20. ;
 291   //Calibration of PHOS
 292   fOutputContainer->Add(new TH3F("PHOS_mod1_th1","Inv.Mass distr. per channel",64,0.,64,56,0.,56,200,0.,1.)) ;
 293   fOutputContainer->Add(new TH3F("PHOS_mod2_th1","Inv.Mass distr. per channel",64,0.,64,56,0.,56,200,0.,1.)) ;
 294   fOutputContainer->Add(new TH3F("PHOS_mod3_th1","Inv.Mass distr. per channel",64,0.,64,56,0.,56,200,0.,1.)) ;
 295   fOutputContainer->Add(new TH3F("PHOS_mod4_th1","Inv.Mass distr. per channel",64,0.,64,56,0.,56,200,0.,1.)) ;
 296   fOutputContainer->Add(new TH3F("PHOS_mod5_th1","Inv.Mass distr. per channel",64,0.,64,56,0.,56,200,0.,1.)) ;
 297
 298   fOutputContainer->Add(new TH3F("PHOS_mod1_th2","Inv.Mass distr. per channel",64,0.,64,56,0.,56,200,0.,1.)) ;
 299   fOutputContainer->Add(new TH3F("PHOS_mod2_th2","Inv.Mass distr. per channel",64,0.,64,56,0.,56,200,0.,1.)) ;
 300   fOutputContainer->Add(new TH3F("PHOS_mod3_th2","Inv.Mass distr. per channel",64,0.,64,56,0.,56,200,0.,1.)) ;
 301   fOutputContainer->Add(new TH3F("PHOS_mod4_th2","Inv.Mass distr. per channel",64,0.,64,56,0.,56,200,0.,1.)) ;
 302   fOutputContainer->Add(new TH3F("PHOS_mod5_th2","Inv.Mass distr. per channel",64,0.,64,56,0.,56,200,0.,1.)) ;
 303
 304   //Pi0 histograms
 305   //Vary Conversion cuts
 306   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_OnFly","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 307   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_Offline","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 308   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_On_Kink","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 309   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_Off_Kink","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 310   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_On_dEdx","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 311   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_Off_dEdx","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 312   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_On_Prob","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 313   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_Off_Prob","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 314   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_On_R120","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 315   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_Off_R120","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 316   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_On_Z","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 317   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_Off_Z","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 318   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_On_chi","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 319   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_Off_chi","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 320   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_On_Eta","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 321   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_Off_Eta","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 322   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_On_Wcut","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 323   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_Off_Wcut","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 324
 325   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_OnFly","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 326   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_Offline","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 327   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_On_Kink","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 328   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_Off_Kink","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 329   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_On_dEdx","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 330   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_Off_dEdx","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 331   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_On_Prob","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 332   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_Off_Prob","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 333   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_On_R120","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 334   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_Off_R120","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 335   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_On_Z","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 336   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_Off_Z","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 337   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_On_chi","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 338   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_Off_chi","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 339   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_On_Eta","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 340   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_Off_Eta","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 341   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_On_Wcut","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 342   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_Off_Wcut","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 343
 344   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_OnFly","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 345   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_Offline","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 346   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_On_Kink","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 347   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_Off_Kink","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 348   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_On_dEdx","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 349   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_Off_dEdx","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 350   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_On_Prob","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 351   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_Off_Prob","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 352   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_On_R120","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 353   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_Off_R120","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 354   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_On_Z","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 355   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_Off_Z","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 356   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_On_chi","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 357   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_Off_chi","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 358   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_On_Eta","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 359   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_Off_Eta","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 360   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_On_Wcut","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 361   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_Off_Wcut","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 362
 363   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_OnFly","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 364   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_Offline","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 365   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_On_Kink","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 366   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_Off_Kink","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 367   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_On_dEdx","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 368   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_Off_dEdx","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 369   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_On_Prob","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 370   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_Off_Prob","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 371   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_On_R120","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 372   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_Off_R120","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 373   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_On_Z","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 374   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_Off_Z","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 375   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_On_chi","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 376   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_Off_chi","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 377   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_On_Eta","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 378   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_Off_Eta","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 379   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_On_Wcut","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 380   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_Off_Wcut","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 381
 382   //PHOS PID variations
 383   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_all","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 384   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_Disp","Mass vs pt, disp cut",400,0.,1.,npt,0.,ptmax)) ;
 385   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_TOF","Mass vs pt, TOF cut",400,0.,1.,npt,0.,ptmax)) ;
 386   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_Neutral","Mass vs pt, Neutral cut",400,0.,1.,npt,0.,ptmax)) ;
 387   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_DispNeutral","Mass vs pt, Disp & neutral cut",400,0.,1.,npt,0.,ptmax)) ;
 388
 389   //PHOS PID module-by-module
 390   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_mod1_all","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 391   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_mod2_all","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 392   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_mod3_all","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 393   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_mod4_all","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 394   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_mod5_all","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 395   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_mod1_Disp","Mass vs pt, disp cut",400,0.,1.,npt,0.,ptmax)) ;
 396   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_mod2_Disp","Mass vs pt, disp cut",400,0.,1.,npt,0.,ptmax)) ;
 397   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_mod3_Disp","Mass vs pt, disp cut",400,0.,1.,npt,0.,ptmax)) ;
 398   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_mod4_Disp","Mass vs pt, disp cut",400,0.,1.,npt,0.,ptmax)) ;
 399   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_mod5_Disp","Mass vs pt, disp cut",400,0.,1.,npt,0.,ptmax)) ;
 400   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_mod1_TOF","Mass vs pt, TOF cut",400,0.,1.,npt,0.,ptmax)) ;
 401   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_mod2_TOF","Mass vs pt, TOF cut",400,0.,1.,npt,0.,ptmax)) ;
 402   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_mod3_TOF","Mass vs pt, TOF cut",400,0.,1.,npt,0.,ptmax)) ;
 403   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_mod4_TOF","Mass vs pt, TOF cut",400,0.,1.,npt,0.,ptmax)) ;
 404   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_mod5_TOF","Mass vs pt, TOF cut",400,0.,1.,npt,0.,ptmax)) ;
 405   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_mod1_Neutral","Mass vs pt, Neutral cut",400,0.,1.,npt,0.,ptmax)) ;
 406   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_mod2_Neutral","Mass vs pt, Neutral cut",400,0.,1.,npt,0.,ptmax)) ;
 407   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_mod3_Neutral","Mass vs pt, Neutral cut",400,0.,1.,npt,0.,ptmax)) ;
 408   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_mod4_Neutral","Mass vs pt, Neutral cut",400,0.,1.,npt,0.,ptmax)) ;
 409   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_mod5_Neutral","Mass vs pt, Neutral cut",400,0.,1.,npt,0.,ptmax)) ;
 410   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_mod1_DispNeutral","Mass vs pt, Disp & neutral cut",400,0.,1.,npt,0.,ptmax)) ;
 411   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_mod2_DispNeutral","Mass vs pt, Disp & neutral cut",400,0.,1.,npt,0.,ptmax)) ;
 412   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_mod3_DispNeutral","Mass vs pt, Disp & neutral cut",400,0.,1.,npt,0.,ptmax)) ;
 413   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_mod4_DispNeutral","Mass vs pt, Disp & neutral cut",400,0.,1.,npt,0.,ptmax)) ;
 414   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_mod5_DispNeutral","Mass vs pt, Disp & neutral cut",400,0.,1.,npt,0.,ptmax)) ;
 415
 416   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_all","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 417   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_Disp","Mass vs pt, disp cut",400,0.,1.,npt,0.,ptmax)) ;
 418   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_TOF","Mass vs pt, TOF cut",400,0.,1.,npt,0.,ptmax)) ;
 419   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_Neutral","Mass vs pt, Neutral cut",400,0.,1.,npt,0.,ptmax)) ;
 420   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_DispNeutral","Mass vs pt, Disp & neutral cut",400,0.,1.,npt,0.,ptmax)) ;
 421
 422   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_mod1_single","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 423   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_mod2_single","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 424   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_mod3_single","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 425   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_mod4_single","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 426   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_mod5_single","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 427
 428   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_mod1_all","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 429   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_mod2_all","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 430   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_mod3_all","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 431   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_mod4_all","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 432   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_mod5_all","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 433   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_mod1_Disp","Mass vs pt, disp cut",400,0.,1.,npt,0.,ptmax)) ;
 434   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_mod2_Disp","Mass vs pt, disp cut",400,0.,1.,npt,0.,ptmax)) ;
 435   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_mod3_Disp","Mass vs pt, disp cut",400,0.,1.,npt,0.,ptmax)) ;
 436   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_mod4_Disp","Mass vs pt, disp cut",400,0.,1.,npt,0.,ptmax)) ;
 437   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_mod5_Disp","Mass vs pt, disp cut",400,0.,1.,npt,0.,ptmax)) ;
 438   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_mod1_TOF","Mass vs pt, TOF cut",400,0.,1.,npt,0.,ptmax)) ;
 439   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_mod2_TOF","Mass vs pt, TOF cut",400,0.,1.,npt,0.,ptmax)) ;
 440   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_mod3_TOF","Mass vs pt, TOF cut",400,0.,1.,npt,0.,ptmax)) ;
 441   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_mod4_TOF","Mass vs pt, TOF cut",400,0.,1.,npt,0.,ptmax)) ;
 442   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_mod5_TOF","Mass vs pt, TOF cut",400,0.,1.,npt,0.,ptmax)) ;
 443   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_mod1_Neutral","Mass vs pt, Neutral cut",400,0.,1.,npt,0.,ptmax)) ;
 444   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_mod2_Neutral","Mass vs pt, Neutral cut",400,0.,1.,npt,0.,ptmax)) ;
 445   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_mod3_Neutral","Mass vs pt, Neutral cut",400,0.,1.,npt,0.,ptmax)) ;
 446   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_mod4_Neutral","Mass vs pt, Neutral cut",400,0.,1.,npt,0.,ptmax)) ;
 447   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_mod5_Neutral","Mass vs pt, Neutral cut",400,0.,1.,npt,0.,ptmax)) ;
 448   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_mod1_DispNeutral","Mass vs pt, Disp & neutral cut",400,0.,1.,npt,0.,ptmax)) ;
 449   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_mod2_DispNeutral","Mass vs pt, Disp & neutral cut",400,0.,1.,npt,0.,ptmax)) ;
 450   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_mod3_DispNeutral","Mass vs pt, Disp & neutral cut",400,0.,1.,npt,0.,ptmax)) ;
 451   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_mod4_DispNeutral","Mass vs pt, Disp & neutral cut",400,0.,1.,npt,0.,ptmax)) ;
 452   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_mod5_DispNeutral","Mass vs pt, Disp & neutral cut",400,0.,1.,npt,0.,ptmax)) ;
 453
 454   //Single photon spectrum
 455   //Conversion
 456   fOutputContainer->Add(new TH1F("Single_conv","Single photon spectrum",npt,0.,ptmax)) ;
 457   fOutputContainer->Add(new TH1F("Single_conv_OnFly","Single photon spectrum",npt,0.,ptmax)) ;
 458   fOutputContainer->Add(new TH1F("Single_conv_Offline","Single photon spectrum",npt,0.,ptmax)) ;
 459   fOutputContainer->Add(new TH1F("Single_conv_On_Kink","Single photon spectrum",npt,0.,ptmax)) ;
 460   fOutputContainer->Add(new TH1F("Single_conv_Off_Kink","Single photon spectrum",npt,0.,ptmax)) ;
 461   fOutputContainer->Add(new TH1F("Single_conv_On_dEdx","Single photon spectrum",npt,0.,ptmax)) ;
 462   fOutputContainer->Add(new TH1F("Single_conv_Off_dEdx","Single photon spectrum",npt,0.,ptmax)) ;
 463   fOutputContainer->Add(new TH1F("Single_conv_On_Prob","Single photon spectrum",npt,0.,ptmax)) ;
 464   fOutputContainer->Add(new TH1F("Single_conv_Off_Prob","Single photon spectrum",npt,0.,ptmax)) ;
 465   fOutputContainer->Add(new TH1F("Single_conv_On_R120","Single photon spectrum",npt,0.,ptmax)) ;
 466   fOutputContainer->Add(new TH1F("Single_conv_Off_R120","Single photon spectrum",npt,0.,ptmax)) ;
 467   fOutputContainer->Add(new TH1F("Single_conv_On_Z","Single photon spectrum",npt,0.,ptmax)) ;
 468   fOutputContainer->Add(new TH1F("Single_conv_Off_Z","Single photon spectrum",npt,0.,ptmax)) ;
 469   fOutputContainer->Add(new TH1F("Single_conv_On_chi","Single photon spectrum",npt,0.,ptmax)) ;
 470   fOutputContainer->Add(new TH1F("Single_conv_Off_chi","Single photon spectrum",npt,0.,ptmax)) ;
 471   fOutputContainer->Add(new TH1F("Single_conv_On_Eta","Single photon spectrum",npt,0.,ptmax)) ;
 472   fOutputContainer->Add(new TH1F("Single_conv_Off_Eta","Single photon spectrum",npt,0.,ptmax)) ;
 473   fOutputContainer->Add(new TH1F("Single_conv_On_Wcut","Single photon spectrum",npt,0.,ptmax)) ;
 474   fOutputContainer->Add(new TH1F("Single_conv_Off_Wcut","Single photon spectrum",npt,0.,ptmax)) ;
 475
 476   //PHOS
 477   fOutputContainer->Add(new TH1F("PHOS_single_mod1_all","Single photon spectrum",npt,0.,ptmax)) ;
 478   fOutputContainer->Add(new TH1F("PHOS_single_mod2_all","Single photon spectrum",npt,0.,ptmax)) ;
 479   fOutputContainer->Add(new TH1F("PHOS_single_mod3_all","Single photon spectrum",npt,0.,ptmax)) ;
 480   fOutputContainer->Add(new TH1F("PHOS_single_mod4_all","Single photon spectrum",npt,0.,ptmax)) ;
 481   fOutputContainer->Add(new TH1F("PHOS_single_mod5_all","Single photon spectrum",npt,0.,ptmax)) ;
 482
 483   fOutputContainer->Add(new TH1F("PHOS_single_mod1_disp","Single photon spectrum",npt,0.,ptmax)) ;
 484   fOutputContainer->Add(new TH1F("PHOS_single_mod2_disp","Single photon spectrum",npt,0.,ptmax)) ;
 485   fOutputContainer->Add(new TH1F("PHOS_single_mod3_disp","Single photon spectrum",npt,0.,ptmax)) ;
 486   fOutputContainer->Add(new TH1F("PHOS_single_mod4_disp","Single photon spectrum",npt,0.,ptmax)) ;
 487   fOutputContainer->Add(new TH1F("PHOS_single_mod5_disp","Single photon spectrum",npt,0.,ptmax)) ;
 488
 489   fOutputContainer->Add(new TH1F("PHOS_single_mod1_neutral","Single photon spectrum",npt,0.,ptmax)) ;
 490   fOutputContainer->Add(new TH1F("PHOS_single_mod2_neutral","Single photon spectrum",npt,0.,ptmax)) ;
 491   fOutputContainer->Add(new TH1F("PHOS_single_mod3_neutral","Single photon spectrum",npt,0.,ptmax)) ;
 492   fOutputContainer->Add(new TH1F("PHOS_single_mod4_neutral","Single photon spectrum",npt,0.,ptmax)) ;
 493   fOutputContainer->Add(new TH1F("PHOS_single_mod5_neutral","Single photon spectrum",npt,0.,ptmax)) ;
 494
 495   fOutputContainer->Add(new TH1F("PHOS_single_mod1_dispneutral","Single photon spectrum",npt,0.,ptmax)) ;
 496   fOutputContainer->Add(new TH1F("PHOS_single_mod2_dispneutral","Single photon spectrum",npt,0.,ptmax)) ;
 497   fOutputContainer->Add(new TH1F("PHOS_single_mod3_dispneutral","Single photon spectrum",npt,0.,ptmax)) ;
 498   fOutputContainer->Add(new TH1F("PHOS_single_mod4_dispneutral","Single photon spectrum",npt,0.,ptmax)) ;
 499   fOutputContainer->Add(new TH1F("PHOS_single_mod5_dispneutral","Single photon spectrum",npt,0.,ptmax)) ;
 500
 501   fOutputContainer->Add(new TH1F("PHOS_single_mod1_dist1","Single photon spectrum",npt,0.,ptmax)) ;
 502   fOutputContainer->Add(new TH1F("PHOS_single_mod2_dist1","Single photon spectrum",npt,0.,ptmax)) ;
 503   fOutputContainer->Add(new TH1F("PHOS_single_mod3_dist1","Single photon spectrum",npt,0.,ptmax)) ;
 504   fOutputContainer->Add(new TH1F("PHOS_single_mod4_dist1","Single photon spectrum",npt,0.,ptmax)) ;
 505   fOutputContainer->Add(new TH1F("PHOS_single_mod5_dist1","Single photon spectrum",npt,0.,ptmax)) ;
 506
 507   fOutputContainer->Add(new TH1F("PHOS_single_mod1_dist2","Single photon spectrum",npt,0.,ptmax)) ;
 508   fOutputContainer->Add(new TH1F("PHOS_single_mod2_dist2","Single photon spectrum",npt,0.,ptmax)) ;
 509   fOutputContainer->Add(new TH1F("PHOS_single_mod3_dist2","Single photon spectrum",npt,0.,ptmax)) ;
 510   fOutputContainer->Add(new TH1F("PHOS_single_mod4_dist2","Single photon spectrum",npt,0.,ptmax)) ;
 511   fOutputContainer->Add(new TH1F("PHOS_single_mod5_dist2","Single photon spectrum",npt,0.,ptmax)) ;
 512
 513   fOutputContainer->Add(new TH3F("EMCAL_Ep_0","PHOS E/p ratio",24,0.,24,48,0.,48,200,0.,2.)) ;
 514   fOutputContainer->Add(new TH3F("EMCAL_Ep_1","PHOS E/p ratio",24,0.,24,48,0.,48,200,0.,2.)) ;
 515   fOutputContainer->Add(new TH3F("EMCAL_Ep_2","PHOS E/p ratio",24,0.,24,48,0.,48,200,0.,2.)) ;
 516   fOutputContainer->Add(new TH3F("EMCAL_Ep_3","PHOS E/p ratio",24,0.,24,48,0.,48,200,0.,2.)) ;
 517   fOutputContainer->Add(new TH3F("EMCAL_EpPE","PHOS E/p vs E vs p",200,0.,2.,npt,0.,ptmax,npt,0.,ptmax)) ;
 518
 519   fOutputContainer->Add(new TH3F("EMCAL_beyond","EMCAL clusters not in EMCAL",200,-300.,300.,200,-500.,0.,200,-100.,100.)) ;
 520
 521   fOutputContainer->Add(new TH3F("EMCAL_mod0_th1","Inv.Mass distr. per channel",24,0.,24,48,0.,48,200,0.,1.)) ;
 522   fOutputContainer->Add(new TH3F("EMCAL_mod1_th1","Inv.Mass distr. per channel",24,0.,24,48,0.,48,200,0.,1.)) ;
 523   fOutputContainer->Add(new TH3F("EMCAL_mod2_th1","Inv.Mass distr. per channel",24,0.,24,48,0.,48,200,0.,1.)) ;
 524   fOutputContainer->Add(new TH3F("EMCAL_mod3_th1","Inv.Mass distr. per channel",24,0.,24,48,0.,48,200,0.,1.)) ;
 525
 526   fOutputContainer->Add(new TH3F("EMCAL_mod0_th2","Inv.Mass distr. per channel",24,0.,24,48,0.,48,200,0.,1.)) ;
 527   fOutputContainer->Add(new TH3F("EMCAL_mod1_th2","Inv.Mass distr. per channel",24,0.,24,48,0.,48,200,0.,1.)) ;
 528   fOutputContainer->Add(new TH3F("EMCAL_mod2_th2","Inv.Mass distr. per channel",24,0.,24,48,0.,48,200,0.,1.)) ;
 529   fOutputContainer->Add(new TH3F("EMCAL_mod3_th2","Inv.Mass distr. per channel",24,0.,24,48,0.,48,200,0.,1.)) ;
 530
 531   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_mod0_single","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 532   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_mod1_single","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 533   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_mod2_single","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 534   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_mod3_single","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 535
 536   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_mod0_all","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 537   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_mod1_all","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 538   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_mod2_all","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 539   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_mod3_all","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 540   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_mod4_all","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 541   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_mod5_all","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 542   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_mod0_Disp","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 543   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_mod1_Disp","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 544   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_mod2_Disp","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 545   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_mod3_Disp","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 546   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_mod4_Disp","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 547   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_mod5_Disp","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 548   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_mod0_TOF","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 549   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_mod1_TOF","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 550   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_mod2_TOF","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 551   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_mod3_TOF","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 552   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_mod4_TOF","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 553   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_mod5_TOF","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 554   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_mod0_Neutral","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 555   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_mod1_Neutral","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 556   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_mod2_Neutral","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 557   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_mod3_Neutral","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 558   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_mod4_Neutral","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 559   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_mod5_Neutral","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 560   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_mod0_DispNeutral","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 561   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_mod1_DispNeutral","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 562   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_mod2_DispNeutral","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 563   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_mod3_DispNeutral","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 564   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_mod4_DispNeutral","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 565   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_mod5_DispNeutral","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 566
 567   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_mod0_all","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 568   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_mod1_all","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 569   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_mod2_all","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 570   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_mod3_all","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 571   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_mod4_all","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 572   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_mod5_all","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 573   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_mod0_Disp","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 574   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_mod1_Disp","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 575   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_mod2_Disp","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 576   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_mod3_Disp","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 577   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_mod4_Disp","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 578   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_mod5_Disp","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 579   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_mod0_TOF","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 580   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_mod1_TOF","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 581   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_mod2_TOF","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 582   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_mod3_TOF","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 583   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_mod4_TOF","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 584   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_mod5_TOF","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 585   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_mod0_Neutral","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 586   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_mod1_Neutral","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 587   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_mod2_Neutral","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 588   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_mod3_Neutral","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 589   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_mod4_Neutral","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 590   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_mod5_Neutral","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 591   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_mod0_DispNeutral","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 592   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_mod1_DispNeutral","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 593   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_mod2_DispNeutral","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 594   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_mod3_DispNeutral","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 595   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_mod4_DispNeutral","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 596   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_mod5_DispNeutral","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 597
 598   //MC info
 599   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_unitEta","Primary #pi^{0}",npt,0.,ptmax)) ;
 600   fOutputContainer->Add(new TH1F("hMC_CaloConv_allPi0","Primary #pi^{0}",npt,0.,ptmax)) ;
 601   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0PHOSacc","#pi^{0} decayed in PHOS acc",npt,0.,ptmax)) ;
 602   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0EMCALacc","#pi^{0} decayed in EMCAL acc",npt,0.,ptmax)) ;
 603   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_PHOS_conv","#pi^{0} decayed in PHOS acc asnd conv. photon",npt,0.,ptmax)) ;
 604   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_EMCAL_conv","#pi^{0} decayed in EMCAL acc asnd conv. photon",npt,0.,ptmax)) ;
 605   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_bothphot_conv","#pi^{0} both photons converted",npt,0.,ptmax)) ;
 606   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0__convPhotInCalo","#pi^{0} photon in calo converted",npt,0.,ptmax)) ;
 607   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0_PHOSacc","#pi^{0} photon converted and V0 found",npt,0.,ptmax)) ;
 608   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0_EMCALacc","#pi^{0} photon converted and V0 found",npt,0.,ptmax)) ;
 609   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0_PHOSclu","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
 610   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0_EMCALclu_ptRec","#pi^{0} V0 and cluster in EMCAL found (rec pt)",npt,0.,ptmax)) ;
 611   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0_PHOSclu_ptRec","#pi^{0} V0 and cluster in PHOS found(rec pt)",npt,0.,ptmax)) ;
 612   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0_EMCALclu","#pi^{0} V0 and cluster in EMCAL found",npt,0.,ptmax)) ;
 613   fOutputContainer->Add(new TH2F("hMC_CaloConv_pi0_v0_PHOSclu_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax)) ;
 614   fOutputContainer->Add(new TH2F("hMC_CaloConv_pi0_v0_EMCALclu_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax)) ;
 615
 616   fOutputContainer->Add(new TH1F("hMC_CaloConv_phot","Primary photons",npt,0.,ptmax)) ;
 617   fOutputContainer->Add(new TH1F("hMC_CaloConv_gammaPHOSacc","Photons in PHOS acc",npt,0.,ptmax)) ;
 618   fOutputContainer->Add(new TH1F("hMC_CaloConv_gammaEMCALacc","Photons in EMCAL acc",npt,0.,ptmax)) ;
 619   fOutputContainer->Add(new TH1F("hMC_CaloConv_gamma_conv","Converted photons",npt,0.,ptmax)) ;
 620   fOutputContainer->Add(new TH1F("hMC_CaloConv_gamma_v0","Converted photons with V0",npt,0.,ptmax)) ;
 621   fOutputContainer->Add(new TH2F("hMC_CaloConv_gamma_v0_devsE","Converted photons with V0",200,-1.,1.,npt,0.,ptmax)) ;
 622   fOutputContainer->Add(new TH1F("hMC_CaloConv_gamma_PHOSclu","Photons with cluster in PHOS",npt,0.,ptmax)) ;
 623   fOutputContainer->Add(new TH1F("hMC_CaloConv_gamma_EMCALclu","Photons with cluster in EMCAL",npt,0.,ptmax)) ;
 624   fOutputContainer->Add(new TH1F("hMC_CaloConv_gamma_PHOSclu_recE","Photons with cluster in PHOS",npt,0.,ptmax)) ;
 625   fOutputContainer->Add(new TH1F("hMC_CaloConv_gamma_EMCALclu_recE","Photons with cluster in EMCAL",npt,0.,ptmax)) ;
 626   fOutputContainer->Add(new TH2F("hMC_CaloConv_gamma_PHOSclu_devsE","Photons with cluster in PHOS",200,-1.,1.,npt,0.,ptmax)) ;
 627   fOutputContainer->Add(new TH2F("hMC_CaloConv_gamma_EMCALclu_devsE","Photons with cluster in EMCAL",200,-1.,1.,npt,0.,ptmax)) ;
 628
 629
 630   fOutputContainer->Add(new TH3F("All_chi2_eta_pt","MC chi2 vs eta vs phi",100,0.,100.,200,-2.,2.,npt,0.,ptmax)) ;
 631   fOutputContainer->Add(new TH2F("All_w_vs_m","MC w vs m",300,0.,TMath::Pi(),400,0.,1.)) ;
 632   fOutputContainer->Add(new TH3F("MC_V0_pt_eta_phi","MC pt vs eta vs phi",npt,0.,ptmax,200,-2.,2.,200,0.,TMath::TwoPi())) ;
 633   fOutputContainer->Add(new TH3F("MC_V0_m_eta_pt","MC m vs eta vs phi",400,0.,1.,200,-2.,2.,npt,0.,ptmax)) ;
 634   fOutputContainer->Add(new TH3F("MC_V0_chi2_eta_pt","MC chi2 vs eta vs phi",100,0.,100.,200,-2.,2.,npt,0.,ptmax)) ;
 635   fOutputContainer->Add(new TH2F("MC_V0_w_vs_m","MC w vs m",300,0.,TMath::Pi(),400,0.,1.)) ;
 636
 637   fOutputContainer->SetName(GetName());
 638 }
 639 //______________________________________________________________________
 640  void AliAnalysisTaskCaloConv::InitGeometry()
 641 {
 642   //If not done yet, create Geometry for PHOS and EMCAL
 643   //and read misalignment matrixes from ESD/AOD (AOD not implemented yet)
 644   //
 645
 646    if(fPHOSgeom && fEMCALgeom){ //already initialized
 647      return ;
 648    }
 649
 650    AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent()) ;
 651    if(!esd )
 652      AliFatal("Can not read geometry matrixes from ESD/AOD: NO ESD") ;
 653    if(!fPHOSgeom){//reading PHOS matrixes
 654      fPHOSgeom = new AliPHOSGeoUtils("IHEP","");
 655      for(Int_t mod=0; mod<5; mod++){
 656        if(esd){
 657          const TGeoHMatrix* m=esd->GetPHOSMatrix(mod) ;
 658          if(m)
 659            fPHOSgeom->SetMisalMatrix(m, mod) ;
 660        }
 661      }
 662    }
 663    if(!fEMCALgeom){
 664      fEMCALgeom = new AliEMCALGeoUtils("EMCAL_FIRSTYEAR");
 665      for(Int_t mod=0; mod < (fEMCALgeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){ //<---Gustavo, could you check???
 666        if(esd){
 667          const TGeoHMatrix* m=esd->GetEMCALMatrix(mod) ;
 668          if(m)
 669            fEMCALgeom->SetMisalMatrix(m, mod) ;
 670        }
 671      }
 672    }
 673 }
 674 //________________________________________________________________
 675 void AliAnalysisTaskCaloConv::SelectPHOSPhotons(){
 676
 677   // Loop over all CaloClusters
 678   if(fPHOSEvent)
 679     fPHOSEvent->Clear() ;
 680   else
 681     fPHOSEvent = new TClonesArray("TLorentzVector",10) ;
 682   Int_t inPHOS = 0;
 683   TLorentzVector pi0 ;
 684
 685   //vertex
 686   Double_t vtx[3];
 687   vtx[0] = fESDEvent->GetPrimaryVertex()->GetX();
 688   vtx[1] = fESDEvent->GetPrimaryVertex()->GetY();
 689   vtx[2] = fESDEvent->GetPrimaryVertex()->GetZ();
 690   for (Int_t i=0; i<fESDEvent->GetNumberOfCaloClusters(); i++) {
 691     AliESDCaloCluster * clu = fESDEvent->GetCaloCluster(i);
 692     if(!clu->IsPHOS())
 693       continue ;
 694     TLorentzVector p ;
 695     clu ->GetMomentum(p ,vtx);
 696     Bool_t isNeutral = kTRUE ;
 697     Bool_t isDispOK = kTRUE ;
 698     Bool_t isTOFOK = kTRUE ;
 699     Int_t iMod,iX,iZ ;
 700
 701     isNeutral = clu->GetEmcCpvDistance()>5. ;  //To be improved
 702     isDispOK = kFALSE ;
 703     Double_t l0=clu->GetM02(),l1=clu->GetM20() ;
 704     if(l1>= 0   && l0>= 0   && l1 < 0.1 && l0 < 0.1) isDispOK=kFALSE ;
 705     if(l1>= 0   && l0 > 0.5 && l1 < 0.1 && l0 < 1.5) isDispOK=kTRUE ;
 706     if(l1>= 0   && l0 > 2.0 && l1 < 0.1 && l0 < 2.7) isDispOK=kFALSE ;
 707     if(l1>= 0   && l0 > 2.7 && l1 < 0.1 && l0 < 4.0) isDispOK=kFALSE ;
 708     if(l1 > 0.1 && l1 < 0.7 && l0 > 0.7 && l0 < 2.1) isDispOK=kTRUE ;
 709     if(l1 > 0.1 && l1 < 0.3 && l0 > 3.0 && l0 < 5.0) isDispOK=kFALSE  ;
 710     if(l1 > 0.3 && l1 < 0.7 && l0 > 2.5 && l0 < 4.0) isDispOK=kFALSE ;
 711     if(l1 > 0.7 && l1 < 1.3 && l0 > 1.0 && l0 < 1.6) isDispOK=kTRUE ;
 712     if(l1 > 0.7 && l1 < 1.3 && l0 > 1.6 && l0 < 3.5) isDispOK=kTRUE ;
 713     if(l1 > 1.3 && l1 < 3.5 && l0 > 1.3 && l0 < 3.5) isDispOK=kTRUE ;
 714
 715     Float_t xyz[3] = {0,0,0};
 716     clu->GetPosition(xyz);   //Global position in ALICE system
 717     TVector3 global(xyz) ;
 718     Int_t relid[4] ;
 719     if(!fPHOSgeom->GlobalPos2RelId(global,relid)){
 720       FillHistogram("PHOS_beyond",xyz[0],xyz[1],xyz[2]) ;
 721       printf("PHOS_beyond: x=%f, y=%f, z=%f \n",xyz[0],xyz[1],xyz[2]) ;
 722       continue ;
 723     }
 724     iMod=relid[0] ;
 725     iX=relid[2];
 726     iZ=relid[3] ;
 727
 728     p.SetBit(kCaloPIDdisp,isDispOK) ;
 729     p.SetBit(kCaloPIDtof,isTOFOK) ;
 730     p.SetBit(kCaloPIDneutral,isNeutral) ;
 731     p.SetBit(BIT(16+iMod),kTRUE) ;
 732     new((*fPHOSEvent)[inPHOS]) TLorentzVector(p) ;
 733     inPHOS++ ;
 734
 735
 736     //Single photon spectra
 737     Double_t pt= p.Pt() ;
 738     TString skey="PHOS_single_"; skey+="mod" ; skey+=iMod ; skey+="_all" ;
 739     FillHistogram(skey,pt) ;
 740     if(isDispOK){
 741       skey="PHOS_single_"; skey+="mod" ; skey+=iMod ; skey+="_disp" ;
 742       FillHistogram(skey,pt) ;
 743     }
 744     if(isNeutral){
 745       skey="PHOS_single_"; skey+="mod" ; skey+=iMod ; skey+="_neutral" ;
 746       FillHistogram(skey,pt) ;
 747     }
 748     if(isNeutral && isDispOK){
 749       skey="PHOS_single_"; skey+="mod" ; skey+=iMod ; skey+="_dispneutral" ;
 750       FillHistogram(skey,pt) ;
 751     }
 752     //Distance to bad channel
 753     if(clu->GetDistanceToBadChannel()>2.2){
 754       skey="PHOS_single_"; skey+="mod" ; skey+=iMod ; skey+="_dist1" ;
 755       FillHistogram(skey,pt) ;
 756     }
 757     if(clu->GetDistanceToBadChannel()>4.4){
 758       skey="PHOS_single_"; skey+="mod" ; skey+=iMod ; skey+="_dist2" ;
 759       FillHistogram(skey,pt) ;
 760     }
 761
 762     //Fill histogams for calibration
 763     if(clu->E()<fPi0Thresh1 ) continue;
 764     for(Int_t iconv=0;iconv<fConvEvent->GetEntriesFast();iconv++){
 765       TLorentzVector *gammaConv=static_cast<TLorentzVector*>(fConvEvent->At(iconv)) ;
 766       pi0=*gammaConv+p;
 767       skey="PHOS_"; skey+="mod" ; skey+=iMod ; skey+="_th1" ;
 768       FillHistogram(skey,iX-0.5, iZ-0.5,pi0.M());
 769       if(isNeutral){
 770         skey="PHOS_"; skey+="mod"; skey+=iMod ; skey+="_th2" ;
 771         FillHistogram(skey,iX-0.5, iZ-0.5,pi0.M());
 772       }
 773     }
 774   }
 775 }
 776 //____________________________________________________________
 777 void AliAnalysisTaskCaloConv::SelectEMCALPhotons(){
 778
 779   // Loop over all CaloClusters
 780   if(fEMCALEvent)
 781     fEMCALEvent->Clear() ;
 782   else
 783     fEMCALEvent = new TClonesArray("TLorentzVector",10) ;
 784   Int_t inEMCAL = 0 ; //, inEMCALRecal=0;
 785   TLorentzVector pi0 ;
 786
 787   //vertex
 788   Double_t vtx[3];
 789   vtx[0] = fESDEvent->GetPrimaryVertex()->GetX();
 790   vtx[1] = fESDEvent->GetPrimaryVertex()->GetY();
 791   vtx[2] = fESDEvent->GetPrimaryVertex()->GetZ();
 792   Int_t nEMCAL=0 ; //For QA
 793   for (Int_t i=0; i<fESDEvent->GetNumberOfCaloClusters(); i++) {
 794     AliESDCaloCluster * clu = fESDEvent->GetCaloCluster(i);
 795     if(clu->IsPHOS())
 796       continue ;
 797     TLorentzVector p ;
 798     TLorentzVector pRecal ;
 799     clu ->GetMomentum(p ,vtx);
 800     Bool_t isNeutral = kTRUE ;
 801     Bool_t isDispOK = kTRUE ;
 802     Bool_t isTOFOK = kTRUE ;
 803     Int_t iMod,iX,iZ ;
 804
 805     if(clu->E()>0.1 ){
 806       nEMCAL++ ;
 807       isNeutral = clu->GetEmcCpvDistance()>10. ;  //To be improved
 808       if(clu->GetTOF()>550.e-9 && clu->GetTOF()<750.e-9)
 809         isTOFOK=kTRUE ;
 810       else
 811         isTOFOK=kFALSE ;
 812       Float_t phi = p.Phi();
 813       if(phi < 0) phi+=TMath::TwoPi();
 814       Int_t absId = -1;
 815       fEMCALgeom->GetAbsCellIdFromEtaPhi(p.Eta(),phi, absId);
 816       iMod=fEMCALgeom->GetSuperModuleNumber(absId) ;
 817       Int_t imod = -1, iTower = -1, iIphi = -1, iIeta = -1, iphi = -1, ieta = -1;
 818       fEMCALgeom->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
 819       fEMCALgeom->GetCellPhiEtaIndexInSModule(imod,iTower, iIphi, iIeta,iphi,ieta);
 820       if(imod<0 || imod>5){
 821         printf("EMCAL: Beyond the geometry!\n") ;
 822         printf("phi=%f, eta=%f, absId=%d, SM=%d \n",p.Eta(),phi, absId, imod) ;
 823         continue ;
 824       }
 825
 826       iX=iphi+1 ;
 827       iZ=ieta+1 ;
 828       p.SetBit(kCaloPIDdisp,isDispOK) ;
 829       p.SetBit(kCaloPIDtof,isTOFOK) ;
 830       p.SetBit(kCaloPIDneutral,isNeutral) ;
 831       p.SetBit(BIT(17+imod),kTRUE) ;
 832       new((*fEMCALEvent)[inEMCAL]) TLorentzVector(p) ;
 833       inEMCAL++ ;
 834
 835       //Fill histograms for recalibration
 836       if(clu->E()<fPi0Thresh1) continue ;
 837       for(Int_t iconv=0;iconv<fConvEvent->GetEntriesFast();iconv++){
 838         TLorentzVector *gammaConv=static_cast<TLorentzVector*>(fConvEvent->At(iconv)) ;
 839         pi0=*gammaConv+p;
 840         TString skey="EMCAL_"; skey+="mod" ; skey+=iMod ; skey+="_th1" ;
 841         FillHistogram(skey,iX-0.5, iZ-0.5,pi0.M());
 842         if(clu->E()>fPi0Thresh2){
 843           skey="EMCAL_"; skey+="mod"; skey+=iMod ; skey+="_th2" ;
 844           FillHistogram(skey,iX-0.5, iZ-0.5,pi0.M());
 845         }
 846       }
 847     }
 848   }
 849 }
 850 //______________________________________________________________________
 851 void AliAnalysisTaskCaloConv::SelectConvPhotons(){
 852   //Fill list of conversion photons
 853   //that is scan v0s and select photon-like
 854
 855   //set some constants
 856   const Double_t cutSigmaMass=0.0001;  //Constraint on photon mass
 857   const Bool_t useImprovedVertex=kTRUE ; //Use verted with converted photon?
 858 //  const Double_t zrSlope = TMath::Tan(2*TMath::ATan(TMath::Exp(-fetaCut)));
 859   const Double_t zrSlope = TMath::Tan(2*TMath::ATan(TMath::Exp(-1.2)));
 860   const Double_t zOffset = 7.;
 861
 862   if(!fConvEvent)
 863     fConvEvent = new TClonesArray("TLorentzVector",10) ;
 864   else
 865     fConvEvent->Clear() ;
 866
 867   //No primary vertex in event => scip
 868   if(fESDEvent->GetPrimaryVertex()->GetNContributors()<=0) {
 869     return;
 870   }
 871
 872   Int_t inConv=0 ;
 873   for(Int_t iv0=0; iv0<fESDEvent->GetNumberOfV0s();iv0++){
 874     AliESDv0 * v0 = fESDEvent->GetV0(iv0) ;
 875
 876     AliESDtrack * pos = fESDEvent->GetTrack(v0->GetPindex()) ;
 877     AliESDtrack * neg = fESDEvent->GetTrack(v0->GetNindex()) ;
 878     const AliExternalTrackParam * paramPos = v0->GetParamP() ;
 879     const AliExternalTrackParam * paramNeg = v0->GetParamN() ;
 880     if(pos->GetSign() <0){//change tracks
 881       pos=neg ;
 882       neg=fESDEvent->GetTrack(v0->GetPindex()) ;
 883       paramPos=paramNeg ;
 884       paramNeg=v0->GetParamP() ;
 885     }
 886     AliKFParticle negKF(*paramNeg,11);
 887     AliKFParticle posKF(*paramPos,-11);
 888     AliKFParticle photKF(negKF,posKF) ;
 889     photKF.SetMassConstraint(0,cutSigmaMass);
 890     TLorentzVector photLV;
 891     photLV.SetXYZM(negKF.Px()+posKF.Px(),negKF.Py()+posKF.Py(),negKF.Pz()+negKF.Pz(),0.) ;
 892
 893     if(useImprovedVertex){
 894       AliKFVertex primaryVertexImproved(*(fESDEvent->GetPrimaryVertex()));
 895       primaryVertexImproved+=photKF;
 896       photKF.SetProductionVertex(primaryVertexImproved);
 897     }
 898
 899     Double_t m, width ;
 900     photKF.GetMass(m,width);
 901
 902     //Parameters for correction function
 903     Double_t a[3]={photLV.Pt(),photLV.Eta(),m} ;
 904     fConvCFCont->Fill(a,0) ;
 905
 906     //select V0 finder
 907     Bool_t isOnFly=kTRUE ;
 908     //select V0Finder
 909     if (v0->GetOnFlyStatus()){
 910       fConvCFCont->Fill(a,1) ;
 911     }
 912     else{
 913       isOnFly=kFALSE ;
 914       fConvCFCont->Fill(a,2) ;
 915     }
 916
 917     //remove like sign pairs
 918     if(pos->GetSign() == neg->GetSign()){
 919       continue ;
 920     }
 921     if(isOnFly)
 922       fConvCFCont->Fill(a,3) ;
 923     else
 924       fConvCFCont->Fill(a,4) ;
 925
 926
 927     if( !(pos->GetStatus() & AliESDtrack::kTPCrefit) ||
 928         !(neg->GetStatus() & AliESDtrack::kTPCrefit) ){
 929       continue;
 930     }
 931     if(isOnFly)
 932       fConvCFCont->Fill(a,5) ;
 933     else
 934       fConvCFCont->Fill(a,6) ;
 935
 936     Bool_t isKink=kFALSE ;
 937     if( neg->GetKinkIndex(0) > 0 ||
 938         pos->GetKinkIndex(0) > 0) {
 939       isKink=kTRUE;
 940     }
 941     if(!isKink){
 942       if(isOnFly)
 943         fConvCFCont->Fill(a,7) ;
 944       else
 945         fConvCFCont->Fill(a,8) ;
 946     }
 947
 948     Bool_t isdEdx=kTRUE;
 949     if( fESDpid->NumberOfSigmasTPC(pos,AliPID::kElectron)<fnSigmaBelowElectronLine ||
 950         fESDpid->NumberOfSigmasTPC(pos,AliPID::kElectron)>fnSigmaAboveElectronLine ||
 951         fESDpid->NumberOfSigmasTPC(neg,AliPID::kElectron)<fnSigmaBelowElectronLine ||
 952         fESDpid->NumberOfSigmasTPC(neg,AliPID::kElectron)>fnSigmaAboveElectronLine ){
 953          isdEdx=kFALSE;
 954     }
 955     if( pos->P()>fpnSigmaAbovePionLine){
 956       if(fESDpid->NumberOfSigmasTPC(pos,AliPID::kElectron)>fnSigmaBelowElectronLine &&
 957          fESDpid->NumberOfSigmasTPC(pos,AliPID::kElectron)<fnSigmaAboveElectronLine&&
 958          fESDpid->NumberOfSigmasTPC(pos,AliPID::kPion)<fnSigmaAbovePionLine){
 959           isdEdx=kFALSE;
 960       }
 961     }
 962     if( neg->P()>fpnSigmaAbovePionLine){
 963       if(fESDpid->NumberOfSigmasTPC(neg,AliPID::kElectron)>fnSigmaBelowElectronLine &&
 964          fESDpid->NumberOfSigmasTPC(neg,AliPID::kElectron)<fnSigmaAboveElectronLine&&
 965          fESDpid->NumberOfSigmasTPC(neg,AliPID::kPion)<fnSigmaAbovePionLine){
 966           isdEdx=kFALSE;
 967       }
 968     }
 969
 970     //Check the pid probability
 971     Bool_t isProb=kTRUE ;
 972     Double_t posProbArray[10];
 973     Double_t negProbArray[10];
 974     neg->GetTPCpid(negProbArray);
 975     pos->GetTPCpid(posProbArray);
 976     if(negProbArray[AliPID::kElectron]<fprobCut || posProbArray[AliPID::kElectron]<fprobCut){
 977       isProb=kFALSE ;
 978     }
 979     if(!isKink && isProb){
 980       if(isOnFly)
 981         fConvCFCont->Fill(a,9) ;
 982       else
 983         fConvCFCont->Fill(a,10) ;
 984     }
 985
 986     Double_t v0x,v0y,v0z;
 987     v0->GetXYZ(v0x,v0y,v0z) ;
 988     Double_t r=TMath::Sqrt(v0x*v0x + v0y*v0y) ;
 989     if(r>fmaxR){ // cuts on distance from collision point
 990       continue;
 991     }
 992     Bool_t isStrictR=kFALSE ;
 993     if(r<120.)
 994       isStrictR=kTRUE ;
 995     if(isOnFly)
 996       fConvCFCont->Fill(a,11) ;
 997     else
 998       fConvCFCont->Fill(a,12) ;
 999
1000
1001     if((TMath::Abs(v0z)*zrSlope)-zOffset > r ){ // cuts out regions where we do not reconstruct
1002       continue;
1003     }
1004     if(isOnFly)
1005       fConvCFCont->Fill(a,13) ;
1006     else
1007       fConvCFCont->Fill(a,14) ;
1008
1009     if(TMath::Abs(v0z) > fmaxZ ){ // cuts out regions where we do not reconstruct
1010       continue;
1011     }
1012     Bool_t isStrictZ=kFALSE ;
1013     const Double_t strictZcut=0.9 ;
1014     if((TMath::Abs(v0z)*zrSlope)-zOffset < strictZcut*r || TMath::Abs(v0z) < strictZcut*fmaxZ)
1015       isStrictZ=kTRUE ;
1016
1017     if(isOnFly)
1018       fConvCFCont->Fill(a,15) ;
1019     else
1020       fConvCFCont->Fill(a,16) ;
1021
1022
1023     if(photKF.GetNDF()<=0){
1024       continue;
1025     }
1026     if(isOnFly)
1027       fConvCFCont->Fill(a,17) ;
1028     else
1029       fConvCFCont->Fill(a,18) ;
1030
1031     Double_t chi2V0 = photKF.GetChi2()/photKF.GetNDF();
1032     FillHistogram("All_chi2_eta_pt",chi2V0,photLV.Eta(),photLV.Pt()) ;
1033
1034     if(chi2V0 > fchi2CutConversion || chi2V0 <=0){
1035       continue;
1036     }
1037     Bool_t isStrictChi=kFALSE ;
1038     if(chi2V0 < 0.7*fchi2CutConversion && chi2V0 >0){
1039       isStrictChi=kTRUE;
1040     }
1041
1042     if(isOnFly)
1043       fConvCFCont->Fill(a,19) ;
1044     else
1045       fConvCFCont->Fill(a,20) ;
1046
1047     const Double_t wideEtaCut=1.2 ;
1048     if(TMath::Abs(photLV.Eta())> wideEtaCut){
1049       continue;
1050     }
1051     Bool_t isWideEta=kTRUE ;
1052     if(TMath::Abs(photLV.Eta())< fetaCut){
1053       isWideEta=kFALSE;
1054     }
1055
1056
1057     if(photLV.Pt()<fptCut){
1058       continue;
1059     }
1060     if(isOnFly)
1061       fConvCFCont->Fill(a,21) ;
1062     else
1063       fConvCFCont->Fill(a,22) ;
1064
1065
1066     Double_t w=PlanarityAngle(paramPos,paramNeg) ;
1067     Bool_t isPlanarityCut = (0.08-0.22*w > m || 0.15*(w-2.4)>m) ;
1068     FillHistogram("All_w_vs_m",w,m) ;
1069
1070     photLV.SetBit(kConvOnFly,isOnFly) ;
1071     photLV.SetBit(kConvKink,isKink) ;
1072     photLV.SetBit(kConvdEdx,isdEdx) ;
1073     photLV.SetBit(kConvProb,isProb) ;
1074     photLV.SetBit(kConvR,isStrictR) ;
1075     photLV.SetBit(kConvZR,isStrictZ) ;
1076     photLV.SetBit(kConvNDF,isStrictChi) ;
1077     photLV.SetBit(kConvEta,isWideEta) ;
1078     photLV.SetBit(kConvPlan,isPlanarityCut) ;
1079
1080     new((*fConvEvent)[inConv]) TLorentzVector(photLV) ;
1081     inConv++ ;
1082
1083     //Single photon spectrum
1084     Double_t pt=photLV.Pt() ;
1085     if(isOnFly){
1086       //Default
1087       if(!isKink && isdEdx && isProb && !isWideEta)
1088         FillHistogram("Single_conv_OnFly",pt) ;
1089       if(isdEdx && isProb && !isWideEta)
1090         FillHistogram("Single_conv_On_Kink",pt) ;
1091       if(!isKink && isProb && !isWideEta)
1092         FillHistogram("Single_conv_On_dEdx",pt) ;
1093       if(!isKink && isdEdx && !isWideEta)
1094         FillHistogram("Single_conv_On_Prob",pt) ;
1095       if(!isKink && isdEdx && isProb && !isWideEta && isStrictR)
1096         FillHistogram("Single_conv_On_R120",pt) ;
1097       if(!isKink && isdEdx && isProb && !isWideEta && isStrictZ)
1098         FillHistogram("Single_conv_On_Z",pt) ;
1099       if(!isKink && isdEdx && isProb && !isWideEta && isStrictChi)
1100         FillHistogram("Single_conv_On_chi",pt) ;
1101       if(!isKink && isdEdx && isProb)
1102         FillHistogram("Single_conv_On_Eta",pt) ;
1103       if(!isKink && isdEdx && isProb && !isWideEta && isPlanarityCut)
1104         FillHistogram("Single_conv_On_Wcut",pt) ;
1105     }
1106     else{
1107       if(!isKink && isdEdx && isProb && !isWideEta)
1108         FillHistogram("Single_conv_Offline",pt) ;
1109       if(isdEdx && isProb && !isWideEta)
1110         FillHistogram("Single_conv_Off_Kink",pt) ;
1111       if(!isKink && isProb && !isWideEta)
1112         FillHistogram("Single_conv_Off_dEdx",pt) ;
1113       if(!isKink && isdEdx && !isWideEta)
1114         FillHistogram("Single_conv_Off_Prob",pt) ;
1115       if(!isKink && isdEdx && isProb && !isWideEta && isStrictR)
1116         FillHistogram("Single_conv_Off_R120",pt) ;
1117       if(!isKink && isdEdx && isProb && !isWideEta && isStrictZ)
1118         FillHistogram("Single_conv_Off_Z",pt) ;
1119       if(!isKink && isdEdx && isProb && !isWideEta && isStrictChi)
1120         FillHistogram("Single_conv_Off_chi",pt) ;
1121       if(!isKink && isdEdx && isProb)
1122         FillHistogram("Single_conv_Off_Eta",pt) ;
1123       if(!isKink && isdEdx && isProb && !isWideEta && isPlanarityCut)
1124         FillHistogram("Single_conv_Off_Wcut",pt) ;
1125     }
1126
1127     //Fill MC information
1128     if(fStack){
1129       TParticle * negativeMC = fStack->Particle(TMath::Abs(neg->GetLabel()));
1130       TParticle * positiveMC = fStack->Particle(TMath::Abs(pos->GetLabel()));
1131
1132       if(negativeMC && positiveMC){
1133         if(negativeMC->GetMother(0) != positiveMC->GetMother(0))
1134           continue ;
1135       }
1136
1137       if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){
1138         continue;
1139       }
1140       if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){
1141         continue;
1142       }
1143
1144       TParticle * v0Gamma = fStack->Particle(negativeMC->GetMother(0));
1145
1146       if(negativeMC->GetUniqueID() != 5 || positiveMC->GetUniqueID() !=5){ // id 5 is conversion
1147         continue;
1148       }
1149       if(v0Gamma->GetPdgCode() == 22){
1150         FillHistogram("MC_V0_pt_eta_phi",v0Gamma->Pt(),v0Gamma->Eta(),v0Gamma->Phi()) ;
1151         FillHistogram("MC_V0_m_eta_pt",m,v0Gamma->Eta(),v0Gamma->Pt()) ;
1152         FillHistogram("MC_V0_chi2_eta_pt",chi2V0,v0Gamma->Eta(),v0Gamma->Pt()) ;
1153         FillHistogram("MC_V0_w_vs_m",w,m) ;
1154       }
1155     }
1156   }
1157 }
1158 //______________________________________________________________________
1159 void AliAnalysisTaskCaloConv::FillRealMixed(){
1160   // Fills real (same event) and Mixed (different events) inv.mass dsitributions
1161   // Moves current event to the list of stored events at the end
1162
1163   Double_t vtx[3];
1164   vtx[0] = fESDEvent->GetPrimaryVertex()->GetX();
1165   vtx[1] = fESDEvent->GetPrimaryVertex()->GetY();
1166   vtx[2] = fESDEvent->GetPrimaryVertex()->GetZ();
1167
1168   //Vtx class z-bin
1169   Int_t zvtx = (Int_t)((vtx[2]+10.)/2.) ;
1170   if(zvtx<0)zvtx=0 ;
1171   if(zvtx>9)zvtx=9 ;
1172
1173   Double_t run = fESDEvent->GetRunNumber()+0.5;
1174   TString trigClasses = fESDEvent->GetFiredTriggerClasses();
1175   if(trigClasses.Contains("CINT1B-ABCE-NOPF-ALL"))
1176     FillHistogram("hRunTrigger",run,0.5) ;
1177   else
1178     FillHistogram("hRunTrigger",run,1.5) ;
1179
1180   FillHistogram("hEvents",0.5) ;
1181   FillHistogram("hVtxBin",zvtx-0.5) ;
1182   FillHistogram("hRunEvents",run) ;
1183
1184   //check if containers for mixed events exist
1185   //create new if necessary
1186   if(!fPHOSEvents[zvtx]) fPHOSEvents[zvtx]=new TList() ;
1187   if(!fEMCALEvents[zvtx]) fEMCALEvents[zvtx]=new TList() ;
1188   if(!fConvEvents[zvtx]) fConvEvents[zvtx]=new TList() ;
1189
1190   Int_t nPHOS=fPHOSEvent->GetEntriesFast() ;
1191   Int_t nEMCAL=fEMCALEvent->GetEntriesFast() ;
1192   Int_t nConv = fConvEvent->GetEntriesFast() ;
1193   //Some QA histograms
1194   FillHistogram("hRunConvs",run,double(nConv)) ;
1195   FillHistogram("hRunPHOS", run,double(nPHOS)) ;
1196   FillHistogram("hRunEMCAL",run,double(nEMCAL)) ;
1197
1198   //Fill Real distributions
1199   for(Int_t iPHOS=0; iPHOS<nPHOS;iPHOS++){
1200     TLorentzVector * cal = static_cast<TLorentzVector*>(fPHOSEvent->At(iPHOS)) ;
1201     for(Int_t iConv = 0; iConv<nConv; iConv++){
1202       TLorentzVector * cnv = static_cast<TLorentzVector*>(fConvEvent->At(iConv)) ;
1203       TLorentzVector pi=*cal + *cnv ;
1204       if(cnv->TestBit(kConvOnFly) && !cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta)){
1205         //Non-linearity check
1206         FillHistogram("PHOS_Re_mvsPt_all",pi.M(),pi.Pt()) ;
1207         if(cal->TestBit(kCaloPIDdisp))
1208           FillHistogram("PHOS_Re_mvsPt_Disp",pi.M(),pi.Pt()) ;
1209         if(cal->TestBit(kCaloPIDtof))
1210           FillHistogram("PHOS_Re_mvsPt_TOF",pi.M(),pi.Pt()) ;
1211         if(cal->TestBit(kCaloPIDneutral))
1212           FillHistogram("PHOS_Re_mvsPt_Neutral",pi.M(),pi.Pt()) ;
1213         if(cal->TestBit(kCaloPIDneutral) && cal->TestBit(kCaloPIDdisp))
1214           FillHistogram("PHOS_Re_mvsPt_DispNeutral",pi.M(),pi.Pt()) ;
1215       }
1216       //Vary Conversion cuts
1217       if(cnv->TestBit(kConvOnFly)){
1218         //Default
1219         if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta))
1220           FillHistogram("PHOS_Re_mvsPt_OnFly",pi.M(),pi.Pt()) ;
1221         if(cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta))
1222           FillHistogram("PHOS_Re_mvsPt_On_Kink",pi.M(),pi.Pt()) ;
1223         if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta))
1224           FillHistogram("PHOS_Re_mvsPt_On_dEdx",pi.M(),pi.Pt()) ;
1225         if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && !cnv->TestBit(kConvEta))
1226           FillHistogram("PHOS_Re_mvsPt_On_Prob",pi.M(),pi.Pt()) ;
1227         if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvR))
1228           FillHistogram("PHOS_Re_mvsPt_On_R120",pi.M(),pi.Pt()) ;
1229         if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvZR))
1230           FillHistogram("PHOS_Re_mvsPt_On_Z",pi.M(),pi.Pt()) ;
1231         if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvNDF))
1232           FillHistogram("PHOS_Re_mvsPt_On_chi",pi.M(),pi.Pt()) ;
1233         if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb))
1234           FillHistogram("PHOS_Re_mvsPt_On_Eta",pi.M(),pi.Pt()) ;
1235         if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvPlan))
1236           FillHistogram("PHOS_Re_mvsPt_On_Wcut",pi.M(),pi.Pt()) ;
1237       }
1238       else{
1239         //Default
1240         if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta))
1241           FillHistogram("PHOS_Re_mvsPt_Offline",pi.M(),pi.Pt()) ;
1242         if(cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta))
1243           FillHistogram("PHOS_Re_mvsPt_Off_Kink",pi.M(),pi.Pt()) ;
1244           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta))
1245           FillHistogram("PHOS_Re_mvsPt_Off_dEdx",pi.M(),pi.Pt()) ;
1246         if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && !cnv->TestBit(kConvEta))
1247           FillHistogram("PHOS_Re_mvsPt_Off_Prob",pi.M(),pi.Pt()) ;
1248         if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvR))
1249           FillHistogram("PHOS_Re_mvsPt_Off_R120",pi.M(),pi.Pt()) ;
1250         if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvZR))
1251           FillHistogram("PHOS_Re_mvsPt_Off_Z",pi.M(),pi.Pt()) ;
1252         if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvNDF))
1253           FillHistogram("PHOS_Re_mvsPt_Off_chi",pi.M(),pi.Pt()) ;
1254         if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb))
1255           FillHistogram("PHOS_Re_mvsPt_Off_Eta",pi.M(),pi.Pt()) ;
1256         if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvPlan))
1257           FillHistogram("PHOS_Re_mvsPt_Off_Wcut",pi.M(),pi.Pt()) ;
1258       }
1259     }
1260   }
1261   //PHOS module-dependent histograms
1262   for(Int_t iPHOS=0; iPHOS<nPHOS;iPHOS++){
1263     TLorentzVector * cal = static_cast<TLorentzVector*>(fPHOSEvent->At(iPHOS)) ;
1264     Int_t mod=1;
1265     while(!cal->TestBit(BIT(16+mod)) && mod<5)mod++ ;
1266     TString base("PHOS_Re_mvsPt_mod") ; base+=mod ;
1267     TString full ;
1268     for(Int_t iConv = 0; iConv<nConv; iConv++){
1269       TLorentzVector * cnv = static_cast<TLorentzVector*>(fConvEvent->At(iConv)) ;
1270       TLorentzVector pi=*cal + *cnv ;
1271       full=base ; full+="_single" ;
1272       if(cnv->TestBit(kConvOnFly) && !cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta)){
1273         FillHistogram(full,pi.M(),cal->Pt()) ;
1274         full=base ; full+="_all" ;
1275         FillHistogram(full,pi.M(),pi.Pt()) ;
1276         if(cal->TestBit(kCaloPIDdisp)){
1277           full=base ; full+="_Disp" ;
1278           FillHistogram(full,pi.M(),pi.Pt()) ;
1279         }
1280         if(cal->TestBit(kCaloPIDtof)){
1281           full=base ; full+="_TOF" ;
1282           FillHistogram(full,pi.M(),pi.Pt()) ;
1283         }
1284           if(cal->TestBit(kCaloPIDneutral)){
1285           full=base ; full+="_Neutral" ;
1286           FillHistogram(full,pi.M(),pi.Pt()) ;
1287         }
1288         if(cal->TestBit(kCaloPIDneutral) && cal->TestBit(kCaloPIDdisp)){
1289           full=base ; full+="_DispNeutral" ;
1290             FillHistogram(full,pi.M(),pi.Pt()) ;
1291         }
1292       }
1293     }
1294   }
1295   for(Int_t iEMCAL=0; iEMCAL<nEMCAL;iEMCAL++){
1296     TLorentzVector * cal = static_cast<TLorentzVector*>(fEMCALEvent->At(iEMCAL)) ;
1297     Int_t mod=0;
1298     while(!cal->TestBit(BIT(17+mod)) && mod<6)mod++ ;
1299     TString base("EMCAL_Re_mvsPt_mod") ; base+=mod ;
1300     TString full ;
1301     for(Int_t iConv = 0; iConv<nConv; iConv++){
1302       TLorentzVector * cnv = static_cast<TLorentzVector*>(fConvEvent->At(iConv)) ;
1303       TLorentzVector pi=*cal + *cnv ;
1304       full=base+"_single" ;
1305       if(cnv->TestBit(kConvOnFly) && !cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta)){
1306         FillHistogram(full,pi.M(),cal->Pt()) ;
1307         full=base+"_all" ;
1308         FillHistogram(full,pi.M(),pi.Pt()) ;
1309         if(cal->TestBit(kCaloPIDdisp)){
1310           full=base+"_Disp" ;
1311           FillHistogram(full,pi.M(),pi.Pt()) ;
1312         }
1313         if(cal->TestBit(kCaloPIDtof)){
1314           full=base+"_TOF" ;
1315           FillHistogram(full,pi.M(),pi.Pt()) ;
1316         }
1317         if(cal->TestBit(kCaloPIDneutral)){
1318           full=base+"_Neutral" ;
1319           FillHistogram(full,pi.M(),pi.Pt()) ;
1320         }
1321         if(cal->TestBit(kCaloPIDneutral) && cal->TestBit(kCaloPIDdisp)){
1322           full=base+"_DispNeutral" ;
1323           FillHistogram(full,pi.M(),pi.Pt()) ;
1324         }
1325       }
1326       //Vary Conversion cuts
1327       if(cnv->TestBit(kConvOnFly)){
1328         //Default
1329         if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta))
1330           FillHistogram("EMCAL_Re_mvsPt_OnFly",pi.M(),pi.Pt()) ;
1331         if(cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta))
1332           FillHistogram("EMCAL_Re_mvsPt_On_Kink",pi.M(),pi.Pt()) ;
1333         if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta))
1334           FillHistogram("EMCAL_Re_mvsPt_On_dEdx",pi.M(),pi.Pt()) ;
1335         if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && !cnv->TestBit(kConvEta))
1336           FillHistogram("EMCAL_Re_mvsPt_On_Prob",pi.M(),pi.Pt()) ;
1337         if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvR))
1338           FillHistogram("EMCAL_Re_mvsPt_On_R120",pi.M(),pi.Pt()) ;
1339         if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvZR))
1340           FillHistogram("EMCAL_Re_mvsPt_On_Z",pi.M(),pi.Pt()) ;
1341         if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvNDF))
1342           FillHistogram("EMCAL_Re_mvsPt_On_chi",pi.M(),pi.Pt()) ;
1343         if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb))
1344           FillHistogram("EMCAL_Re_mvsPt_On_Eta",pi.M(),pi.Pt()) ;
1345         if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvPlan))
1346           FillHistogram("EMCAL_Re_mvsPt_On_Wcut",pi.M(),pi.Pt()) ;
1347       }
1348       else{
1349         //Default
1350         if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta))
1351           FillHistogram("EMCAL_Re_mvsPt_Offline",pi.M(),pi.Pt()) ;
1352         if(cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta))
1353           FillHistogram("EMCAL_Re_mvsPt_Off_Kink",pi.M(),pi.Pt()) ;
1354           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta))
1355           FillHistogram("EMCAL_Re_mvsPt_Off_dEdx",pi.M(),pi.Pt()) ;
1356         if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && !cnv->TestBit(kConvEta))
1357           FillHistogram("EMCAL_Re_mvsPt_Off_Prob",pi.M(),pi.Pt()) ;
1358         if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvR))
1359           FillHistogram("EMCAL_Re_mvsPt_Off_R120",pi.M(),pi.Pt()) ;
1360         if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvZR))
1361           FillHistogram("EMCAL_Re_mvsPt_Off_Z",pi.M(),pi.Pt()) ;
1362         if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvNDF))
1363           FillHistogram("EMCAL_Re_mvsPt_Off_chi",pi.M(),pi.Pt()) ;
1364         if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb))
1365           FillHistogram("EMCAL_Re_mvsPt_Off_Eta",pi.M(),pi.Pt()) ;
1366         if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvPlan))
1367           FillHistogram("EMCAL_Re_mvsPt_Off_Wcut",pi.M(),pi.Pt()) ;
1368       }
1369     }
1370   }
1371   //Now fill mixed
1372   TList * prevPHOS = fPHOSEvents[zvtx] ;
1373   TList * prevEMCAL = fEMCALEvents[zvtx] ;
1374   TList * prevConv = fConvEvents[zvtx] ;
1375
1376   //PHOS
1377   for(Int_t iPHOS=0; iPHOS<nPHOS;iPHOS++){
1378     TLorentzVector * cal = static_cast<TLorentzVector*>(fPHOSEvent->At(iPHOS)) ;
1379     for(Int_t ev=0; ev<prevConv->GetSize();ev++){
1380       TClonesArray * mixConv = static_cast<TClonesArray*>(prevConv->At(ev)) ;
1381       for(Int_t iConv = 0; iConv<mixConv->GetEntriesFast(); iConv++){
1382         TLorentzVector * cnv = static_cast<TLorentzVector*>(mixConv->At(iConv)) ;
1383         TLorentzVector pi=*cal + *cnv ;
1384         if(cnv->TestBit(kConvOnFly) && !cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta)){
1385           FillHistogram("PHOS_Mi_mvsPt_all",pi.M(),pi.Pt()) ;
1386           if(cal->TestBit(kCaloPIDdisp))
1387             FillHistogram("PHOS_Mi_mvsPt_Disp",pi.M(),pi.Pt()) ;
1388           if(cal->TestBit(kCaloPIDtof))
1389             FillHistogram("PHOS_Mi_mvsPt_TOF",pi.M(),pi.Pt()) ;
1390             if(cal->TestBit(kCaloPIDneutral))
1391             FillHistogram("PHOS_Mi_mvsPt_Neutral",pi.M(),pi.Pt()) ;
1392           if(cal->TestBit(kCaloPIDneutral) && cal->TestBit(kCaloPIDdisp))
1393             FillHistogram("PHOS_Mi_mvsPt_DispNeutral",pi.M(),pi.Pt()) ;
1394         }
1395         //Vary Conversion cuts
1396         if(cnv->TestBit(kConvOnFly)){
1397           //Default
1398           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta))
1399             FillHistogram("PHOS_Mi_mvsPt_OnFly",pi.M(),pi.Pt()) ;
1400           if(cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta))
1401             FillHistogram("PHOS_Mi_mvsPt_On_Kink",pi.M(),pi.Pt()) ;
1402           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta))
1403             FillHistogram("PHOS_Mi_mvsPt_On_dEdx",pi.M(),pi.Pt()) ;
1404           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && !cnv->TestBit(kConvEta))
1405             FillHistogram("PHOS_Mi_mvsPt_On_Prob",pi.M(),pi.Pt()) ;
1406           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvR))
1407             FillHistogram("PHOS_Mi_mvsPt_On_R120",pi.M(),pi.Pt()) ;
1408           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvZR))
1409             FillHistogram("PHOS_Mi_mvsPt_On_Z",pi.M(),pi.Pt()) ;
1410           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvNDF))
1411             FillHistogram("PHOS_Mi_mvsPt_On_chi",pi.M(),pi.Pt()) ;
1412           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb))
1413             FillHistogram("PHOS_Mi_mvsPt_On_Eta",pi.M(),pi.Pt()) ;
1414           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvPlan))
1415             FillHistogram("PHOS_Mi_mvsPt_On_Wcut",pi.M(),pi.Pt()) ;
1416         }
1417         else{
1418           //Default
1419           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta))
1420             FillHistogram("PHOS_Mi_mvsPt_Offline",pi.M(),pi.Pt()) ;
1421           if(cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta))
1422             FillHistogram("PHOS_Mi_mvsPt_Off_Kink",pi.M(),pi.Pt()) ;
1423             if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta))
1424             FillHistogram("PHOS_Mi_mvsPt_Off_dEdx",pi.M(),pi.Pt()) ;
1425           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && !cnv->TestBit(kConvEta))
1426             FillHistogram("PHOS_Mi_mvsPt_Off_Prob",pi.M(),pi.Pt()) ;
1427           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvR))
1428             FillHistogram("PHOS_Mi_mvsPt_Off_R120",pi.M(),pi.Pt()) ;
1429           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvZR))
1430             FillHistogram("PHOS_Mi_mvsPt_Off_Z",pi.M(),pi.Pt()) ;
1431           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvNDF))
1432             FillHistogram("PHOS_Mi_mvsPt_Off_chi",pi.M(),pi.Pt()) ;
1433           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb))
1434             FillHistogram("PHOS_Mi_mvsPt_Off_Eta",pi.M(),pi.Pt()) ;
1435           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvPlan))
1436             FillHistogram("PHOS_Mi_mvsPt_Off_Wcut",pi.M(),pi.Pt()) ;
1437         }
1438       }
1439     }
1440   }
1441   for(Int_t iConv = 0; iConv<nConv; iConv++){
1442     TLorentzVector * cnv = static_cast<TLorentzVector*>(fConvEvent->At(iConv)) ;
1443     for(Int_t ev=0; ev<prevPHOS->GetSize();ev++){
1444       TClonesArray * mixPHOS = static_cast<TClonesArray*>(prevPHOS->At(ev)) ;
1445       for(Int_t iPHOS=0; iPHOS<mixPHOS->GetEntriesFast();iPHOS++){
1446         TLorentzVector * cal = static_cast<TLorentzVector*>(mixPHOS->At(iPHOS)) ;
1447         TLorentzVector pi=*cal + *cnv ;
1448         if(cnv->TestBit(kConvOnFly) && !cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta)){
1449           FillHistogram("PHOS_Mi_mvsPt_all",pi.M(),pi.Pt()) ;
1450           if(cal->TestBit(kCaloPIDdisp))
1451             FillHistogram("PHOS_Mi_mvsPt_Disp",pi.M(),pi.Pt()) ;
1452           if(cal->TestBit(kCaloPIDtof))
1453             FillHistogram("PHOS_Mi_mvsPt_TOF",pi.M(),pi.Pt()) ;
1454           if(cal->TestBit(kCaloPIDneutral))
1455             FillHistogram("PHOS_Mi_mvsPt_Neutral",pi.M(),pi.Pt()) ;
1456           if(cal->TestBit(kCaloPIDneutral) && cal->TestBit(kCaloPIDdisp))
1457             FillHistogram("PHOS_Mi_mvsPt_DispNeutral",pi.M(),pi.Pt()) ;
1458         }
1459         //Vary Conversion cuts
1460         if(cnv->TestBit(kConvOnFly)){
1461           //Default
1462           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta))
1463             FillHistogram("PHOS_Mi_mvsPt_OnFly",pi.M(),pi.Pt()) ;
1464           if(cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta))
1465             FillHistogram("PHOS_Mi_mvsPt_On_Kink",pi.M(),pi.Pt()) ;
1466           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta))
1467             FillHistogram("PHOS_Mi_mvsPt_On_dEdx",pi.M(),pi.Pt()) ;
1468           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && !cnv->TestBit(kConvEta))
1469             FillHistogram("PHOS_Mi_mvsPt_On_Prob",pi.M(),pi.Pt()) ;
1470           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvR))
1471             FillHistogram("PHOS_Mi_mvsPt_On_R120",pi.M(),pi.Pt()) ;
1472           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvZR))
1473             FillHistogram("PHOS_Mi_mvsPt_On_Z",pi.M(),pi.Pt()) ;
1474           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvNDF))
1475             FillHistogram("PHOS_Mi_mvsPt_On_chi",pi.M(),pi.Pt()) ;
1476           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb))
1477             FillHistogram("PHOS_Mi_mvsPt_On_Eta",pi.M(),pi.Pt()) ;
1478           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvPlan))
1479             FillHistogram("PHOS_Mi_mvsPt_On_Wcut",pi.M(),pi.Pt()) ;
1480         }
1481         else{
1482           //Default
1483           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta))
1484             FillHistogram("PHOS_Mi_mvsPt_Offline",pi.M(),pi.Pt()) ;
1485           if(cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta))
1486             FillHistogram("PHOS_Mi_mvsPt_Off_Kink",pi.M(),pi.Pt()) ;
1487             if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta))
1488             FillHistogram("PHOS_Mi_mvsPt_Off_dEdx",pi.M(),pi.Pt()) ;
1489           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && !cnv->TestBit(kConvEta))
1490             FillHistogram("PHOS_Mi_mvsPt_Off_Prob",pi.M(),pi.Pt()) ;
1491           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvR))
1492             FillHistogram("PHOS_Mi_mvsPt_Off_R120",pi.M(),pi.Pt()) ;
1493           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvZR))
1494             FillHistogram("PHOS_Mi_mvsPt_Off_Z",pi.M(),pi.Pt()) ;
1495           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvNDF))
1496             FillHistogram("PHOS_Mi_mvsPt_Off_chi",pi.M(),pi.Pt()) ;
1497           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb))
1498             FillHistogram("PHOS_Mi_mvsPt_Off_Eta",pi.M(),pi.Pt()) ;
1499           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvPlan))
1500             FillHistogram("PHOS_Mi_mvsPt_Off_Wcut",pi.M(),pi.Pt()) ;
1501         }
1502       }
1503     }
1504   }
1505
1506   //PHOS module dependent
1507   for(Int_t iPHOS=0; iPHOS<nPHOS;iPHOS++){
1508     TLorentzVector * cal = static_cast<TLorentzVector*>(fPHOSEvent->At(iPHOS)) ;
1509     Int_t mod=1;
1510     while(!cal->TestBit(BIT(16+mod)) && mod<5)mod++ ;
1511     TString base("PHOS_Mi_mvsPt_mod") ; base+=mod ;
1512     TString full ;
1513     for(Int_t ev=0; ev<prevConv->GetSize();ev++){
1514       TClonesArray * mixConv = static_cast<TClonesArray*>(prevConv->At(ev)) ;
1515       for(Int_t iConv = 0; iConv<mixConv->GetEntriesFast(); iConv++){
1516         TLorentzVector * cnv = static_cast<TLorentzVector*>(mixConv->At(iConv)) ;
1517         TLorentzVector pi=*cal + *cnv ;
1518         full=base+"_all" ;
1519         if(cnv->TestBit(kConvOnFly) && !cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta)){
1520           FillHistogram(full,pi.M(),pi.Pt()) ;
1521           if(cal->TestBit(kCaloPIDdisp)){
1522             full=base+"_Disp" ;
1523             FillHistogram(full,pi.M(),pi.Pt()) ;
1524           }
1525           if(cal->TestBit(kCaloPIDtof)){
1526             full=base+"_TOF" ;
1527             FillHistogram(full,pi.M(),pi.Pt()) ;
1528           }
1529           if(cal->TestBit(kCaloPIDneutral)){
1530             full=base+"_Neutral" ;
1531             FillHistogram(full,pi.M(),pi.Pt()) ;
1532           }
1533           if(cal->TestBit(kCaloPIDneutral) && cal->TestBit(kCaloPIDdisp)){
1534             full=base+"_DispNeutral" ;
1535             FillHistogram(full,pi.M(),pi.Pt()) ;
1536           }
1537         }
1538       }
1539     }
1540   }
1541   for(Int_t iConv = 0; iConv<nConv; iConv++){
1542     TLorentzVector * cnv = static_cast<TLorentzVector*>(fConvEvent->At(iConv)) ;
1543     for(Int_t ev=0; ev<prevPHOS->GetSize();ev++){
1544       TClonesArray * mixPHOS = static_cast<TClonesArray*>(prevPHOS->At(ev)) ;
1545       for(Int_t iPHOS=0; iPHOS<mixPHOS->GetEntriesFast();iPHOS++){
1546         TLorentzVector * cal = static_cast<TLorentzVector*>(mixPHOS->At(iPHOS)) ;
1547         Int_t mod=1;
1548         while(!cal->TestBit(BIT(16+mod)) && mod<5)mod++ ;
1549         TString base("PHOS_Mi_mvsPt_mod") ; base+=mod ;
1550         TString full ;
1551         TLorentzVector pi=*cal + *cnv ;
1552         full=base+"_all" ;
1553         if(cnv->TestBit(kConvOnFly) && !cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta)){
1554           FillHistogram(full,pi.M(),pi.Pt()) ;
1555           if(cal->TestBit(kCaloPIDdisp)){
1556             full=base+"_Disp" ;
1557             FillHistogram(full,pi.M(),pi.Pt()) ;
1558           }
1559           if(cal->TestBit(kCaloPIDtof)){
1560             full=base+"_TOF" ;
1561             FillHistogram(full,pi.M(),pi.Pt()) ;
1562           }
1563           if(cal->TestBit(kCaloPIDneutral)){
1564             full=base+"_Neutral" ;
1565             FillHistogram(full,pi.M(),pi.Pt()) ;
1566           }
1567           if(cal->TestBit(kCaloPIDneutral) && cal->TestBit(kCaloPIDdisp)){
1568             full=base+"_DispNeutral" ;
1569             FillHistogram(full,pi.M(),pi.Pt()) ;
1570           }
1571         }
1572       }
1573     }
1574   }
1575
1576
1577   //EMCAL
1578   for(Int_t iEMCAL=0; iEMCAL<nEMCAL;iEMCAL++){
1579     TLorentzVector * cal = static_cast<TLorentzVector*>(fEMCALEvent->At(iEMCAL)) ;
1580     Int_t mod=0;
1581     while(!cal->TestBit(BIT(17+mod)) && mod<6)mod++ ;
1582     TString base("EMCAL_Mi_mvsPt_mod") ; base+=mod ;
1583     TString full ;
1584     for(Int_t ev=0; ev<prevConv->GetSize();ev++){
1585       TClonesArray * mixConv = static_cast<TClonesArray*>(prevConv->At(ev)) ;
1586       for(Int_t iConv = 0; iConv<mixConv->GetEntriesFast(); iConv++){
1587         TLorentzVector * cnv = static_cast<TLorentzVector*>(mixConv->At(iConv)) ;
1588         TLorentzVector pi=*cal + *cnv ;
1589         full=base+"_all" ;
1590         if(cnv->TestBit(kConvOnFly) && !cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta)){
1591           FillHistogram(full,pi.M(),pi.Pt()) ;
1592           if(cal->TestBit(kCaloPIDdisp)){
1593             full=base+"_Disp" ;
1594             FillHistogram(full,pi.M(),pi.Pt()) ;
1595           }
1596           if(cal->TestBit(kCaloPIDtof)){
1597             full=base+"_TOF" ;
1598             FillHistogram(full,pi.M(),pi.Pt()) ;
1599           }
1600           if(cal->TestBit(kCaloPIDneutral)){
1601             full=base+"_Neutral" ;
1602             FillHistogram(full,pi.M(),pi.Pt()) ;
1603           }
1604           if(cal->TestBit(kCaloPIDneutral) && cal->TestBit(kCaloPIDdisp)){
1605             full=base+"_Neutral" ;
1606             FillHistogram(full,pi.M(),pi.Pt()) ;
1607           }
1608         }
1609         if(cnv->TestBit(kConvOnFly)){
1610           //Default
1611           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta))
1612             FillHistogram("EMCAL_Mi_mvsPt_OnFly",pi.M(),pi.Pt()) ;
1613           if(cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta))
1614             FillHistogram("EMCAL_Mi_mvsPt_On_Kink",pi.M(),pi.Pt()) ;
1615           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta))
1616             FillHistogram("EMCAL_Mi_mvsPt_On_dEdx",pi.M(),pi.Pt()) ;
1617           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && !cnv->TestBit(kConvEta))
1618             FillHistogram("EMCAL_Mi_mvsPt_On_Prob",pi.M(),pi.Pt()) ;
1619           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvR))
1620             FillHistogram("EMCAL_Mi_mvsPt_On_R120",pi.M(),pi.Pt()) ;
1621           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvZR))
1622             FillHistogram("EMCAL_Mi_mvsPt_On_Z",pi.M(),pi.Pt()) ;
1623           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvNDF))
1624             FillHistogram("EMCAL_Mi_mvsPt_On_chi",pi.M(),pi.Pt()) ;
1625           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb))
1626             FillHistogram("EMCAL_Mi_mvsPt_On_Eta",pi.M(),pi.Pt()) ;
1627           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvPlan))
1628             FillHistogram("EMCAL_Mi_mvsPt_On_Wcut",pi.M(),pi.Pt()) ;
1629         }
1630         else{
1631           //Default
1632           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta))
1633             FillHistogram("EMCAL_Mi_mvsPt_Offline",pi.M(),pi.Pt()) ;
1634           if(cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta))
1635             FillHistogram("EMCAL_Mi_mvsPt_Off_Kink",pi.M(),pi.Pt()) ;
1636           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta))
1637             FillHistogram("EMCAL_Mi_mvsPt_Off_dEdx",pi.M(),pi.Pt()) ;
1638           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && !cnv->TestBit(kConvEta))
1639             FillHistogram("EMCAL_Mi_mvsPt_Off_Prob",pi.M(),pi.Pt()) ;
1640           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvR))
1641             FillHistogram("EMCAL_Mi_mvsPt_Off_R120",pi.M(),pi.Pt()) ;
1642           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvZR))
1643             FillHistogram("EMCAL_Mi_mvsPt_Off_Z",pi.M(),pi.Pt()) ;
1644           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvNDF))
1645             FillHistogram("EMCAL_Mi_mvsPt_Off_chi",pi.M(),pi.Pt()) ;
1646           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb))
1647             FillHistogram("EMCAL_Mi_mvsPt_Off_Eta",pi.M(),pi.Pt()) ;
1648           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvPlan))
1649             FillHistogram("EMCAL_Mi_mvsPt_Off_Wcut",pi.M(),pi.Pt()) ;
1650         }
1651       }
1652     }
1653   }
1654   for(Int_t iConv = 0; iConv<nConv; iConv++){
1655     TLorentzVector * cnv = static_cast<TLorentzVector*>(fConvEvent->At(iConv)) ;
1656     for(Int_t ev=0; ev<prevEMCAL->GetSize();ev++){
1657       TClonesArray * mixEMCAL = static_cast<TClonesArray*>(prevEMCAL->At(ev)) ;
1658       for(Int_t iEMCAL=0; iEMCAL<mixEMCAL->GetEntriesFast();iEMCAL++){
1659         TLorentzVector * cal = static_cast<TLorentzVector*>(mixEMCAL->At(iEMCAL)) ;
1660         Int_t mod=0;
1661         while(!cal->TestBit(BIT(17+mod)) && mod<6)mod++ ;
1662         TString base("EMCAL_Mi_mvsPt_mod") ; base+=mod ;
1663         TString full ;
1664         TLorentzVector pi=*cal + *cnv ;
1665         full=base+"_all" ;
1666         if(cnv->TestBit(kConvOnFly) && !cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta)){
1667           FillHistogram(full,pi.M(),pi.Pt()) ;
1668           if(cal->TestBit(kCaloPIDdisp)){
1669             full=base+"_Disp" ;
1670             FillHistogram(full,pi.M(),pi.Pt()) ;
1671           }
1672           if(cal->TestBit(kCaloPIDtof)){
1673             full=base+"_TOF" ;
1674             FillHistogram(full,pi.M(),pi.Pt()) ;
1675           }
1676           if(cal->TestBit(kCaloPIDneutral)){
1677             full=base+"_Neutral" ;
1678               FillHistogram(full,pi.M(),pi.Pt()) ;
1679           }
1680           if(cal->TestBit(kCaloPIDneutral) && cal->TestBit(kCaloPIDdisp)){
1681             full=base+"_DispNeutral" ;
1682             FillHistogram(full,pi.M(),pi.Pt()) ;
1683           }
1684         }
1685         if(cnv->TestBit(kConvOnFly)){
1686           //Default
1687           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta))
1688             FillHistogram("EMCAL_Mi_mvsPt_OnFly",pi.M(),pi.Pt()) ;
1689           if(cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta))
1690             FillHistogram("EMCAL_Mi_mvsPt_On_Kink",pi.M(),pi.Pt()) ;
1691           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta))
1692             FillHistogram("EMCAL_Mi_mvsPt_On_dEdx",pi.M(),pi.Pt()) ;
1693           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && !cnv->TestBit(kConvEta))
1694             FillHistogram("EMCAL_Mi_mvsPt_On_Prob",pi.M(),pi.Pt()) ;
1695           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvR))
1696             FillHistogram("EMCAL_Mi_mvsPt_On_R120",pi.M(),pi.Pt()) ;
1697           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvZR))
1698             FillHistogram("EMCAL_Mi_mvsPt_On_Z",pi.M(),pi.Pt()) ;
1699           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvNDF))
1700             FillHistogram("EMCAL_Mi_mvsPt_On_chi",pi.M(),pi.Pt()) ;
1701           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb))
1702             FillHistogram("EMCAL_Mi_mvsPt_On_Eta",pi.M(),pi.Pt()) ;
1703           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvPlan))
1704             FillHistogram("EMCAL_Mi_mvsPt_On_Wcut",pi.M(),pi.Pt()) ;
1705         }
1706         else{
1707           //Default
1708           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta))
1709             FillHistogram("EMCAL_Mi_mvsPt_Offline",pi.M(),pi.Pt()) ;
1710           if(cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta))
1711             FillHistogram("EMCAL_Mi_mvsPt_Off_Kink",pi.M(),pi.Pt()) ;
1712           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta))
1713             FillHistogram("EMCAL_Mi_mvsPt_Off_dEdx",pi.M(),pi.Pt()) ;
1714           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && !cnv->TestBit(kConvEta))
1715             FillHistogram("EMCAL_Mi_mvsPt_Off_Prob",pi.M(),pi.Pt()) ;
1716           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvR))
1717             FillHistogram("EMCAL_Mi_mvsPt_Off_R120",pi.M(),pi.Pt()) ;
1718           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvZR))
1719             FillHistogram("EMCAL_Mi_mvsPt_Off_Z",pi.M(),pi.Pt()) ;
1720           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvNDF))
1721             FillHistogram("EMCAL_Mi_mvsPt_Off_chi",pi.M(),pi.Pt()) ;
1722           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb))
1723             FillHistogram("EMCAL_Mi_mvsPt_Off_Eta",pi.M(),pi.Pt()) ;
1724           if(!cnv->TestBit(kConvKink) && cnv->TestBit(kConvdEdx) && cnv->TestBit(kConvProb) && !cnv->TestBit(kConvEta) && cnv->TestBit(kConvPlan))
1725             FillHistogram("EMCAL_Mi_mvsPt_Off_Wcut",pi.M(),pi.Pt()) ;
1726         }
1727       }
1728     }
1729   }
1730
1731
1732
1733   //Now we either add current events to stack or remove
1734   //Here we have some difficulty: conversion and calorimeter photons have different spectra.
1735   //So to correctly reproduce combinatorial background we have to preserve average number of
1736   //photons of both kinds per event. Therefore we should not reject empty PHOS/EMCAL events
1737   //though it will cost some memory. Reject only those events where no photons anywhere
1738
1739   if((fPHOSEvent->GetEntriesFast()>0 || fEMCALEvent->GetEntriesFast()>0) && fConvEvent->GetEntriesFast()>0){
1740     prevPHOS->AddFirst(fPHOSEvent) ;
1741     fPHOSEvent=0;
1742     prevEMCAL->AddFirst(fEMCALEvent) ;
1743     fEMCALEvent=0 ;
1744     prevConv->AddFirst(fConvEvent) ;
1745     fConvEvent=0 ;
1746     if(prevPHOS->GetSize()>100){//Remove redundant events
1747       TClonesArray * tmp = static_cast<TClonesArray*>(prevPHOS->Last()) ;
1748       prevPHOS->RemoveLast() ;
1749       delete tmp ;
1750       tmp = static_cast<TClonesArray*>(prevEMCAL->Last()) ;
1751       prevEMCAL->RemoveLast() ;
1752       delete tmp ;
1753       tmp = static_cast<TClonesArray*>(prevConv->Last()) ;
1754       prevConv->RemoveLast() ;
1755       delete tmp ;
1756     }
1757   }
1758
1759 }
1760 //___________________________________________________________________________
1761 void AliAnalysisTaskCaloConv::ProcessMC(){
1762
1763   //fill histograms for efficiensy etc. calculation
1764   if(!fStack) return ;
1765
1766   const Double_t rcut = 1. ; //cut for primary particles
1767   Double_t vtx[3];
1768   vtx[0] = fESDEvent->GetPrimaryVertex()->GetX();
1769   vtx[1] = fESDEvent->GetPrimaryVertex()->GetY();
1770   vtx[2] = fESDEvent->GetPrimaryVertex()->GetZ();
1771
1772   //---------First pi0s-----------------------------
1773   for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++) {
1774     TParticle* particle = (TParticle *)fStack->Particle(iTracks);
1775     if(particle->GetPdgCode() != 111)
1776       continue ;
1777
1778     if(particle->R() >rcut)
1779       continue ;
1780
1781    Double_t pt = particle->Pt() ;
1782    //Total number of pi0 with creation radius <1 cm
1783    FillHistogram("hMC_CaloConv_allPi0",pt) ;
1784    if(TMath::Abs(particle->Eta())<1.)
1785      FillHistogram("hMC_CaloConv_pi0_unitEta",pt) ;
1786
1787    //Check if one of photons converted
1788    if(particle->GetNDaughters()!=2)
1789      continue ; //Do not account Dalitz decays
1790
1791    TParticle * gamma1 = fStack->Particle(particle->GetFirstDaughter());
1792    TParticle * gamma2 = fStack->Particle(particle->GetLastDaughter());
1793    //Number of pi0s decayed into acceptance
1794    Bool_t inAcc1 = (TMath::Abs(gamma1->Eta())<0.9) ;
1795    Bool_t inAcc2 = (TMath::Abs(gamma2->Eta())<0.9) ;
1796    Int_t mod ;
1797    Double_t x,z ;
1798    Bool_t hitPHOS1 = fPHOSgeom->ImpactOnEmc(gamma1, mod, z,x) ;
1799    Bool_t hitPHOS2 = fPHOSgeom->ImpactOnEmc(gamma2, mod, z,x) ;
1800    Bool_t hitEMCAL1= fEMCALgeom->Impact(gamma1) ;
1801    Bool_t hitEMCAL2= fEMCALgeom->Impact(gamma2) ;
1802
1803    Bool_t goodPair=kFALSE ;
1804    if((inAcc1 && hitPHOS2) || (inAcc2 && hitPHOS1)){
1805      FillHistogram("hMC_CaloConv_pi0PHOSacc",pt) ;
1806      goodPair=kTRUE ;
1807    }
1808    if((inAcc1 && hitEMCAL2) || (inAcc2 && hitEMCAL1)){
1809      FillHistogram("hMC_CaloConv_pi0EMCALacc",pt) ;
1810      goodPair=kTRUE ;
1811    }
1812    if(!goodPair){
1813      continue ;
1814    }
1815
1816    Bool_t converted1 = kFALSE ;
1817    if(gamma1->GetNDaughters()==2){
1818      TParticle * e1=fStack->Particle(gamma1->GetFirstDaughter()) ;
1819      TParticle * e2=fStack->Particle(gamma1->GetLastDaughter()) ;
1820      if(TMath::Abs(e1->GetPdgCode())==11 && TMath::Abs(e2->GetPdgCode())==11){ //conversion
1821        if(e1->R()<180.)
1822          converted1 = kTRUE ;
1823      }
1824    }
1825    Bool_t converted2 = kFALSE ;
1826    if(gamma2->GetNDaughters()==2){
1827      TParticle * e1=fStack->Particle(gamma2->GetFirstDaughter()) ;
1828      TParticle * e2=fStack->Particle(gamma2->GetLastDaughter()) ;
1829      if(TMath::Abs(e1->GetPdgCode())==11 && TMath::Abs(e2->GetPdgCode())==11){ //conversion
1830        if(e1->R()<180.)
1831          converted2 = kTRUE ;
1832      }
1833    }
1834
1835    //Number of pi0s with one photon converted
1836    if((converted1 && !converted2 && hitPHOS2) || (!converted1 && hitPHOS1 && converted2))
1837      FillHistogram("hMC_CaloConv_pi0_PHOS_conv",pt) ;
1838
1839    if((converted1 && !converted2 && hitEMCAL2) || (!converted1 && hitEMCAL1 && converted2))
1840      FillHistogram("hMC_CaloConv_pi0_EMCAL_conv",pt) ;
1841
1842    //Both converted
1843    if(converted1 && converted2) {
1844      FillHistogram("hMC_CaloConv_pi0_bothphot_conv",pt) ;
1845      continue ;
1846    }
1847
1848    //photon pointing calorimeter converted
1849    if((converted1 && hitPHOS1 && !hitEMCAL2) || (converted2 && hitPHOS2 && !hitEMCAL1) ||
1850       (converted1 && hitEMCAL1 && !hitPHOS2) || (converted2 && hitEMCAL2 && !hitPHOS1)){
1851      FillHistogram("hMC_CaloConv_pi0__convPhotInCalo",pt) ;
1852       continue ;
1853    }
1854
1855    //Converted pi0 with v0 and photon PHOS or EMCAL
1856    Bool_t foundV01=kFALSE, foundV02=kFALSE ;
1857    TLorentzVector pConv ;
1858    for(Int_t iv0=0; iv0<fESDEvent->GetNumberOfV0s();iv0++){
1859      AliESDv0 * v0 = fESDEvent->GetV0(iv0) ;
1860
1861      TParticle * negativeMC = fStack->Particle(TMath::Abs(fESDEvent->GetTrack(v0->GetNindex())->GetLabel()));
1862      TParticle * positiveMC = fStack->Particle(TMath::Abs(fESDEvent->GetTrack(v0->GetPindex())->GetLabel()));
1863
1864      if(negativeMC && positiveMC){
1865        if(negativeMC->GetMother(0) != positiveMC->GetMother(0))
1866          continue ;
1867      }
1868
1869      if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){
1870        continue;
1871      }
1872      if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){
1873        continue;
1874      }
1875
1876      TParticle * v0Gamma = fStack->Particle(negativeMC->GetMother(0));
1877      Bool_t same = (v0Gamma == gamma1) ;
1878      TParticle * tmp = v0Gamma ;
1879      while(!same && tmp->GetFirstMother()>=0){
1880        tmp = fStack->Particle(tmp->GetFirstMother());
1881        same = (tmp == gamma1) ;
1882      }
1883      if(same){
1884        foundV01 = kTRUE ;
1885        const AliExternalTrackParam * paramPos = v0->GetParamP() ;
1886        const AliExternalTrackParam * paramNeg = v0->GetParamN() ;
1887        AliKFParticle negKF(*paramNeg,11);
1888        AliKFParticle posKF(*paramPos,-11);
1889        pConv.SetXYZM(negKF.Px()+posKF.Px(),negKF.Py()+posKF.Py(),negKF.Pz()+negKF.Pz(),0.) ;
1890        break ;
1891      }
1892      same = (v0Gamma == gamma2) ;
1893      tmp = v0Gamma ;
1894      while(!same && tmp->GetFirstMother()>=0){
1895        tmp = fStack->Particle(tmp->GetFirstMother());
1896        same = (tmp == gamma2) ;
1897      }
1898      if(same){
1899        foundV02 = kTRUE ;
1900        const AliExternalTrackParam * paramPos = v0->GetParamP() ;
1901        const AliExternalTrackParam * paramNeg = v0->GetParamN() ;
1902        AliKFParticle negKF(*paramNeg,11);
1903        AliKFParticle posKF(*paramPos,-11);
1904        pConv.SetXYZM(negKF.Px()+posKF.Px(),negKF.Py()+posKF.Py(),negKF.Pz()+negKF.Pz(),0.) ;
1905        break ;
1906      }
1907    }
1908
1909    goodPair=kFALSE ;
1910    if((foundV01 && hitPHOS2) || (foundV02 && hitPHOS1)){
1911      FillHistogram("hMC_CaloConv_pi0_v0_PHOSacc",pt) ;
1912      goodPair=kTRUE;
1913    }
1914    if((foundV01 && hitEMCAL2) || (foundV02 && hitEMCAL1)){
1915      FillHistogram("hMC_CaloConv_pi0_v0_EMCALacc",pt) ;
1916      goodPair=kTRUE;
1917    }
1918    if(!goodPair){
1919      continue ;
1920    }
1921
1922    //Converted pi0 with v0 and cluster in PHOS/EMCAL
1923    Bool_t cluInPHOS = kFALSE,cluInEMCAL=kFALSE ;
1924    TLorentzVector pCalo ;
1925    for (Int_t i=0; i<fESDEvent->GetNumberOfCaloClusters(); i++) {
1926      AliESDCaloCluster * clu = fESDEvent->GetCaloCluster(i);
1927      Int_t iprim = clu->GetLabel() ; //# of particle hit PHOS/EMCAL
1928      Bool_t matched = kFALSE ;
1929      while(iprim>=0 ) {
1930        if(iprim==particle->GetFirstDaughter() || iprim==particle->GetLastDaughter()){
1931          matched=kTRUE ;
1932          break ;
1933        }
1934        else{
1935          iprim=fStack->Particle(iprim)->GetFirstMother() ;
1936        }
1937      }
1938      if(!matched)
1939        continue ;
1940      if(clu->IsPHOS() && (hitPHOS1 || hitPHOS2)){
1941        cluInPHOS=kTRUE ;
1942        clu->GetMomentum(pCalo ,vtx);
1943        break ;
1944      }
1945      if(!clu->IsPHOS() && (hitEMCAL1 || hitEMCAL2)){
1946        cluInEMCAL=kTRUE ;
1947        clu->GetMomentum(pCalo ,vtx);
1948        break ;
1949      }
1950    }
1951
1952    if(cluInPHOS){
1953      FillHistogram("hMC_CaloConv_pi0_v0_PHOSclu",pt) ;
1954      Double_t m=(pCalo+pConv).M() ;
1955      Double_t ptm=(pCalo+pConv).Pt() ;
1956      FillHistogram("hMC_CaloConv_pi0_v0_PHOSclu_ptRec",ptm) ;
1957      FillHistogram("hMC_CaloConv_pi0_v0_PHOSclu_mvsPt",m,ptm) ;
1958    }
1959    if(cluInEMCAL){
1960      FillHistogram("hMC_CaloConv_pi0_v0_EMCALclu",pt) ;
1961      Double_t m=(pCalo+pConv).M() ;
1962      Double_t ptm=(pCalo+pConv).Pt() ;
1963      FillHistogram("hMC_CaloConv_pi0_v0_EMCALclu_ptRec",ptm) ;
1964      FillHistogram("hMC_CaloConv_pi0_v0_EMCALclu_mvsPt",m,ptm) ;
1965    }
1966   }
1967
1968
1969
1970    //------------- now photons ----------------
1971    for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++) {
1972      TParticle* particle = (TParticle *)fStack->Particle(iTracks);
1973      if(particle->GetPdgCode() != 22)
1974        continue ;
1975
1976      if(particle->R() >rcut)
1977        continue ;
1978
1979      if(TMath::Abs(particle->Eta())>0.9)
1980        continue ;
1981
1982      Double_t pt = particle->Pt() ;
1983      //Total number of pi0 with creation radius <1 cm
1984      FillHistogram("hMC_CaloConv_phot",pt) ;
1985
1986      Int_t mod ;
1987      Double_t x,z ;
1988      Bool_t hitPHOS = fPHOSgeom->ImpactOnEmc(particle, mod, z,x) ;
1989      Bool_t hitEMCAL= fEMCALgeom->Impact(particle) ;
1990
1991      //Photons in PHOS/EMCAL acceptance
1992      if(hitPHOS)
1993        FillHistogram("hMC_CaloConv_gammaPHOSacc",pt) ;
1994      if(hitEMCAL)
1995        FillHistogram("hMC_CaloConv_gammaEMCALacc",pt) ;
1996
1997      //number of photons converted
1998      Bool_t converted = kFALSE ;
1999      if(particle->GetNDaughters()==2){
2000        TParticle * e1=fStack->Particle(particle->GetFirstDaughter()) ;
2001        TParticle * e2=fStack->Particle(particle->GetLastDaughter()) ;
2002        if(TMath::Abs(e1->GetPdgCode())==11 && TMath::Abs(e2->GetPdgCode())==11){ //conversion
2003          if(e1->R()<180.)
2004            converted = kTRUE ;
2005        }
2006      }
2007      if(converted)
2008        FillHistogram("hMC_CaloConv_gamma_conv",pt) ;
2009
2010      //Converted photons with V0
2011      TLorentzVector pConv ;
2012      Bool_t foundV0=kFALSE ;
2013      for(Int_t iv0=0; iv0<fESDEvent->GetNumberOfV0s();iv0++){
2014        AliESDv0 * v0 = fESDEvent->GetV0(iv0) ;
2015
2016        TParticle * negativeMC = fStack->Particle(TMath::Abs(fESDEvent->GetTrack(v0->GetNindex())->GetLabel()));
2017        TParticle * positiveMC = fStack->Particle(TMath::Abs(fESDEvent->GetTrack(v0->GetPindex())->GetLabel()));
2018
2019        if(negativeMC && positiveMC){
2020          if(negativeMC->GetMother(0) != positiveMC->GetMother(0))
2021            continue ;
2022        }
2023
2024        if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){
2025          continue;
2026        }
2027        if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){
2028          continue;
2029        }
2030
2031        TParticle * v0Gamma = fStack->Particle(negativeMC->GetMother(0));
2032        Bool_t same = (v0Gamma == particle) ;
2033        TParticle * tmp = v0Gamma ;
2034        while(!same && tmp->GetFirstMother()>=0){
2035          tmp = fStack->Particle(tmp->GetFirstMother());
2036          same = (tmp == particle) ;
2037        }
2038        if(same){
2039          foundV0 = kTRUE ;
2040          const AliExternalTrackParam * paramPos = v0->GetParamP() ;
2041          const AliExternalTrackParam * paramNeg = v0->GetParamN() ;
2042          AliKFParticle negKF(*paramNeg,11);
2043          AliKFParticle posKF(*paramPos,-11);
2044          pConv.SetXYZM(negKF.Px()+posKF.Px(),negKF.Py()+posKF.Py(),negKF.Pz()+negKF.Pz(),0.) ;
2045          break ;
2046        }
2047      }
2048      if(foundV0){
2049        FillHistogram("hMC_CaloConv_gamma_v0",pt) ;
2050        FillHistogram("hMC_CaloConv_gamma_v0_devsE",(particle->Energy()-pConv.E())/particle->Energy(),particle->Energy()) ;
2051      }
2052
2053       //Registered in PHOS/EMCAL
2054      Bool_t cluInPHOS = kFALSE,cluInEMCAL=kFALSE ;
2055      TLorentzVector pCalo ;
2056      for (Int_t i=0; i<fESDEvent->GetNumberOfCaloClusters(); i++) {
2057        AliESDCaloCluster * clu = fESDEvent->GetCaloCluster(i);
2058        Int_t iprim = clu->GetLabel() ; //# of particle hit PHOS/EMCAL
2059        Bool_t matched = kFALSE ;
2060        while(iprim>=0 ) {
2061          if(iprim==iTracks){
2062            matched=kTRUE ;
2063            break ;
2064          }
2065          else{
2066            iprim=fStack->Particle(iprim)->GetFirstMother() ;
2067          }
2068        }
2069        if(!matched)
2070          continue ;
2071        if(clu->IsPHOS() && hitPHOS){
2072          cluInPHOS=kTRUE ;
2073          clu->GetMomentum(pCalo ,vtx);
2074          break ;
2075        }
2076        if(!clu->IsPHOS() && hitEMCAL){
2077          cluInEMCAL=kTRUE ;
2078          clu->GetMomentum(pCalo ,vtx);
2079          break ;
2080        }
2081      }
2082
2083      if(cluInPHOS){
2084        FillHistogram("hMC_CaloConv_gamma_PHOSclu",pt) ;
2085        FillHistogram("hMC_CaloConv_gamma_PHOSclu_recE",pCalo.E()) ;
2086        FillHistogram("hMC_CaloConv_gamma_PHOSclu_devsE",(particle->Energy()-pCalo.E())/particle->Energy(),particle->Energy()) ;
2087      }
2088      if(cluInEMCAL){
2089        FillHistogram("hMC_CaloConv_gamma_EMCALclu",pt) ;
2090        FillHistogram("hMC_CaloConv_gamma_EMCALclu_recE",pCalo.E()) ;
2091        FillHistogram("hMC_CaloConv_gamma_EMCALclu_devsE",(particle->Energy()-pCalo.E())/particle->Energy(),particle->Energy()) ;
2092      }
2093   }
2094 }
2095 //_____________________________________________________________________________
2096 void AliAnalysisTaskCaloConv::FillHistogram(const char * key,Double_t x)const{
2097   TH1F * tmp = dynamic_cast<TH1F*>(fOutputContainer->FindObject(key)) ;
2098   if(!tmp)
2099     AliError(Form("can not find histogram <%s> ",key)) ;
2100   tmp->Fill(x) ;
2101 }
2102 //_____________________________________________________________________________
2103 void AliAnalysisTaskCaloConv::FillHistogram(const char * key,Double_t x,Double_t y)const{
2104   TObject * tmp = fOutputContainer->FindObject(key) ;
2105   if(tmp->IsA() == TClass::GetClass("TH1F")){
2106     ((TH1F*)tmp)->Fill(x,y) ;
2107     return ;
2108   }
2109   if(tmp->IsA() == TClass::GetClass("TH2F")){
2110     ((TH2F*)tmp)->Fill(x,y) ;
2111     return ;
2112   }
2113   AliError(Form("Calling FillHistogram with 2 parameters for histo <%s> of type %s",key,tmp->IsA())) ;
2114 }
2115
2116
2117 //_____________________________________________________________________________
2118 void AliAnalysisTaskCaloConv::FillHistogram(const char * key,Double_t x,Double_t y, Double_t z) const{
2119   //Fills 1D histograms with key
2120   TObject * tmp = fOutputContainer->FindObject(key) ;
2121   if(tmp->IsA() == TClass::GetClass("TH2F")){
2122     ((TH2F*)tmp)->Fill(x,y,z) ;
2123     return ;
2124   }
2125   if(tmp->IsA() == TClass::GetClass("TH3F")){
2126     ((TH3F*)tmp)->Fill(x,y,z) ;
2127     return ;
2128   }
2129 }
2130 //______________________________________________________________________________
2131 Double_t AliAnalysisTaskCaloConv::PlanarityAngle(const AliExternalTrackParam * paramPos,const AliExternalTrackParam * paramNeg)const {
2132   //calculate angle between e+e- plain and perpendicular to MF
2133   //We need sign of MagField to calculate orienation
2134
2135   TVector3 u(paramPos->Px()+paramNeg->Px(),paramPos->Py()+paramNeg->Py(),paramPos->Pz()+paramNeg->Pz()) ;
2136   u.Unit() ;
2137   TVector3 vPos(paramPos->Px(),paramPos->Py(),paramPos->Pz()) ;
2138   TVector3 vNeg(paramNeg->Px(),paramNeg->Py(),paramNeg->Pz()) ;
2139   TVector3 v=vPos.Cross(vNeg) ;
2140   TVector3 w = u.Cross(v);
2141   TVector3 z(0,0,1.);
2142   TVector3 ua=u.Cross(z);
2143   Double_t wa = w.Angle(ua);
2144   Double_t mfield=fESDEvent->GetMagneticField() ;
2145   if(mfield>0.)
2146     return wa;      //normal field
2147   else
2148     return TMath::Pi()-wa ; //reverse field
2149
2150 }
2151
2152
2153
2154
2155