PWG4/AliAnalysisTaskGammaConversion.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: Ana Marin, Kathrin Koch, Kenneth Aamodt                        *
   5  * Version 1.0                                                            *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 // root
  17 #include <TROOT.h>
  18 #include <TSystem.h>
  19 #include <TInterpreter.h>
  20 #include <TChain.h>
  21 #include <TFile.h>
  22 #include <Riostream.h>
  23
  24 // analysis
  25 //#include "AliAnalysisTaskSE.h"
  26 #include "AliAnalysisTaskGammaConversion.h"
  27 #include "AliAnalysisManager.h"
  28 #include "AliESDInputHandler.h"
  29 #include "AliMCEventHandler.h"
  30 #include "AliMCEvent.h"
  31 #include "AliESDEvent.h"
  32 #include "AliAODEvent.h"
  33 #include "AliAODHandler.h"
  34 #include "AliStack.h"
  35 #include "AliLog.h"
  36 //#include "TLorentzVector.h"
  37 #include "AliKFVertex.h"
  38
  39 ClassImp(AliAnalysisTaskGammaConversion)
  40
  41
  42 AliAnalysisTaskGammaConversion::AliAnalysisTaskGammaConversion():
  43   AliAnalysisTaskSE(),
  44   fV0Reader(NULL),
  45   fStack(NULL),
  46   fOutputContainer(NULL),
  47   fHistograms(NULL),
  48   fDoMCTruth(kFALSE),
  49   fMCAllGammas(),
  50   fMCPi0s(),
  51   fMCEtas(),
  52   fMCGammaChi_c(),
  53   fKFReconstructedGammas(),
  54   fElectronMass(-1),
  55   fGammaMass(-1),
  56   fPi0Mass(-1),
  57   fEtaMass(-1),
  58   fGammaWidth(-1),
  59   fPi0Width(-1),
  60   fEtaWidth(-1),
  61   fCalculateBackground(kFALSE)
  62 {
  63   // Default constructor
  64   // Common I/O in slot 0
  65   DefineInput (0, TChain::Class());
  66   DefineOutput(0, TTree::Class());
  67
  68   // Your private output
  69   DefineOutput(1, TList::Class());
  70 }
  71
  72 AliAnalysisTaskGammaConversion::AliAnalysisTaskGammaConversion(const char* name):
  73   AliAnalysisTaskSE(name),
  74   fV0Reader(NULL),
  75   fStack(NULL),
  76   fOutputContainer(0x0),
  77   fHistograms(NULL),
  78   fDoMCTruth(kFALSE),
  79   fMCAllGammas(),
  80   fMCPi0s(),
  81   fMCEtas(),
  82   fMCGammaChi_c(),
  83   fKFReconstructedGammas(),
  84   fElectronMass(-1),
  85   fGammaMass(-1),
  86   fPi0Mass(-1),
  87   fEtaMass(-1),
  88   fGammaWidth(-1),
  89   fPi0Width(-1),
  90   fEtaWidth(-1),
  91   fCalculateBackground(kFALSE)
  92 {
  93   // Common I/O in slot 0
  94   DefineInput (0, TChain::Class());
  95   DefineOutput(0, TTree::Class());
  96
  97   // Your private output
  98   DefineOutput(1, TList::Class());
  99 }
 100
 101 AliAnalysisTaskGammaConversion::~AliAnalysisTaskGammaConversion()
 102 {
 103   // Remove all pointers
 104
 105   if(fOutputContainer){
 106     fOutputContainer->Clear() ;
 107     delete fOutputContainer ;
 108   }
 109   if(fHistograms){
 110     delete fHistograms;
 111   }
 112   if(fV0Reader){
 113     delete fV0Reader;
 114   }
 115 }
 116
 117
 118 void AliAnalysisTaskGammaConversion::Init()
 119 {
 120   // Initialization
 121 }
 122
 123
 124 void AliAnalysisTaskGammaConversion::Exec(Option_t */*option*/)
 125 {
 126   // Execute analysis for current event
 127   //
 128
 129   ConnectInputData("");
 130
 131   //clear vectors
 132   fMCAllGammas.clear();
 133   fMCPi0s.clear();
 134   fMCEtas.clear();
 135   fMCGammaChi_c.clear();
 136   /*
 137   for(UInt_t ikf=0;ikf<fKFReconstructedGammas.size();ikf++){
 138     delete fKFReconstructedGammas[ikf];
 139     fKFReconstructedGammas[ikf]=NULL;
 140   }
 141   */
 142   fKFReconstructedGammas.clear();
 143
 144   //Clear the data in the v0Reader
 145   fV0Reader->UpdateEventByEventData();
 146   // Process the v0 information
 147   if(fDoMCTruth){
 148     ProcessMCData();
 149   }
 150
 151   // Process the v0 information
 152   ProcessV0s();
 153
 154   if(fCalculateBackground){//calculate background if flag is set
 155     CalculateBackground();
 156   }
 157
 158   // Process reconstructed gammas
 159   ProcessGammasForNeutralMesonAnalysis();
 160
 161   PostData(1, fOutputContainer);
 162
 163 }
 164
 165 void AliAnalysisTaskGammaConversion::ConnectInputData(Option_t */*option*/){
 166
 167   if(fV0Reader == NULL){
 168     // Write warning here cuts and so on are default if this ever happens
 169   }
 170   fV0Reader->Initialize();
 171 }
 172
 173 void AliAnalysisTaskGammaConversion::ProcessMCData(){
 174
 175   fStack = fV0Reader->GetMCStack();
 176
 177   for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++) {
 178     TParticle* particle = (TParticle *)fStack->Particle(iTracks);
 179
 180     if (!particle) {
 181       //print warning here
 182       continue;
 183     }
 184
 185     if(particle->Pt()<fV0Reader->GetPtCut()){
 186       continue;
 187     }
 188
 189     if(TMath::Abs(particle->Eta())> fV0Reader->GetEtaCut()){
 190       continue;
 191     }
 192
 193     if(particle->R()>fV0Reader->GetMaxRCut()){ // cuts on distance from collision point
 194       continue;
 195     }
 196
 197     Double_t tmpPhi=particle->Phi();
 198     if(particle->Phi()> TMath::Pi()){
 199       tmpPhi = particle->Phi()-(2*TMath::Pi());
 200     }
 201
 202
 203     //process the gammas
 204     if (particle->GetPdgCode()== 22){
 205       fMCAllGammas.push_back(particle);
 206       if(particle->GetMother(0)>-1){ //Means we have a mother
 207         if( fStack->Particle(particle->GetMother(0))->GetPdgCode() != 22 ){//Checks for a non gamma mother.
 208           if(fHistograms->fMC_Gamma_Energy){fHistograms->fMC_Gamma_Energy->Fill(particle->Energy());}
 209           if(fHistograms->fMC_Gamma_Pt){fHistograms->fMC_Gamma_Pt->Fill(particle->Pt());}
 210           if(fHistograms->fMC_Gamma_Eta){fHistograms->fMC_Gamma_Eta->Fill(particle->Eta());}
 211
 212           /*      Double_t tmpPhi=particle->Phi();
 213           if(particle->Phi()> TMath::Pi()){
 214             tmpPhi = particle->Phi()-(2*TMath::Pi());
 215             }*/
 216           if(fHistograms->fMC_Gamma_Phi){fHistograms->fMC_Gamma_Phi->Fill(tmpPhi);}
 217
 218           //adding the conversion points from all gammas with e+e- daughters
 219           if(particle->GetNDaughters() == 2){
 220             TParticle* daughter0 = (TParticle*)fStack->Particle(particle->GetFirstDaughter());
 221             TParticle* daughter1 = (TParticle*)fStack->Particle(particle->GetLastDaughter());
 222
 223             if(daughter0->R()>fV0Reader->GetMaxRCut() || daughter1->R()>fV0Reader->GetMaxRCut()){
 224               continue;
 225             }
 226
 227             if(daughter0->GetPdgCode() == -11 && daughter1->GetPdgCode() == 11 ||
 228                daughter0->GetPdgCode() == 11 && daughter1->GetPdgCode() == -11){
 229
 230               // begin Mapping
 231               Int_t rBin    = fHistograms->GetRBin(daughter0->R());
 232               Int_t phiBin  = fHistograms->GetPhiBin(daughter0->Phi());
 233
 234               if(fHistograms->fMC_Mapping[phiBin][rBin] != NULL){fHistograms->fMC_Mapping[phiBin][rBin]->Fill(daughter0->Vz(), particle->Eta());}
 235               if(fHistograms->fMC_Mapping_Phi[phiBin] != NULL){fHistograms->fMC_Mapping_Phi[phiBin]->Fill(daughter0->Vz(), particle->Eta());}
 236               if(fHistograms->fMC_Mapping_R[rBin] != NULL){fHistograms->fMC_Mapping_R[rBin]->Fill(daughter0->Vz(), particle->Eta());}
 237               //end mapping
 238
 239
 240               if(fHistograms->fMC_EP_R){fHistograms->fMC_EP_R->Fill(daughter0->R());}
 241               if(fHistograms->fMC_EP_Z_R){fHistograms->fMC_EP_Z_R->Fill(daughter0->Vz(),daughter0->R());}
 242               if(fHistograms->fMC_EP_X_Y){fHistograms->fMC_EP_X_Y->Fill(daughter0->Vx(),daughter0->Vy());}
 243               if(fHistograms->fMC_EP_OpeningAngle){fHistograms->fMC_EP_OpeningAngle->Fill(GetMCOpeningAngle(daughter0, daughter1));}
 244
 245             }
 246           }
 247         }
 248         if( fStack->Particle(particle->GetMother(0))->GetPdgCode()==10441 ||//chi_c0
 249             fStack->Particle(particle->GetMother(0))->GetPdgCode()==20443 ||//psi_2S
 250             fStack->Particle(particle->GetMother(0))->GetPdgCode()==445  //chi_c2
 251         ){
 252           fMCGammaChi_c.push_back(particle);
 253          }
 254       }
 255       else{//means we have a primary particle
 256         if(fHistograms->fMC_DirectGamma_Energy){fHistograms->fMC_DirectGamma_Energy->Fill(particle->Energy());}
 257         if(fHistograms->fMC_DirectGamma_Pt){fHistograms->fMC_DirectGamma_Pt->Fill(particle->Pt());}
 258         if(fHistograms->fMC_DirectGamma_Eta){fHistograms->fMC_DirectGamma_Eta->Fill(particle->Eta());}
 259         /*
 260         Double_t tmpPhi=particle->Phi();
 261         if(particle->Phi()> TMath::Pi()){
 262           tmpPhi = particle->Phi()-(2*TMath::Pi());
 263         }
 264         */
 265         if(fHistograms->fMC_DirectGamma_Phi){fHistograms->fMC_DirectGamma_Phi->Fill(tmpPhi);}
 266
 267         //adding the conversion points from all gammas with e+e- daughters
 268         if(particle->GetNDaughters() == 2){
 269           TParticle* daughter0 = (TParticle*)fStack->Particle(particle->GetFirstDaughter());
 270           TParticle* daughter1 = (TParticle*)fStack->Particle(particle->GetLastDaughter());
 271           if(daughter0->GetPdgCode() == -11 && daughter1->GetPdgCode() == 11 ||
 272              daughter0->GetPdgCode() == 11 && daughter1->GetPdgCode() == -11){
 273
 274             if(fHistograms->fMC_EP_R){fHistograms->fMC_EP_R->Fill(daughter0->R());}
 275             if(fHistograms->fMC_EP_Z_R){fHistograms->fMC_EP_Z_R->Fill(daughter0->Vz(),daughter0->R());}
 276             if(fHistograms->fMC_EP_X_Y){fHistograms->fMC_EP_X_Y->Fill(daughter0->Vx(),daughter0->Vy());}
 277             if(fHistograms->fMC_EP_OpeningAngle){fHistograms->fMC_EP_OpeningAngle->Fill(GetMCOpeningAngle(daughter0,daughter1));}
 278
 279           }
 280         }
 281
 282       }
 283     }
 284     else if (TMath::Abs(particle->GetPdgCode())== 11){ // Means we have an electron or a positron
 285       if(particle->GetMother(0)>-1){ // means we have a mother
 286         if( fStack->Particle(particle->GetMother(0))->GetPdgCode()==22 ){ // Means we have a gamma mother
 287           if(particle->GetPdgCode() == 11){//electron
 288             if(fHistograms->fMC_E_Energy){fHistograms->fMC_E_Energy->Fill(particle->Energy());}
 289             if(fHistograms->fMC_E_Pt){fHistograms->fMC_E_Pt->Fill(particle->Pt());}
 290             if(fHistograms->fMC_E_Eta){fHistograms->fMC_E_Eta->Fill(particle->Eta());}
 291             /*      Double_t tmpPhi=particle->Phi();
 292             if(particle->Phi()> TMath::Pi()){
 293               tmpPhi = particle->Phi()-(2*TMath::Pi());
 294               }*/
 295             if(fHistograms->fMC_E_Phi){fHistograms->fMC_E_Phi->Fill(tmpPhi);}
 296           }
 297           if(particle->GetPdgCode() == -11){//positron
 298             if(fHistograms->fMC_P_Energy){fHistograms->fMC_P_Energy->Fill(particle->Energy());}
 299             if(fHistograms->fMC_P_Pt){fHistograms->fMC_P_Pt->Fill(particle->Pt());}
 300             if(fHistograms->fMC_P_Eta){fHistograms->fMC_P_Eta->Fill(particle->Eta());}
 301             /*
 302             Double_t tmpPhi=particle->Phi();
 303             if(particle->Phi()> TMath::Pi()){
 304               tmpPhi = particle->Phi()-(2*TMath::Pi());
 305             }
 306             */
 307             if(fHistograms->fMC_P_Phi){fHistograms->fMC_P_Phi->Fill(tmpPhi);}
 308           }
 309         }
 310       }
 311     }
 312     else if(particle->GetNDaughters() == 2){
 313
 314       TParticle* daughter0 = (TParticle*)fStack->Particle(particle->GetFirstDaughter());
 315       TParticle* daughter1 = (TParticle*)fStack->Particle(particle->GetLastDaughter());
 316       if(daughter0->GetPdgCode() == 22 && daughter1->GetPdgCode() == 22){//check for gamma gamma daughters
 317
 318         if(particle->GetPdgCode()==111){//Pi0
 319           if( iTracks >= fStack->GetNprimary()){
 320
 321             if(fHistograms->fMC_Pi0Secondaries_Eta){fHistograms->fMC_Pi0Secondaries_Eta->Fill(particle->Eta());}
 322             /*
 323               Double_t tmpPhi=particle->Phi();
 324               if(particle->Phi()> TMath::Pi()){
 325               tmpPhi = particle->Phi()-(2*TMath::Pi());
 326               }
 327             */
 328             if(fHistograms->fMC_Pi0Secondaries_Phi){fHistograms->fMC_Pi0Secondaries_Phi->Fill(tmpPhi);}
 329             if(fHistograms->fMC_Pi0Secondaries_Pt){fHistograms->fMC_Pi0Secondaries_Pt->Fill(particle->Pt());}
 330             if(fHistograms->fMC_Pi0Secondaries_Energy){fHistograms->fMC_Pi0Secondaries_Energy->Fill(particle->Energy());}
 331             if(fHistograms->fMC_Pi0Secondaries_R){fHistograms->fMC_Pi0Secondaries_R->Fill(particle->R());}
 332             if(fHistograms->fMC_Pi0Secondaries_Z_R){fHistograms->fMC_Pi0Secondaries_Z_R->Fill(particle->Vz(),particle->R());}
 333             if(fHistograms->fMC_Pi0Secondaries_OpeningAngleGamma){fHistograms->fMC_Pi0Secondaries_OpeningAngleGamma->Fill(GetMCOpeningAngle(daughter0,daughter1));}
 334             if(fHistograms->fMC_Pi0Secondaries_X_Y){fHistograms->fMC_Pi0Secondaries_X_Y->Fill(particle->Vx(),particle->Vy());}//only fill from one daughter to avoid multiple filling
 335           }
 336           else{
 337             if(fHistograms->fMC_Pi0_Eta){fHistograms->fMC_Pi0_Eta->Fill(particle->Eta());}
 338             /*
 339             Double_t tmpPhi=particle->Phi();
 340             if(particle->Phi()> TMath::Pi()){
 341               tmpPhi = particle->Phi()-(2*TMath::Pi());
 342             }
 343             */
 344             if(fHistograms->fMC_Pi0_Phi){fHistograms->fMC_Pi0_Phi->Fill(tmpPhi);}
 345             if(fHistograms->fMC_Pi0_Pt){fHistograms->fMC_Pi0_Pt->Fill(particle->Pt());}
 346             if(fHistograms->fMC_Pi0_Energy){fHistograms->fMC_Pi0_Energy->Fill(particle->Energy());}
 347             if(fHistograms->fMC_Pi0_R){fHistograms->fMC_Pi0_R->Fill(particle->R());}
 348             if(fHistograms->fMC_Pi0_Z_R){fHistograms->fMC_Pi0_Z_R->Fill(particle->Vz(),particle->R());}
 349             if(fHistograms->fMC_Pi0_OpeningAngleGamma){fHistograms->fMC_Pi0_OpeningAngleGamma->Fill(GetMCOpeningAngle(daughter0,daughter1));}
 350             if(fHistograms->fMC_Pi0_X_Y){fHistograms->fMC_Pi0_X_Y->Fill(particle->Vx(),particle->Vy());}//only fill from one daughter to avoid multiple filling
 351           }
 352         }
 353         else if(particle->GetPdgCode()==221){//Eta
 354           if(fHistograms->fMC_Eta_Eta){fHistograms->fMC_Eta_Eta->Fill(particle->Eta());}
 355           /*
 356             Double_t tmpPhi=particle->Phi();
 357           if(particle->Phi()> TMath::Pi()){
 358             tmpPhi = particle->Phi()-(2*TMath::Pi());
 359           }
 360           */
 361           if(fHistograms->fMC_Eta_Phi){fHistograms->fMC_Eta_Phi->Fill(tmpPhi);}
 362           if(fHistograms->fMC_Eta_Pt){fHistograms->fMC_Eta_Pt->Fill(particle->Pt());}
 363           if(fHistograms->fMC_Eta_Energy){fHistograms->fMC_Eta_Energy->Fill(particle->Energy());}
 364           if(fHistograms->fMC_Eta_R){fHistograms->fMC_Eta_R->Fill(particle->R());}
 365           if(fHistograms->fMC_Eta_Z_R){fHistograms->fMC_Eta_Z_R->Fill(particle->Vz(),particle->R());}
 366           if(fHistograms->fMC_Eta_OpeningAngleGamma){fHistograms->fMC_Eta_OpeningAngleGamma->Fill(GetMCOpeningAngle(daughter0,daughter1));}
 367           if(fHistograms->fMC_Eta_X_Y){fHistograms->fMC_Eta_X_Y->Fill(particle->Vx(),particle->Vy());}//only fill from one daughter to avoid multiple filling
 368         }
 369
 370         //the match data should be filled no matter which mother the gamma-gamma comes from
 371         if(fHistograms->fMC_Match_Gamma_R){fHistograms->fMC_Match_Gamma_R->Fill(particle->R());}
 372         if(fHistograms->fMC_Match_Gamma_Z_R){fHistograms->fMC_Match_Gamma_Z_R->Fill(particle->Vz(),particle->R());}
 373         if(fHistograms->fMC_Match_Gamma_X_Y){fHistograms->fMC_Match_Gamma_X_Y->Fill(particle->Vx(),particle->Vy());}
 374         if(fHistograms->fMC_Match_Gamma_Mass){fHistograms->fMC_Match_Gamma_Mass->Fill(particle->GetCalcMass());}
 375         if(fHistograms->fMC_Match_Gamma_OpeningAngle){fHistograms->fMC_Match_Gamma_OpeningAngle->Fill(GetMCOpeningAngle(daughter0,daughter1));}
 376         if(fHistograms->fMC_Match_Gamma_Energy){fHistograms->fMC_Match_Gamma_Energy->Fill(particle->Energy());}
 377         if(fHistograms->fMC_Match_Gamma_Pt){fHistograms->fMC_Match_Gamma_Pt->Fill(particle->Pt());}
 378         if(fHistograms->fMC_Match_Gamma_Eta){fHistograms->fMC_Match_Gamma_Eta->Fill(particle->Eta());}
 379         /*
 380         Double_t tmpPhi=particle->Phi();
 381         if(particle->Phi()> TMath::Pi()){
 382           tmpPhi = particle->Phi()-(2*TMath::Pi());
 383         }
 384         */
 385         if(fHistograms->fMC_Match_Gamma_Phi){fHistograms->fMC_Match_Gamma_Phi->Fill(tmpPhi);}
 386       }
 387     }
 388   }
 389 }
 390
 391 void AliAnalysisTaskGammaConversion::ProcessV0s(){
 392   Int_t nSurvivingV0s=0;
 393   while(fV0Reader->NextV0()){
 394     nSurvivingV0s++;
 395     //-------------------------- filling v0 information -------------------------------------
 396     if(fHistograms->fESD_EP_OpeningAngle){fHistograms->fESD_EP_OpeningAngle->Fill(fV0Reader->GetOpeningAngle());}
 397     if(fHistograms->fESD_EP_R){fHistograms->fESD_EP_R->Fill(fV0Reader->GetXYRadius());}
 398     if(fHistograms->fESD_EP_Z_R){fHistograms->fESD_EP_Z_R->Fill(fV0Reader->GetZ(),fV0Reader->GetXYRadius());}
 399     if(fHistograms->fESD_EP_X_Y){fHistograms->fESD_EP_X_Y->Fill(fV0Reader->GetX(),fV0Reader->GetY());}
 400
 401
 402     if(fHistograms->fESD_E_Energy){fHistograms->fESD_E_Energy->Fill(fV0Reader->GetNegativeTrackEnergy());}
 403     if(fHistograms->fESD_E_Pt){fHistograms->fESD_E_Pt->Fill(fV0Reader->GetNegativeTrackPt());}
 404     if(fHistograms->fESD_E_Eta){fHistograms->fESD_E_Eta->Fill(fV0Reader->GetNegativeTrackEta());}
 405     if(fHistograms->fESD_E_Phi){fHistograms->fESD_E_Phi->Fill(fV0Reader->GetNegativeTrackPhi());}
 406
 407     if(fHistograms->fESD_P_Energy){fHistograms->fESD_P_Energy->Fill(fV0Reader->GetPositiveTrackEnergy());}
 408     if(fHistograms->fESD_P_Pt){fHistograms->fESD_P_Pt->Fill(fV0Reader->GetPositiveTrackPt());}
 409     if(fHistograms->fESD_P_Eta){fHistograms->fESD_P_Eta->Fill(fV0Reader->GetPositiveTrackEta());}
 410     if(fHistograms->fESD_P_Phi){fHistograms->fESD_P_Phi->Fill(fV0Reader->GetPositiveTrackPhi());}
 411
 412     if(fHistograms->fESD_Gamma_Energy){fHistograms->fESD_Gamma_Energy->Fill(fV0Reader->GetMotherCandidateEnergy());}
 413     if(fHistograms->fESD_Gamma_Pt){fHistograms->fESD_Gamma_Pt->Fill(fV0Reader->GetMotherCandidatePt());}
 414     if(fHistograms->fESD_Gamma_Eta){fHistograms->fESD_Gamma_Eta->Fill(fV0Reader->GetMotherCandidateEta());}
 415     if(fHistograms->fESD_Gamma_Phi){fHistograms->fESD_Gamma_Phi->Fill(fV0Reader->GetMotherCandidatePhi());}
 416
 417
 418     // begin mapping
 419     Int_t rBin    = fHistograms->GetRBin(fV0Reader->GetXYRadius());
 420     Int_t phiBin  = fHistograms->GetPhiBin(fV0Reader->GetNegativeTrackPhi());
 421     Double_t motherCandidateEta= fV0Reader->GetMotherCandidateEta();
 422     if(fHistograms->fESD_Mapping[phiBin][rBin] != NULL){fHistograms->fESD_Mapping[phiBin][rBin]->Fill(fV0Reader->GetZ(),motherCandidateEta);}
 423     if(fHistograms->fESD_Mapping_Phi[phiBin] != NULL){fHistograms->fESD_Mapping_Phi[phiBin]->Fill(fV0Reader->GetZ(),motherCandidateEta);}
 424     if(fHistograms->fESD_Mapping_R[rBin] != NULL){fHistograms->fESD_Mapping_R[rBin]->Fill(fV0Reader->GetZ(),motherCandidateEta);}
 425     // end mapping
 426
 427
 428     fKFReconstructedGammas.push_back(*fV0Reader->GetMotherCandidateKFCombination());
 429
 430     //----------------------------------- checking for "real" conversions (MC match) --------------------------------------
 431     if(fDoMCTruth){
 432       if(fV0Reader->HasSameMCMother() == kFALSE){
 433         continue;
 434       }
 435
 436       TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
 437       TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
 438       if(negativeMC->GetPdgCode()!=11 || positiveMC->GetPdgCode()!=-11){
 439         continue;
 440       }
 441
 442       if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
 443         if(fHistograms->fESD_Match_Gamma_X_Y){fHistograms->fESD_Match_Gamma_X_Y->Fill(fV0Reader->GetX(),fV0Reader->GetY());}
 444         if(fHistograms->fESD_Match_Gamma_OpeningAngle){fHistograms->fESD_Match_Gamma_OpeningAngle->Fill(fV0Reader->GetOpeningAngle());}
 445         if(fHistograms->fESD_Match_Gamma_Pt){fHistograms->fESD_Match_Gamma_Pt->Fill(fV0Reader->GetMotherCandidatePt());}
 446         if(fHistograms->fESD_Match_Gamma_Energy){fHistograms->fESD_Match_Gamma_Energy->Fill(fV0Reader->GetMotherCandidateEnergy());}
 447         if(fHistograms->fESD_Match_Gamma_Eta){fHistograms->fESD_Match_Gamma_Eta->Fill(fV0Reader->GetMotherCandidateEta());}
 448
 449         if(fHistograms->fESD_Match_Gamma_Phi){fHistograms->fESD_Match_Gamma_Phi->Fill(fV0Reader->GetMotherCandidatePhi());}
 450         if(fHistograms->fESD_Match_Gamma_Mass){fHistograms->fESD_Match_Gamma_Mass->Fill(fV0Reader->GetMotherCandidateMass());}
 451         if(fHistograms->fESD_Match_Gamma_Width){fHistograms->fESD_Match_Gamma_Width->Fill(fV0Reader->GetMotherCandidateWidth());}
 452         if(fHistograms->fESD_Match_Gamma_Chi2){fHistograms->fESD_Match_Gamma_Chi2->Fill(fV0Reader->GetMotherCandidateChi2());}
 453         if(fHistograms->fESD_Match_Gamma_NDF){fHistograms->fESD_Match_Gamma_NDF->Fill(fV0Reader->GetMotherCandidateNDF());}
 454         if(fHistograms->fESD_Match_Gamma_R){fHistograms->fESD_Match_Gamma_R->Fill(fV0Reader->GetXYRadius());}
 455         if(fHistograms->fESD_Match_Gamma_Z_R){fHistograms->fESD_Match_Gamma_Z_R->Fill(fV0Reader->GetZ(),fV0Reader->GetXYRadius());}
 456
 457         //resolution
 458         Double_t mc_pt   = fV0Reader->GetMotherMCParticle()->Pt();
 459         Double_t esd_pt  = fV0Reader->GetMotherCandidatePt();
 460         Double_t res_dPt = ((esd_pt - mc_pt)/mc_pt)*100;
 461         if(fHistograms->fResolution_dPt != NULL){fHistograms->fResolution_dPt->Fill(mc_pt,res_dPt);}
 462         if(fHistograms->fResolution_MC_Pt != NULL){fHistograms->fResolution_MC_Pt->Fill(mc_pt);}
 463         if(fHistograms->fResolution_ESD_Pt != NULL){fHistograms->fResolution_ESD_Pt->Fill(esd_pt);}
 464
 465         Double_t res_dZ = ((fV0Reader->GetZ() -fV0Reader->GetNegativeMCParticle()->Vz())/fV0Reader->GetNegativeMCParticle()->Vz())*100;
 466         if(fHistograms->fResolution_dZ != NULL){fHistograms->fResolution_dZ->Fill(fV0Reader->GetNegativeMCParticle()->Vz(),res_dZ);}
 467         if(fHistograms->fResolution_MC_Z != NULL){fHistograms->fResolution_MC_Z->Fill(fV0Reader->GetNegativeMCParticle()->Vz());}
 468         if(fHistograms->fResolution_ESD_Z != NULL){fHistograms->fResolution_ESD_Z->Fill(fV0Reader->GetZ());}
 469
 470         Double_t res_dR = ((fV0Reader->GetXYRadius() - fV0Reader->GetNegativeMCParticle()->R())/fV0Reader->GetNegativeMCParticle()->R())*100;
 471         if(fHistograms->fResolution_dR != NULL){fHistograms->fResolution_dR->Fill(fV0Reader->GetNegativeMCParticle()->R(),res_dR);}
 472         if(fHistograms->fResolution_MC_R != NULL){fHistograms->fResolution_MC_R->Fill(fV0Reader->GetNegativeMCParticle()->R());}
 473         if(fHistograms->fResolution_ESD_R != NULL){fHistograms->fResolution_ESD_R->Fill(fV0Reader->GetXYRadius());}
 474         if(fHistograms->fResolution_dR_dPt != NULL){fHistograms->fResolution_dR_dPt->Fill(res_dR,res_dPt);}
 475       }
 476     }
 477   }
 478   if(fHistograms->fNumberOfSurvivingV0s != NULL){fHistograms->fNumberOfSurvivingV0s->Fill(nSurvivingV0s);}
 479   if(fHistograms->fNumberOfV0s != NULL){fHistograms->fNumberOfV0s->Fill(fV0Reader->GetNumberOfV0s());}
 480 }
 481
 482 void AliAnalysisTaskGammaConversion::ProcessGammasForNeutralMesonAnalysis(){
 483
 484   for(UInt_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammas.size();firstGammaIndex++){
 485     for(UInt_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fKFReconstructedGammas.size();secondGammaIndex++){
 486       AliKFParticle * twoGammaDecayCandidateDaughter0 = &fKFReconstructedGammas[firstGammaIndex];
 487       AliKFParticle * twoGammaDecayCandidateDaughter1 = &fKFReconstructedGammas[secondGammaIndex];
 488
 489       // Pi0's
 490       AliKFParticle *pi0Candidate = new AliKFParticle(*twoGammaDecayCandidateDaughter0,*twoGammaDecayCandidateDaughter1);
 491
 492       //      pi0Candidate->SetMassConstraint(fPi0Mass,fPi0Width);
 493
 494       Double_t massPi0 =0.;
 495       Double_t widthPi0 = 0.;
 496       Double_t chi2Pi0 =10000.;
 497       pi0Candidate->GetMass(massPi0,widthPi0);
 498       if(pi0Candidate->GetNDF()>0){
 499         chi2Pi0 = pi0Candidate->GetChi2()/pi0Candidate->GetNDF();
 500         //      if(chi2Pi0>0 && chi2Pi0<fChi2Cut){//TODO  find this out
 501         if(chi2Pi0>0 && chi2Pi0<10000){//TODO  find this out see line above
 502
 503           /*
 504           TVector3 vertexDaughter0(twoGammaDecayCandidateDaughter0->Px(),twoGammaDecayCandidateDaughter0->Py(),twoGammaDecayCandidateDaughter0->Pz());
 505           TVector3 vertexDaughter1(twoGammaDecayCandidateDaughter1->Px(),twoGammaDecayCandidateDaughter1->Py(),twoGammaDecayCandidateDaughter1->Pz());
 506           TVector3 vertexPi0Candidate(pi0Candidate->Px(),pi0Candidate->Py(),pi0Candidate->Pz());
 507
 508           Double_t openingAnglePi0 = vertexDaughter0.Angle(vertexDaughter1);
 509
 510           Double_t radiusPi0 = sqrt( vertexPi0Candidate.x()*vertexPi0Candidate.x() + vertexPi0Candidate.y()*vertexPi0Candidate.y() );
 511           */
 512
 513           TVector3 vectorPi0Candidate(pi0Candidate->Px(),pi0Candidate->Py(),pi0Candidate->Pz());
 514
 515           Double_t openingAnglePi0 = twoGammaDecayCandidateDaughter0->GetAngle(*twoGammaDecayCandidateDaughter1);
 516
 517           //Double_t vtx000[3] = {0,0,0};
 518           //   NOT SURE IF THIS DOES WHAT WE WANT: remember to ask Sergey  Double_t radiusPi0 = pi0Candidate->GetDistanceFromVertex(vtx000);
 519           //Calculating by hand the radius
 520           Double_t tmpX= pi0Candidate->GetX();
 521           Double_t tmpY= pi0Candidate->GetY();
 522
 523           Double_t radiusPi0 = TMath::Sqrt(tmpX*tmpX + tmpY*tmpY);
 524
 525           if(fHistograms->fESD_Pi0_OpeningAngleGamma){fHistograms->fESD_Pi0_OpeningAngleGamma->Fill(openingAnglePi0);}
 526           if(fHistograms->fESD_Pi0_Energy){fHistograms->fESD_Pi0_Energy->Fill(pi0Candidate->GetE());}
 527           if(fHistograms->fESD_Pi0_Pt){fHistograms->fESD_Pi0_Pt->Fill(sqrt(pi0Candidate->GetPx()*pi0Candidate->GetPx()+pi0Candidate->GetPy()*pi0Candidate->GetPy()));}
 528           if(fHistograms->fESD_Pi0_Eta){fHistograms->fESD_Pi0_Eta->Fill(vectorPi0Candidate.Eta());}
 529           if(fHistograms->fESD_Pi0_Phi){fHistograms->fESD_Pi0_Phi->Fill(vectorPi0Candidate.Phi());}
 530           if(fHistograms->fESD_Pi0_Mass){fHistograms->fESD_Pi0_Mass->Fill(massPi0);}
 531           if(fHistograms->fESD_Pi0_R){fHistograms->fESD_Pi0_R->Fill(radiusPi0);}
 532           if(fHistograms->fESD_Pi0_Z_R){fHistograms->fESD_Pi0_Z_R->Fill(tmpY,radiusPi0);}
 533           if(fHistograms->fESD_Pi0_X_Y){fHistograms->fESD_Pi0_X_Y->Fill(tmpX,tmpY);}
 534         }
 535       }
 536       delete pi0Candidate;
 537
 538       //Eta's
 539       AliKFParticle *etaCandidate = new AliKFParticle(*twoGammaDecayCandidateDaughter0,*twoGammaDecayCandidateDaughter1);
 540
 541       //      etaCandidate->SetMassConstraint(fEtaMass,fEtaWidth);
 542       Double_t massEta =0.;
 543       Double_t widthEta = 0.;
 544       Double_t chi2Eta =10000.;
 545       etaCandidate->GetMass(massEta,widthEta);
 546       if(etaCandidate->GetNDF()>0){
 547         chi2Eta = etaCandidate->GetChi2()/etaCandidate->GetNDF();
 548         //      if(chi2Eta>0 && chi2Eta<fChi2Cut){
 549         if(chi2Eta>0 && chi2Eta<100000){// TODO find this out se line above
 550           /*
 551           TVector3 vertexDaughter0(twoGammaDecayCandidateDaughter0->Px(),twoGammaDecayCandidateDaughter0->Py(),twoGammaDecayCandidateDaughter0->Pz());
 552           TVector3 vertexDaughter1(twoGammaDecayCandidateDaughter1->Px(),twoGammaDecayCandidateDaughter1->Py(),twoGammaDecayCandidateDaughter1->Pz());
 553           TVector3 vertexEtaCandidate(etaCandidate->Px(),etaCandidate->Py(),etaCandidate->Pz());
 554           Double_t openingAngleEta = vertexDaughter0.Angle(vertexDaughter1);
 555
 556           Double_t radiusEta = sqrt( etaCandidate->GetX()*etaCandidate->GetX() + vertexEtaCandidate.y()*vertexEtaCandidate.y() );
 557           */
 558
 559           //Calculating by hand the radius
 560           Double_t tmpX= etaCandidate->GetX();
 561           Double_t tmpY= etaCandidate->GetY();
 562
 563           Double_t radiusEta = TMath::Sqrt(tmpX*tmpX + tmpY*tmpY);
 564
 565           Double_t openingAngleEta = twoGammaDecayCandidateDaughter0->GetAngle(*twoGammaDecayCandidateDaughter1);
 566
 567           TVector3 vectorEtaCandidate(etaCandidate->Px(),etaCandidate->Py(),etaCandidate->Pz());
 568
 569           if(fHistograms->fESD_Eta_OpeningAngleGamma){fHistograms->fESD_Eta_OpeningAngleGamma->Fill(openingAngleEta);}
 570           if(fHistograms->fESD_Eta_Energy){fHistograms->fESD_Eta_Energy->Fill(etaCandidate->GetE());}
 571           if(fHistograms->fESD_Eta_Pt){fHistograms->fESD_Eta_Pt->Fill(sqrt(etaCandidate->GetPx()*etaCandidate->GetPx()+etaCandidate->GetPy()*etaCandidate->GetPy()));}
 572           if(fHistograms->fESD_Eta_Eta){fHistograms->fESD_Eta_Eta->Fill(vectorEtaCandidate.Eta());}
 573           if(fHistograms->fESD_Eta_Phi){fHistograms->fESD_Eta_Phi->Fill(vectorEtaCandidate.Phi());}
 574           if(fHistograms->fESD_Eta_Mass){fHistograms->fESD_Eta_Mass->Fill(massEta);}
 575           if(fHistograms->fESD_Eta_R){fHistograms->fESD_Eta_R->Fill(radiusEta);}
 576           if(fHistograms->fESD_Eta_Z_R){fHistograms->fESD_Eta_Z_R->Fill(etaCandidate->GetZ(),radiusEta);}
 577           if(fHistograms->fESD_Eta_X_Y){fHistograms->fESD_Eta_X_Y->Fill(tmpX,tmpY);}
 578         }
 579       }
 580       delete etaCandidate;
 581     }
 582   }
 583
 584 }
 585
 586 void AliAnalysisTaskGammaConversion::CalculateBackground(){
 587
 588   for(UInt_t iCurrent=0;iCurrent<fV0Reader->fCurrentEventGoodV0s.size();iCurrent++){
 589     AliKFParticle * currentEventGoodV0 = &fV0Reader->fCurrentEventGoodV0s.at(iCurrent);
 590     for(UInt_t iPrevious=0;iPrevious<fV0Reader->fCurrentEventGoodV0s.size();iPrevious++){
 591       AliKFParticle * previousEventGoodV0 = &fV0Reader->fPreviousEventGoodV0s.at(iPrevious);
 592
 593       AliKFParticle *backgroundCandidate = new AliKFParticle(*currentEventGoodV0,*previousEventGoodV0);
 594
 595       Double_t massBG =0.;
 596       Double_t widthBG = 0.;
 597       Double_t chi2BG =10000.;
 598       backgroundCandidate->GetMass(massBG,widthBG);
 599       if(backgroundCandidate->GetNDF()>0){
 600         chi2BG = backgroundCandidate->GetChi2()/backgroundCandidate->GetNDF();
 601         //      if(chi2Pi0>0 && chi2Pi0<fChi2Cut){//TODO  find this out
 602         if(chi2BG>0 && chi2BG<fV0Reader->GetChi2Cut()){//TODO  find this out see line above
 603
 604           TVector3 vectorBGCandidate(backgroundCandidate->Px(),backgroundCandidate->Py(),backgroundCandidate->Pz());
 605
 606           Double_t openingAngleBG = currentEventGoodV0->GetAngle(*previousEventGoodV0);
 607
 608           //Calculating by hand the radius (find a better way)
 609           Double_t tmpX= backgroundCandidate->GetX();
 610           Double_t tmpY= backgroundCandidate->GetY();
 611
 612           Double_t radiusBG = TMath::Sqrt(tmpX*tmpX + tmpY*tmpY);
 613
 614           if(fHistograms->fESD_Background_OpeningAngleGamma){fHistograms->fESD_Background_OpeningAngleGamma->Fill(openingAngleBG);}
 615           if(fHistograms->fESD_Background_Energy){fHistograms->fESD_Background_Energy->Fill(backgroundCandidate->GetE());}
 616           if(fHistograms->fESD_Background_Pt){fHistograms->fESD_Background_Pt->Fill(sqrt(backgroundCandidate->GetPx()*backgroundCandidate->GetPx()+backgroundCandidate->GetPy()*backgroundCandidate->GetPy()));}
 617           if(fHistograms->fESD_Background_Eta){fHistograms->fESD_Background_Eta->Fill(vectorBGCandidate.Eta());}
 618           if(fHistograms->fESD_Background_Phi){fHistograms->fESD_Background_Phi->Fill(vectorBGCandidate.Phi());}
 619           if(fHistograms->fESD_Background_Mass){fHistograms->fESD_Background_Mass->Fill(massBG);}
 620           if(fHistograms->fESD_Background_R){fHistograms->fESD_Background_R->Fill(radiusBG);}
 621           if(fHistograms->fESD_Background_Z_R){fHistograms->fESD_Background_Z_R->Fill(tmpY,radiusBG);}
 622           if(fHistograms->fESD_Background_X_Y){fHistograms->fESD_Background_X_Y->Fill(tmpX,tmpY);}
 623         }
 624       }
 625       delete backgroundCandidate;
 626     }
 627   }
 628 }
 629
 630 void AliAnalysisTaskGammaConversion::Terminate(Option_t */*option*/)
 631 {
 632   // Terminate analysis
 633   //
 634   AliDebug(1,"Do nothing in Terminate");
 635 }
 636
 637 void AliAnalysisTaskGammaConversion::UserCreateOutputObjects()
 638 {
 639   // Create the output container
 640
 641   fOutputContainer = fHistograms->GetOutputContainer();
 642   fOutputContainer->SetName(GetName()) ;
 643 }
 644
 645 Double_t AliAnalysisTaskGammaConversion::GetMCOpeningAngle(TParticle* daughter0,TParticle* daughter1){
 646   //helper function
 647   TVector3 V3D0(daughter0->Px(),daughter0->Py(),daughter0->Pz());
 648   TVector3 V3D1(daughter1->Px(),daughter1->Py(),daughter1->Pz());
 649   return V3D0.Angle(V3D1);
 650 }
 651