PWGGA/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 #include "TDirectory.h"
  28 #include "TLorentzVector.h"
  29 #include "TRandom.h"
  30
  31 // analysis
  32 #include "AliAnalysisManager.h"
  33 #include "AliESDInputHandler.h"
  34 #include "AliAnalysisTaskCaloConv.h"
  35 #include "AliStack.h"
  36 #include "AliLog.h"
  37 #include "AliESDEvent.h"
  38 #include "AliESDpid.h"
  39 #include "AliESDtrackCuts.h"
  40 #include "AliESDtrackCuts.h"
  41 #include "AliCFContainer.h"   // for CF
  42 #include "AliESDCaloCluster.h"
  43 #include "AliPHOSGeoUtils.h"
  44 #include "AliEMCALGeometry.h"
  45 #include "AliCFContainer.h"
  46 #include "AliMCEventHandler.h"
  47 #include "AliMCEvent.h"
  48 #include "AliESDv0.h"
  49 #include "AliKFParticle.h"
  50 #include "AliKFVertex.h"
  51
  52 class Riostream;
  53 class TFile;
  54
  55 ClassImp(AliAnalysisTaskCaloConv)
  56
  57
  58 AliAnalysisTaskCaloConv::AliAnalysisTaskCaloConv():
  59 AliAnalysisTaskSE(),
  60   fESDEvent(NULL),
  61   fESDpid(NULL),
  62   fESDtrackCuts(NULL),
  63   fStack(NULL),
  64   fOutputContainer(NULL),
  65   fCFOutputContainer(NULL),
  66   fConvCFCont(0x0),
  67   fPHOSCFCont(0x0),
  68   fEMCALCFCont(0x0),
  69   fPi0CFCont(0x0),
  70   fCentr(0.),
  71   fTriggerCINT1B(kFALSE),
  72   fToUseCF(kFALSE),
  73   fMinOpeningAngleGhostCut(0.),
  74   fPHOSgeom(0x0),
  75   fEMCALgeom(0x0),
  76   fPi0Thresh1(0.5),
  77   fPi0Thresh2(1.),
  78   fBadDistCutPHOS(3.3),
  79   fBadDistCutEMCAL(6.),
  80   fGammaV0s(),
  81   fGammaPHOS(),
  82   fGammaEMCAL(),
  83   fConvEvent(NULL) ,
  84   fPHOSEvent(NULL),
  85   fEMCALEvent(NULL),
  86   fnSigmaAboveElectronLine(5.),
  87   fnSigmaBelowElectronLine(-3.),
  88   fnSigmaAbovePionLine(0.),
  89   fpnSigmaAbovePionLine(1.),
  90   fprobCut(0.),
  91   fmaxR(180.),
  92   fmaxZ(240.),
  93   fetaCut(0.9),
  94   fptCut(0.02),
  95   fchi2CutConversion(30.)
  96 {
  97   // Default constructor
  98   Int_t nBin=10 ;
  99   for(Int_t i=0;i<nBin;i++){
 100     fPHOSEvents[i]=0 ;
 101     fEMCALEvents[i]=0;
 102     fConvEvents[i]=0;
 103   }
 104   char key[55] ;
 105   for(Int_t i=0; i<6; i++){
 106     snprintf(key,55,"PHOS_BadMap_mod%d",i) ;
 107     fPHOSBadMap[i]=new TH2I(key,"Bad Modules map",64,0.,64.,56,0.,56.) ;
 108   }
 109   for(Int_t i=0; i<10; i++){
 110     snprintf(key,55,"EMCAL_BadMap_mod%d",i) ;
 111     fEMCALBadMap[i] = new TH2I(key,"Bad Modules map",24,0.,24.,48,0.,48.) ;
 112   }
 113 }
 114 AliAnalysisTaskCaloConv::AliAnalysisTaskCaloConv(const char* name):
 115   AliAnalysisTaskSE(name),
 116   fESDEvent(NULL),
 117   fESDpid(NULL),
 118   fESDtrackCuts(NULL),
 119   fStack(NULL),
 120   fOutputContainer(NULL),
 121   fCFOutputContainer(NULL),
 122   fConvCFCont(0x0),
 123   fPHOSCFCont(0x0),
 124   fEMCALCFCont(0x0),
 125   fPi0CFCont(0x0),
 126   fCentr(0.),
 127   fTriggerCINT1B(kFALSE),
 128   fToUseCF(kFALSE),
 129   fMinOpeningAngleGhostCut(0.),
 130   fPHOSgeom(0x0),
 131   fEMCALgeom(0x0),
 132   fPi0Thresh1(0.5),
 133   fPi0Thresh2(1.),
 134   fBadDistCutPHOS(3.3),
 135   fBadDistCutEMCAL(6.),
 136   fGammaV0s(),
 137   fGammaPHOS(),
 138   fGammaEMCAL(),
 139   fConvEvent(NULL) ,
 140   fPHOSEvent(NULL),
 141   fEMCALEvent(NULL),
 142   fnSigmaAboveElectronLine(5.),
 143   fnSigmaBelowElectronLine(-3.),
 144   fnSigmaAbovePionLine(0.),
 145   fpnSigmaAbovePionLine(1.),
 146   fprobCut(0.),
 147   fmaxR(180.),
 148   fmaxZ(240.),
 149   fetaCut(0.9),
 150   fptCut(0.02),
 151   fchi2CutConversion(30.)
 152 {
 153   // Common I/O in slot 0
 154   DefineInput (0, TChain::Class());
 155   DefineOutput(0, TTree::Class());
 156
 157   // Your private output
 158   DefineOutput(1, TList::Class());
 159   DefineOutput(2, TList::Class());  // for CF
 160
 161   Int_t nBin=10 ;
 162   for(Int_t i=0;i<nBin;i++){
 163     fPHOSEvents[i]=0 ;
 164     fEMCALEvents[i]=0;
 165     fConvEvents[i]=0;
 166   }
 167   char key[55] ;
 168   for(Int_t i=0; i<6; i++){
 169     snprintf(key,55,"PHOS_BadMap_mod%d",i) ;
 170     fPHOSBadMap[i]=new TH2I(key,"Bad Modules map",64,0.,64.,56,0.,56.) ;
 171   }
 172   for(Int_t i=0; i<10; i++){
 173     snprintf(key,55,"EMCAL_BadMap_mod%d",i) ;
 174     fEMCALBadMap[i] = new TH2I(key,"Bad Modules map",24,0.,24.,48,0.,48.) ;
 175   }
 176 //  fESDpid = new AliESDpid;
 177 }
 178 //_____________________________________________________
 179 AliAnalysisTaskCaloConv::~AliAnalysisTaskCaloConv()
 180 {
 181   // Remove all pointers
 182
 183   if (AliAnalysisManager::GetAnalysisManager()->GetAnalysisType() !=
 184       AliAnalysisManager::kProofAnalysis) {
 185
 186     if(fOutputContainer){
 187       fOutputContainer->Clear() ;
 188       delete fOutputContainer ;
 189     }
 190     if(fCFOutputContainer){
 191       fCFOutputContainer->Clear() ;
 192       delete fCFOutputContainer ;
 193     }
 194
 195     if(fPHOSgeom){
 196       delete fPHOSgeom ;
 197       fPHOSgeom=0x0 ;
 198     }
 199
 200     if(fEMCALgeom){
 201       delete fEMCALgeom ;
 202       fEMCALgeom=0x0;
 203     }
 204
 205     for(Int_t ivtx=0; ivtx<10; ivtx++){
 206       if(fPHOSEvents[ivtx]){
 207         delete fPHOSEvents[ivtx] ;
 208           fPHOSEvents[ivtx]=0x0 ;
 209     }
 210       if(fEMCALEvents[ivtx]){
 211         delete fEMCALEvents[ivtx] ;
 212         fEMCALEvents[ivtx]=0x0 ;
 213       }
 214       if(fConvEvents[ivtx]){
 215         delete fConvEvents[ivtx] ;
 216         fConvEvents[ivtx]=0x0 ;
 217       }
 218     }
 219     for(Int_t i=0; i<6; i++)
 220       if(fPHOSBadMap[i]){
 221         delete fPHOSBadMap[i] ;
 222           fPHOSBadMap[i]=0 ;
 223       }
 224     for(Int_t i=0; i<10; i++)
 225       if(fEMCALBadMap[i]){
 226        delete fEMCALBadMap[i];
 227        fEMCALBadMap[i]=0 ;
 228     }
 229   }
 230 }
 231 //_____________________________________________________
 232 void AliAnalysisTaskCaloConv::Init()
 233 {
 234   // Initialization
 235   // AliLog::SetGlobalLogLevel(AliLog::kError);
 236 }
 237 //_____________________________________________________
 238 void AliAnalysisTaskCaloConv::UserExec(Option_t */*option*/)
 239 {
 240   // Execute analysis for current event
 241   // First select conversion and calorimeter photons
 242   // then construct inv. mass distributions
 243   //First try to find Stack information.
 244   if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()){
 245     if(static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent())
 246       fStack = static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent()->Stack();
 247   }
 248
 249
 250   AliESDInputHandler *esdHandler=dynamic_cast<AliESDInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 251   if(!fESDpid){
 252     if( esdHandler && esdHandler->GetESDpid()){
 253       fESDpid=new AliESDpid(*(esdHandler->GetESDpid())) ;
 254     }
 255     else {
 256       fESDpid=new AliESDpid;
 257       Double_t alephParameters[5];
 258       if(fStack){// simulation
 259         alephParameters[0] = 2.15898e+00/50.;
 260         alephParameters[1] = 1.75295e+01;
 261         alephParameters[2] = 3.40030e-09;
 262         alephParameters[3] = 1.96178e+00;
 263         alephParameters[4] = 3.91720e+00;
 264         fESDpid->GetTOFResponse().SetTimeResolution(80.);
 265       }
 266       else{// data
 267         alephParameters[0] = 0.0283086;
 268         alephParameters[1] = 2.63394e+01;
 269         alephParameters[2] = 5.04114e-11;
 270         alephParameters[3] = 2.12543e+00;
 271         alephParameters[4] = 4.88663e+00;
 272         fESDpid->GetTOFResponse().SetTimeResolution(130.);
 273         fESDpid->GetTPCResponse().SetMip(47.9);
 274       }
 275
 276       fESDpid->GetTPCResponse().SetBetheBlochParameters(
 277         alephParameters[0],alephParameters[1],alephParameters[2],
 278         alephParameters[3],alephParameters[4]);
 279       fESDpid->GetTPCResponse().SetSigma(3.79301e-03, 2.21280e+04);
 280     }
 281   }
 282
 283   if(!fESDtrackCuts){
 284 //    fESDtrackCuts= AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
 285     fESDtrackCuts = new AliESDtrackCuts;
 286
 287     // TPC
 288     fESDtrackCuts->SetMinNClustersTPC(70);
 289     fESDtrackCuts->SetMaxChi2PerClusterTPC(4);
 290     fESDtrackCuts->SetAcceptKinkDaughters(kFALSE);
 291     fESDtrackCuts->SetRequireTPCRefit(kTRUE);
 292     // ITS
 293     fESDtrackCuts->SetRequireITSRefit(kTRUE);
 294     fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 295                                            AliESDtrackCuts::kAny);
 296 //    if(selPrimaries) {
 297       // 7*(0.0026+0.0050/pt^1.01)
 298       fESDtrackCuts->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
 299 //    }
 300     fESDtrackCuts->SetMaxDCAToVertexZ(2);
 301     fESDtrackCuts->SetDCAToVertex2D(kFALSE);
 302     fESDtrackCuts->SetRequireSigmaToVertex(kFALSE);
 303
 304
 305
 306 /*
 307     fESDtrackCuts = new AliESDtrackCuts("AliESDtrackCuts","AliESDtrackCuts");
 308
 309     fESDtrackCuts->SetAcceptKinkDaughters(kFALSE);
 310     fESDtrackCuts->SetMinNClustersTPC(70);
 311     fESDtrackCuts->SetMaxChi2PerClusterTPC(4);
 312     fESDtrackCuts->SetRequireTPCRefit(kTRUE);
 313     fESDtrackCuts->SetRequireITSRefit(kTRUE);
 314 //    fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny); //TEMPORARY <-> REMOVE
 315 */
 316   }
 317
 318
 319   fESDEvent=(AliESDEvent*)InputEvent();
 320   FillHistogram("hEventsTrig",0.5) ;
 321   Bool_t isSelected = esdHandler && ((esdHandler->IsEventSelected()& AliVEvent::kMB) == AliVEvent::kMB);
 322   if(!isSelected){
 323     printf("Not selected !!!!! \n") ;
 324     PostData(1, fOutputContainer);
 325     return ;
 326   }
 327   FillHistogram("hEventsTrig",1.5) ;
 328
 329   //Take Only events with proper trigger
 330   //No trigger in MC data => no check
 331 //  if(!fStack && !fESDEvent->IsTriggerClassFired("CINT1B-ABCE-NOPF-ALL")){
 332   //for LHC10e
 333   if(!fStack && !fESDEvent->IsTriggerClassFired("CINT1-B-NOPF-ALLNOTRD")){
 334     PostData(1, fOutputContainer);
 335     return ;
 336   }
 337   FillHistogram("hEventsTrig",2.5) ;
 338
 339   //checks if we have a prim vertex
 340   if(fESDEvent->GetPrimaryVertex()->GetNContributors()<=0) {
 341     PostData(1, fOutputContainer);
 342     return ;
 343   }
 344   FillHistogram("hEventsTrig",3.5) ;
 345
 346   if(TMath::Abs(fESDEvent->GetPrimaryVertex()->GetZ())>10.){
 347     PostData(1, fOutputContainer);
 348     return ;
 349   }
 350   FillHistogram("hEventsTrig",4.5) ;
 351
 352
 353   //Calculate charged multiplicity
 354   Int_t trackCounter = 0;
 355   for (Int_t i=0;i<fESDEvent->GetNumberOfTracks();++i) {
 356     AliESDtrack *track = new AliESDtrack(*fESDEvent->GetTrack(i)) ;
 357     if(fESDtrackCuts->AcceptTrack(track) &&  TMath::Abs(track->Eta())< 0.9)
 358       trackCounter++;
 359     delete track;
 360   }
 361   fCentr=trackCounter+0.5 ;
 362   FillHistogram("hMult",fCentr) ;
 363
 364   //Init geometry if not done yet
 365   InitGeometry();
 366
 367   //Select conversion and calorimeter photons
 368   //Conversion photons should go first since they are used in calibration in SelectCALOPhotons()
 369   SelectConvPhotons() ;
 370   SelectPHOSPhotons() ;
 371   SelectEMCALPhotons() ;
 372   //Fill MC histograms if MC is present
 373   ProcessMC();
 374   FillRealMixed() ;
 375
 376   PostData(1, fOutputContainer);
 377   if(fToUseCF)
 378     PostData(2, fCFOutputContainer);  // for CF
 379
 380
 381 }
 382 //____________________________________________________________
 383 void AliAnalysisTaskCaloConv::ConnectInputData(Option_t *option){
 384   // see header file for documentation
 385
 386   AliAnalysisTaskSE::ConnectInputData(option);
 387
 388 }
 389 //____________________________________________________________
 390 void AliAnalysisTaskCaloConv::UserCreateOutputObjects()
 391 {
 392   //UserCreateOutputObjects
 393   if(fDebug)gDirectory->Print() ;
 394   // Create the output container
 395   if(fOutputContainer != NULL){
 396     delete fOutputContainer;
 397     fOutputContainer = NULL;
 398   }
 399   fOutputContainer = new TList();
 400   fOutputContainer->SetOwner(kTRUE);
 401
 402   if(fCFOutputContainer != NULL){
 403     delete fCFOutputContainer;
 404     fCFOutputContainer = NULL;
 405   }
 406   //===========Correction Framework ======================
 407   if(fToUseCF){
 408     fCFOutputContainer = new TList();
 409     fCFOutputContainer->SetOwner(kTRUE);
 410
 411     //bins: pt,eta,mass
 412     Int_t iBin[3]={500,40,100};
 413     fConvCFCont = new AliCFContainer("ConvContainer","container for converted photons", 23,3,iBin);
 414     fConvCFCont->SetBinLimits(0,0.,50.);
 415     fConvCFCont->SetBinLimits(1,-2.,2.) ;
 416     fConvCFCont->SetBinLimits(2,0.,1.);
 417     fCFOutputContainer->Add(fConvCFCont) ;
 418
 419     fPHOSCFCont = new AliCFContainer("PHOSContainer","container for PHOS photons", 10,2,iBin);
 420     fPHOSCFCont->SetBinLimits(0,0.,50.);
 421     fPHOSCFCont->SetBinLimits(1,-2.,2.) ;
 422     fCFOutputContainer->Add(fPHOSCFCont) ;
 423
 424     fEMCALCFCont = new AliCFContainer("EMCALContainer","container for EMCAL photons", 10,2,iBin);
 425     fEMCALCFCont->SetBinLimits(0,0.,50.);
 426     fEMCALCFCont->SetBinLimits(1,-2.,2.) ;
 427     fCFOutputContainer->Add(fEMCALCFCont) ;
 428
 429     fPi0CFCont = new AliCFContainer("Pi0Container","container for EMCAL photons", 10,2,iBin);
 430     fPi0CFCont->SetBinLimits(0,0.,50.);
 431     fPi0CFCont->SetBinLimits(1,-2.,2.) ;
 432     fCFOutputContainer->Add(fPi0CFCont) ;
 433
 434   }
 435   //========================================================
 436
 437   //Adding the histograms to the output container
 438   Int_t firstRun= 125000 ;
 439   Int_t lastRun = 135000 ;
 440   Int_t nRuns =lastRun-firstRun+1 ;
 441
 442   //Run QA histigrams
 443   fOutputContainer->Add(new TH2F("hRunTrigger","Triggers fired",nRuns,float(firstRun),float(lastRun),2,0.,2.)) ;
 444   fOutputContainer->Add(new TH1F("hRunEvents","Events per run",nRuns,float(firstRun),float(lastRun))) ;
 445   fOutputContainer->Add(new TH1F("hRunConvs","Conversion photons per run",nRuns,float(firstRun),float(lastRun))) ;
 446   fOutputContainer->Add(new TH1F("hRunPHOS","PHOS photons per run",nRuns,float(firstRun),float(lastRun))) ;
 447   fOutputContainer->Add(new TH1F("hRunEMCAL","EMCAL photons per run",nRuns,float(firstRun),float(lastRun))) ;
 448   fOutputContainer->Add(new TH1F("hVtxBin","Vtx distribution",10,0.,10.)) ;
 449   fOutputContainer->Add(new TH1F("hEvents","Events processed",1,0.,1.)) ;
 450   fOutputContainer->Add(new TH1F("hEventsTrig","Events processed",10,0.,10.)) ;
 451
 452   fOutputContainer->Add(new TH2F("hQA_PHOS_mod1_soft","number of clusters per cell",64,0.,64.,56,0.,56.)) ;
 453   fOutputContainer->Add(new TH2F("hQA_PHOS_mod2_soft","number of clusters per cell",64,0.,64.,56,0.,56.)) ;
 454   fOutputContainer->Add(new TH2F("hQA_PHOS_mod3_soft","number of clusters per cell",64,0.,64.,56,0.,56.)) ;
 455   fOutputContainer->Add(new TH2F("hQA_PHOS_mod4_soft","number of clusters per cell",64,0.,64.,56,0.,56.)) ;
 456   fOutputContainer->Add(new TH2F("hQA_PHOS_mod5_soft","number of clusters per cell",64,0.,64.,56,0.,56.)) ;
 457   fOutputContainer->Add(new TH2F("hQA_PHOS_mod1_hard","number of clusters per cell",64,0.,64.,56,0.,56.)) ;
 458   fOutputContainer->Add(new TH2F("hQA_PHOS_mod2_hard","number of clusters per cell",64,0.,64.,56,0.,56.)) ;
 459   fOutputContainer->Add(new TH2F("hQA_PHOS_mod3_hard","number of clusters per cell",64,0.,64.,56,0.,56.)) ;
 460   fOutputContainer->Add(new TH2F("hQA_PHOS_mod4_hard","number of clusters per cell",64,0.,64.,56,0.,56.)) ;
 461   fOutputContainer->Add(new TH2F("hQA_PHOS_mod5_hard","number of clusters per cell",64,0.,64.,56,0.,56.)) ;
 462
 463   fOutputContainer->Add(new TH2F("hQA_EMCAL_SM0_soft","number of clusters per cell",24,0.,24,48,0.,48.)) ;
 464   fOutputContainer->Add(new TH2F("hQA_EMCAL_SM1_soft","number of clusters per cell",24,0.,24,48,0.,48.)) ;
 465   fOutputContainer->Add(new TH2F("hQA_EMCAL_SM2_soft","number of clusters per cell",24,0.,24,48,0.,48.)) ;
 466   fOutputContainer->Add(new TH2F("hQA_EMCAL_SM3_soft","number of clusters per cell",24,0.,24,48,0.,48.)) ;
 467   fOutputContainer->Add(new TH2F("hQA_EMCAL_SM4_soft","number of clusters per cell",24,0.,24,48,0.,48.)) ;
 468   fOutputContainer->Add(new TH2F("hQA_EMCAL_SM5_soft","number of clusters per cell",24,0.,24,48,0.,48.)) ;
 469   fOutputContainer->Add(new TH2F("hQA_EMCAL_SM6_soft","number of clusters per cell",24,0.,24,48,0.,48.)) ;
 470   fOutputContainer->Add(new TH2F("hQA_EMCAL_SM7_soft","number of clusters per cell",24,0.,24,48,0.,48.)) ;
 471   fOutputContainer->Add(new TH2F("hQA_EMCAL_SM8_soft","number of clusters per cell",24,0.,24,48,0.,48.)) ;
 472   fOutputContainer->Add(new TH2F("hQA_EMCAL_SM9_soft","number of clusters per cell",24,0.,24,48,0.,48.)) ;
 473   fOutputContainer->Add(new TH2F("hQA_EMCAL_SM0_hard","number of clusters per cell",24,0.,24,48,0.,48.)) ;
 474   fOutputContainer->Add(new TH2F("hQA_EMCAL_SM1_hard","number of clusters per cell",24,0.,24,48,0.,48.)) ;
 475   fOutputContainer->Add(new TH2F("hQA_EMCAL_SM2_hard","number of clusters per cell",24,0.,24,48,0.,48.)) ;
 476   fOutputContainer->Add(new TH2F("hQA_EMCAL_SM3_hard","number of clusters per cell",24,0.,24,48,0.,48.)) ;
 477   fOutputContainer->Add(new TH2F("hQA_EMCAL_SM4_hard","number of clusters per cell",24,0.,24,48,0.,48.)) ;
 478   fOutputContainer->Add(new TH2F("hQA_EMCAL_SM5_hard","number of clusters per cell",24,0.,24,48,0.,48.)) ;
 479   fOutputContainer->Add(new TH2F("hQA_EMCAL_SM6_hard","number of clusters per cell",24,0.,24,48,0.,48.)) ;
 480   fOutputContainer->Add(new TH2F("hQA_EMCAL_SM7_hard","number of clusters per cell",24,0.,24,48,0.,48.)) ;
 481   fOutputContainer->Add(new TH2F("hQA_EMCAL_SM8_hard","number of clusters per cell",24,0.,24,48,0.,48.)) ;
 482   fOutputContainer->Add(new TH2F("hQA_EMCAL_SM9_hard","number of clusters per cell",24,0.,24,48,0.,48.)) ;
 483
 484   fOutputContainer->Add(new TH2F("hQA_ConvPhiEta","Number of V0s phi eta",100,0.,TMath::TwoPi(),40,-1.5,1.5)) ;
 485
 486   fOutputContainer->Add(new TH2F("hdEdx","dEdx of acceptaed electrons",1000,0.,10.,150,0.,150.)) ;
 487
 488   fOutputContainer->Add(new TH1F("hMult","Multiplicity",200,0.,200.)) ;
 489
 490   Int_t npt=200 ;
 491   Double_t ptmax=20. ;
 492   //Calibration of PHOS
 493   fOutputContainer->Add(new TH3F("PHOS_mod1_th1","Inv.Mass distr. per channel",64,0.,64,56,0.,56,100,0.,0.5)) ;
 494   fOutputContainer->Add(new TH3F("PHOS_mod2_th1","Inv.Mass distr. per channel",64,0.,64,56,0.,56,100,0.,0.5)) ;
 495   fOutputContainer->Add(new TH3F("PHOS_mod3_th1","Inv.Mass distr. per channel",64,0.,64,56,0.,56,100,0.,0.5)) ;
 496   fOutputContainer->Add(new TH3F("PHOS_mod4_th1","Inv.Mass distr. per channel",64,0.,64,56,0.,56,100,0.,0.5)) ;
 497   fOutputContainer->Add(new TH3F("PHOS_mod5_th1","Inv.Mass distr. per channel",64,0.,64,56,0.,56,100,0.,0.5)) ;
 498
 499   fOutputContainer->Add(new TH3F("PHOS_mod1_th2","Inv.Mass distr. per channel",64,0.,64,56,0.,56,100,0.,0.5)) ;
 500   fOutputContainer->Add(new TH3F("PHOS_mod2_th2","Inv.Mass distr. per channel",64,0.,64,56,0.,56,100,0.,0.5)) ;
 501   fOutputContainer->Add(new TH3F("PHOS_mod3_th2","Inv.Mass distr. per channel",64,0.,64,56,0.,56,100,0.,0.5)) ;
 502   fOutputContainer->Add(new TH3F("PHOS_mod4_th2","Inv.Mass distr. per channel",64,0.,64,56,0.,56,100,0.,0.5)) ;
 503   fOutputContainer->Add(new TH3F("PHOS_mod5_th2","Inv.Mass distr. per channel",64,0.,64,56,0.,56,100,0.,0.5)) ;
 504
 505   //Pi0 histograms
 506   //Vary Conversion cuts
 507   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_OnFly","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 508   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_Offline","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 509   fOutputContainer->Add(new TH3F("PHOS_Re_mvsPt_OnFly_mult","Mass vs pt",400,0.,1.,npt,0.,ptmax,150,0.,150.)) ;
 510   fOutputContainer->Add(new TH3F("PHOS_Re_mvsPt_Offline_mult","Mass vs pt",400,0.,1.,npt,0.,ptmax,150,0.,150.)) ;
 511   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_On_dEdx","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 512   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_Off_dEdx","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 513 //  fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_On_Prob","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 514 //  fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_Off_Prob","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 515   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_On_R120","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 516   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_Off_R120","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 517   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_On_Z","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 518   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_Off_Z","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 519   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_On_chi","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 520   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_Off_chi","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 521   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_On_Eta","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 522   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_Off_Eta","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 523   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_On_Wcut","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 524   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_On_Wcut_Neu","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 525   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_Off_Wcut","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 526   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_Off_Wcut_Neu","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 527   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_On_ArmQt","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 528   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_Off_ArmQt","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 529
 530   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_OnFly","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 531   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_Offline","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 532   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_On_dEdx","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 533   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_Off_dEdx","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 534 //  fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_On_Prob","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 535 //  fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_Off_Prob","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 536   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_On_R120","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 537   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_Off_R120","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 538   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_On_Z","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 539   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_Off_Z","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 540   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_On_chi","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 541   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_Off_chi","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 542   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_On_Eta","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 543   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_Off_Eta","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 544   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_On_Wcut","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 545   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_Off_Wcut","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 546   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_On_ArmQt","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 547   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_Off_ArmQt","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 548
 549   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_OnFly","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 550   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_Offline","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 551   fOutputContainer->Add(new TH3F("EMCAL_Re_mvsPt_OnFly_mult","Mass vs pt",400,0.,1.,npt,0.,ptmax,30,0.,60.)) ;
 552   fOutputContainer->Add(new TH3F("EMCAL_Re_mvsPt_Offline_mult","Mass vs pt",400,0.,1.,npt,0.,ptmax,30,0.,60.)) ;
 553   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_On_ArmQt","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 554   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_Off_ArmQt","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 555   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_On_dEdx","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 556   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_Off_dEdx","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 557 //  fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_On_Prob","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 558 //  fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_Off_Prob","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 559   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_On_R120","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 560   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_Off_R120","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 561   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_On_Z","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 562   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_Off_Z","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 563   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_On_chi","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 564   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_Off_chi","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 565   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_On_Eta","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 566   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_Off_Eta","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 567   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_On_Wcut","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 568   fOutputContainer->Add(new TH2F("EMCAL_Re_mvsPt_Off_Wcut","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 569
 570   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_OnFly","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 571   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_Offline","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 572   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_On_ArmQt","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 573   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_Off_ArmQt","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 574   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_On_dEdx","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 575   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_Off_dEdx","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 576 //  fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_On_Prob","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 577 //  fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_Off_Prob","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 578   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_On_R120","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 579   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_Off_R120","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 580   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_On_Z","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 581   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_Off_Z","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 582   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_On_chi","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 583   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_Off_chi","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 584   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_On_Eta","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 585   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_Off_Eta","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 586   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_On_Wcut","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 587   fOutputContainer->Add(new TH2F("EMCAL_Mi_mvsPt_Off_Wcut","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 588
 589   //PHOS PID variations
 590   fOutputContainer->Add(new TH3F("PHOS_Re_mvsPt_alpha","Mass vs pt vs PHOS E",400,0.,1.,npt,0.,ptmax,10,0.,1.)) ;
 591   fOutputContainer->Add(new TH3F("PHOS_Re_mvsPt_E","Mass vs pt vs PHOS E",400,0.,1.,npt,0.,ptmax,100,0.,10.)) ;
 592   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_all","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 593   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_all_dist","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 594   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_Disp","Mass vs pt, disp cut",400,0.,1.,npt,0.,ptmax)) ;
 595   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_TOF","Mass vs pt, TOF cut",400,0.,1.,npt,0.,ptmax)) ;
 596   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_Neutral","Mass vs pt, Neutral cut",400,0.,1.,npt,0.,ptmax)) ;
 597   fOutputContainer->Add(new TH2F("PHOS_Re_mvsPt_DispNeutral","Mass vs pt, Disp & neutral cut",400,0.,1.,npt,0.,ptmax)) ;
 598
 599   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_all","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 600   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_all_dist","Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 601   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_Disp","Mass vs pt, disp cut",400,0.,1.,npt,0.,ptmax)) ;
 602   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_TOF","Mass vs pt, TOF cut",400,0.,1.,npt,0.,ptmax)) ;
 603   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_Neutral","Mass vs pt, Neutral cut",400,0.,1.,npt,0.,ptmax)) ;
 604   fOutputContainer->Add(new TH2F("PHOS_Mi_mvsPt_DispNeutral","Mass vs pt, Disp & neutral cut",400,0.,1.,npt,0.,ptmax)) ;
 605
 606   char key[155] ;
 607   for(Int_t mod=1; mod<=5;mod++){
 608     snprintf(key,155,"PHOS_Re_mvsPt_mod%d_single",mod) ;
 609     fOutputContainer->Add(new TH2F(key,"Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 610     snprintf(key,155,"PHOS_Re_mvsPt_mod%d_all",mod) ;
 611     fOutputContainer->Add(new TH2F(key,"Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 612   }
 613
 614   //Single photon spectrum
 615   //Conversion
 616   fOutputContainer->Add(new TH1F("Single_conv_OnFly","Single photon spectrum",npt,0.,ptmax)) ;
 617   fOutputContainer->Add(new TH1F("Single_conv_Offline","Single photon spectrum",npt,0.,ptmax)) ;
 618   fOutputContainer->Add(new TH1F("Single_conv_On_ArmQt","Single photon spectrum",npt,0.,ptmax)) ;
 619   fOutputContainer->Add(new TH1F("Single_conv_Off_ArmQt","Single photon spectrum",npt,0.,ptmax)) ;
 620   fOutputContainer->Add(new TH1F("Single_conv_On_dEdx","Single photon spectrum",npt,0.,ptmax)) ;
 621   fOutputContainer->Add(new TH1F("Single_conv_Off_dEdx","Single photon spectrum",npt,0.,ptmax)) ;
 622 //  fOutputContainer->Add(new TH1F("Single_conv_On_Prob","Single photon spectrum",npt,0.,ptmax)) ;
 623 //  fOutputContainer->Add(new TH1F("Single_conv_Off_Prob","Single photon spectrum",npt,0.,ptmax)) ;
 624   fOutputContainer->Add(new TH1F("Single_conv_On_R120","Single photon spectrum",npt,0.,ptmax)) ;
 625   fOutputContainer->Add(new TH1F("Single_conv_Off_R120","Single photon spectrum",npt,0.,ptmax)) ;
 626   fOutputContainer->Add(new TH1F("Single_conv_On_Z","Single photon spectrum",npt,0.,ptmax)) ;
 627   fOutputContainer->Add(new TH1F("Single_conv_Off_Z","Single photon spectrum",npt,0.,ptmax)) ;
 628   fOutputContainer->Add(new TH1F("Single_conv_On_chi","Single photon spectrum",npt,0.,ptmax)) ;
 629   fOutputContainer->Add(new TH1F("Single_conv_Off_chi","Single photon spectrum",npt,0.,ptmax)) ;
 630   fOutputContainer->Add(new TH1F("Single_conv_On_Eta","Single photon spectrum",npt,0.,ptmax)) ;
 631   fOutputContainer->Add(new TH1F("Single_conv_Off_Eta","Single photon spectrum",npt,0.,ptmax)) ;
 632   fOutputContainer->Add(new TH1F("Single_conv_On_Wcut","Single photon spectrum",npt,0.,ptmax)) ;
 633   fOutputContainer->Add(new TH1F("Single_conv_Off_Wcut","Single photon spectrum",npt,0.,ptmax)) ;
 634
 635   //PHOS
 636   fOutputContainer->Add(new TH2F("PHOS_single_all_mult","Single photon spectrum",npt,0.,ptmax,150,0.,150.)) ;
 637   fOutputContainer->Add(new TH2F("PHOS_single_disp_mult","Single photon spectrum",npt,0.,ptmax,150,0.,150.)) ;
 638   fOutputContainer->Add(new TH2F("PHOS_single_neu_mult","Single photon spectrum",npt,0.,ptmax,150,0.,150.)) ;
 639
 640   for(Int_t mod=1; mod<=5;mod++){
 641     snprintf(key,155,"PHOS_single_mod%d_all",mod) ;
 642     fOutputContainer->Add(new TH1F(key,"Single photon spectrum",npt,0.,ptmax)) ;
 643     snprintf(key,155,"PHOS_single_mod%d_disp",mod) ;
 644     fOutputContainer->Add(new TH1F(key,"Single photon spectrum",npt,0.,ptmax)) ;
 645     snprintf(key,155,"PHOS_single_mod%d_neutral",mod) ;
 646     fOutputContainer->Add(new TH1F(key,"Single photon spectrum",npt,0.,ptmax)) ;
 647     snprintf(key,155,"PHOS_single_mod%d_dispneutral",mod) ;
 648     fOutputContainer->Add(new TH1F(key,"Single photon spectrum",npt,0.,ptmax)) ;
 649     snprintf(key,155,"PHOS_single_mod%d_dist1",mod) ;
 650     fOutputContainer->Add(new TH1F(key,"Single photon spectrum",npt,0.,ptmax)) ;
 651     snprintf(key,155,"PHOS_single_mod%d_dist2",mod) ;
 652     fOutputContainer->Add(new TH1F(key,"Single photon spectrum",npt,0.,ptmax)) ;
 653   }
 654
 655   for(Int_t mod=0; mod<4;mod++){
 656     snprintf(key,155,"EMCAL_mod%d_th1",mod) ;
 657     fOutputContainer->Add(new TH3F(key,"Inv.Mass distr. per channel",24,0.,24,48,0.,48,100,0.,0.5)) ;
 658     snprintf(key,155,"EMCAL_mod%d_th2",mod) ;
 659     fOutputContainer->Add(new TH3F(key,"Inv.Mass distr. per channel",24,0.,24,48,0.,48,100,0.,0.5)) ;
 660
 661     snprintf(key,155,"EMCAL_Re_mvsPt_mod%d_single",mod) ;
 662     fOutputContainer->Add(new TH2F(key,"Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 663
 664     snprintf(key,155,"EMCAL_Re_mvsPt_mod%d_all",mod) ;
 665     fOutputContainer->Add(new TH2F(key,"Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 666     snprintf(key,155,"EMCAL_Re_mvsPt_mod%d_Disp",mod) ;
 667     fOutputContainer->Add(new TH2F(key,"Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 668     snprintf(key,155,"EMCAL_Re_mvsPt_mod%d_TOF",mod) ;
 669     fOutputContainer->Add(new TH2F(key,"Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 670     snprintf(key,155,"EMCAL_Re_mvsPt_mod%d_Neutral",mod) ;
 671     fOutputContainer->Add(new TH2F(key,"Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 672     snprintf(key,155,"EMCAL_Re_mvsPt_mod%d_DispNeutral",mod) ;
 673     fOutputContainer->Add(new TH2F(key,"Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 674
 675     snprintf(key,155,"EMCAL_Mi_mvsPt_mod%d_all",mod) ;
 676     fOutputContainer->Add(new TH2F(key,"Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 677     snprintf(key,155,"EMCAL_Mi_mvsPt_mod%d_Disp",mod) ;
 678     fOutputContainer->Add(new TH2F(key,"Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 679     snprintf(key,155,"EMCAL_Mi_mvsPt_mod%d_TOF",mod) ;
 680     fOutputContainer->Add(new TH2F(key,"Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 681     snprintf(key,155,"EMCAL_Mi_mvsPt_mod%d_Neutral",mod) ;
 682     fOutputContainer->Add(new TH2F(key,"Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 683     snprintf(key,155,"EMCAL_Mi_mvsPt_mod%d_DispNeutral",mod) ;
 684     fOutputContainer->Add(new TH2F(key,"Mass vs pt",400,0.,1.,npt,0.,ptmax)) ;
 685   }
 686
 687   //MC info
 688   fOutputContainer->Add(new TH2F("hMC_CaloConv_pi0_unitEta","Primary #pi^{0}",npt,0.,ptmax,150,0.,150.)) ;
 689   fOutputContainer->Add(new TH2F("hMC_CaloConv_eta_unitEta","Primary #pi^{0}",npt,0.,ptmax,150,0.,150.)) ;
 690   fOutputContainer->Add(new TH1F("hMC_CaloConv_allpi0","Primary #pi^{0}",npt,0.,ptmax)) ;
 691   fOutputContainer->Add(new TH1F("hMC_CaloConv_alleta","Primary #pi^{0}",npt,0.,ptmax)) ;
 692   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0PHOSacc","#pi^{0} decayed in PHOS acc",npt,0.,ptmax)) ;
 693   fOutputContainer->Add(new TH1F("hMC_CaloConv_etaPHOSacc","#pi^{0} decayed in PHOS acc",npt,0.,ptmax)) ;
 694   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0EMCALacc","#pi^{0} decayed in EMCAL acc",npt,0.,ptmax)) ;
 695   fOutputContainer->Add(new TH1F("hMC_CaloConv_etaEMCALacc","#pi^{0} decayed in EMCAL acc",npt,0.,ptmax)) ;
 696   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_PHOS_conv","#pi^{0} decayed in PHOS acc asnd conv. photon",npt,0.,ptmax)) ;
 697   fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_PHOS_conv","#pi^{0} decayed in PHOS acc asnd conv. photon",npt,0.,ptmax)) ;
 698   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_EMCAL_conv","#pi^{0} decayed in EMCAL acc asnd conv. photon",npt,0.,ptmax)) ;
 699   fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_EMCAL_conv","#pi^{0} decayed in EMCAL acc asnd conv. photon",npt,0.,ptmax)) ;
 700   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_bothphot_conv","#pi^{0} both photons converted",npt,0.,ptmax)) ;
 701   fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_bothphot_conv","#pi^{0} both photons converted",npt,0.,ptmax)) ;
 702   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_convPhotInCalo","#pi^{0} photon in calo converted",npt,0.,ptmax)) ;
 703   fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_convPhotInCalo","#pi^{0} photon in calo converted",npt,0.,ptmax)) ;
 704   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0onfly_PHOSacc","#pi^{0} photon converted and V0 found",npt,0.,ptmax)) ;
 705   fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0onfly_PHOSacc","#pi^{0} photon converted and V0 found",npt,0.,ptmax)) ;
 706   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0offline_PHOSacc","#pi^{0} photon converted and V0 found",npt,0.,ptmax)) ;
 707   fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0offline_PHOSacc","#pi^{0} photon converted and V0 found",npt,0.,ptmax)) ;
 708   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0onfly_EMCALacc","#pi^{0} photon converted and V0 found",npt,0.,ptmax)) ;
 709   fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0onfly_EMCALacc","#pi^{0} photon converted and V0 found",npt,0.,ptmax)) ;
 710   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0offline_EMCALacc","#pi^{0} photon converted and V0 found",npt,0.,ptmax)) ;
 711   fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0offline_EMCALacc","#pi^{0} photon converted and V0 found",npt,0.,ptmax)) ;
 712   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0onfly_PHOSclu","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
 713   fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0onfly_PHOSclu","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
 714   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0onfly_PHOSclu_pid","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
 715   fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0onfly_PHOSclu_pid","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
 716   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0offline_PHOSclu","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
 717   fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0offline_PHOSclu","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
 718   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0offline_PHOSclu_pid","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
 719   fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0offline_PHOSclu_pid","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
 720   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0onfly_PHOSclu_good","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
 721   fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0onfly_PHOSclu_good","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
 722   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0offline_PHOSclu_good","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
 723   fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0offline_PHOSclu_good","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
 724   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0onfly_PHOSclu_mod1","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
 725   fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0onfly_PHOSclu_mod1","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
 726   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0onfly_PHOSclu_mod2","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
 727   fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0onfly_PHOSclu_mod2","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
 728   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0onfly_PHOSclu_mod3","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
 729   fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0onfly_PHOSclu_mod3","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
 730   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0onfly_PHOSclu_mod4","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
 731   fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0onfly_PHOSclu_mod4","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
 732   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0onfly_PHOSclu_mod5","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
 733   fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0onfly_PHOSclu_mod5","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
 734   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0offline_PHOSclu_mod1","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
 735   fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0offline_PHOSclu_mod1","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
 736   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0offline_PHOSclu_mod2","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
 737   fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0offline_PHOSclu_mod2","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
 738   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0offline_PHOSclu_mod3","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
 739   fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0offline_PHOSclu_mod3","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
 740   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0offline_PHOSclu_mod4","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
 741   fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0offline_PHOSclu_mod4","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
 742   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0offline_PHOSclu_mod5","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
 743   fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0offline_PHOSclu_mod5","#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
 744
 745   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0on_PHOSclu_ptRec","#pi^{0} V0 and cluster in PHOS found(rec pt)",npt,0.,ptmax)) ;
 746   fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0on_PHOSclu_ptRec","#pi^{0} V0 and cluster in PHOS found(rec pt)",npt,0.,ptmax)) ;
 747   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0off_PHOSclu_ptRec","#pi^{0} V0 and cluster in PHOS found(rec pt)",npt,0.,ptmax)) ;
 748   fOutputContainer->Add(new TH1F("hMC_CaloConv_eta_v0off_PHOSclu_ptRec","#pi^{0} V0 and cluster in PHOS found(rec pt)",npt,0.,ptmax)) ;
 749   fOutputContainer->Add(new TH2F("hMC_CaloConv_pi0_v0on_PHOSclu_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax)) ;
 750   fOutputContainer->Add(new TH2F("hMC_CaloConv_eta_v0on_PHOSclu_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax)) ;
 751   fOutputContainer->Add(new TH2F("hMC_CaloConv_pi0_v0off_PHOSclu_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax)) ;
 752   fOutputContainer->Add(new TH2F("hMC_CaloConv_eta_v0off_PHOSclu_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax)) ;
 753
 754   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0off_EMCALclu_ptRec","#pi^{0} V0 and cluster in EMCAL found (rec pt)",npt,0.,ptmax)) ;
 755   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0on_EMCALclu_ptRec","#pi^{0} V0 and cluster in EMCAL found (rec pt)",npt,0.,ptmax)) ;
 756   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0onfly_EMCALclu","#pi^{0} V0 and cluster in EMCAL found",npt,0.,ptmax)) ;
 757   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0onfly_EMCALclu_pid","#pi^{0} V0 and cluster in EMCAL found",npt,0.,ptmax)) ;
 758   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0offline_EMCALclu","#pi^{0} V0 and cluster in EMCAL found",npt,0.,ptmax)) ;
 759   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0offline_EMCALclu_pid","#pi^{0} V0 and cluster in EMCAL found",npt,0.,ptmax)) ;
 760   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0onfly_EMCALclu_good","#pi^{0} V0 and cluster in EMCAL found",npt,0.,ptmax)) ;
 761   fOutputContainer->Add(new TH1F("hMC_CaloConv_pi0_v0offline_EMCALclu_good","#pi^{0} V0 and cluster in EMCAL found",npt,0.,ptmax)) ;
 762   fOutputContainer->Add(new TH2F("hMC_CaloConv_pi0_v0on_EMCALclu_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax)) ;
 763   fOutputContainer->Add(new TH2F("hMC_CaloConv_pi0_v0off_EMCALclu_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax)) ;
 764
 765   fOutputContainer->Add(new TH3F("hMC_Resid_PHOS_Phot_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax,10,0.,1.)) ;
 766   fOutputContainer->Add(new TH3F("hMC_Resid_PHOS_Pi0_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax,10,0.,1.)) ;
 767   fOutputContainer->Add(new TH3F("hMC_Resid_PHOS_eta_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax,10,0.,1.)) ;
 768   fOutputContainer->Add(new TH3F("hMC_Resid_PHOS_K_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax,10,0.,1.)) ;
 769   fOutputContainer->Add(new TH3F("hMC_Resid_PHOS_pi_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax,10,0.,1.)) ;
 770   fOutputContainer->Add(new TH3F("hMC_Resid_PHOS_pbar_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax,10,0.,1.)) ;
 771   fOutputContainer->Add(new TH3F("hMC_Resid_PHOS_other_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax,10,0.,1.)) ;
 772   fOutputContainer->Add(new TH3F("hMC_Resid_EMCAL_Phot_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax,10,0.,1.)) ;
 773   fOutputContainer->Add(new TH3F("hMC_Resid_EMCAL_Pi0_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax,10,0.,1.)) ;
 774   fOutputContainer->Add(new TH3F("hMC_Resid_EMCAL_eta_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax,10,0.,1.)) ;
 775   fOutputContainer->Add(new TH3F("hMC_Resid_EMCAL_K_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax,10,0.,1.)) ;
 776   fOutputContainer->Add(new TH3F("hMC_Resid_EMCAL_pi_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax,10,0.,1.)) ;
 777   fOutputContainer->Add(new TH3F("hMC_Resid_EMCAL_pbar_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax,10,0.,1.)) ;
 778   fOutputContainer->Add(new TH3F("hMC_Resid_EMCAL_other_mvsPt","m vs pt for rec pi0s",400,0.,1.,npt,0.,ptmax,10,0.,1.)) ;
 779
 780
 781   fOutputContainer->Add(new TH1F("hMC_CaloConv_phot","Primary photons",npt,0.,ptmax)) ;
 782   fOutputContainer->Add(new TH1F("hMC_CaloConv_gammaPHOSacc","Photons in PHOS acc",npt,0.,ptmax)) ;
 783   fOutputContainer->Add(new TH1F("hMC_CaloConv_gammaEMCALacc","Photons in EMCAL acc",npt,0.,ptmax)) ;
 784   fOutputContainer->Add(new TH1F("hMC_CaloConv_gamma_conv","Converted photons",npt,0.,ptmax)) ;
 785   fOutputContainer->Add(new TH1F("hMC_CaloConv_gamma_v0","Converted photons with V0",npt,0.,ptmax)) ;
 786   fOutputContainer->Add(new TH2F("hMC_CaloConv_gamma_v0_devsE","Converted photons with V0",200,-1.,1.,npt,0.,ptmax)) ;
 787   fOutputContainer->Add(new TH1F("hMC_CaloConv_gamma_PHOSclu","Photons with cluster in PHOS",npt,0.,ptmax)) ;
 788   fOutputContainer->Add(new TH1F("hMC_CaloConv_gamma_PHOSclu_dist1","Photons with cluster in PHOS",npt,0.,ptmax)) ;
 789   fOutputContainer->Add(new TH1F("hMC_CaloConv_gamma_PHOSclu_dist2","Photons with cluster in PHOS",npt,0.,ptmax)) ;
 790   fOutputContainer->Add(new TH1F("hMC_CaloConv_gamma_EMCALclu","Photons with cluster in EMCAL",npt,0.,ptmax)) ;
 791   fOutputContainer->Add(new TH1F("hMC_CaloConv_gamma_EMCALclu_dist1","Photons with cluster in EMCAL",npt,0.,ptmax)) ;
 792   fOutputContainer->Add(new TH1F("hMC_CaloConv_gamma_EMCALclu_dist2","Photons with cluster in EMCAL",npt,0.,ptmax)) ;
 793   fOutputContainer->Add(new TH1F("hMC_CaloConv_gamma_PHOSclu_recE","Photons with cluster in PHOS",npt,0.,ptmax)) ;
 794   fOutputContainer->Add(new TH1F("hMC_CaloConv_gamma_EMCALclu_recE","Photons with cluster in EMCAL",npt,0.,ptmax)) ;
 795   fOutputContainer->Add(new TH2F("hMC_CaloConv_gamma_PHOSclu_devsE","Photons with cluster in PHOS",200,-1.,1.,npt,0.,ptmax)) ;
 796   fOutputContainer->Add(new TH2F("hMC_CaloConv_gamma_EMCALclu_devsE","Photons with cluster in EMCAL",200,-1.,1.,npt,0.,ptmax)) ;
 797
 798
 799   //Non-linearity test
 800   char keym[55] ;
 801   for(Int_t iw=0;iw<10;iw++){ //resolution
 802     for(Int_t in=0;in<10;in++){
 803       snprintf(keym,55,"hMC_nonlinearity_w%d_n%d",iw,in) ;
 804       fOutputContainer->Add(new TH2F(keym,"m vs pt, nonlinearity test" ,200,0.,0.5,npt,0.,ptmax)) ;
 805       snprintf(keym,55,"hMC_nonlinearity_ConvPHOS_w%d_n%d",iw,in) ;
 806       fOutputContainer->Add(new TH2F(keym,"m vs pt, nonlinearity test" ,200,0.,0.5,npt,0.,ptmax)) ;
 807       snprintf(keym,55,"hMC_nonlinearity_EMCAL_w%d_n%d",iw,in) ;
 808       fOutputContainer->Add(new TH2F(keym,"m vs pt, nonlinearity test" ,200,0.,0.5,npt,0.,ptmax)) ;
 809       snprintf(keym,55,"hMC_CaloConv_pi0_v0onfly_PHOSclu_w%d_n%d",iw,in) ;
 810       fOutputContainer->Add(new TH1F(keym,"#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
 811       snprintf(keym,55,"hMC_CaloConv_pi0_v0onfly_ConvPHOSclu_w%d_n%d",iw,in) ;
 812       fOutputContainer->Add(new TH1F(keym,"#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
 813       snprintf(keym,55,"hMC_CaloConv_pi0_v0onfly_EMCALclu_w%d_n%d",iw,in) ;
 814       fOutputContainer->Add(new TH1F(keym,"#pi^{0} V0 and cluster in PHOS found",npt,0.,ptmax)) ;
 815     }
 816   }
 817
 818
 819   fOutputContainer->Add(new TH3F("All_chi2_eta_pt","MC chi2 vs eta vs phi",100,0.,100.,200,-2.,2.,npt,0.,ptmax)) ;
 820   fOutputContainer->Add(new TH2F("All_w_vs_m","MC w vs m",300,0.,TMath::Pi(),400,0.,1.)) ;
 821   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())) ;
 822   fOutputContainer->Add(new TH3F("MC_V0_m_eta_pt","MC m vs eta vs phi",400,0.,1.,200,-2.,2.,npt,0.,ptmax)) ;
 823   fOutputContainer->Add(new TH3F("MC_V0_chi2_eta_pt","MC chi2 vs eta vs phi",100,0.,100.,200,-2.,2.,npt,0.,ptmax)) ;
 824   fOutputContainer->Add(new TH2F("MC_V0_w_vs_m","MC w vs m",300,0.,TMath::Pi(),400,0.,1.)) ;
 825
 826   fOutputContainer->SetName(GetName());
 827 }
 828 //______________________________________________________________________
 829  void AliAnalysisTaskCaloConv::InitGeometry()
 830 {
 831   //If not done yet, create Geometry for PHOS and EMCAL
 832   //and read misalignment matrixes from ESD/AOD (AOD not implemented yet)
 833   //
 834
 835    if(fPHOSgeom && fEMCALgeom){ //already initialized
 836      return ;
 837    }
 838
 839    AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent()) ;
 840    if(!esd )
 841      AliFatal("Can not read geometry matrixes from ESD/AOD: NO ESD") ;
 842    if(!fPHOSgeom){//reading PHOS matrixes
 843      fPHOSgeom = new AliPHOSGeoUtils("IHEP","");
 844      for(Int_t mod=0; mod<5; mod++){
 845        if(esd){
 846          const TGeoHMatrix* m=esd->GetPHOSMatrix(mod) ;
 847          if(m)
 848            fPHOSgeom->SetMisalMatrix(m, mod) ;
 849        }
 850      }
 851    }
 852    if(!fEMCALgeom){
 853      fEMCALgeom = AliEMCALGeometry::GetInstance("EMCAL_FIRSTYEARV1");
 854      for(Int_t mod=0; mod < (fEMCALgeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){ //<---Gustavo, could you check???
 855        if(esd){
 856          const TGeoHMatrix* m=esd->GetEMCALMatrix(mod) ;
 857          if(m)
 858            fEMCALgeom->SetMisalMatrix(m, mod) ;
 859        }
 860      }
 861    }
 862 }
 863 //________________________________________________________________
 864 void AliAnalysisTaskCaloConv::SelectPHOSPhotons(){
 865   //SelectPHOSPhotons
 866   // Loop over all CaloClusters
 867   if(fPHOSEvent)
 868     fPHOSEvent->Clear() ;
 869   else
 870     fPHOSEvent = new TClonesArray("TLorentzVector",10) ;
 871   Int_t inPHOS = 0;
 872   TLorentzVector pi0 ;
 873
 874   //vertex
 875   Double_t vtx[3];
 876   vtx[0] = fESDEvent->GetPrimaryVertex()->GetX();
 877   vtx[1] = fESDEvent->GetPrimaryVertex()->GetY();
 878   vtx[2] = fESDEvent->GetPrimaryVertex()->GetZ();
 879   for (Int_t i=0; i<fESDEvent->GetNumberOfCaloClusters(); i++) {
 880     AliESDCaloCluster * clu = fESDEvent->GetCaloCluster(i);
 881     if(!clu->IsPHOS())
 882       continue ;
 883     TLorentzVector p ;
 884     clu ->GetMomentum(p ,vtx);
 885     if(p.Energy()<0.25){
 886       continue ;
 887     }
 888     if(clu->GetNCells()<=2){
 889       continue ;
 890     }
 891
 892     Bool_t isNeutral = kTRUE ;
 893     Bool_t isDispOK = kTRUE ;
 894     Bool_t isTOFOK = kTRUE ;
 895     Int_t iMod,iX,iZ ;
 896
 897     isNeutral = clu->GetEmcCpvDistance()>5. ;  //To be improved
 898     isDispOK = kFALSE ;
 899     Double_t l0=clu->GetM02(),l1=clu->GetM20() ;
 900     if(l1>= 0   && l0>= 0   && l1 < 0.1 && l0 < 0.1) isDispOK=kFALSE ;
 901     if(l1>= 0   && l0 > 0.5 && l1 < 0.1 && l0 < 1.5) isDispOK=kTRUE ;
 902     if(l1>= 0   && l0 > 2.0 && l1 < 0.1 && l0 < 2.7) isDispOK=kFALSE ;
 903     if(l1>= 0   && l0 > 2.7 && l1 < 0.1 && l0 < 4.0) isDispOK=kFALSE ;
 904     if(l1 > 0.1 && l1 < 0.7 && l0 > 0.7 && l0 < 2.1) isDispOK=kTRUE ;
 905     if(l1 > 0.1 && l1 < 0.3 && l0 > 3.0 && l0 < 5.0) isDispOK=kFALSE  ;
 906     if(l1 > 0.3 && l1 < 0.7 && l0 > 2.5 && l0 < 4.0) isDispOK=kFALSE ;
 907     if(l1 > 0.7 && l1 < 1.3 && l0 > 1.0 && l0 < 1.6) isDispOK=kTRUE ;
 908     if(l1 > 0.7 && l1 < 1.3 && l0 > 1.6 && l0 < 3.5) isDispOK=kTRUE ;
 909     if(l1 > 1.3 && l1 < 3.5 && l0 > 1.3 && l0 < 3.5) isDispOK=kTRUE ;
 910
 911     Float_t xyz[3] = {0,0,0};
 912     clu->GetPosition(xyz);   //Global position in ALICE system
 913     TVector3 global(xyz) ;
 914     Int_t relid[4] ;
 915     if(!fPHOSgeom->GlobalPos2RelId(global,relid)){
 916       printf("PHOS_beyond: x=%f, y=%f, z=%f \n",xyz[0],xyz[1],xyz[2]) ;
 917       continue ;
 918     }
 919     iMod=relid[0] ;
 920     iX=relid[2];
 921     iZ=relid[3] ;
 922     if(!IsGoodChannel("PHOS",iMod,iX,iZ))
 923       continue ;
 924
 925     Bool_t closeToBad=(clu->GetDistanceToBadChannel()>fBadDistCutPHOS) ;
 926
 927     p.SetBit(kCaloPIDdisp,isDispOK) ;
 928     p.SetBit(kCaloPIDtof,isTOFOK) ;
 929     p.SetBit(kCaloPIDneutral,isNeutral) ;
 930     p.SetBit(BIT(17+iMod),kTRUE) ;
 931     p.SetBit(kCaloDistBad,closeToBad) ;
 932     new((*fPHOSEvent)[inPHOS]) TLorentzVector(p) ;
 933     fGammaPHOS[inPHOS] = i ;
 934     inPHOS++ ;
 935
 936
 937     //Single photon spectra
 938     Double_t pt= p.Pt() ;
 939     TString skey="PHOS_single_"; skey+="mod" ; skey+=iMod ; skey+="_all" ;
 940     FillHistogram(skey,pt) ;
 941     FillHistogram("PHOS_single_all_mult",pt,fCentr) ;
 942     if(isDispOK){
 943       skey="PHOS_single_"; skey+="mod" ; skey+=iMod ; skey+="_disp" ;
 944       FillHistogram("PHOS_single_disp_mult",pt,fCentr) ;
 945       FillHistogram(skey,pt) ;
 946     }
 947     if(isNeutral){
 948       skey="PHOS_single_"; skey+="mod" ; skey+=iMod ; skey+="_neutral" ;
 949       FillHistogram(skey,pt) ;
 950       FillHistogram("PHOS_single_neu_mult",pt,fCentr) ;
 951     }
 952     if(isNeutral && isDispOK){
 953       skey="PHOS_single_"; skey+="mod" ; skey+=iMod ; skey+="_dispneutral" ;
 954       FillHistogram(skey,pt) ;
 955     }
 956     //Distance to bad channel
 957     if(clu->GetDistanceToBadChannel()>fBadDistCutPHOS){
 958       skey="PHOS_single_"; skey+="mod" ; skey+=iMod ; skey+="_dist1" ;
 959       FillHistogram(skey,pt) ;
 960     }
 961     if(clu->GetDistanceToBadChannel()>2.*fBadDistCutPHOS){
 962       skey="PHOS_single_"; skey+="mod" ; skey+=iMod ; skey+="_dist2" ;
 963       FillHistogram(skey,pt) ;
 964     }
 965
 966     //Fill QA
 967     if(clu->E()>0.5){
 968       skey="hQA_PHOS_mod"; skey+=iMod; skey+="_soft" ;
 969       FillHistogram(skey,iX-0.5, iZ-0.5,1.) ;
 970       if(clu->E()>1.5){
 971         skey="hQA_PHOS_mod"; skey+=iMod; skey+="_hard" ;
 972         FillHistogram(skey,iX-0.5, iZ-0.5,1.) ;
 973       }
 974     }
 975
 976     //Fill histogams for calibration
 977     if(clu->E()<fPi0Thresh1 ) continue;
 978     for(Int_t iconv=0;iconv<fConvEvent->GetEntriesFast();iconv++){
 979       TLorentzVector *gammaConv=static_cast<TLorentzVector*>(fConvEvent->At(iconv)) ;
 980       pi0=*gammaConv+p;
 981       skey="PHOS_"; skey+="mod" ; skey+=iMod ; skey+="_th1" ;
 982       FillHistogram(skey,iX-0.5, iZ-0.5,pi0.M());
 983       if(isNeutral){
 984         skey="PHOS_"; skey+="mod"; skey+=iMod ; skey+="_th2" ;
 985         FillHistogram(skey,iX-0.5, iZ-0.5,pi0.M());
 986       }
 987     }
 988   }
 989 }
 990 //____________________________________________________________
 991 void AliAnalysisTaskCaloConv::SelectEMCALPhotons(){
 992   //SelectEMCALPhotons
 993   // Loop over all CaloClusters
 994   if(fEMCALEvent)
 995     fEMCALEvent->Clear() ;
 996   else
 997     fEMCALEvent = new TClonesArray("TLorentzVector",10) ;
 998   Int_t inEMCAL = 0 ; //, inEMCALRecal=0;
 999   TLorentzVector pi0 ;
1000
1001   //vertex
1002   Double_t vtx[3]={0.,0.,0.};
1003   vtx[0] = fESDEvent->GetPrimaryVertex()->GetX();
1004   vtx[1] = fESDEvent->GetPrimaryVertex()->GetY();
1005   vtx[2] = fESDEvent->GetPrimaryVertex()->GetZ();
1006   Int_t nEMCAL=0 ; //For QA
1007   for (Int_t i=0; i<fESDEvent->GetNumberOfCaloClusters(); i++) {
1008     AliESDCaloCluster * clu = fESDEvent->GetCaloCluster(i);
1009     if(clu->IsPHOS())
1010       continue ;
1011     TLorentzVector p ;
1012     TLorentzVector pRecal ;
1013     clu ->GetMomentum(p ,vtx);
1014     Bool_t isNeutral = kTRUE ;
1015     Bool_t isDispOK = kTRUE ;
1016     Bool_t isTOFOK = kTRUE ;
1017     Int_t iMod,iX,iZ ;
1018
1019     if(clu->E()>0.1 ){
1020       nEMCAL++ ;
1021       isNeutral = clu->GetEmcCpvDistance()>10. ;  //To be improved
1022       if(clu->GetTOF()>550.e-9 && clu->GetTOF()<750.e-9)
1023         isTOFOK=kTRUE ;
1024       else
1025         isTOFOK=kFALSE ;
1026       Float_t phi = p.Phi();
1027       if(phi < 0) phi+=TMath::TwoPi();
1028       Int_t absId = -1;
1029       fEMCALgeom->GetAbsCellIdFromEtaPhi(p.Eta(),phi, absId);
1030       iMod=fEMCALgeom->GetSuperModuleNumber(absId) ;
1031       Int_t imod = -1, iTower = -1, iIphi = -1, iIeta = -1, iphi = -1, ieta = -1;
1032       fEMCALgeom->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
1033       fEMCALgeom->GetCellPhiEtaIndexInSModule(imod,iTower, iIphi, iIeta,iphi,ieta);
1034       if(imod<0 || imod>5){
1035         printf("EMCAL: Beyond the geometry!\n") ;
1036         printf("phi=%f, eta=%f, absId=%d, SM=%d \n",p.Eta(),phi, absId, imod) ;
1037         continue ;
1038       }
1039
1040       iX=iphi+1 ;
1041       iZ=ieta+1 ;
1042       if(!IsGoodChannel("EMCAL",iMod,iX,iZ))
1043         continue ;
1044       p.SetBit(kCaloPIDdisp,isDispOK) ;
1045       p.SetBit(kCaloPIDtof,isTOFOK) ;
1046       p.SetBit(kCaloPIDneutral,isNeutral) ;
1047       p.SetBit(BIT(17+imod),kTRUE) ;
1048       new((*fEMCALEvent)[inEMCAL]) TLorentzVector(p) ;
1049       fGammaEMCAL[inEMCAL] = i ;
1050       inEMCAL++ ;
1051
1052
1053       //Fill QA histograms
1054       if(clu->E()>0.5 && iMod>=0){ //Sometimes modules is negative not found??
1055         TString skey="hQA_EMCAL_SM";skey+=iMod ; skey+="_soft" ;
1056         FillHistogram(skey,iX-0.5, iZ-0.5,1.) ;
1057         if(clu->E()>1.5){
1058           skey="hQA_EMCAL_SM";skey+=iMod ; skey+="_hard" ;
1059           FillHistogram(skey,iX-0.5, iZ-0.5,1.) ;
1060         }
1061       }
1062
1063       //Fill histograms for recalibration
1064       if(clu->E()<fPi0Thresh1) continue ;
1065       for(Int_t iconv=0;iconv<fConvEvent->GetEntriesFast();iconv++){
1066         TLorentzVector *gammaConv=static_cast<TLorentzVector*>(fConvEvent->At(iconv)) ;
1067         pi0=*gammaConv+p;
1068         TString skey="EMCAL_"; skey+="mod" ; skey+=iMod ; skey+="_th1" ;
1069         FillHistogram(skey,iX-0.5, iZ-0.5,pi0.M());
1070         if(clu->E()>fPi0Thresh2){
1071           skey="EMCAL_"; skey+="mod"; skey+=iMod ; skey+="_th2" ;
1072           FillHistogram(skey,iX-0.5, iZ-0.5,pi0.M());
1073         }
1074       }
1075     }
1076   }
1077 }
1078 //______________________________________________________________________
1079 void AliAnalysisTaskCaloConv::SelectConvPhotons(){
1080   //Fill list of conversion photons
1081   //that is scan v0s and select photon-like
1082
1083   //set some constants
1084   const Double_t cutSigmaMass=0.0001;  //Constraint on photon mass
1085   const Bool_t useImprovedVertex=kTRUE ; //Use verted with converted photon?
1086 //  const Double_t zrSlope = TMath::Tan(2*TMath::ATan(TMath::Exp(-fetaCut)));
1087   const Double_t zrSlope12 = TMath::Tan(2*TMath::ATan(TMath::Exp(-1.2)));
1088   const Double_t zrSlope09 = TMath::Tan(2*TMath::ATan(TMath::Exp(-0.9)));
1089   const Double_t zOffset = 7.;
1090
1091   if(!fConvEvent)
1092     fConvEvent = new TClonesArray("TLorentzVector",10) ;
1093   else
1094     fConvEvent->Clear() ;
1095
1096   //No primary vertex in event => scip
1097   if(fESDEvent->GetPrimaryVertex()->GetNContributors()<=0) {
1098     return;
1099   }
1100
1101   Int_t inConv=0 ;
1102   for(Int_t iv0=0; iv0<fESDEvent->GetNumberOfV0s();iv0++){
1103     AliESDv0 * v0 = fESDEvent->GetV0(iv0) ;
1104
1105     AliESDtrack * pos = fESDEvent->GetTrack(v0->GetPindex()) ;
1106     AliESDtrack * neg = fESDEvent->GetTrack(v0->GetNindex()) ;
1107     const AliExternalTrackParam * paramPos = v0->GetParamP() ;
1108     const AliExternalTrackParam * paramNeg = v0->GetParamN() ;
1109     if(pos->GetSign() <0){//change tracks
1110       pos=neg ;
1111       neg=fESDEvent->GetTrack(v0->GetPindex()) ;
1112       paramPos=paramNeg ;
1113       paramNeg=v0->GetParamP() ;
1114     }
1115     AliKFParticle negKF(*paramNeg,11);
1116     AliKFParticle posKF(*paramPos,-11);
1117     AliKFParticle photKF(negKF,posKF) ;
1118     photKF.SetMassConstraint(0,cutSigmaMass);
1119
1120     if(useImprovedVertex){
1121       AliKFVertex primaryVertexImproved(*(fESDEvent->GetPrimaryVertex()));
1122       //if Vtx do created
1123       if(primaryVertexImproved.GetNContributors()>1){
1124         primaryVertexImproved+=photKF;
1125         photKF.SetProductionVertex(primaryVertexImproved);
1126       }
1127     }
1128     Double_t m=0., width=0. ;
1129     photKF.GetMass(m,width);
1130
1131     TLorentzVector photLV;
1132 //    photLV.SetXYZM(negKF.Px()+posKF.Px(),negKF.Py()+posKF.Py(),negKF.Pz()+negKF.Pz(),0.) ;
1133 //    photLV.SetXYZT(photKF.GetPx(),photKF.GetPy(),photKF.GetPz(),photKF.GetE()) ;
1134     photLV.SetXYZM(photKF.GetPx(),photKF.GetPy(),photKF.GetPz(),0.) ;  //Produces slightly better pi0 width
1135
1136     //Parameters for correction function
1137     Double_t a[3]={photLV.Pt(),photLV.Eta(),m} ;
1138     if(fToUseCF)
1139       fConvCFCont->Fill(a,0) ;
1140
1141     //select V0 finder
1142     Bool_t isOnFly=kTRUE ;
1143     //select V0Finder
1144     if (v0->GetOnFlyStatus()){
1145       if(fToUseCF)
1146         fConvCFCont->Fill(a,1) ;
1147     }
1148     else{
1149       isOnFly=kFALSE ;
1150       if(fToUseCF)
1151         fConvCFCont->Fill(a,2) ;
1152     }
1153
1154     //Number of TPC clusters
1155     if(neg->GetNcls(1) <2 || pos->GetNcls(1) <2){
1156       continue ;
1157     }
1158
1159     //remove like sign pairs
1160     if(pos->GetSign() == neg->GetSign()){
1161       continue ;
1162     }
1163     if(fToUseCF){
1164       if(isOnFly)
1165         fConvCFCont->Fill(a,3) ;
1166       else
1167         fConvCFCont->Fill(a,4) ;
1168     }
1169
1170     if( !(pos->GetStatus() & AliESDtrack::kTPCrefit) ||
1171         !(neg->GetStatus() & AliESDtrack::kTPCrefit) ){
1172       continue;
1173     }
1174     if(fToUseCF){
1175       if(isOnFly)
1176         fConvCFCont->Fill(a,5) ;
1177       else
1178         fConvCFCont->Fill(a,6) ;
1179     }
1180
1181     if( neg->GetKinkIndex(0) > 0 ||
1182         pos->GetKinkIndex(0) > 0) {
1183       continue ;
1184     }
1185
1186     //First rough PID
1187     if( fESDpid->NumberOfSigmasTPC(pos,AliPID::kElectron)<fnSigmaBelowElectronLine ||
1188         fESDpid->NumberOfSigmasTPC(pos,AliPID::kElectron)>fnSigmaAboveElectronLine ||
1189         fESDpid->NumberOfSigmasTPC(neg,AliPID::kElectron)<fnSigmaBelowElectronLine ||
1190         fESDpid->NumberOfSigmasTPC(neg,AliPID::kElectron)>fnSigmaAboveElectronLine ){
1191         continue ;
1192     }
1193     const Double_t minPnSigmaAbovePionLine = 1. ;
1194     const Double_t maxPnSigmaAbovePionLine = 3. ;
1195     const Double_t nSigmaAbovePionLine = 0 ;
1196     if(pos->P()>minPnSigmaAbovePionLine && pos->P()<maxPnSigmaAbovePionLine ){
1197       if(fESDpid->NumberOfSigmasTPC(pos,AliPID::kPion)<nSigmaAbovePionLine){
1198           continue ;
1199         }
1200     }
1201     if(neg->P()>minPnSigmaAbovePionLine && neg->P()<maxPnSigmaAbovePionLine){
1202       if(fESDpid->NumberOfSigmasTPC(neg,AliPID::kPion)<nSigmaAbovePionLine){
1203           continue ;
1204       }
1205     }
1206     //Strict dEdx
1207     Bool_t isdEdx=kTRUE;
1208     if(pos->P()>minPnSigmaAbovePionLine && pos->P()<maxPnSigmaAbovePionLine ){
1209       if(fESDpid->NumberOfSigmasTPC(pos,AliPID::kPion)<2.){
1210         isdEdx=kFALSE;
1211       }
1212     }
1213     if(neg->P()>minPnSigmaAbovePionLine && neg->P()<maxPnSigmaAbovePionLine){
1214       if(fESDpid->NumberOfSigmasTPC(neg,AliPID::kPion)<2.){
1215         isdEdx=kFALSE;
1216       }
1217     }
1218
1219
1220     //Kaon rejection
1221     const Double_t minPKaonRejection=1.5 ;
1222     const Double_t sigmaAroundLine=1. ;
1223     if(neg->P()<minPKaonRejection ){
1224       if(TMath::Abs(fESDpid->NumberOfSigmasTPC(neg,AliPID::kKaon))<sigmaAroundLine){
1225         isdEdx=kFALSE;
1226       }
1227     }
1228     if(pos->P()<minPKaonRejection ){
1229       if(TMath::Abs(fESDpid->NumberOfSigmasTPC(pos,AliPID::kKaon))<sigmaAroundLine){
1230         isdEdx=kFALSE;
1231       }
1232     }
1233
1234     //Proton rejection
1235     const Double_t minPProtonRejection=2. ;
1236     if(neg->P()<minPProtonRejection){
1237       if(TMath::Abs(fESDpid->NumberOfSigmasTPC(neg,AliPID::kProton))<sigmaAroundLine){
1238         isdEdx=kFALSE;
1239       }
1240     }
1241     if(pos->P()<minPProtonRejection ){
1242       if(TMath::Abs(fESDpid->NumberOfSigmasTPC(pos,AliPID::kProton))<sigmaAroundLine){
1243         isdEdx=kFALSE;
1244       }
1245     }
1246
1247     const Double_t minPPionRejection=0.5 ;
1248     if(neg->P()<minPPionRejection ){
1249       if(TMath::Abs(fESDpid->NumberOfSigmasTPC(neg,AliPID::kPion))<sigmaAroundLine){
1250         isdEdx=kFALSE;
1251       }
1252     }
1253     if(pos->P()<minPPionRejection ){
1254       if( TMath::Abs(fESDpid->NumberOfSigmasTPC(pos,AliPID::kPion))<sigmaAroundLine){
1255         isdEdx=kFALSE;
1256       }
1257     }
1258
1259
1260     if(isdEdx){
1261       FillHistogram("hdEdx",paramPos->GetP(),pos->GetTPCsignal()) ;
1262       FillHistogram("hdEdx",paramNeg->GetP(),neg->GetTPCsignal()) ;
1263     }
1264
1265
1266     //Check the pid probability
1267     Bool_t isProb=kTRUE ;
1268     Double_t posProbArray[10];
1269     Double_t negProbArray[10];
1270     neg->GetTPCpid(negProbArray);
1271     pos->GetTPCpid(posProbArray);
1272     if(negProbArray[AliPID::kElectron]<fprobCut || posProbArray[AliPID::kElectron]<fprobCut){
1273       isProb=kFALSE ;
1274     }
1275     if(fToUseCF){
1276       if(isOnFly)
1277         fConvCFCont->Fill(a,9) ;
1278       else
1279         fConvCFCont->Fill(a,10) ;
1280     }
1281
1282     Double_t v0x=0.,v0y=0.,v0z=0.;
1283     v0->GetXYZ(v0x,v0y,v0z) ;
1284     Double_t r=TMath::Sqrt(v0x*v0x + v0y*v0y) ;
1285     //Remove Dalitz
1286     const Double_t rMin=2.8 ;
1287     if(r<rMin)
1288       continue ;
1289     if(r>fmaxR){ // cuts on distance from collision point
1290       continue;
1291     }
1292     Bool_t isStrictR=kFALSE ;
1293     if(r<120.)
1294       isStrictR=kTRUE ;
1295     if(fToUseCF){
1296       if(isOnFly)
1297         fConvCFCont->Fill(a,11) ;
1298       else
1299         fConvCFCont->Fill(a,12) ;
1300     }
1301
1302
1303     if((TMath::Abs(v0z)*zrSlope12)-zOffset > r ){ // cuts out regions where we do not reconstruct
1304       continue;
1305     }
1306     if(fToUseCF){
1307       if(isOnFly)
1308         fConvCFCont->Fill(a,13) ;
1309       else
1310         fConvCFCont->Fill(a,14) ;
1311     }
1312
1313     if(TMath::Abs(v0z) > fmaxZ ){ // cuts out regions where we do not reconstruct
1314       continue;
1315     }
1316     Bool_t isStrictZ=kFALSE ;
1317     if((TMath::Abs(v0z)*zrSlope09)-zOffset < r )
1318       isStrictZ=kTRUE ;
1319
1320     if(fToUseCF){
1321       if(isOnFly)
1322         fConvCFCont->Fill(a,15) ;
1323       else
1324         fConvCFCont->Fill(a,16) ;
1325     }
1326
1327     if(photKF.GetNDF()<=0){
1328       continue;
1329     }
1330     if(fToUseCF){
1331       if(isOnFly)
1332         fConvCFCont->Fill(a,17) ;
1333       else
1334         fConvCFCont->Fill(a,18) ;
1335     }
1336
1337     Double_t chi2V0 = photKF.GetChi2()/photKF.GetNDF();
1338     FillHistogram("All_chi2_eta_pt",chi2V0,photLV.Eta(),photLV.Pt()) ;
1339
1340     if(chi2V0 > fchi2CutConversion || chi2V0 <=0){
1341       continue;
1342     }
1343     Bool_t isStrictChi=kFALSE ;
1344     if(chi2V0 < 0.7*fchi2CutConversion && chi2V0 >0){
1345       isStrictChi=kTRUE;
1346     }
1347
1348     if(fToUseCF){
1349       if(isOnFly)
1350         fConvCFCont->Fill(a,19) ;
1351       else
1352         fConvCFCont->Fill(a,20) ;
1353     }
1354
1355     const Double_t wideEtaCut=1.2 ;
1356     if(TMath::Abs(photLV.Eta())> wideEtaCut){
1357       continue;
1358     }
1359     if(TMath::Abs(paramPos->Eta())> wideEtaCut ||
1360        TMath::Abs(paramNeg->Eta())> wideEtaCut ){
1361       continue ;
1362     }
1363
1364     Bool_t isWideEta=kTRUE ;
1365     if(TMath::Abs(photLV.Eta())< fetaCut &&
1366        TMath::Abs(paramPos->Eta())<fetaCut  &&
1367        TMath::Abs(paramNeg->Eta()) < fetaCut){
1368       isWideEta=kFALSE;
1369     }
1370
1371
1372     if(photLV.Pt()<fptCut){
1373       continue;
1374     }
1375     if(fToUseCF){
1376       if(isOnFly)
1377         fConvCFCont->Fill(a,21) ;
1378       else
1379         fConvCFCont->Fill(a,22) ;
1380     }
1381
1382
1383     //Just QA plot
1384     if(photLV.Pt()>0.5){
1385        Double_t phi=photLV.Phi() ;
1386        while(phi<0.)phi+=TMath::TwoPi() ;
1387        while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
1388        FillHistogram("hQA_ConvPhiEta",phi,photLV.Eta()) ;
1389     }
1390
1391     Double_t w=PlanarityAngle(paramPos,paramNeg) ;
1392     Bool_t isPlanarityCut = (0.08-0.22*w > m || 0.15*(w-2.4)>m) ;
1393     FillHistogram("All_w_vs_m",w,m) ;
1394
1395     const Double_t armenterosAlphaCut=0.05 ;
1396     Double_t armenterosQtAlfa[2]={0.,0.}  ;
1397     GetArmenterosQtAlfa(&posKF, &negKF, &photKF, armenterosQtAlfa ) ;
1398     Bool_t isArmQt=(armenterosQtAlfa[1]<armenterosAlphaCut) ;
1399
1400     photLV.SetBit(kConvOnFly,isOnFly) ;
1401     photLV.SetBit(kConvArmQt,isArmQt) ;
1402     photLV.SetBit(kConvdEdx,isdEdx) ;
1403     photLV.SetBit(kConvProb,isProb) ;
1404     photLV.SetBit(kConvR,isStrictR) ;
1405     photLV.SetBit(kConvZR,isStrictZ) ;
1406     photLV.SetBit(kConvNDF,isStrictChi) ;
1407     photLV.SetBit(kConvEta,isWideEta) ;
1408     photLV.SetBit(kConvPlan,isPlanarityCut) ;
1409
1410     new((*fConvEvent)[inConv]) TLorentzVector(photLV) ;
1411     fGammaV0s[inConv] = iv0 ;
1412     inConv++ ;
1413
1414     //Single photon spectrum
1415     Double_t pt=photLV.Pt() ;
1416     if(isOnFly){
1417       //Default
1418       if(!isWideEta)
1419         FillHistogram("Single_conv_OnFly",pt) ;
1420       if(isdEdx && !isWideEta)
1421         FillHistogram("Single_conv_On_dEdx",pt) ;
1422       if(!isWideEta && isStrictR)
1423         FillHistogram("Single_conv_On_R120",pt) ;
1424       if( !isWideEta && isStrictZ)
1425         FillHistogram("Single_conv_On_Z",pt) ;
1426       if(!isWideEta && isStrictChi)
1427         FillHistogram("Single_conv_On_chi",pt) ;
1428       if(1)
1429         FillHistogram("Single_conv_On_Eta",pt) ;
1430       if(!isWideEta && isPlanarityCut)
1431         FillHistogram("Single_conv_On_Wcut",pt) ;
1432       if(!isWideEta && isArmQt)
1433         FillHistogram("Single_conv_On_ArmQt",pt) ;
1434     }
1435     else{
1436       if(!isWideEta)
1437         FillHistogram("Single_conv_Offline",pt) ;
1438       if(isdEdx && !isWideEta)
1439         FillHistogram("Single_conv_Off_dEdx",pt) ;
1440       if(!isWideEta && isStrictR)
1441         FillHistogram("Single_conv_Off_R120",pt) ;
1442       if(!isWideEta && isStrictZ)
1443         FillHistogram("Single_conv_Off_Z",pt) ;
1444       if(!isWideEta && isStrictChi)
1445         FillHistogram("Single_conv_Off_chi",pt) ;
1446       if(1)
1447         FillHistogram("Single_conv_Off_Eta",pt) ;
1448       if(!isWideEta && isPlanarityCut)
1449         FillHistogram("Single_conv_Off_Wcut",pt) ;
1450       if(!isWideEta && isArmQt)
1451         FillHistogram("Single_conv_Off_ArmQt",pt) ;
1452     }
1453
1454     //Fill MC information
1455     if(fStack){
1456       TParticle * negativeMC = fStack->Particle(TMath::Abs(neg->GetLabel()));
1457       TParticle * positiveMC = fStack->Particle(TMath::Abs(pos->GetLabel()));
1458
1459       if(!negativeMC || !positiveMC)
1460           continue ;
1461
1462       if(negativeMC->GetMother(0) != positiveMC->GetMother(0))
1463         continue ;
1464
1465       if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){
1466         continue;
1467       }
1468       if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){
1469         continue;
1470       }
1471
1472       TParticle * v0Gamma = fStack->Particle(negativeMC->GetMother(0));
1473
1474       if(negativeMC->GetUniqueID() != 5 || positiveMC->GetUniqueID() !=5){ // id 5 is conversion
1475         continue;
1476       }
1477       if(v0Gamma->GetPdgCode() == 22){
1478         FillHistogram("MC_V0_pt_eta_phi",v0Gamma->Pt(),v0Gamma->Eta(),v0Gamma->Phi()) ;
1479         FillHistogram("MC_V0_m_eta_pt",m,v0Gamma->Eta(),v0Gamma->Pt()) ;
1480         FillHistogram("MC_V0_chi2_eta_pt",chi2V0,v0Gamma->Eta(),v0Gamma->Pt()) ;
1481         FillHistogram("MC_V0_w_vs_m",w,m) ;
1482       }
1483     }
1484   }
1485 }
1486 //______________________________________________________________________
1487 void AliAnalysisTaskCaloConv::FillRealMixed(){
1488   // Fills real (same event) and Mixed (different events) inv.mass dsitributions
1489   // Moves current event to the list of stored events at the end
1490
1491   Double_t vtx[3];
1492   vtx[0] = fESDEvent->GetPrimaryVertex()->GetX();
1493   vtx[1] = fESDEvent->GetPrimaryVertex()->GetY();
1494   vtx[2] = fESDEvent->GetPrimaryVertex()->GetZ();
1495
1496   //Vtx class z-bin
1497   Int_t zvtx = (Int_t)((vtx[2]+10.)/2.) ;
1498   if(zvtx<0)zvtx=0 ;
1499   if(zvtx>9)zvtx=9 ;
1500
1501   Double_t run = fESDEvent->GetRunNumber()+0.5;
1502   TString trigClasses = fESDEvent->GetFiredTriggerClasses();
1503   if(trigClasses.Contains("CINT1B-ABCE-NOPF-ALL"))
1504     FillHistogram("hRunTrigger",run,0.5) ;
1505   else
1506     FillHistogram("hRunTrigger",run,1.5) ;
1507
1508   FillHistogram("hEvents",0.5) ;
1509   FillHistogram("hVtxBin",zvtx-0.5) ;
1510   FillHistogram("hRunEvents",run) ;
1511
1512   //check if containers for mixed events exist
1513   //create new if necessary
1514   if(!fPHOSEvents[zvtx]) fPHOSEvents[zvtx]=new TList() ;
1515   if(!fEMCALEvents[zvtx]) fEMCALEvents[zvtx]=new TList() ;
1516   if(!fConvEvents[zvtx]) fConvEvents[zvtx]=new TList() ;
1517
1518   Int_t nPHOS=fPHOSEvent->GetEntriesFast() ;
1519   Int_t nEMCAL=fEMCALEvent->GetEntriesFast() ;
1520   Int_t nConv = fConvEvent->GetEntriesFast() ;
1521   //Some QA histograms
1522   //Calculate number of good converion photons
1523   Int_t nConvGood=0 ;
1524   for(Int_t iConv = 0; iConv<nConv; iConv++){
1525     TLorentzVector * cnv = static_cast<TLorentzVector*>(fConvEvent->At(iConv)) ;
1526     if(cnv->TestBit(kConvOnFly) && !cnv->TestBit(kConvEta)){
1527       nConvGood++ ;
1528     }
1529   }
1530   FillHistogram("hRunConvs",run,double(nConvGood)) ;
1531   FillHistogram("hRunPHOS", run,double(nPHOS)) ;
1532   FillHistogram("hRunEMCAL",run,double(nEMCAL)) ;
1533
1534   //Fill Real distributions
1535   for(Int_t iPHOS=0; iPHOS<nPHOS;iPHOS++){
1536     TLorentzVector * cal = static_cast<TLorentzVector*>(fPHOSEvent->At(iPHOS)) ;
1537     for(Int_t iConv = 0; iConv<nConv; iConv++){
1538       TLorentzVector * cnv = static_cast<TLorentzVector*>(fConvEvent->At(iConv)) ;
1539       TLorentzVector pi=*cal + *cnv ;
1540       Double_t alpha=TMath::Abs(cal->Energy()-cnv->Energy())/(cal->Energy()+cnv->Energy()) ;
1541       if(cnv->TestBit(kConvOnFly) && !cnv->TestBit(kConvEta)){
1542         FillHistogram("PHOS_Re_mvsPt_all",pi.M(),pi.Pt()) ;
1543         char keym[55] ;
1544         for(Int_t iw=0;iw<10;iw++){ //resolution
1545           for(Int_t in=0;in<10;in++){
1546             snprintf(keym,55,"hMC_nonlinearity_w%d_n%d",iw,in) ;
1547             Double_t mMod=0.,ptMod=0. ;
1548             Recalibrate(mMod, ptMod, cal, cnv, iw, in) ;
1549             FillHistogram(keym,mMod,ptMod) ;
1550             snprintf(keym,55,"hMC_nonlinearity_ConvPHOS_w%d_n%d",iw,in) ;
1551             RecalibrateConvPHOS(mMod, ptMod, cal, cnv, iw, in) ;
1552             FillHistogram(keym,mMod,ptMod) ;
1553
1554           }
1555         }
1556         if(cal->TestBit(kCaloDistBad))
1557           FillHistogram("PHOS_Re_mvsPt_all_dist",pi.M(),pi.Pt()) ;
1558         FillHistogram("PHOS_Re_mvsPt_E",pi.M(),pi.Pt(),cal->Energy()) ;
1559         FillHistogram("PHOS_Re_mvsPt_alpha",pi.M(),pi.Pt(),alpha) ;
1560         if(cal->TestBit(kCaloPIDdisp))
1561           FillHistogram("PHOS_Re_mvsPt_Disp",pi.M(),pi.Pt()) ;
1562         if(cal->TestBit(kCaloPIDtof))
1563           FillHistogram("PHOS_Re_mvsPt_TOF",pi.M(),pi.Pt()) ;
1564         if(cal->TestBit(kCaloPIDneutral))
1565           FillHistogram("PHOS_Re_mvsPt_Neutral",pi.M(),pi.Pt()) ;
1566         if(cal->TestBit(kCaloPIDneutral) && cal->TestBit(kCaloPIDdisp))
1567           FillHistogram("PHOS_Re_mvsPt_DispNeutral",pi.M(),pi.Pt()) ;
1568       }
1569       //Vary Conversion cuts
1570       if(cnv->TestBit(kConvOnFly)){
1571         //Default
1572         if(!cnv->TestBit(kConvEta)){
1573           FillHistogram("PHOS_Re_mvsPt_OnFly",pi.M(),pi.Pt()) ;
1574           FillHistogram("PHOS_Re_mvsPt_OnFly_mult",pi.M(),pi.Pt(),fCentr) ;
1575         }
1576         if(cnv->TestBit(kConvdEdx) && !cnv->TestBit(kConvEta))
1577           FillHistogram("PHOS_Re_mvsPt_On_dEdx",pi.M(),pi.Pt()) ;
1578         if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvR))
1579           FillHistogram("PHOS_Re_mvsPt_On_R120",pi.M(),pi.Pt()) ;
1580         if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvZR))
1581           FillHistogram("PHOS_Re_mvsPt_On_Z",pi.M(),pi.Pt()) ;
1582         if( !cnv->TestBit(kConvEta) && cnv->TestBit(kConvNDF))
1583           FillHistogram("PHOS_Re_mvsPt_On_chi",pi.M(),pi.Pt()) ;
1584         if(1)
1585           FillHistogram("PHOS_Re_mvsPt_On_Eta",pi.M(),pi.Pt()) ;
1586         if( !cnv->TestBit(kConvEta) && cnv->TestBit(kConvPlan)){
1587           FillHistogram("PHOS_Re_mvsPt_On_Wcut",pi.M(),pi.Pt()) ;
1588           if(cal->TestBit(kCaloPIDneutral))
1589             FillHistogram("PHOS_Re_mvsPt_On_Wcut_Neu",pi.M(),pi.Pt()) ;
1590         }
1591         if( !cnv->TestBit(kConvEta) && cnv->TestBit(kConvArmQt))
1592           FillHistogram("PHOS_Re_mvsPt_On_ArmQt",pi.M(),pi.Pt()) ;
1593       }
1594       else{
1595         //Default
1596         if(!cnv->TestBit(kConvEta)){
1597           FillHistogram("PHOS_Re_mvsPt_Offline",pi.M(),pi.Pt()) ;
1598           FillHistogram("PHOS_Re_mvsPt_Offline_mult",pi.M(),pi.Pt(),fCentr) ;
1599         }
1600         if(cnv->TestBit(kConvdEdx) && !cnv->TestBit(kConvEta))
1601           FillHistogram("PHOS_Re_mvsPt_Off_dEdx",pi.M(),pi.Pt()) ;
1602         if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvR))
1603           FillHistogram("PHOS_Re_mvsPt_Off_R120",pi.M(),pi.Pt()) ;
1604         if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvZR))
1605           FillHistogram("PHOS_Re_mvsPt_Off_Z",pi.M(),pi.Pt()) ;
1606         if( !cnv->TestBit(kConvEta) && cnv->TestBit(kConvNDF))
1607           FillHistogram("PHOS_Re_mvsPt_Off_chi",pi.M(),pi.Pt()) ;
1608         if(1)
1609           FillHistogram("PHOS_Re_mvsPt_Off_Eta",pi.M(),pi.Pt()) ;
1610         if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvPlan)){
1611           FillHistogram("PHOS_Re_mvsPt_Off_Wcut",pi.M(),pi.Pt()) ;
1612           if(cal->TestBit(kCaloPIDneutral))
1613             FillHistogram("PHOS_Re_mvsPt_Off_Wcut_Neu",pi.M(),pi.Pt()) ;
1614         }
1615         if( !cnv->TestBit(kConvEta) && cnv->TestBit(kConvArmQt))
1616           FillHistogram("PHOS_Re_mvsPt_Off_ArmQt",pi.M(),pi.Pt()) ;
1617       }
1618     }
1619   }
1620   //PHOS module-dependent histograms
1621   for(Int_t iPHOS=0; iPHOS<nPHOS;iPHOS++){
1622     TLorentzVector * cal = static_cast<TLorentzVector*>(fPHOSEvent->At(iPHOS)) ;
1623     Int_t mod=1;
1624     while(!cal->TestBit(BIT(17+mod)) && mod<5)mod++ ;
1625     TString base("PHOS_Re_mvsPt_mod") ; base+=mod ;
1626     TString full ;
1627     for(Int_t iConv = 0; iConv<nConv; iConv++){
1628       TLorentzVector * cnv = static_cast<TLorentzVector*>(fConvEvent->At(iConv)) ;
1629       TLorentzVector pi=*cal + *cnv ;
1630       full=base ; full+="_single" ;
1631       if(cnv->TestBit(kConvOnFly) && !cnv->TestBit(kConvEta)){
1632         FillHistogram(full,pi.M(),cal->Pt()) ;
1633         full=base ; full+="_all" ;
1634         FillHistogram(full,pi.M(),pi.Pt()) ;
1635       }
1636     }
1637   }
1638   for(Int_t iEMCAL=0; iEMCAL<nEMCAL;iEMCAL++){
1639     TLorentzVector * cal = static_cast<TLorentzVector*>(fEMCALEvent->At(iEMCAL)) ;
1640     Int_t mod=0;
1641     while(!cal->TestBit(BIT(17+mod)) && mod<6)mod++ ;
1642     TString base("EMCAL_Re_mvsPt_mod") ; base+=mod ;
1643     TString full ;
1644     for(Int_t iConv = 0; iConv<nConv; iConv++){
1645       TLorentzVector * cnv = static_cast<TLorentzVector*>(fConvEvent->At(iConv)) ;
1646       TLorentzVector pi=*cal + *cnv ;
1647 //      Double_t alpha=TMath::Abs(cal->Energy()-cnv->Energy())/(cal->Energy()+cnv->Energy()) ;
1648       full=base+"_single" ;
1649       if(cnv->TestBit(kConvOnFly) && !cnv->TestBit(kConvEta)){
1650         FillHistogram(full,pi.M(),cal->Pt()) ;
1651         full=base+"_all" ;
1652         FillHistogram(full,pi.M(),pi.Pt()) ;
1653         char keym[55] ;
1654         for(Int_t iw=0;iw<10;iw++){ //resolution
1655           for(Int_t in=0;in<10;in++){
1656             snprintf(keym,55,"hMC_nonlinearity_EMCAL_w%d_n%d",iw,in) ;
1657             Double_t mMod=0.,ptMod=0. ;
1658             RecalibrateEMCAL(mMod, ptMod, cal, cnv, iw, in) ;
1659             FillHistogram(keym,mMod,ptMod) ;
1660           }
1661         }
1662       }
1663       //Vary Conversion cuts
1664       if(cnv->TestBit(kConvOnFly)){
1665         //Default
1666         if(!cnv->TestBit(kConvEta)){
1667           FillHistogram("EMCAL_Re_mvsPt_OnFly",pi.M(),pi.Pt()) ;
1668           FillHistogram("EMCAL_Re_mvsPt_OnFly_mult",pi.M(),pi.Pt(),fCentr) ;
1669         }
1670         if(cnv->TestBit(kConvdEdx) && !cnv->TestBit(kConvEta))
1671           FillHistogram("EMCAL_Re_mvsPt_On_dEdx",pi.M(),pi.Pt()) ;
1672         if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvR))
1673           FillHistogram("EMCAL_Re_mvsPt_On_R120",pi.M(),pi.Pt()) ;
1674         if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvZR))
1675           FillHistogram("EMCAL_Re_mvsPt_On_Z",pi.M(),pi.Pt()) ;
1676         if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvNDF))
1677           FillHistogram("EMCAL_Re_mvsPt_On_chi",pi.M(),pi.Pt()) ;
1678         if(1)
1679           FillHistogram("EMCAL_Re_mvsPt_On_Eta",pi.M(),pi.Pt()) ;
1680         if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvPlan))
1681           FillHistogram("EMCAL_Re_mvsPt_On_Wcut",pi.M(),pi.Pt()) ;
1682         if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvArmQt))
1683           FillHistogram("EMCAL_Re_mvsPt_On_ArmQt",pi.M(),pi.Pt()) ;
1684       }
1685       else{
1686         //Default
1687         if(!cnv->TestBit(kConvEta)){
1688           FillHistogram("EMCAL_Re_mvsPt_Offline",pi.M(),pi.Pt()) ;
1689           FillHistogram("EMCAL_Re_mvsPt_Offline_mult",pi.M(),pi.Pt(),fCentr) ;
1690         }
1691         if(cnv->TestBit(kConvdEdx) && !cnv->TestBit(kConvEta))
1692           FillHistogram("EMCAL_Re_mvsPt_Off_dEdx",pi.M(),pi.Pt()) ;
1693         if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvR))
1694           FillHistogram("EMCAL_Re_mvsPt_Off_R120",pi.M(),pi.Pt()) ;
1695         if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvZR))
1696           FillHistogram("EMCAL_Re_mvsPt_Off_Z",pi.M(),pi.Pt()) ;
1697         if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvNDF))
1698           FillHistogram("EMCAL_Re_mvsPt_Off_chi",pi.M(),pi.Pt()) ;
1699         if(1)
1700           FillHistogram("EMCAL_Re_mvsPt_Off_Eta",pi.M(),pi.Pt()) ;
1701         if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvPlan))
1702           FillHistogram("EMCAL_Re_mvsPt_Off_Wcut",pi.M(),pi.Pt()) ;
1703         if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvArmQt))
1704           FillHistogram("EMCAL_Re_mvsPt_Off_ArmQt",pi.M(),pi.Pt()) ;
1705       }
1706     }
1707   }
1708   //Now fill mixed
1709   TList * prevPHOS = fPHOSEvents[zvtx] ;
1710   TList * prevEMCAL = fEMCALEvents[zvtx] ;
1711   TList * prevConv = fConvEvents[zvtx] ;
1712
1713   //PHOS
1714   for(Int_t iPHOS=0; iPHOS<nPHOS;iPHOS++){
1715     TLorentzVector * cal = static_cast<TLorentzVector*>(fPHOSEvent->At(iPHOS)) ;
1716     for(Int_t ev=0; ev<prevConv->GetSize();ev++){
1717       TClonesArray * mixConv = static_cast<TClonesArray*>(prevConv->At(ev)) ;
1718       for(Int_t iConv = 0; iConv<mixConv->GetEntriesFast(); iConv++){
1719         TLorentzVector * cnv = static_cast<TLorentzVector*>(mixConv->At(iConv)) ;
1720         TLorentzVector pi=*cal + *cnv ;
1721         if(cnv->TestBit(kConvOnFly) && !cnv->TestBit(kConvEta)){
1722           FillHistogram("PHOS_Mi_mvsPt_all",pi.M(),pi.Pt()) ;
1723           if(cal->TestBit(kCaloPIDdisp))
1724             FillHistogram("PHOS_Mi_mvsPt_Disp",pi.M(),pi.Pt()) ;
1725           if(cal->TestBit(kCaloPIDtof))
1726             FillHistogram("PHOS_Mi_mvsPt_TOF",pi.M(),pi.Pt()) ;
1727             if(cal->TestBit(kCaloPIDneutral))
1728             FillHistogram("PHOS_Mi_mvsPt_Neutral",pi.M(),pi.Pt()) ;
1729           if(cal->TestBit(kCaloPIDneutral) && cal->TestBit(kCaloPIDdisp))
1730             FillHistogram("PHOS_Mi_mvsPt_DispNeutral",pi.M(),pi.Pt()) ;
1731         }
1732         //Vary Conversion cuts
1733         if(cnv->TestBit(kConvOnFly)){
1734           //Default
1735           if(!cnv->TestBit(kConvEta))
1736             FillHistogram("PHOS_Mi_mvsPt_OnFly",pi.M(),pi.Pt()) ;
1737           if( cnv->TestBit(kConvdEdx) && !cnv->TestBit(kConvEta))
1738             FillHistogram("PHOS_Mi_mvsPt_On_dEdx",pi.M(),pi.Pt()) ;
1739           if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvR))
1740             FillHistogram("PHOS_Mi_mvsPt_On_R120",pi.M(),pi.Pt()) ;
1741           if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvZR))
1742             FillHistogram("PHOS_Mi_mvsPt_On_Z",pi.M(),pi.Pt()) ;
1743           if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvNDF))
1744             FillHistogram("PHOS_Mi_mvsPt_On_chi",pi.M(),pi.Pt()) ;
1745           if(1)
1746             FillHistogram("PHOS_Mi_mvsPt_On_Eta",pi.M(),pi.Pt()) ;
1747           if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvPlan))
1748             FillHistogram("PHOS_Mi_mvsPt_On_Wcut",pi.M(),pi.Pt()) ;
1749           if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvArmQt))
1750             FillHistogram("PHOS_Mi_mvsPt_On_ArmQt",pi.M(),pi.Pt()) ;
1751         }
1752         else{
1753           //Default
1754           if(!cnv->TestBit(kConvEta))
1755             FillHistogram("PHOS_Mi_mvsPt_Offline",pi.M(),pi.Pt()) ;
1756           if(cnv->TestBit(kConvdEdx) && !cnv->TestBit(kConvEta))
1757             FillHistogram("PHOS_Mi_mvsPt_Off_dEdx",pi.M(),pi.Pt()) ;
1758           if( !cnv->TestBit(kConvEta) && cnv->TestBit(kConvR))
1759             FillHistogram("PHOS_Mi_mvsPt_Off_R120",pi.M(),pi.Pt()) ;
1760           if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvZR))
1761             FillHistogram("PHOS_Mi_mvsPt_Off_Z",pi.M(),pi.Pt()) ;
1762           if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvNDF))
1763             FillHistogram("PHOS_Mi_mvsPt_Off_chi",pi.M(),pi.Pt()) ;
1764           if(1)
1765             FillHistogram("PHOS_Mi_mvsPt_Off_Eta",pi.M(),pi.Pt()) ;
1766           if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvPlan))
1767             FillHistogram("PHOS_Mi_mvsPt_Off_Wcut",pi.M(),pi.Pt()) ;
1768           if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvArmQt))
1769             FillHistogram("PHOS_Mi_mvsPt_Off_ArmQt",pi.M(),pi.Pt()) ;
1770         }
1771       }
1772     }
1773   }
1774   for(Int_t iConv = 0; iConv<nConv; iConv++){
1775     TLorentzVector * cnv = static_cast<TLorentzVector*>(fConvEvent->At(iConv)) ;
1776     for(Int_t ev=0; ev<prevPHOS->GetSize();ev++){
1777       TClonesArray * mixPHOS = static_cast<TClonesArray*>(prevPHOS->At(ev)) ;
1778       for(Int_t iPHOS=0; iPHOS<mixPHOS->GetEntriesFast();iPHOS++){
1779         TLorentzVector * cal = static_cast<TLorentzVector*>(mixPHOS->At(iPHOS)) ;
1780         TLorentzVector pi=*cal + *cnv ;
1781         if(cnv->TestBit(kConvOnFly) && !cnv->TestBit(kConvEta)){
1782           FillHistogram("PHOS_Mi_mvsPt_all",pi.M(),pi.Pt()) ;
1783           if(cal->TestBit(kCaloPIDdisp))
1784             FillHistogram("PHOS_Mi_mvsPt_Disp",pi.M(),pi.Pt()) ;
1785           if(cal->TestBit(kCaloPIDtof))
1786             FillHistogram("PHOS_Mi_mvsPt_TOF",pi.M(),pi.Pt()) ;
1787           if(cal->TestBit(kCaloPIDneutral))
1788             FillHistogram("PHOS_Mi_mvsPt_Neutral",pi.M(),pi.Pt()) ;
1789           if(cal->TestBit(kCaloPIDneutral) && cal->TestBit(kCaloPIDdisp))
1790             FillHistogram("PHOS_Mi_mvsPt_DispNeutral",pi.M(),pi.Pt()) ;
1791         }
1792         //Vary Conversion cuts
1793         if(cnv->TestBit(kConvOnFly)){
1794           //Default
1795           if(!cnv->TestBit(kConvEta))
1796             FillHistogram("PHOS_Mi_mvsPt_OnFly",pi.M(),pi.Pt()) ;
1797           if(cnv->TestBit(kConvdEdx) && !cnv->TestBit(kConvEta))
1798             FillHistogram("PHOS_Mi_mvsPt_On_dEdx",pi.M(),pi.Pt()) ;
1799           if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvR))
1800             FillHistogram("PHOS_Mi_mvsPt_On_R120",pi.M(),pi.Pt()) ;
1801           if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvZR))
1802             FillHistogram("PHOS_Mi_mvsPt_On_Z",pi.M(),pi.Pt()) ;
1803           if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvNDF))
1804             FillHistogram("PHOS_Mi_mvsPt_On_chi",pi.M(),pi.Pt()) ;
1805           if(1)
1806             FillHistogram("PHOS_Mi_mvsPt_On_Eta",pi.M(),pi.Pt()) ;
1807           if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvPlan))
1808             FillHistogram("PHOS_Mi_mvsPt_On_Wcut",pi.M(),pi.Pt()) ;
1809           if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvArmQt))
1810             FillHistogram("PHOS_Mi_mvsPt_On_ArmQt",pi.M(),pi.Pt()) ;
1811         }
1812         else{
1813           //Default
1814           if(!cnv->TestBit(kConvEta))
1815             FillHistogram("PHOS_Mi_mvsPt_Offline",pi.M(),pi.Pt()) ;
1816           if(cnv->TestBit(kConvdEdx) && !cnv->TestBit(kConvEta))
1817             FillHistogram("PHOS_Mi_mvsPt_Off_dEdx",pi.M(),pi.Pt()) ;
1818           if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvR))
1819             FillHistogram("PHOS_Mi_mvsPt_Off_R120",pi.M(),pi.Pt()) ;
1820           if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvZR))
1821             FillHistogram("PHOS_Mi_mvsPt_Off_Z",pi.M(),pi.Pt()) ;
1822           if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvNDF))
1823             FillHistogram("PHOS_Mi_mvsPt_Off_chi",pi.M(),pi.Pt()) ;
1824           if(1)
1825             FillHistogram("PHOS_Mi_mvsPt_Off_Eta",pi.M(),pi.Pt()) ;
1826           if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvPlan))
1827             FillHistogram("PHOS_Mi_mvsPt_Off_Wcut",pi.M(),pi.Pt()) ;
1828           if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvArmQt))
1829             FillHistogram("PHOS_Mi_mvsPt_Off_ArmQt",pi.M(),pi.Pt()) ;
1830         }
1831       }
1832     }
1833   }
1834
1835 /*
1836   //PHOS module dependent
1837   for(Int_t iPHOS=0; iPHOS<nPHOS;iPHOS++){
1838     TLorentzVector * cal = static_cast<TLorentzVector*>(fPHOSEvent->At(iPHOS)) ;
1839     Int_t mod=1;
1840     while(!cal->TestBit(BIT(17+mod)) && mod<5)mod++ ;
1841     TString base("PHOS_Mi_mvsPt_mod") ; base+=mod ;
1842     TString full ;
1843     for(Int_t ev=0; ev<prevConv->GetSize();ev++){
1844       TClonesArray * mixConv = static_cast<TClonesArray*>(prevConv->At(ev)) ;
1845       for(Int_t iConv = 0; iConv<mixConv->GetEntriesFast(); iConv++){
1846         TLorentzVector * cnv = static_cast<TLorentzVector*>(mixConv->At(iConv)) ;
1847         TLorentzVector pi=*cal + *cnv ;
1848         full=base+"_all" ;
1849         if(cnv->TestBit(kConvOnFly) && !cnv->TestBit(kConvEta)){
1850           FillHistogram(full,pi.M(),pi.Pt()) ;
1851         }
1852       }
1853     }
1854   }
1855   for(Int_t iConv = 0; iConv<nConv; iConv++){
1856     TLorentzVector * cnv = static_cast<TLorentzVector*>(fConvEvent->At(iConv)) ;
1857     for(Int_t ev=0; ev<prevPHOS->GetSize();ev++){
1858       TClonesArray * mixPHOS = static_cast<TClonesArray*>(prevPHOS->At(ev)) ;
1859       for(Int_t iPHOS=0; iPHOS<mixPHOS->GetEntriesFast();iPHOS++){
1860         TLorentzVector * cal = static_cast<TLorentzVector*>(mixPHOS->At(iPHOS)) ;
1861         Int_t mod=1;
1862         while(!cal->TestBit(BIT(17+mod)) && mod<5)mod++ ;
1863         TString base("PHOS_Mi_mvsPt_mod") ; base+=mod ;
1864         TString full ;
1865         TLorentzVector pi=*cal + *cnv ;
1866         full=base+"_all" ;
1867         if(cnv->TestBit(kConvOnFly) && !cnv->TestBit(kConvEta)){
1868           FillHistogram(full,pi.M(),pi.Pt()) ;
1869         }
1870       }
1871     }
1872   }
1873 */
1874
1875   //EMCAL
1876   for(Int_t iEMCAL=0; iEMCAL<nEMCAL;iEMCAL++){
1877     TLorentzVector * cal = static_cast<TLorentzVector*>(fEMCALEvent->At(iEMCAL)) ;
1878     Int_t mod=0;
1879     while(!cal->TestBit(BIT(17+mod)) && mod<6)mod++ ;
1880     TString base("EMCAL_Mi_mvsPt_mod") ; base+=mod ;
1881     TString full ;
1882     for(Int_t ev=0; ev<prevConv->GetSize();ev++){
1883       TClonesArray * mixConv = static_cast<TClonesArray*>(prevConv->At(ev)) ;
1884       for(Int_t iConv = 0; iConv<mixConv->GetEntriesFast(); iConv++){
1885         TLorentzVector * cnv = static_cast<TLorentzVector*>(mixConv->At(iConv)) ;
1886         TLorentzVector pi=*cal + *cnv ;
1887         full=base+"_all" ;
1888         if(cnv->TestBit(kConvOnFly) && !cnv->TestBit(kConvEta)){
1889           FillHistogram(full,pi.M(),pi.Pt()) ;
1890         }
1891         if(cnv->TestBit(kConvOnFly)){
1892           //Default
1893           if(!cnv->TestBit(kConvEta))
1894             FillHistogram("EMCAL_Mi_mvsPt_OnFly",pi.M(),pi.Pt()) ;
1895           if(cnv->TestBit(kConvdEdx) && !cnv->TestBit(kConvEta))
1896             FillHistogram("EMCAL_Mi_mvsPt_On_dEdx",pi.M(),pi.Pt()) ;
1897           if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvR))
1898             FillHistogram("EMCAL_Mi_mvsPt_On_R120",pi.M(),pi.Pt()) ;
1899           if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvZR))
1900             FillHistogram("EMCAL_Mi_mvsPt_On_Z",pi.M(),pi.Pt()) ;
1901           if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvNDF))
1902             FillHistogram("EMCAL_Mi_mvsPt_On_chi",pi.M(),pi.Pt()) ;
1903           if(1)
1904             FillHistogram("EMCAL_Mi_mvsPt_On_Eta",pi.M(),pi.Pt()) ;
1905           if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvPlan))
1906             FillHistogram("EMCAL_Mi_mvsPt_On_Wcut",pi.M(),pi.Pt()) ;
1907           if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvArmQt))
1908             FillHistogram("EMCAL_Mi_mvsPt_On_ArmQt",pi.M(),pi.Pt()) ;
1909         }
1910         else{
1911           //Default
1912           if(!cnv->TestBit(kConvEta))
1913             FillHistogram("EMCAL_Mi_mvsPt_Offline",pi.M(),pi.Pt()) ;
1914           if( cnv->TestBit(kConvdEdx) && !cnv->TestBit(kConvEta))
1915             FillHistogram("EMCAL_Mi_mvsPt_Off_dEdx",pi.M(),pi.Pt()) ;
1916           if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvR))
1917             FillHistogram("EMCAL_Mi_mvsPt_Off_R120",pi.M(),pi.Pt()) ;
1918           if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvZR))
1919             FillHistogram("EMCAL_Mi_mvsPt_Off_Z",pi.M(),pi.Pt()) ;
1920           if( !cnv->TestBit(kConvEta) && cnv->TestBit(kConvNDF))
1921             FillHistogram("EMCAL_Mi_mvsPt_Off_chi",pi.M(),pi.Pt()) ;
1922           if(1)
1923             FillHistogram("EMCAL_Mi_mvsPt_Off_Eta",pi.M(),pi.Pt()) ;
1924           if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvPlan))
1925             FillHistogram("EMCAL_Mi_mvsPt_Off_Wcut",pi.M(),pi.Pt()) ;
1926           if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvArmQt))
1927             FillHistogram("EMCAL_Mi_mvsPt_Off_ArmQt",pi.M(),pi.Pt()) ;
1928         }
1929       }
1930     }
1931   }
1932 /*
1933   for(Int_t iConv = 0; iConv<nConv; iConv++){
1934     TLorentzVector * cnv = static_cast<TLorentzVector*>(fConvEvent->At(iConv)) ;
1935     for(Int_t ev=0; ev<prevEMCAL->GetSize();ev++){
1936       TClonesArray * mixEMCAL = static_cast<TClonesArray*>(prevEMCAL->At(ev)) ;
1937       for(Int_t iEMCAL=0; iEMCAL<mixEMCAL->GetEntriesFast();iEMCAL++){
1938         TLorentzVector * cal = static_cast<TLorentzVector*>(mixEMCAL->At(iEMCAL)) ;
1939         Int_t mod=0;
1940         while(!cal->TestBit(BIT(17+mod)) && mod<6)mod++ ;
1941         TString base("EMCAL_Mi_mvsPt_mod") ; base+=mod ;
1942         TString full ;
1943         TLorentzVector pi=*cal + *cnv ;
1944         full=base+"_all" ;
1945         if(cnv->TestBit(kConvOnFly) && !cnv->TestBit(kConvEta)){
1946           FillHistogram(full,pi.M(),pi.Pt()) ;
1947           if(cal->TestBit(kCaloPIDdisp)){
1948             full=base+"_Disp" ;
1949             FillHistogram(full,pi.M(),pi.Pt()) ;
1950           }
1951           if(cal->TestBit(kCaloPIDtof)){
1952             full=base+"_TOF" ;
1953             FillHistogram(full,pi.M(),pi.Pt()) ;
1954           }
1955           if(cal->TestBit(kCaloPIDneutral)){
1956             full=base+"_Neutral" ;
1957               FillHistogram(full,pi.M(),pi.Pt()) ;
1958           }
1959           if(cal->TestBit(kCaloPIDneutral) && cal->TestBit(kCaloPIDdisp)){
1960             full=base+"_DispNeutral" ;
1961             FillHistogram(full,pi.M(),pi.Pt()) ;
1962           }
1963         }
1964         if(cnv->TestBit(kConvOnFly)){
1965           //Default
1966           if(!cnv->TestBit(kConvEta))
1967             FillHistogram("EMCAL_Mi_mvsPt_OnFly",pi.M(),pi.Pt()) ;
1968           if(cnv->TestBit(kConvdEdx) && !cnv->TestBit(kConvEta))
1969             FillHistogram("EMCAL_Mi_mvsPt_On_dEdx",pi.M(),pi.Pt()) ;
1970           if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvR))
1971             FillHistogram("EMCAL_Mi_mvsPt_On_R120",pi.M(),pi.Pt()) ;
1972           if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvZR))
1973             FillHistogram("EMCAL_Mi_mvsPt_On_Z",pi.M(),pi.Pt()) ;
1974           if( !cnv->TestBit(kConvEta) && cnv->TestBit(kConvNDF))
1975             FillHistogram("EMCAL_Mi_mvsPt_On_chi",pi.M(),pi.Pt()) ;
1976           if(1)
1977             FillHistogram("EMCAL_Mi_mvsPt_On_Eta",pi.M(),pi.Pt()) ;
1978           if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvPlan))
1979             FillHistogram("EMCAL_Mi_mvsPt_On_Wcut",pi.M(),pi.Pt()) ;
1980           if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvArmQt))
1981             FillHistogram("EMCAL_Mi_mvsPt_On_ArmQt",pi.M(),pi.Pt()) ;
1982         }
1983         else{
1984           //Default
1985           if(!cnv->TestBit(kConvEta))
1986             FillHistogram("EMCAL_Mi_mvsPt_Offline",pi.M(),pi.Pt()) ;
1987           if(cnv->TestBit(kConvdEdx) && !cnv->TestBit(kConvEta))
1988             FillHistogram("EMCAL_Mi_mvsPt_Off_dEdx",pi.M(),pi.Pt()) ;
1989           if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvR))
1990             FillHistogram("EMCAL_Mi_mvsPt_Off_R120",pi.M(),pi.Pt()) ;
1991           if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvZR))
1992             FillHistogram("EMCAL_Mi_mvsPt_Off_Z",pi.M(),pi.Pt()) ;
1993           if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvNDF))
1994             FillHistogram("EMCAL_Mi_mvsPt_Off_chi",pi.M(),pi.Pt()) ;
1995           if(1)
1996             FillHistogram("EMCAL_Mi_mvsPt_Off_Eta",pi.M(),pi.Pt()) ;
1997           if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvPlan))
1998             FillHistogram("EMCAL_Mi_mvsPt_Off_Wcut",pi.M(),pi.Pt()) ;
1999           if(!cnv->TestBit(kConvEta) && cnv->TestBit(kConvArmQt))
2000             FillHistogram("EMCAL_Mi_mvsPt_Off_ArmQt",pi.M(),pi.Pt()) ;
2001         }
2002       }
2003     }
2004   }
2005 */
2006
2007
2008
2009   //Now we either add current events to stack or remove
2010   //Here we have some difficulty: conversion and calorimeter photons have different spectra.
2011   //So to correctly reproduce combinatorial background we have to preserve average number of
2012   //photons of both kinds per event. Therefore we should not reject empty PHOS/EMCAL events
2013   //though it will cost some memory. Reject only those events where no photons anywhere
2014
2015   if((fPHOSEvent->GetEntriesFast()>0 || fEMCALEvent->GetEntriesFast()>0) && fConvEvent->GetEntriesFast()>0){
2016     prevPHOS->AddFirst(fPHOSEvent) ;
2017     fPHOSEvent=0;
2018     prevEMCAL->AddFirst(fEMCALEvent) ;
2019     fEMCALEvent=0 ;
2020     prevConv->AddFirst(fConvEvent) ;
2021     fConvEvent=0 ;
2022     if(prevPHOS->GetSize()>100){//Remove redundant events
2023       TClonesArray * tmp = static_cast<TClonesArray*>(prevPHOS->Last()) ;
2024       prevPHOS->RemoveLast() ;
2025       delete tmp ;
2026       tmp = static_cast<TClonesArray*>(prevEMCAL->Last()) ;
2027       prevEMCAL->RemoveLast() ;
2028       delete tmp ;
2029       tmp = static_cast<TClonesArray*>(prevConv->Last()) ;
2030       prevConv->RemoveLast() ;
2031       delete tmp ;
2032     }
2033   }
2034
2035 }
2036 //___________________________________________________________________________
2037 void AliAnalysisTaskCaloConv::ProcessMC(){
2038   //ProcessMC
2039   //fill histograms for efficiensy etc. calculation
2040   if(!fStack) return ;
2041
2042   const Double_t rcut = 1. ; //cut for primary particles
2043   Double_t vtx[3];
2044   vtx[0] = fESDEvent->GetPrimaryVertex()->GetX();
2045   vtx[1] = fESDEvent->GetPrimaryVertex()->GetY();
2046   vtx[2] = fESDEvent->GetPrimaryVertex()->GetZ();
2047
2048   Int_t nPHOS=fPHOSEvent->GetEntriesFast() ;
2049   Int_t nEMCAL=fEMCALEvent->GetEntriesFast() ;
2050   Int_t nConv = fConvEvent->GetEntriesFast() ;
2051
2052   //---------First pi0/eta-----------------------------
2053   char partName[10] ;
2054   char hkey[55] ;
2055   for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++) {
2056     TParticle* particle = (TParticle *)fStack->Particle(iTracks);
2057     if(particle->GetPdgCode() == 111)
2058       snprintf(partName,10,"pi0") ;
2059     else
2060       if(particle->GetPdgCode() == 221)
2061         snprintf(partName,10,"eta") ;
2062       else
2063         continue ;
2064
2065     //Primary particle
2066     if(particle->R() >rcut)
2067       continue ;
2068
2069     Double_t pt = particle->Pt() ;
2070     //Total number of pi0 with creation radius <1 cm
2071     snprintf(hkey,55,"hMC_CaloConv_all%s",partName) ;
2072     FillHistogram(hkey,pt) ;
2073     if(TMath::Abs(particle->Y())<1.){
2074       snprintf(hkey,55,"hMC_CaloConv_%s_unitEta",partName) ;
2075       FillHistogram(hkey,pt,fCentr) ;
2076     }
2077
2078     //Check if one of photons converted
2079     if(particle->GetNDaughters()!=2)
2080      continue ; //Do not account Dalitz decays
2081
2082     TParticle * gamma1 = fStack->Particle(particle->GetFirstDaughter());
2083     TParticle * gamma2 = fStack->Particle(particle->GetLastDaughter());
2084     //Number of pi0s decayed into acceptance
2085     Bool_t inAcc1 = (TMath::Abs(gamma1->Eta())<0.9) ;
2086     Bool_t inAcc2 = (TMath::Abs(gamma2->Eta())<0.9) ;
2087     Int_t mod1,mod2 ;
2088     Double_t x=0.,z=0. ;
2089     Bool_t hitPHOS1 = fPHOSgeom->ImpactOnEmc(gamma1, mod1, z,x) ;
2090     Bool_t hitPHOS2 = fPHOSgeom->ImpactOnEmc(gamma2, mod2, z,x) ;
2091     Bool_t hitEMCAL1= fEMCALgeom->Impact(gamma1) ;
2092     Bool_t hitEMCAL2= fEMCALgeom->Impact(gamma2) ;
2093
2094     Bool_t goodPair=kFALSE ;
2095     if((inAcc1 && hitPHOS2) || (inAcc2 && hitPHOS1)){
2096       snprintf(hkey,55,"hMC_CaloConv_%sPHOSacc",partName) ;
2097       FillHistogram(hkey,pt) ;
2098       goodPair=kTRUE ;
2099     }
2100     if((inAcc1 && hitEMCAL2) || (inAcc2 && hitEMCAL1)){
2101       snprintf(hkey,55,"hMC_CaloConv_%sEMCALacc",partName) ;
2102       FillHistogram(hkey,pt) ;
2103        goodPair=kTRUE ;
2104     }
2105     if(!goodPair){
2106       continue ;
2107     }
2108
2109     Bool_t converted1 = kFALSE ;
2110     if(gamma1->GetNDaughters()==2){
2111       TParticle * e1=fStack->Particle(gamma1->GetFirstDaughter()) ;
2112       TParticle * e2=fStack->Particle(gamma1->GetLastDaughter()) ;
2113       if(TMath::Abs(e1->GetPdgCode())==11 && TMath::Abs(e2->GetPdgCode())==11){ //conversion
2114         if(e1->R()<180.)
2115           converted1 = kTRUE ;
2116       }
2117     }
2118     Bool_t converted2 = kFALSE ;
2119     if(gamma2->GetNDaughters()==2){
2120       TParticle * e1=fStack->Particle(gamma2->GetFirstDaughter()) ;
2121       TParticle * e2=fStack->Particle(gamma2->GetLastDaughter()) ;
2122       if(TMath::Abs(e1->GetPdgCode())==11 && TMath::Abs(e2->GetPdgCode())==11){ //conversion
2123         if(e1->R()<180.)
2124           converted2 = kTRUE ;
2125       }
2126     }
2127
2128     //Number of pi0s with one photon converted
2129     if((converted1 && !converted2 && hitPHOS2) || (!converted1 && hitPHOS1 && converted2)) {
2130       snprintf(hkey,55,"hMC_CaloConv_%s_PHOS_conv",partName) ;
2131       FillHistogram(hkey,pt) ;
2132     }
2133
2134     if((converted1 && !converted2 && hitEMCAL2) || (!converted1 && hitEMCAL1 && converted2)) {
2135       snprintf(hkey,55,"hMC_CaloConv_%s_EMCAL_conv",partName) ;
2136       FillHistogram(hkey,pt) ;
2137     }
2138
2139     //Both converted
2140     if(converted1 && converted2) {
2141       snprintf(hkey,55,"hMC_CaloConv_%s_bothphot_conv",partName) ;
2142       FillHistogram(hkey,pt) ;
2143         continue ;
2144     }
2145
2146     //photon pointing calorimeter converted
2147     if((converted1 && hitPHOS1 && !hitEMCAL2) || (converted2 && hitPHOS2 && !hitEMCAL1) ||
2148        (converted1 && hitEMCAL1 && !hitPHOS2) || (converted2 && hitEMCAL2 && !hitPHOS1)){
2149       snprintf(hkey,55,"hMC_CaloConv_%s_convPhotInCalo",partName) ;
2150       FillHistogram(hkey,pt) ;
2151        continue ;
2152     }
2153
2154     //Converted pi0 with v0 and photon PHOS or EMCAL
2155     Bool_t foundV01onfly=kFALSE, foundV01offline=kFALSE, foundV02onfly=kFALSE, foundV02offline=kFALSE ;
2156     Bool_t foundV01onflyPID=kFALSE, foundV01offlinePID=kFALSE, foundV02onflyPID=kFALSE, foundV02offlinePID=kFALSE ;
2157     TLorentzVector pConvOn,pConvOff ;
2158     for(Int_t iv0=0; iv0<fESDEvent->GetNumberOfV0s();iv0++){
2159       AliESDv0 * v0 = fESDEvent->GetV0(iv0) ;
2160
2161       TParticle * negativeMC = fStack->Particle(TMath::Abs(fESDEvent->GetTrack(v0->GetNindex())->GetLabel()));
2162       TParticle * positiveMC = fStack->Particle(TMath::Abs(fESDEvent->GetTrack(v0->GetPindex())->GetLabel()));
2163
2164       if(negativeMC->GetMother(0) != positiveMC->GetMother(0))
2165         continue ;
2166
2167       if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){
2168         continue;
2169        }
2170       if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){
2171         continue;
2172       }
2173
2174       TParticle * v0Gamma = fStack->Particle(negativeMC->GetMother(0));
2175       Bool_t same = (v0Gamma == gamma1) ;
2176       TParticle * tmp = v0Gamma ;
2177       while(!same && tmp->GetFirstMother()>=0){
2178         tmp = fStack->Particle(tmp->GetFirstMother());
2179        same = (tmp == gamma1) ;
2180      }
2181      if(same){
2182        if(v0->GetOnFlyStatus())
2183          foundV01onfly = kTRUE ;
2184        else
2185          foundV01offline= kTRUE ;
2186        for(Int_t iconv=0; iconv<nConv;iconv++){
2187          if(fGammaV0s[iconv] == iv0){
2188            TLorentzVector * cnv = static_cast<TLorentzVector*>(fConvEvent->At(iconv)) ;
2189            //default cuts
2190            if(!cnv->TestBit(kConvEta)){
2191              if(v0->GetOnFlyStatus()){
2192                pConvOn= *cnv ;
2193                foundV01onflyPID = kTRUE ;
2194              }
2195              else{
2196                pConvOff= *cnv ;
2197                foundV01offlinePID = kTRUE ;
2198              }
2199            }
2200            break ;
2201          }
2202        }
2203        continue ;
2204      }
2205      same = (v0Gamma == gamma2) ;
2206      tmp = v0Gamma ;
2207      while(!same && tmp->GetFirstMother()>=0){
2208        tmp = fStack->Particle(tmp->GetFirstMother());
2209        same = (tmp == gamma2) ;
2210      }
2211      if(same){
2212        if(v0->GetOnFlyStatus())
2213          foundV02onfly = kTRUE ;
2214        else
2215          foundV02offline = kTRUE ;
2216        for(Int_t iconv=0; iconv<nConv;iconv++){
2217          if(fGammaV0s[iconv] == iv0){
2218            TLorentzVector * cnv = static_cast<TLorentzVector*>(fConvEvent->At(iconv)) ;
2219            //default cuts
2220            if(!cnv->TestBit(kConvEta)){
2221              if(v0->GetOnFlyStatus()){
2222                pConvOn= *cnv ;
2223                foundV02onflyPID = kTRUE ;
2224              }
2225              else{
2226                pConvOff= *cnv ;
2227                foundV02offlinePID = kTRUE ;
2228              }
2229            }
2230            break ;
2231          }
2232        }
2233      }
2234    }
2235
2236    goodPair=kFALSE ;
2237    if((foundV01onfly && hitPHOS2) || (foundV02onfly && hitPHOS1)){
2238      snprintf(hkey,55,"hMC_CaloConv_%s_v0onfly_PHOSacc",partName) ;
2239      FillHistogram(hkey,pt) ;
2240      goodPair=kTRUE;
2241    }
2242    if((foundV01offline && hitPHOS2) || (foundV02offline && hitPHOS1)){
2243      snprintf(hkey,55,"hMC_CaloConv_%s_v0offline_PHOSacc",partName) ;
2244      FillHistogram(hkey,pt) ;
2245      goodPair=kTRUE;
2246    }
2247    if((foundV01onfly && hitEMCAL2) || (foundV02onfly && hitEMCAL1)){
2248      snprintf(hkey,55,"hMC_CaloConv_%s_v0onfly_EMCALacc",partName) ;
2249      FillHistogram(hkey,pt) ;
2250      goodPair=kTRUE;
2251    }
2252    if((foundV01offline && hitEMCAL2) || (foundV02offline && hitEMCAL1)){
2253      snprintf(hkey,55,"hMC_CaloConv_%s_v0offline_EMCALacc",partName) ;
2254      FillHistogram(hkey,pt) ;
2255      goodPair=kTRUE;
2256    }
2257    if(!goodPair){
2258      continue ;
2259    }
2260
2261    //Converted pi0 with v0 and cluster in PHOS/EMCAL
2262    Bool_t cluInPHOS = kFALSE,cluInEMCAL=kFALSE ;
2263    Bool_t cluInPHOSpid = kFALSE,cluInEMCALpid=kFALSE ;
2264    Bool_t closeToBad= kFALSE ;
2265    TLorentzVector pCalo ;
2266    for (Int_t i=0; i<fESDEvent->GetNumberOfCaloClusters(); i++) {
2267      AliESDCaloCluster * clu = fESDEvent->GetCaloCluster(i);
2268      Int_t iprim = clu->GetLabel() ; //# of particle hit PHOS/EMCAL
2269      Bool_t matched = kFALSE ;
2270      while(iprim>=0 ) {
2271        if(iprim==particle->GetFirstDaughter() || iprim==particle->GetLastDaughter()){
2272          matched=kTRUE ;
2273          break ;
2274        }
2275        else{
2276          iprim=fStack->Particle(iprim)->GetFirstMother() ;
2277        }
2278      }
2279      if(!matched)
2280        continue ;
2281      if(clu->IsPHOS() && (hitPHOS1 || hitPHOS2)){
2282        cluInPHOS=kTRUE ;
2283        //Check if cluster passed PID
2284        for(Int_t inPHOS=0; inPHOS<nPHOS;inPHOS++){
2285          if(fGammaPHOS[inPHOS] == i){
2286            cluInPHOSpid=kTRUE ;
2287            break ;
2288          }
2289        }
2290        clu->GetMomentum(pCalo ,vtx);
2291        if(clu->GetDistanceToBadChannel()<fBadDistCutPHOS)
2292          closeToBad=kTRUE ;
2293        break ;
2294      }
2295      if(!clu->IsPHOS() && (hitEMCAL1 || hitEMCAL2)){
2296        cluInEMCAL=kTRUE ;
2297        //Check if cluster passed PID
2298        for(Int_t inEMCAL=0; inEMCAL<nEMCAL;inEMCAL++){
2299          if(fGammaEMCAL[inEMCAL] == i){
2300            cluInPHOSpid=kTRUE ;
2301            break ;
2302          }
2303        }
2304        clu->GetMomentum(pCalo ,vtx);
2305        if(clu->GetDistanceToBadChannel()<fBadDistCutEMCAL)
2306          closeToBad=kTRUE ;
2307        break ;
2308      }
2309    }
2310
2311    if(cluInPHOS){
2312      //OnFly
2313      if(foundV01onfly ||foundV02onfly){
2314        snprintf(hkey,55,"hMC_CaloConv_%s_v0onfly_PHOSclu",partName) ;
2315        FillHistogram(hkey,pt) ;
2316
2317        if((foundV01onflyPID ||foundV02onflyPID) && cluInPHOSpid){
2318        snprintf(hkey,55,"hMC_CaloConv_%s_v0onfly_PHOSclu_pid",partName) ;
2319          FillHistogram(hkey,pt) ;
2320          for(Int_t iw=0;iw<10;iw++){ //resolution
2321            for(Int_t in=0;in<10;in++){
2322              char keym[55] ;
2323              snprintf(keym,55,"hMC_CaloConv_%s_v0onfly_PHOSclu_w%d_n%d",partName,iw,in) ;
2324              Double_t mMod=0.,ptMod=0. ;
2325              Recalibrate(mMod, ptMod, &pCalo, &pConvOn, iw, in) ;
2326              FillHistogram(keym,ptMod) ;
2327              snprintf(keym,55,"hMC_CaloConv_%s_v0onfly_ConvPHOSclu_w%d_n%d",partName,iw,in) ;
2328              RecalibrateConvPHOS(mMod, ptMod, &pCalo, &pConvOn, iw, in) ;
2329              FillHistogram(keym,ptMod) ;
2330            }
2331          }
2332          if(!closeToBad){
2333            snprintf(hkey,55,"hMC_CaloConv_%s_v0onfly_PHOSclu_good",partName) ;
2334            FillHistogram(hkey,pt) ;
2335          }
2336          Double_t m=(pCalo+pConvOn).M() ;
2337          Double_t ptm=(pCalo+pConvOn).Pt() ;
2338          snprintf(hkey,55,"hMC_CaloConv_%s_v0on_PHOSclu_ptRec",partName) ;
2339          FillHistogram(hkey,ptm) ;
2340          snprintf(hkey,55,"hMC_CaloConv_%s_v0on_PHOSclu_mvsPt",partName) ;
2341          FillHistogram(hkey,m,ptm) ;
2342        }
2343      }
2344
2345      //Offline
2346      if(foundV01offline ||foundV02offline){
2347        snprintf(hkey,55,"hMC_CaloConv_%s_v0offline_PHOSclu",partName) ;
2348        FillHistogram(hkey,pt) ;
2349        if((foundV01offlinePID ||foundV02offlinePID) && cluInPHOSpid){
2350          snprintf(hkey,55,"hMC_CaloConv_%s_v0offline_PHOSclu_pid",partName) ;
2351          FillHistogram(hkey,pt) ;
2352          Double_t m=(pCalo+pConvOff).M() ;
2353          Double_t ptm=(pCalo+pConvOff).Pt() ;
2354          snprintf(hkey,55,"hMC_CaloConv_%s_v0off_PHOSclu_ptRec",partName) ;
2355          FillHistogram(hkey,ptm) ;
2356          snprintf(hkey,55,"hMC_CaloConv_%s_v0off_PHOSclu_mvsPt",partName) ;
2357          FillHistogram(hkey,m,ptm) ;
2358          if(!closeToBad){
2359            snprintf(hkey,55,"hMC_CaloConv_%s_v0offline_PHOSclu_good",partName) ;
2360            FillHistogram(hkey,pt) ;
2361          }
2362        }
2363      }
2364
2365      if((foundV01onflyPID ||foundV02onflyPID) && cluInPHOSpid){
2366        TString base("hMC_CaloConv_") ; base+=partName; base+="_v0onfly_PHOSclu_mod" ;
2367        if(hitPHOS1)
2368          base+=mod1 ;
2369        else
2370          base+=mod2 ;
2371        FillHistogram(base.Data(),pt) ;
2372      }
2373
2374      if((foundV01offlinePID ||foundV02offlinePID) && cluInPHOSpid){
2375        TString base("hMC_CaloConv_") ; base+=partName; base+="_v0offline_PHOSclu_mod" ;
2376        if(hitPHOS1)
2377          base+=mod1 ;
2378        else
2379          base+=mod2 ;
2380        FillHistogram(base.Data(),pt) ;
2381      }
2382    }
2383    if(cluInEMCAL && strcmp(partName,"pi0")==0){
2384      //OnFly
2385      if(foundV01onfly ||foundV02onfly){
2386        FillHistogram("hMC_CaloConv_pi0_v0onfly_EMCALclu",pt) ;
2387
2388        if((foundV01onflyPID ||foundV02onflyPID) && cluInEMCALpid){
2389          FillHistogram("hMC_CaloConv_pi0_v0onfly_EMCALclu_pid",pt) ;
2390          for(Int_t iw=0;iw<10;iw++){ //resolution
2391            for(Int_t in=0;in<10;in++){
2392              char keym[55] ;
2393              snprintf(keym,55,"hMC_CaloConv_pi0_v0onfly_EMCALclu_w%d_n%d",iw,in) ;
2394              Double_t mMod=0.,ptMod=0. ;
2395              RecalibrateEMCAL(mMod, ptMod, &pCalo, &pConvOn, iw, in) ;
2396              FillHistogram(keym,ptMod) ;
2397            }
2398          }
2399          if(!closeToBad)
2400            FillHistogram("hMC_CaloConv_pi0_v0onfly_EMCALclu_good",pt) ;
2401          Double_t m=(pCalo+pConvOn).M() ;
2402          Double_t ptm=(pCalo+pConvOn).Pt() ;
2403          FillHistogram("hMC_CaloConv_pi0_v0on_EMCALclu_ptRec",ptm) ;
2404          FillHistogram("hMC_CaloConv_pi0_v0on_EMCALclu_mvsPt",m,ptm) ;
2405        }
2406      }
2407
2408      //Offline
2409      if(foundV01offline ||foundV02offline){
2410        FillHistogram("hMC_CaloConv_pi0_v0offline_EMCALclu",pt) ;
2411        if((foundV01offlinePID ||foundV02offlinePID) && cluInEMCALpid){
2412          FillHistogram("hMC_CaloConv_pi0_v0offline_EMCALclu_pid",pt) ;
2413          Double_t m=(pCalo+pConvOff).M() ;
2414          Double_t ptm=(pCalo+pConvOff).Pt() ;
2415          FillHistogram("hMC_CaloConv_pi0_v0off_EMCALclu_ptRec",ptm) ;
2416          FillHistogram("hMC_CaloConv_pi0_v0off_EMCALclu_mvsPt",m,ptm) ;
2417          if(!closeToBad)
2418            FillHistogram("hMC_CaloConv_pi0_v0offline_EMCALclu_good",pt) ;
2419        }
2420      }
2421    }
2422   }
2423
2424   //Construct Inv mass distributions for residual correlations
2425   if(fPHOSEvent && fConvEvent){
2426     for(Int_t iPHOS=0; iPHOS<fPHOSEvent->GetEntriesFast();iPHOS++){
2427       TLorentzVector * cal = static_cast<TLorentzVector*>(fPHOSEvent->At(iPHOS)) ;
2428       Int_t iclu=fGammaPHOS[iPHOS] ;
2429       AliESDCaloCluster * clu = fESDEvent->GetCaloCluster(iclu);
2430       Int_t iprimPHOS = clu->GetLabel() ; //# of particle hit PHOS/EMCAL
2431       for(Int_t iConv = 0; iConv<fConvEvent->GetEntriesFast(); iConv++){
2432         TLorentzVector * cnv = static_cast<TLorentzVector*>(fConvEvent->At(iConv)) ;
2433         if(!cnv->TestBit(kConvOnFly) || cnv->TestBit(kConvEta))
2434           continue;
2435
2436         Int_t iv0=fGammaV0s[iConv] ;
2437         AliESDv0 * v0 = fESDEvent->GetV0(iv0) ;
2438         Int_t iprimNeg = TMath::Abs(fESDEvent->GetTrack(v0->GetNindex())->GetLabel()) ;
2439         Int_t iprimPos = TMath::Abs(fESDEvent->GetTrack(v0->GetPindex())->GetLabel()) ;
2440
2441         //Check if there was a common ancistor
2442         Bool_t found = kFALSE ;
2443         Int_t curPHOS=iprimPHOS ;
2444         Int_t commonA=-1 ;
2445         while(!found && curPHOS>-1){
2446           Int_t curNeg=iprimNeg ;
2447           while(!found && curNeg>-1){
2448             if(curNeg==curPHOS){
2449               found=kTRUE ;
2450               commonA=curPHOS ;
2451             }
2452             else{
2453               curNeg=fStack->Particle(curNeg)->GetFirstMother() ;
2454             }
2455           }
2456           curPHOS=fStack->Particle(curPHOS)->GetFirstMother() ;
2457         }
2458         found = kFALSE ;
2459         curPHOS=iprimPHOS ;
2460         Int_t commonB=-1 ;
2461         while(!found && curPHOS>-1){
2462           Int_t curPos=iprimPos ;
2463           while(!found && curPos>-1){
2464             if(curPos==curPHOS){
2465               found=kTRUE ;
2466               commonB=curPHOS ;
2467             }
2468             else{
2469               curPos=fStack->Particle(curPos)->GetFirstMother() ;
2470             }
2471           }
2472           curPHOS=fStack->Particle(curPHOS)->GetFirstMother() ;
2473         }
2474         if(commonA != commonB){
2475            //Strange
2476            AliInfo(Form("CommonA=%d, commonB=%d",commonA,commonB)) ;
2477         }
2478         if(commonA>-1){//There was common particles
2479           Int_t pdg = fStack->Particle(commonA)->GetPdgCode() ;
2480           TLorentzVector pi=*cal + *cnv ;
2481           Double_t m=pi.M() ;
2482           Double_t pt=pi.Pt() ;
2483           Double_t alpha=TMath::Abs(cal->Energy()-cnv->Energy())/(cal->Energy()+cnv->Energy()) ;
2484           switch(pdg){
2485           case 11:
2486           case -11:
2487           case 22: //conversion
2488             FillHistogram("hMC_Resid_PHOS_Phot_mvsPt",m,pt,alpha) ;
2489             break ;
2490           case 111: //pi0
2491             FillHistogram("hMC_Resid_PHOS_Pi0_mvsPt",m,pt,alpha) ;
2492             break ;
2493           case 221: //eta
2494               FillHistogram("hMC_Resid_PHOS_eta_mvsPt",m,pt,alpha) ;
2495             break ;
2496           case 321: //K+
2497           case -321: //K-
2498           case 310:  //K0s
2499           case 130:  //K0L
2500             FillHistogram("hMC_Resid_PHOS_K_mvsPt",m,pt,alpha) ;
2501             break ;
2502           case 211:
2503           case -211:
2504             FillHistogram("hMC_Resid_PHOS_pi_mvsPt",m,pt,alpha) ;
2505             break ;
2506           case -2212:  //pbar
2507           case -2112:  //nbar
2508             FillHistogram("hMC_Resid_PHOS_pbar_mvsPt",m,pt,alpha) ;
2509             break ;
2510           default: //else
2511             FillHistogram("hMC_Resid_PHOS_other_mvsPt",m,pt,alpha) ;
2512             break ;
2513           }
2514         }
2515        }
2516      }
2517    }
2518
2519
2520   if(fEMCALEvent && fConvEvent){
2521     for(Int_t iEMCAL=0; iEMCAL<fEMCALEvent->GetEntriesFast();iEMCAL++){
2522       TLorentzVector * cal = static_cast<TLorentzVector*>(fEMCALEvent->At(iEMCAL)) ;
2523       Int_t iclu=fGammaEMCAL[iEMCAL] ;
2524       AliESDCaloCluster * clu = fESDEvent->GetCaloCluster(iclu);
2525       Int_t iprimEMCAL = clu->GetLabel() ; //# of particle hit EMCAL
2526       for(Int_t iConv = 0; iConv<fConvEvent->GetEntriesFast(); iConv++){
2527         TLorentzVector * cnv = static_cast<TLorentzVector*>(fConvEvent->At(iConv)) ;
2528         if(!cnv->TestBit(kConvOnFly) || cnv->TestBit(kConvEta))
2529           continue;
2530         Int_t iv0=fGammaV0s[iConv] ;
2531         AliESDv0 * v0 = fESDEvent->GetV0(iv0) ;
2532         Int_t iprimNeg = TMath::Abs(fESDEvent->GetTrack(v0->GetNindex())->GetLabel()) ;
2533         Int_t iprimPos = TMath::Abs(fESDEvent->GetTrack(v0->GetPindex())->GetLabel()) ;
2534
2535         //Check if there was a common ancistor
2536         Bool_t found = kFALSE ;
2537         Int_t curEMCAL=iprimEMCAL ;
2538         Int_t commonA=-1 ;
2539         while(!found && curEMCAL>-1){
2540           Int_t curNeg=iprimNeg ;
2541           while(!found && curNeg>-1){
2542             if(curNeg==curEMCAL){
2543               found=kTRUE ;
2544               commonA=curEMCAL ;
2545             }
2546             else{
2547               curNeg=fStack->Particle(curNeg)->GetFirstMother() ;
2548             }
2549           }
2550           curEMCAL=fStack->Particle(curEMCAL)->GetFirstMother() ;
2551         }
2552         found = kFALSE ;
2553         curEMCAL=iprimEMCAL ;
2554         Int_t commonB=-1 ;
2555         while(!found && curEMCAL>-1){
2556           Int_t curPos=iprimPos ;
2557           while(!found && curPos>-1){
2558             if(curPos==curEMCAL){
2559               found=kTRUE ;
2560               commonB=curEMCAL ;
2561             }
2562             else{
2563               curPos=fStack->Particle(curPos)->GetFirstMother() ;
2564             }
2565           }
2566           curEMCAL=fStack->Particle(curEMCAL)->GetFirstMother() ;
2567         }
2568
2569         if(commonA != commonB){
2570            //Strange
2571            AliInfo(Form("CommonA=%d, commonB=%d",commonA,commonB)) ;
2572         }
2573         if(commonA>-1){//There was common particles
2574           Int_t pdg = fStack->Particle(commonA)->GetPdgCode() ;
2575           TLorentzVector pi=*cal + *cnv ;
2576           Double_t m=pi.M() ;
2577           Double_t pt=pi.Pt() ;
2578           Double_t alpha=TMath::Abs(cal->Energy()-cnv->Energy())/(cal->Energy()+cnv->Energy()) ;
2579           switch(pdg){
2580           case 11:
2581           case -11:
2582           case 22: //conversion
2583             FillHistogram("hMC_Resid_EMCAL_Phot_mvsPt",m,pt,alpha) ;
2584             break ;
2585           case 111: //pi0
2586             FillHistogram("hMC_Resid_EMCAL_Pi0_mvsPt",m,pt,alpha) ;
2587             break ;
2588           case 221: //eta
2589             FillHistogram("hMC_Resid_EMCAL_eta_mvsPt",m,pt,alpha) ;
2590             break ;
2591           case 321: //K+
2592           case -321: //K-
2593           case 310:  //K0s
2594           case 130:  //K0L
2595             FillHistogram("hMC_Resid_EMCAL_K_mvsPt",m,pt,alpha) ;
2596             break ;
2597           case 211:
2598           case -211:
2599             FillHistogram("hMC_Resid_EMCAL_pi_mvsPt",m,pt,alpha) ;
2600             break ;
2601           case -2212:  //pbar
2602           case -2112:  //nbar
2603             FillHistogram("hMC_Resid_EMCAL_pbar_mvsPt",m,pt,alpha) ;
2604             break ;
2605           default: //else
2606             FillHistogram("hMC_Resid_EMCAL_other_mvsPt",m,pt,alpha) ;
2607             break ;
2608           }
2609         }
2610        }
2611      }
2612    }
2613
2614
2615    //------------- now photons ----------------
2616    for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++) {
2617      TParticle* particle = (TParticle *)fStack->Particle(iTracks);
2618      if(particle->GetPdgCode() != 22)
2619        continue ;
2620
2621      if(particle->R() >rcut)
2622        continue ;
2623
2624      if(TMath::Abs(particle->Eta())>0.9)
2625        continue ;
2626
2627      Double_t pt = particle->Pt() ;
2628      //Total number of pi0 with creation radius <1 cm
2629      FillHistogram("hMC_CaloConv_phot",pt) ;
2630
2631      Int_t mod ;
2632      Double_t x=0.,z=0. ;
2633      Bool_t hitPHOS = fPHOSgeom->ImpactOnEmc(particle, mod, z,x) ;
2634      Bool_t hitEMCAL= fEMCALgeom->Impact(particle) ;
2635
2636      //Photons in PHOS/EMCAL acceptance
2637      if(hitPHOS)
2638        FillHistogram("hMC_CaloConv_gammaPHOSacc",pt) ;
2639      if(hitEMCAL)
2640        FillHistogram("hMC_CaloConv_gammaEMCALacc",pt) ;
2641
2642      //number of photons converted
2643      Bool_t converted = kFALSE ;
2644      if(particle->GetNDaughters()==2){
2645        TParticle * e1=fStack->Particle(particle->GetFirstDaughter()) ;
2646        TParticle * e2=fStack->Particle(particle->GetLastDaughter()) ;
2647        if(TMath::Abs(e1->GetPdgCode())==11 && TMath::Abs(e2->GetPdgCode())==11){ //conversion
2648          if(e1->R()<180.)
2649            converted = kTRUE ;
2650        }
2651      }
2652      if(converted)
2653        FillHistogram("hMC_CaloConv_gamma_conv",pt) ;
2654
2655      //Converted photons with V0
2656      TLorentzVector pConv ;
2657      Bool_t foundV0=kFALSE ;
2658      for(Int_t iv0=0; iv0<fESDEvent->GetNumberOfV0s();iv0++){
2659        AliESDv0 * v0 = fESDEvent->GetV0(iv0) ;
2660
2661        TParticle * negativeMC = fStack->Particle(TMath::Abs(fESDEvent->GetTrack(v0->GetNindex())->GetLabel()));
2662        TParticle * positiveMC = fStack->Particle(TMath::Abs(fESDEvent->GetTrack(v0->GetPindex())->GetLabel()));
2663
2664        if(negativeMC->GetMother(0) != positiveMC->GetMother(0))
2665          continue ;
2666
2667        if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){
2668          continue;
2669        }
2670        if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){
2671          continue;
2672        }
2673
2674        TParticle * v0Gamma = fStack->Particle(negativeMC->GetMother(0));
2675        Bool_t same = (v0Gamma == particle) ;
2676        TParticle * tmp = v0Gamma ;
2677        while(!same && tmp->GetFirstMother()>=0){
2678          tmp = fStack->Particle(tmp->GetFirstMother());
2679          same = (tmp == particle) ;
2680        }
2681        if(same){
2682          foundV0 = kTRUE ;
2683          const AliExternalTrackParam * paramPos = v0->GetParamP() ;
2684          const AliExternalTrackParam * paramNeg = v0->GetParamN() ;
2685          AliKFParticle negKF(*paramNeg,11);
2686          AliKFParticle posKF(*paramPos,-11);
2687          pConv.SetXYZM(negKF.Px()+posKF.Px(),negKF.Py()+posKF.Py(),negKF.Pz()+negKF.Pz(),0.) ;
2688          break ;
2689        }
2690      }
2691      if(foundV0){
2692        FillHistogram("hMC_CaloConv_gamma_v0",pt) ;
2693        FillHistogram("hMC_CaloConv_gamma_v0_devsE",(particle->Energy()-pConv.E())/particle->Energy(),particle->Energy()) ;
2694      }
2695
2696       //Registered in PHOS/EMCAL
2697      Bool_t cluInPHOS = kFALSE,cluInEMCAL=kFALSE ;
2698      TLorentzVector pCalo ;
2699      Bool_t dist1=kFALSE, dist2=kFALSE ;
2700      for (Int_t i=0; i<fESDEvent->GetNumberOfCaloClusters(); i++) {
2701        AliESDCaloCluster * clu = fESDEvent->GetCaloCluster(i);
2702        Int_t iprim = clu->GetLabel() ; //# of particle hit PHOS/EMCAL
2703        Bool_t matched = kFALSE ;
2704        while(iprim>=0 ) {
2705          if(iprim==iTracks){
2706            matched=kTRUE ;
2707            break ;
2708          }
2709          else{
2710            iprim=fStack->Particle(iprim)->GetFirstMother() ;
2711          }
2712        }
2713        if(!matched)
2714          continue ;
2715        if(clu->IsPHOS() && hitPHOS){
2716          cluInPHOS=kTRUE ;
2717          clu->GetMomentum(pCalo ,vtx);
2718          if(clu->GetDistanceToBadChannel()<fBadDistCutPHOS)
2719            dist1=kTRUE ;
2720          if(clu->GetDistanceToBadChannel()<2.*fBadDistCutPHOS)
2721            dist2=kTRUE ;
2722          break ;
2723        }
2724        if(!clu->IsPHOS() && hitEMCAL){
2725          cluInEMCAL=kTRUE ;
2726          clu->GetMomentum(pCalo ,vtx);
2727          if(clu->GetDistanceToBadChannel()<fBadDistCutEMCAL)
2728            dist1=kTRUE ;
2729          if(clu->GetDistanceToBadChannel()<2.*fBadDistCutEMCAL)
2730            dist2=kTRUE ;
2731          break ;
2732        }
2733      }
2734
2735      if(cluInPHOS){
2736        FillHistogram("hMC_CaloConv_gamma_PHOSclu",pt) ;
2737        if(!dist1)
2738          FillHistogram("hMC_CaloConv_gamma_PHOSclu_dist1",pt) ;
2739        if(!dist2)
2740          FillHistogram("hMC_CaloConv_gamma_PHOSclu_dist2",pt) ;
2741        FillHistogram("hMC_CaloConv_gamma_PHOSclu_recE",pCalo.E()) ;
2742        FillHistogram("hMC_CaloConv_gamma_PHOSclu_devsE",(particle->Energy()-pCalo.E())/particle->Energy(),particle->Energy()) ;
2743      }
2744      if(cluInEMCAL){
2745        FillHistogram("hMC_CaloConv_gamma_EMCALclu",pt) ;
2746        if(!dist1)
2747          FillHistogram("hMC_CaloConv_gamma_EMCALclu_dist1",pt) ;
2748        if(!dist2)
2749          FillHistogram("hMC_CaloConv_gamma_EMCALclu_dist2",pt) ;
2750        FillHistogram("hMC_CaloConv_gamma_EMCALclu_recE",pCalo.E()) ;
2751        FillHistogram("hMC_CaloConv_gamma_EMCALclu_devsE",(particle->Energy()-pCalo.E())/particle->Energy(),particle->Energy()) ;
2752      }
2753   }
2754 }
2755 //_____________________________________________________________________________
2756 void AliAnalysisTaskCaloConv::FillHistogram(const char * key,Double_t x)const{
2757   //FillHistogram
2758   TH1F * tmp = dynamic_cast<TH1F*>(fOutputContainer->FindObject(key)) ;
2759   if(!tmp){
2760     AliInfo(Form("can not find histogram <%s> ",key)) ;
2761     return ;
2762   }
2763   tmp->Fill(x) ;
2764 }
2765 //_____________________________________________________________________________
2766 void AliAnalysisTaskCaloConv::FillHistogram(const char * key,Double_t x,Double_t y)const{
2767   //FillHistogram
2768   TObject * tmp = fOutputContainer->FindObject(key) ;
2769   if(!tmp){
2770     AliInfo(Form("can not find histogram <%s> ",key)) ;
2771     return ;
2772   }
2773   if(tmp->IsA() == TClass::GetClass("TH1F")){
2774     ((TH1F*)tmp)->Fill(x,y) ;
2775     return ;
2776   }
2777   if(tmp->IsA() == TClass::GetClass("TH2F")){
2778     ((TH2F*)tmp)->Fill(x,y) ;
2779     return ;
2780   }
2781   AliError(Form("Calling FillHistogram with 2 parameters for histo <%s> of type %s",key,tmp->IsA()->GetName())) ;
2782 }
2783
2784 //_____________________________________________________________________________
2785 void AliAnalysisTaskCaloConv::FillHistogram(const char * key,Double_t x,Double_t y, Double_t z) const{
2786   //Fills 1D histograms with key
2787   TObject * tmp = fOutputContainer->FindObject(key) ;
2788   if(!tmp){
2789     AliInfo(Form("can not find histogram <%s> ",key)) ;
2790     return ;
2791   }
2792   if(tmp->IsA() == TClass::GetClass("TH2F")){
2793     ((TH2F*)tmp)->Fill(x,y,z) ;
2794     return ;
2795   }
2796   if(tmp->IsA() == TClass::GetClass("TH3F")){
2797     ((TH3F*)tmp)->Fill(x,y,z) ;
2798     return ;
2799   }
2800 }
2801 //______________________________________________________________________________
2802 Double_t AliAnalysisTaskCaloConv::PlanarityAngle(const AliExternalTrackParam * paramPos,const AliExternalTrackParam * paramNeg)const {
2803   //calculate angle between e+e- plain and perpendicular to MF
2804   //We need sign of MagField to calculate orienation
2805
2806   TVector3 u(paramPos->Px()+paramNeg->Px(),paramPos->Py()+paramNeg->Py(),paramPos->Pz()+paramNeg->Pz()) ;
2807   u.Unit() ;
2808   TVector3 vPos(paramPos->Px(),paramPos->Py(),paramPos->Pz()) ;
2809   TVector3 vNeg(paramNeg->Px(),paramNeg->Py(),paramNeg->Pz()) ;
2810   TVector3 v=vPos.Cross(vNeg) ;
2811   TVector3 w = u.Cross(v);
2812   TVector3 z(0,0,1.);
2813   TVector3 ua=u.Cross(z);
2814   Double_t wa = w.Angle(ua);
2815   Double_t mfield=fESDEvent->GetMagneticField() ;
2816   if(mfield>0.)
2817     return wa;      //normal field
2818   else
2819     return TMath::Pi()-wa ; //reverse field
2820
2821 }
2822 //______________________________________________________________________________
2823 Bool_t AliAnalysisTaskCaloConv::IsGoodChannel(const char * det, Int_t mod, Int_t ix, Int_t iz){
2824 //Check if this channel belogs to the good ones
2825
2826   if(strcmp(det,"PHOS")==0){
2827     if(mod>5 || mod<1){
2828       AliError(Form("No bad map for PHOS module %d ",mod)) ;
2829       return kTRUE ;
2830     }
2831     if(!fPHOSBadMap[mod]){
2832       AliError(Form("No Bad map for PHOS module %d",mod)) ;
2833       return kTRUE ;
2834     }
2835     if(fPHOSBadMap[mod]->GetBinContent(ix,iz)>0)
2836       return kFALSE ;
2837     else
2838       return kTRUE ;
2839   }
2840   else{
2841     if(strcmp(det,"EMCAL")==0){
2842       if(mod>9 || mod<0){
2843         AliError(Form("No bad map for EMCAL module %d ",mod)) ;
2844         return kTRUE ;
2845       }
2846       if(!fEMCALBadMap[mod]){
2847         AliError(Form("No bad map for EMCAL module %d ",mod)) ;
2848         return kTRUE ;
2849       }
2850       if(fEMCALBadMap[mod]->GetBinContent(ix,iz)>0)
2851         return kFALSE ;
2852       else
2853         return kTRUE ;
2854     }
2855     else{
2856       AliError(Form("Can not find bad channels for detector %s ",det)) ;
2857     }
2858   }
2859
2860   return kTRUE ;
2861 }
2862 //______________________________________________________________________________
2863 void AliAnalysisTaskCaloConv::Recalibrate(Double_t &m, Double_t &pt, const TLorentzVector *calo, const TLorentzVector * conv, Int_t iw, Int_t in) {
2864   //Apply decalibration and non-linearity
2865   TLorentzVector calo2(*calo) ;
2866   Double_t en=calo2.E() ;
2867
2868   Double_t sigma=0.06+0.005*iw ; //additional smearing
2869   //Nonlinearity
2870   Double_t a=0.02*(in%6-2.5) ;
2871   Double_t b=0.5+1.*((Int_t)in/6) ;
2872   Double_t enNew=1.-a*TMath::Exp(-en/b) ;
2873   Double_t corr=gRandom->Gaus(enNew,sigma) ;
2874   calo2*=corr ;
2875
2876   m=(calo2+ *conv).M() ;
2877   pt=(calo2+ *conv).Pt() ;
2878
2879 }
2880 //______________________________________________________________________________
2881 void AliAnalysisTaskCaloConv::RecalibrateEMCAL(Double_t &m, Double_t &pt, const TLorentzVector *calo, const TLorentzVector * conv, Int_t iw, Int_t in) {
2882   //Apply decalibration and non-linearity
2883   TLorentzVector calo2(*calo) ;
2884   Double_t en=calo2.E() ;
2885
2886   Double_t sigma=0.04+0.005*iw ; //additional smearing
2887   //Nonlinearity
2888   Double_t a=0.02*(in%6-2.5) ;
2889   Double_t b=0.25+0.5*((Int_t)in/6) ;
2890   Double_t enNew=1.-a*TMath::Exp(-en/b) ;
2891   Double_t corr=gRandom->Gaus(enNew,sigma) ;
2892   calo2*=corr ;
2893
2894   m=(calo2+ *conv).M() ;
2895   pt=(calo2+ *conv).Pt() ;
2896
2897 }
2898 //______________________________________________________________________________
2899 void AliAnalysisTaskCaloConv::RecalibrateConvPHOS(Double_t &m, Double_t &pt, const TLorentzVector *calo, const TLorentzVector * conv, Int_t iw, Int_t in) {
2900   //Apply decalibration and non-linearity
2901
2902   //First default PHOS smearing
2903   TLorentzVector calo2(*calo) ;
2904   Double_t en=calo2.E() ;
2905
2906   Double_t sigma=0.065 ; //additional smearing
2907   //Nonlinearity
2908   Double_t a=0.15 ;
2909   Double_t b=0.45 ;
2910   Double_t enNew=1.+a*TMath::Exp(-en/b) ;
2911   Double_t corr=gRandom->Gaus(enNew,sigma) ;
2912   calo2*=corr ;
2913
2914   //Now conversion photon
2915   TLorentzVector conv2(*conv) ;
2916   //linear offset in z:
2917   Double_t eta=conv2.Eta() ;
2918   Double_t c=1.e-3*iw ;
2919   eta+= c *TMath::Sign(0.9-TMath::Abs(eta),eta) ;
2920
2921   //Smear energy and add nonlinearity
2922   //Nonlinearity
2923   Double_t enConv=conv2.E() ;
2924   Double_t ac=0.02*(in%5) ;
2925   Double_t bc=0.25+0.5*((Int_t)in/5) ;
2926   Double_t enNewc=1.+ac*TMath::Exp(-enConv/bc) ;
2927   corr=gRandom->Gaus(enNewc,0.01) ;
2928   Double_t ptc=conv2.Pt()*corr ;
2929   conv2.SetPtEtaPhiM(ptc,eta,conv2.Phi(),0.) ;
2930
2931   m =(calo2 + conv2).M() ;
2932   pt=(calo2 + conv2).Pt() ;
2933
2934 }
2935 //______________________________________________________________________________
2936 void AliAnalysisTaskCaloConv::GetArmenterosQtAlfa(AliKFParticle* positiveKFParticle, AliKFParticle * negativeKFParticle, AliKFParticle * gammaKFCandidate, Double_t armenterosQtAlfa[2] ){
2937   //see header file for documentation
2938
2939   TVector3 momentumVectorPositiveKF(positiveKFParticle->GetPx(),positiveKFParticle->GetPy(),positiveKFParticle->GetPz());
2940   TVector3 momentumVectorNegativeKF(negativeKFParticle->GetPx(),negativeKFParticle->GetPy(),negativeKFParticle->GetPz());
2941   TVector3 vecV0(gammaKFCandidate->GetPx(),gammaKFCandidate->GetPy(),gammaKFCandidate->GetPz());
2942
2943   Float_t thetaV0pos=TMath::ACos(( momentumVectorPositiveKF* vecV0)/(momentumVectorPositiveKF.Mag() * vecV0.Mag()));
2944   Float_t thetaV0neg=TMath::ACos(( momentumVectorNegativeKF* vecV0)/(momentumVectorNegativeKF.Mag() * vecV0.Mag()));
2945
2946   Float_t alfa =((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)-(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg))/
2947     ((momentumVectorPositiveKF.Mag())*TMath::Cos(thetaV0pos)+(momentumVectorNegativeKF.Mag())*TMath::Cos(thetaV0neg)) ;
2948
2949
2950   Float_t qt = momentumVectorPositiveKF.Mag()*TMath::Sin(thetaV0pos);
2951
2952   armenterosQtAlfa[0]=qt;
2953   armenterosQtAlfa[1]=alfa;
2954
2955 }
2956
2957
2958
2959
2960
2961
2962
2963