PWG4/GammaConv/AliAnalysisTaskGammaConversion.cxx

   1 /**************************************************************************\r
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
   3  *                                                                        *\r
   4  * Author: Ana Marin, Kathrin Koch, Kenneth Aamodt                        *\r
   5  * Version 1.1                                                            *\r
   6  *                                                                        *\r
   7  * Permission to use, copy, modify and distribute this software and its   *\r
   8  * documentation strictly for non-commercial purposes is hereby granted   *\r
   9  * without fee, provided that the above copyright notice appears in all   *\r
  10  * copies and that both the copyright notice and this permission notice   *\r
  11  * appear in the supporting documentation. The authors make no claims     *\r
  12  * about the suitability of this software for any purpose. It is          *\r
  13  * provided "as is" without express or implied warranty.                  *\r
  14  **************************************************************************/\r
  15 \r
  16 ////////////////////////////////////////////////\r
  17 //--------------------------------------------- \r
  18 // Class used to do analysis on conversion pairs\r
  19 //---------------------------------------------\r
  20 ////////////////////////////////////////////////\r
  21 \r
  22 // root\r
  23 #include <TChain.h>\r
  24 \r
  25 // analysis\r
  26 #include "AliAnalysisTaskGammaConversion.h"\r
  27 #include "AliStack.h"\r
  28 #include "AliLog.h"\r
  29 #include <vector>\r
  30 \r
  31 class AliKFVertex;\r
  32 class AliAODHandler;\r
  33 class AliAODEvent;\r
  34 class ALiESDEvent;\r
  35 class AliMCEvent;\r
  36 class AliMCEventHandler;\r
  37 class AliESDInputHandler;\r
  38 class AliAnalysisManager;\r
  39 class Riostream;\r
  40 class TFile;\r
  41 class TInterpreter;\r
  42 class TSystem;\r
  43 class TROOT;\r
  44 \r
  45 ClassImp(AliAnalysisTaskGammaConversion)\r
  46 \r
  47 \r
  48 AliAnalysisTaskGammaConversion::AliAnalysisTaskGammaConversion():\r
  49 AliAnalysisTaskSE(),\r
  50   fV0Reader(NULL),\r
  51   fStack(NULL),\r
  52   fOutputContainer(NULL),\r
  53   fHistograms(NULL),\r
  54   fDoMCTruth(kFALSE),\r
  55   fMCAllGammas(),\r
  56   fMCPi0s(),\r
  57   fMCEtas(),\r
  58   fMCGammaChic(),\r
  59   fKFReconstructedGammas(),\r
  60   fIsTrueReconstructedGammas(),\r
  61   electronv1(),\r
  62   electronv2(),\r
  63   fElectronMass(-1),\r
  64   fGammaMass(-1),\r
  65   fPi0Mass(-1),\r
  66   fEtaMass(-1),\r
  67   fGammaWidth(-1),\r
  68   fPi0Width(-1),\r
  69   fEtaWidth(-1),\r
  70   fMinOpeningAngleGhostCut(0.),\r
  71   fCalculateBackground(kFALSE),\r
  72   fWriteNtuple(kFALSE),\r
  73   fGammaNtuple(NULL),\r
  74   fNeutralMesonNtuple(NULL),\r
  75   fTotalNumberOfAddedNtupleEntries(0)\r
  76 {\r
  77   // Default constructor\r
  78   // Common I/O in slot 0\r
  79   DefineInput (0, TChain::Class());\r
  80   DefineOutput(0, TTree::Class());\r
  81         \r
  82   // Your private output\r
  83   DefineOutput(1, TList::Class());\r
  84 }\r
  85 \r
  86 AliAnalysisTaskGammaConversion::AliAnalysisTaskGammaConversion(const char* name):\r
  87   AliAnalysisTaskSE(name),\r
  88   fV0Reader(NULL),\r
  89   fStack(NULL),\r
  90   fOutputContainer(0x0),\r
  91   fHistograms(NULL),\r
  92   fDoMCTruth(kFALSE),\r
  93   fMCAllGammas(),\r
  94   fMCPi0s(),\r
  95   fMCEtas(),\r
  96   fMCGammaChic(),\r
  97   fKFReconstructedGammas(),\r
  98   fIsTrueReconstructedGammas(),\r
  99   electronv1(),\r
 100   electronv2(),\r
 101   fElectronMass(-1),\r
 102   fGammaMass(-1),\r
 103   fPi0Mass(-1),\r
 104   fEtaMass(-1),\r
 105   fGammaWidth(-1),\r
 106   fPi0Width(-1),\r
 107   fEtaWidth(-1),\r
 108   fMinOpeningAngleGhostCut(0.),\r
 109   fCalculateBackground(kFALSE),\r
 110   fWriteNtuple(kFALSE),\r
 111   fGammaNtuple(NULL),\r
 112   fNeutralMesonNtuple(NULL),\r
 113   fTotalNumberOfAddedNtupleEntries(0)\r
 114 {\r
 115   // Common I/O in slot 0\r
 116   DefineInput (0, TChain::Class());\r
 117   DefineOutput(0, TTree::Class());\r
 118         \r
 119   // Your private output\r
 120   DefineOutput(1, TList::Class());\r
 121 }\r
 122 \r
 123 AliAnalysisTaskGammaConversion::~AliAnalysisTaskGammaConversion() \r
 124 {\r
 125   // Remove all pointers\r
 126         \r
 127   if(fOutputContainer){\r
 128     fOutputContainer->Clear() ; \r
 129     delete fOutputContainer ;\r
 130   }\r
 131   if(fHistograms){\r
 132     delete fHistograms;\r
 133   }\r
 134   if(fV0Reader){\r
 135     delete fV0Reader;\r
 136   }\r
 137 }\r
 138 \r
 139 \r
 140 void AliAnalysisTaskGammaConversion::Init()\r
 141 {\r
 142   // Initialization\r
 143   AliLog::SetGlobalLogLevel(AliLog::kError);\r
 144 }\r
 145 \r
 146 \r
 147 void AliAnalysisTaskGammaConversion::Exec(Option_t */*option*/)\r
 148 {\r
 149   // Execute analysis for current event\r
 150         \r
 151   ConnectInputData("");\r
 152         \r
 153   //clear vectors\r
 154   fMCAllGammas.clear();\r
 155   fMCPi0s.clear();\r
 156   fMCEtas.clear();\r
 157   fMCGammaChic.clear();\r
 158         \r
 159   fKFReconstructedGammas.clear();\r
 160   fIsTrueReconstructedGammas.clear();\r
 161   electronv1.clear();\r
 162   electronv2.clear();\r
 163         \r
 164   //Clear the data in the v0Reader\r
 165   fV0Reader->UpdateEventByEventData();\r
 166 \r
 167   \r
 168   // Process the MC information\r
 169   if(fDoMCTruth){\r
 170     ProcessMCData();\r
 171   }\r
 172   \r
 173   //Process the v0 information with no cuts\r
 174   ProcessV0sNoCut();\r
 175   \r
 176   // Process the v0 information\r
 177   ProcessV0s();\r
 178   \r
 179   //calculate background if flag is set\r
 180   if(fCalculateBackground){\r
 181     CalculateBackground();\r
 182   }\r
 183   \r
 184   // Process reconstructed gammas\r
 185   ProcessGammasForNeutralMesonAnalysis();\r
 186   \r
 187   PostData(1, fOutputContainer);\r
 188         \r
 189 }\r
 190 \r
 191 void AliAnalysisTaskGammaConversion::ConnectInputData(Option_t */*option*/){\r
 192   // see header file for documentation\r
 193         \r
 194   if(fV0Reader == NULL){\r
 195     // Write warning here cuts and so on are default if this ever happens\r
 196   }\r
 197   fV0Reader->Initialize();\r
 198 }\r
 199 \r
 200 \r
 201 \r
 202 void AliAnalysisTaskGammaConversion::ProcessMCData(){\r
 203   // see header file for documentation\r
 204         \r
 205   fStack = fV0Reader->GetMCStack();\r
 206 \r
 207   if(fV0Reader->CheckForPrimaryVertex() == kFALSE){\r
 208     return; // aborts if the primary vertex does not have contributors.\r
 209   }\r
 210 \r
 211   for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++) {\r
 212     TParticle* particle = (TParticle *)fStack->Particle(iTracks);\r
 213                 \r
 214     if (!particle) {\r
 215       //print warning here\r
 216       continue;\r
 217     }\r
 218 \r
 219     if(TMath::Abs(particle->Eta())> fV0Reader->GetEtaCut() )    continue;\r
 220                                         \r
 221     if(particle->R()>fV0Reader->GetMaxRCut())   continue; // cuts on distance from collision point\r
 222                 \r
 223     Double_t tmpPhi=particle->Phi();\r
 224                 \r
 225     if(particle->Phi()> TMath::Pi()){\r
 226       tmpPhi = particle->Phi()-(2*TMath::Pi());\r
 227     }\r
 228                 \r
 229     Double_t rapidity;\r
 230     if(particle->Energy() - particle->Pz() == 0 || particle->Energy() + particle->Pz() == 0){\r
 231       rapidity=0;\r
 232     }\r
 233     else{\r
 234       rapidity = 0.5*(TMath::Log((particle->Energy()+particle->Pz()) / (particle->Energy()-particle->Pz())));\r
 235     }   \r
 236                 \r
 237     //process the gammas\r
 238     if (particle->GetPdgCode() == 22){\r
 239                         \r
 240       if(particle->GetMother(0) >-1 && fStack->Particle(particle->GetMother(0))->GetPdgCode() == 22){\r
 241         continue; // no photon as mothers!\r
 242       }\r
 243 \r
 244       if(particle->GetMother(0) >= fStack->GetNprimary()){\r
 245         continue; // the gamma has a mother, and it is not a primary particle\r
 246       }\r
 247 \r
 248       fMCAllGammas.push_back(particle);\r
 249                         \r
 250       fHistograms->FillHistogram("MC_allGamma_Energy", particle->Energy());\r
 251       fHistograms->FillHistogram("MC_allGamma_Pt", particle->Pt());\r
 252       fHistograms->FillHistogram("MC_allGamma_Eta", particle->Eta());\r
 253       fHistograms->FillHistogram("MC_allGamma_Phi", tmpPhi);\r
 254       fHistograms->FillHistogram("MC_allGamma_Rapid", rapidity);\r
 255                         \r
 256                         \r
 257       if(particle->GetMother(0) < 0){   // direct gamma\r
 258         fHistograms->FillHistogram("MC_allDirectGamma_Energy",particle->Energy());\r
 259         fHistograms->FillHistogram("MC_allDirectGamma_Pt", particle->Pt());\r
 260         fHistograms->FillHistogram("MC_allDirectGamma_Eta", particle->Eta());\r
 261         fHistograms->FillHistogram("MC_allDirectGamma_Phi", tmpPhi);\r
 262         fHistograms->FillHistogram("MC_allDirectGamma_Rapid", rapidity);                                \r
 263       }\r
 264                         \r
 265                         \r
 266       // looking for conversion (electron + positron from pairbuilding (= 5) )\r
 267       TParticle* ePos = NULL;\r
 268       TParticle* eNeg = NULL;\r
 269                         \r
 270       if(particle->GetNDaughters() >= 2){\r
 271         for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){\r
 272           TParticle *tmpDaughter = fStack->Particle(daughterIndex);\r
 273           if(tmpDaughter->GetUniqueID() == 5){\r
 274             if(tmpDaughter->GetPdgCode() == 11){\r
 275               eNeg = tmpDaughter;\r
 276             }\r
 277             else if(tmpDaughter->GetPdgCode() == -11){\r
 278               ePos = tmpDaughter;\r
 279             }\r
 280           }\r
 281         }\r
 282       }\r
 283                         \r
 284                         \r
 285       if(ePos == NULL || eNeg == NULL){ // means we do not have two daughters from pair production\r
 286         continue;\r
 287       }\r
 288                         \r
 289                         \r
 290       Double_t ePosPhi = ePos->Phi();\r
 291       if(ePos->Phi()> TMath::Pi()) ePosPhi = ePos->Phi()-(2*TMath::Pi());\r
 292                         \r
 293       Double_t eNegPhi = eNeg->Phi();\r
 294       if(eNeg->Phi()> TMath::Pi()) eNegPhi = eNeg->Phi()-(2*TMath::Pi());\r
 295                         \r
 296                         \r
 297       if(ePos->Pt()<fV0Reader->GetPtCut() || eNeg->Pt()<fV0Reader->GetPtCut()){\r
 298         continue; // no reconstruction below the Pt cut\r
 299       }\r
 300                                         \r
 301       if(TMath::Abs(ePos->Eta())> fV0Reader->GetEtaCut() || TMath::Abs(eNeg->Eta())> fV0Reader->GetEtaCut()){\r
 302         continue;\r
 303       } \r
 304                                 \r
 305       if(ePos->R()>fV0Reader->GetMaxRCut()){\r
 306         continue; // cuts on distance from collision point\r
 307       }\r
 308       \r
 309       \r
 310       if((TMath::Abs(ePos->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()  > ePos->R()){\r
 311         continue;               // line cut to exclude regions where we do not reconstruct\r
 312       }         \r
 313                 \r
 314       fHistograms->FillHistogram("MC_ConvGamma_Energy", particle->Energy());\r
 315       fHistograms->FillHistogram("MC_ConvGamma_Pt", particle->Pt());\r
 316       fHistograms->FillHistogram("MC_ConvGamma_Eta", particle->Eta());\r
 317       fHistograms->FillHistogram("MC_ConvGamma_Phi", tmpPhi);\r
 318       fHistograms->FillHistogram("MC_ConvGamma_Rapid", rapidity);\r
 319       fHistograms->FillHistogram("MC_ConvGamma_Pt_Eta", particle->Pt(),particle->Eta());\r
 320                         \r
 321       fHistograms->FillHistogram("MC_E_Energy", eNeg->Energy());\r
 322       fHistograms->FillHistogram("MC_E_Pt", eNeg->Pt());\r
 323       fHistograms->FillHistogram("MC_E_Eta", eNeg->Eta());\r
 324       fHistograms->FillHistogram("MC_E_Phi", eNegPhi);\r
 325                         \r
 326       fHistograms->FillHistogram("MC_P_Energy", ePos->Energy());\r
 327       fHistograms->FillHistogram("MC_P_Pt", ePos->Pt());\r
 328       fHistograms->FillHistogram("MC_P_Eta", ePos->Eta());\r
 329       fHistograms->FillHistogram("MC_P_Phi", ePosPhi);\r
 330                         \r
 331                         \r
 332                         \r
 333       //cout << "filled histos for converted gamma, ePos, eNeg" << endl;\r
 334                         \r
 335       // begin Mapping \r
 336       Int_t rBin    = fHistograms->GetRBin(ePos->R());\r
 337       Int_t phiBin  = fHistograms->GetPhiBin(particle->Phi());\r
 338                         \r
 339       TString nameMCMappingPhiR="";\r
 340       nameMCMappingPhiR.Form("MC_Conversion_Mapping-Phi%02d-R%02d",phiBin,rBin);\r
 341       fHistograms->FillHistogram(nameMCMappingPhiR, ePos->Vz(), particle->Eta());\r
 342                         \r
 343       TString nameMCMappingPhi="";\r
 344       nameMCMappingPhi.Form("MC_Conversion_Mapping-Phi%02d",phiBin);\r
 345       fHistograms->FillHistogram(nameMCMappingPhi, particle->Eta());\r
 346                         \r
 347       TString nameMCMappingR="";\r
 348       nameMCMappingR.Form("MC_Conversion_Mapping-R%02d",rBin);\r
 349       fHistograms->FillHistogram(nameMCMappingR, particle->Eta());\r
 350                         \r
 351       TString nameMCMappingPhiInR="";\r
 352       nameMCMappingPhiInR.Form("MC_Conversion_Mapping_Phi_R-%02d",rBin);\r
 353       fHistograms->FillHistogram(nameMCMappingPhiInR, tmpPhi);\r
 354       //end mapping\r
 355                         \r
 356       fHistograms->FillHistogram("MC_Conversion_R",ePos->R());\r
 357       fHistograms->FillHistogram("MC_Conversion_ZR",ePos->Vz(),ePos->R());\r
 358       fHistograms->FillHistogram("MC_Conversion_XY",ePos->Vx(),ePos->Vy());\r
 359       fHistograms->FillHistogram("MC_Conversion_OpeningAngle",GetMCOpeningAngle(ePos, eNeg));\r
 360                         \r
 361       //cout << "mapping is done" << endl;\r
 362                         \r
 363                         \r
 364       if(particle->GetMother(0) < 0){ // no mother = direct gamma, still inside converted\r
 365         fHistograms->FillHistogram("MC_ConvDirectGamma_Energy",particle->Energy());\r
 366         fHistograms->FillHistogram("MC_ConvDirectGamma_Pt", particle->Pt());\r
 367         fHistograms->FillHistogram("MC_ConvDirectGamma_Eta", particle->Eta());\r
 368         fHistograms->FillHistogram("MC_ConvDirectGamma_Phi", tmpPhi);\r
 369         fHistograms->FillHistogram("MC_ConvDirectGamma_Rapid", rapidity);\r
 370                                 \r
 371       } // end direct gamma\r
 372       else{   // mother exits \r
 373         if( fStack->Particle(particle->GetMother(0))->GetPdgCode()==10441 ||//chic0 \r
 374             fStack->Particle(particle->GetMother(0))->GetPdgCode()==20443 ||//psi2S\r
 375             fStack->Particle(particle->GetMother(0))->GetPdgCode()==445  //chic2\r
 376             ){ \r
 377           fMCGammaChic.push_back(particle);\r
 378         }\r
 379       }  // end if mother exits\r
 380     } // end if particle is a photon\r
 381                 \r
 382     if(particle->GetNDaughters() == 2){\r
 383                         \r
 384       TParticle* daughter0 = (TParticle*)fStack->Particle(particle->GetFirstDaughter());\r
 385       TParticle* daughter1 = (TParticle*)fStack->Particle(particle->GetLastDaughter());\r
 386                         \r
 387       if(daughter0->GetPdgCode() != 22 || daughter1->GetPdgCode() != 22) continue; //check for gamma gamma daughters\r
 388                         \r
 389                         \r
 390                         \r
 391       // check for conversions now -> have to pass eta and line cut!\r
 392       Bool_t daughter0Electron = kFALSE;\r
 393       Bool_t daughter0Positron = kFALSE;\r
 394       Bool_t daughter1Electron = kFALSE;\r
 395       Bool_t daughter1Positron = kFALSE;\r
 396                         \r
 397                         \r
 398                         \r
 399       if(daughter0->GetNDaughters() >= 2){\r
 400         for(Int_t TrackIndex=daughter0->GetFirstDaughter();TrackIndex<=daughter0->GetLastDaughter();TrackIndex++){\r
 401           TParticle *tmpDaughter = fStack->Particle(TrackIndex);\r
 402           if(tmpDaughter->GetUniqueID() == 5){\r
 403             if(tmpDaughter->GetPdgCode() == 11){\r
 404               if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){\r
 405                 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()  < tmpDaughter->R() ){\r
 406                   daughter0Electron = kTRUE;\r
 407                 }\r
 408                                                                 \r
 409               }\r
 410             }\r
 411             else if(tmpDaughter->GetPdgCode() == -11){\r
 412               if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){\r
 413                 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()  < tmpDaughter->R() ){\r
 414                   daughter0Positron = kTRUE;\r
 415                 }\r
 416                                                                 \r
 417               }\r
 418                                                         \r
 419             }\r
 420           }\r
 421         }\r
 422       }\r
 423                         \r
 424                         \r
 425         \r
 426       if(daughter1->GetNDaughters() >= 2){\r
 427         for(Int_t TrackIndex=daughter1->GetFirstDaughter();TrackIndex<=daughter1->GetLastDaughter();TrackIndex++){\r
 428           TParticle *tmpDaughter = fStack->Particle(TrackIndex);\r
 429           if(tmpDaughter->GetUniqueID() == 5){\r
 430             if(tmpDaughter->GetPdgCode() == 11){\r
 431               if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){\r
 432                 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()  < tmpDaughter->R() ){\r
 433                   daughter1Electron = kTRUE;\r
 434                 }\r
 435                                                                 \r
 436               }\r
 437             }\r
 438             else if(tmpDaughter->GetPdgCode() == -11){\r
 439               if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){\r
 440                 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()  < tmpDaughter->R() ){\r
 441                   daughter1Positron = kTRUE;\r
 442                 }\r
 443                                                                 \r
 444               }\r
 445                                                         \r
 446             }\r
 447           }\r
 448         }\r
 449       }\r
 450                         \r
 451                                                                                                 \r
 452                         \r
 453                         \r
 454       if(particle->GetPdgCode()==111){     //Pi0\r
 455         if( iTracks >= fStack->GetNprimary()){\r
 456           fHistograms->FillHistogram("MC_Pi0_Secondaries_Eta", particle->Eta());\r
 457           fHistograms->FillHistogram("MC_Pi0_Secondaries_Rapid", rapidity);\r
 458           fHistograms->FillHistogram("MC_Pi0_Secondaries_Phi", tmpPhi);\r
 459           fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt", particle->Pt());\r
 460           fHistograms->FillHistogram("MC_Pi0_Secondaries_Energy", particle->Energy());\r
 461           fHistograms->FillHistogram("MC_Pi0_Secondaries_R", particle->R());\r
 462           fHistograms->FillHistogram("MC_Pi0_Secondaries_ZR", particle->Vz(),particle->R());\r
 463           fHistograms->FillHistogram("MC_Pi0_Secondaries_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
 464           fHistograms->FillHistogram("MC_Pi0_Secondaries_XY", particle->Vx(),particle->Vy());//only fill from one daughter to avoid multiple filling\r
 465                                         \r
 466           if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){\r
 467             fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());\r
 468             fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);\r
 469             if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){\r
 470               fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());\r
 471               fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);\r
 472             }\r
 473           }\r
 474         }\r
 475         else{\r
 476           fHistograms->FillHistogram("MC_Pi0_Eta", particle->Eta());    \r
 477           fHistograms->FillHistogram("MC_Pi0_Rapid", rapidity);\r
 478           fHistograms->FillHistogram("MC_Pi0_Phi", tmpPhi);\r
 479           fHistograms->FillHistogram("MC_Pi0_Pt", particle->Pt());\r
 480           fHistograms->FillHistogram("MC_Pi0_Energy", particle->Energy());\r
 481           fHistograms->FillHistogram("MC_Pi0_R", particle->R());\r
 482           fHistograms->FillHistogram("MC_Pi0_ZR", particle->Vz(),particle->R());\r
 483           fHistograms->FillHistogram("MC_Pi0_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
 484           fHistograms->FillHistogram("MC_Pi0_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling\r
 485                                         \r
 486           if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){\r
 487             fHistograms->FillHistogram("MC_Pi0_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());\r
 488             fHistograms->FillHistogram("MC_Pi0_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);\r
 489             if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){\r
 490               fHistograms->FillHistogram("MC_Pi0_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());\r
 491               fHistograms->FillHistogram("MC_Pi0_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);\r
 492             }\r
 493           }\r
 494         }\r
 495       }\r
 496                         \r
 497       if(particle->GetPdgCode()==221){   //Eta\r
 498         fHistograms->FillHistogram("MC_Eta_Eta", particle->Eta());\r
 499         fHistograms->FillHistogram("MC_Eta_Rapid", rapidity);\r
 500         fHistograms->FillHistogram("MC_Eta_Phi",tmpPhi);\r
 501         fHistograms->FillHistogram("MC_Eta_Pt", particle->Pt());\r
 502         fHistograms->FillHistogram("MC_Eta_Energy", particle->Energy());\r
 503         fHistograms->FillHistogram("MC_Eta_R", particle->R());\r
 504         fHistograms->FillHistogram("MC_Eta_ZR", particle->Vz(),particle->R());\r
 505         fHistograms->FillHistogram("MC_Eta_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
 506         fHistograms->FillHistogram("MC_Eta_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling\r
 507                                 \r
 508         if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){\r
 509           fHistograms->FillHistogram("MC_Eta_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());\r
 510           fHistograms->FillHistogram("MC_Eta_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);\r
 511           if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){\r
 512             fHistograms->FillHistogram("MC_Eta_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());\r
 513             fHistograms->FillHistogram("MC_Eta_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);\r
 514           }\r
 515                                         \r
 516         }\r
 517                                 \r
 518       }\r
 519                         \r
 520       // all motherparticles with 2 gammas as daughters\r
 521       fHistograms->FillHistogram("MC_Mother_R", particle->R());\r
 522       fHistograms->FillHistogram("MC_Mother_ZR", particle->Vz(),particle->R());\r
 523       fHistograms->FillHistogram("MC_Mother_XY", particle->Vx(),particle->Vy());\r
 524       fHistograms->FillHistogram("MC_Mother_Mass", particle->GetCalcMass());\r
 525       fHistograms->FillHistogram("MC_Mother_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
 526       fHistograms->FillHistogram("MC_Mother_Energy", particle->Energy());\r
 527       fHistograms->FillHistogram("MC_Mother_Pt", particle->Pt());\r
 528       fHistograms->FillHistogram("MC_Mother_Eta", particle->Eta());\r
 529       fHistograms->FillHistogram("MC_Mother_Rapid", rapidity);\r
 530       fHistograms->FillHistogram("MC_Mother_Phi",tmpPhi);\r
 531       fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt",particle->GetMass(),particle->Pt());                 \r
 532       if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){\r
 533         fHistograms->FillHistogram("MC_Mother_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());\r
 534         fHistograms->FillHistogram("MC_Mother_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);\r
 535         fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt_withinAcceptance",particle->GetMass(),particle->Pt());                      \r
 536         if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){\r
 537           fHistograms->FillHistogram("MC_Mother_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());\r
 538           fHistograms->FillHistogram("MC_Mother_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);\r
 539           fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt_ConvGamma_withinAcceptance",particle->GetMass(),particle->Pt());                  \r
 540 \r
 541         }\r
 542                                 \r
 543                                 \r
 544       }\r
 545                         \r
 546       //cout << "mother histos are filled" << endl;\r
 547                         \r
 548     } // end if(particle->GetNDaughters() == 2)\r
 549                 \r
 550   }// end for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++)\r
 551         \r
 552   //cout << "right before the end of processMCdata" << endl;\r
 553         \r
 554 } // end ProcessMCData\r
 555 \r
 556 \r
 557 \r
 558 void AliAnalysisTaskGammaConversion::FillNtuple(){\r
 559         \r
 560   if(fGammaNtuple == NULL){\r
 561     return;\r
 562   }\r
 563   Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();\r
 564   for(Int_t i=0;i<numberOfV0s;i++){\r
 565     Float_t values[27] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};\r
 566     AliESDv0 * cV0 = fV0Reader->GetV0(i);\r
 567     Double_t negPID=0;\r
 568     Double_t posPID=0;\r
 569     fV0Reader->GetPIDProbability(negPID,posPID);\r
 570     values[0]=cV0->GetOnFlyStatus();\r
 571     values[1]=fV0Reader->CheckForPrimaryVertex();\r
 572     values[2]=negPID;\r
 573     values[3]=posPID;\r
 574     values[4]=fV0Reader->GetX();\r
 575     values[5]=fV0Reader->GetY();\r
 576     values[6]=fV0Reader->GetZ();\r
 577     values[7]=fV0Reader->GetXYRadius();\r
 578     values[8]=fV0Reader->GetMotherCandidateNDF();\r
 579     values[9]=fV0Reader->GetMotherCandidateChi2();\r
 580     values[10]=fV0Reader->GetMotherCandidateEnergy();\r
 581     values[11]=fV0Reader->GetMotherCandidateEta();\r
 582     values[12]=fV0Reader->GetMotherCandidatePt();\r
 583     values[13]=fV0Reader->GetMotherCandidateMass();\r
 584     values[14]=fV0Reader->GetMotherCandidateWidth();\r
 585     //    values[15]=fV0Reader->GetMotherMCParticle()->Pt();   MOVED TO THE END, HAS TO BE CALLED AFTER HasSameMother NB: still has the same entry in the array\r
 586     values[16]=fV0Reader->GetOpeningAngle();\r
 587     values[17]=fV0Reader->GetNegativeTrackEnergy();\r
 588     values[18]=fV0Reader->GetNegativeTrackPt();\r
 589     values[19]=fV0Reader->GetNegativeTrackEta();\r
 590     values[20]=fV0Reader->GetNegativeTrackPhi();\r
 591     values[21]=fV0Reader->GetPositiveTrackEnergy();\r
 592     values[22]=fV0Reader->GetPositiveTrackPt();\r
 593     values[23]=fV0Reader->GetPositiveTrackEta();\r
 594     values[24]=fV0Reader->GetPositiveTrackPhi();\r
 595     values[25]=fV0Reader->HasSameMCMother();\r
 596     if(values[25] != 0){\r
 597       values[26]=fV0Reader->GetMotherMCParticlePDGCode();\r
 598       values[15]=fV0Reader->GetMotherMCParticle()->Pt();\r
 599     }\r
 600     fTotalNumberOfAddedNtupleEntries++;\r
 601     fGammaNtuple->Fill(values);\r
 602   }\r
 603   fV0Reader->ResetV0IndexNumber();\r
 604         \r
 605 }\r
 606 \r
 607 void AliAnalysisTaskGammaConversion::ProcessV0sNoCut(){\r
 608 \r
 609   Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();\r
 610   for(Int_t i=0;i<numberOfV0s;i++){\r
 611     /*AliESDv0 * cV0 = */fV0Reader->GetV0(i);\r
 612 \r
 613     if(fV0Reader->CheckForPrimaryVertex() == kFALSE){\r
 614       return;\r
 615     }\r
 616     \r
 617     if(fDoMCTruth){\r
 618       \r
 619       if(fV0Reader->HasSameMCMother() == kFALSE){\r
 620         continue;\r
 621       }\r
 622                 \r
 623       TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();\r
 624       TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();\r
 625 \r
 626       if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){\r
 627         continue;\r
 628       }\r
 629       if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){\r
 630         continue;\r
 631       }\r
 632         \r
 633       if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){\r
 634       \r
 635         fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt", fV0Reader->GetMotherCandidatePt());\r
 636         fHistograms->FillHistogram("ESD_NoCutConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());\r
 637         fHistograms->FillHistogram("ESD_NoCutConvGamma_Eta", fV0Reader->GetMotherCandidateEta());                               \r
 638         fHistograms->FillHistogram("ESD_NoCutConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());\r
 639         fHistograms->FillHistogram("ESD_NoCutConvGamma_Mass", fV0Reader->GetMotherCandidateMass());\r
 640         fHistograms->FillHistogram("ESD_NoCutConvGamma_Width", fV0Reader->GetMotherCandidateWidth());\r
 641         fHistograms->FillHistogram("ESD_NoCutConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());\r
 642         fHistograms->FillHistogram("ESD_NoCutConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());\r
 643         fHistograms->FillHistogram("ESD_NoCutConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());\r
 644         fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());\r
 645         \r
 646         fHistograms->FillHistogram("ESD_NoCutConversion_XY", fV0Reader->GetX(),fV0Reader->GetY());\r
 647         fHistograms->FillHistogram("ESD_NoCutConversion_R", fV0Reader->GetXYRadius());\r
 648         fHistograms->FillHistogram("ESD_NoCutConversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());\r
 649         fHistograms->FillHistogram("ESD_NoCutConversion_OpeningAngle", fV0Reader->GetOpeningAngle());\r
 650         \r
 651         /*\r
 652           ESD_NoCutConvGamma_Pt\r
 653           ESD_NoCutConvGamma_Energy\r
 654           ESD_NoCutConvGamma_Eta\r
 655           ESD_NoCutConvGamma_Phi\r
 656           ESD_NoCutConvGamma_Mass\r
 657           ESD_NoCutConvGamma_Width\r
 658           ESD_NoCutConvGamma_Chi2\r
 659           ESD_NoCutConvGamma_NDF\r
 660           ESD_NoCutConvGamma_PtvsEta\r
 661           ESD_NoCutConversion_XY\r
 662           ESD_NoCutConversion_R\r
 663           ESD_NoCutConversion_ZR\r
 664           ESD_NoCutConversion_OpeningAngle\r
 665         */\r
 666       }\r
 667     }\r
 668   }\r
 669   fV0Reader->ResetV0IndexNumber();\r
 670 }\r
 671 \r
 672 void AliAnalysisTaskGammaConversion::ProcessV0s(){\r
 673   // see header file for documentation\r
 674         \r
 675   if(fWriteNtuple == kTRUE){\r
 676     FillNtuple();\r
 677   }\r
 678         \r
 679   Int_t nSurvivingV0s=0;\r
 680   while(fV0Reader->NextV0()){\r
 681     nSurvivingV0s++;\r
 682                 \r
 683                 \r
 684     //-------------------------- filling v0 information -------------------------------------\r
 685     fHistograms->FillHistogram("ESD_Conversion_R", fV0Reader->GetXYRadius());\r
 686     fHistograms->FillHistogram("ESD_Conversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());\r
 687     fHistograms->FillHistogram("ESD_Conversion_XY", fV0Reader->GetX(),fV0Reader->GetY());\r
 688     fHistograms->FillHistogram("ESD_Conversion_OpeningAngle", fV0Reader->GetOpeningAngle());    \r
 689                 \r
 690     fHistograms->FillHistogram("ESD_E_Energy", fV0Reader->GetNegativeTrackEnergy());\r
 691     fHistograms->FillHistogram("ESD_E_Pt", fV0Reader->GetNegativeTrackPt());\r
 692     fHistograms->FillHistogram("ESD_E_Eta", fV0Reader->GetNegativeTrackEta());\r
 693     fHistograms->FillHistogram("ESD_E_Phi", fV0Reader->GetNegativeTrackPhi());\r
 694                 \r
 695     fHistograms->FillHistogram("ESD_P_Energy", fV0Reader->GetPositiveTrackEnergy());\r
 696     fHistograms->FillHistogram("ESD_P_Pt", fV0Reader->GetPositiveTrackPt());\r
 697     fHistograms->FillHistogram("ESD_P_Eta", fV0Reader->GetPositiveTrackEta());\r
 698     fHistograms->FillHistogram("ESD_P_Phi", fV0Reader->GetPositiveTrackPhi());\r
 699                 \r
 700     fHistograms->FillHistogram("ESD_ConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());\r
 701     fHistograms->FillHistogram("ESD_ConvGamma_Pt", fV0Reader->GetMotherCandidatePt());\r
 702     fHistograms->FillHistogram("ESD_ConvGamma_Eta", fV0Reader->GetMotherCandidateEta());\r
 703     fHistograms->FillHistogram("ESD_ConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());\r
 704     fHistograms->FillHistogram("ESD_ConvGamma_Mass", fV0Reader->GetMotherCandidateMass());\r
 705     fHistograms->FillHistogram("ESD_ConvGamma_Width", fV0Reader->GetMotherCandidateWidth());\r
 706     fHistograms->FillHistogram("ESD_ConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());\r
 707     fHistograms->FillHistogram("ESD_ConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());\r
 708     fHistograms->FillHistogram("ESD_ConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());\r
 709     fHistograms->FillHistogram("ESD_ConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());\r
 710                 \r
 711                 \r
 712     // begin mapping\r
 713     Int_t rBin    = fHistograms->GetRBin(fV0Reader->GetXYRadius());\r
 714     Int_t phiBin  = fHistograms->GetPhiBin(fV0Reader->GetNegativeTrackPhi());\r
 715     Double_t motherCandidateEta= fV0Reader->GetMotherCandidateEta();\r
 716                 \r
 717     TString nameESDMappingPhiR="";\r
 718     nameESDMappingPhiR.Form("ESD_Conversion_Mapping-Phi%02d-R%02d",phiBin,rBin);\r
 719     fHistograms->FillHistogram(nameESDMappingPhiR, fV0Reader->GetZ(), motherCandidateEta);\r
 720                 \r
 721     TString nameESDMappingPhi="";\r
 722     nameESDMappingPhi.Form("ESD_Conversion_Mapping-Phi%02d",phiBin);\r
 723     fHistograms->FillHistogram(nameESDMappingPhi, fV0Reader->GetZ(), motherCandidateEta);\r
 724                 \r
 725     TString nameESDMappingR="";\r
 726     nameESDMappingR.Form("ESD_Conversion_Mapping-R%02d",rBin);\r
 727     fHistograms->FillHistogram(nameESDMappingR, fV0Reader->GetZ(), motherCandidateEta);  \r
 728                 \r
 729     TString nameESDMappingPhiInR="";\r
 730     nameESDMappingPhiInR.Form("ESD_Conversion_Mapping_Phi_R-%02d",rBin);\r
 731     fHistograms->FillHistogram(nameESDMappingPhiInR, fV0Reader->GetMotherCandidatePhi());\r
 732     // end mapping\r
 733                 \r
 734     fKFReconstructedGammas.push_back(*fV0Reader->GetMotherCandidateKFCombination());\r
 735     electronv1.push_back(fV0Reader->GetCurrentV0()->GetPindex());\r
 736     electronv2.push_back(fV0Reader->GetCurrentV0()->GetNindex());\r
 737 \r
 738                 \r
 739     //----------------------------------- checking for "real" conversions (MC match) --------------------------------------\r
 740     if(fDoMCTruth){\r
 741                         \r
 742       if(fV0Reader->HasSameMCMother() == kFALSE){\r
 743         fIsTrueReconstructedGammas.push_back(kFALSE);\r
 744         continue;\r
 745       }\r
 746       \r
 747                 \r
 748       TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();\r
 749       TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();\r
 750 \r
 751       if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){\r
 752         fIsTrueReconstructedGammas.push_back(kFALSE);\r
 753         continue;\r
 754       }\r
 755       if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){\r
 756         fIsTrueReconstructedGammas.push_back(kFALSE);\r
 757         continue;\r
 758       }\r
 759         \r
 760       if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){\r
 761         fIsTrueReconstructedGammas.push_back(kTRUE);\r
 762                                 \r
 763         fHistograms->FillHistogram("ESD_TrueConvGamma_Pt", fV0Reader->GetMotherCandidatePt());\r
 764         fHistograms->FillHistogram("ESD_TrueConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());\r
 765         fHistograms->FillHistogram("ESD_TrueConvGamma_Eta", fV0Reader->GetMotherCandidateEta());                                \r
 766         fHistograms->FillHistogram("ESD_TrueConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());\r
 767         fHistograms->FillHistogram("ESD_TrueConvGamma_Mass", fV0Reader->GetMotherCandidateMass());\r
 768         fHistograms->FillHistogram("ESD_TrueConvGamma_Width", fV0Reader->GetMotherCandidateWidth());\r
 769         fHistograms->FillHistogram("ESD_TrueConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());\r
 770         fHistograms->FillHistogram("ESD_TrueConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());\r
 771         fHistograms->FillHistogram("ESD_TrueConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());\r
 772         fHistograms->FillHistogram("ESD_TrueConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());\r
 773         fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLength", /*fV0Reader->GetNegativeTrackLength()*/fV0Reader->GetNegativeNTPCClusters());\r
 774         fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLength", /*fV0Reader->GetPositiveTrackLength()*/fV0Reader->GetPositiveNTPCClusters());\r
 775         fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLengthVSInvMass",/*fV0Reader->GetNegativeTrackLength()*/fV0Reader->GetNegativeNTPCClusters(),fV0Reader->GetMotherCandidateMass());\r
 776         fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLengthVSInvMass",/*fV0Reader->GetPositiveTrackLength()*/fV0Reader->GetPositiveNTPCClusters(),fV0Reader->GetMotherCandidateMass());\r
 777         \r
 778         fHistograms->FillHistogram("ESD_TrueConversion_XY", fV0Reader->GetX(),fV0Reader->GetY());\r
 779         fHistograms->FillHistogram("ESD_TrueConversion_R", fV0Reader->GetXYRadius());\r
 780         fHistograms->FillHistogram("ESD_TrueConversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());\r
 781         fHistograms->FillHistogram("ESD_TrueConversion_OpeningAngle", fV0Reader->GetOpeningAngle());\r
 782 \r
 783         \r
 784         /*\r
 785           fHistograms->FillHistogram("ESD_TrueConversion_XY", fV0Reader->GetX(),fV0Reader->GetY());\r
 786           fHistograms->FillHistogram("ESD_TrueConversion_R", fV0Reader->GetXYRadius());\r
 787           fHistograms->FillHistogram("ESD_TrueConversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());\r
 788           fHistograms->FillHistogram("ESD_TrueConversion_OpeningAngle", fV0Reader->GetOpeningAngle());\r
 789         */\r
 790 \r
 791 \r
 792 \r
 793         //resolution\r
 794         Double_t mcpt   = fV0Reader->GetMotherMCParticle()->Pt();\r
 795         Double_t esdpt  = fV0Reader->GetMotherCandidatePt();\r
 796         Double_t resdPt = 0;\r
 797         if(mcpt != 0){\r
 798           resdPt = ((esdpt - mcpt)/mcpt)*100;\r
 799         }\r
 800                                 \r
 801         fHistograms->FillHistogram("Resolution_dPt", mcpt, resdPt);\r
 802         fHistograms->FillHistogram("Resolution_MC_Pt", mcpt);\r
 803         fHistograms->FillHistogram("Resolution_ESD_Pt", esdpt);\r
 804                                 \r
 805         Double_t resdZ = 0;\r
 806         if(fV0Reader->GetNegativeMCParticle()->Vz() != 0){\r
 807           resdZ = ((fV0Reader->GetZ() -fV0Reader->GetNegativeMCParticle()->Vz())/fV0Reader->GetNegativeMCParticle()->Vz())*100;\r
 808         }\r
 809                                 \r
 810         fHistograms->FillHistogram("Resolution_dZ", fV0Reader->GetNegativeMCParticle()->Vz(), resdZ);\r
 811         fHistograms->FillHistogram("Resolution_MC_Z", fV0Reader->GetNegativeMCParticle()->Vz());\r
 812         fHistograms->FillHistogram("Resolution_ESD_Z", fV0Reader->GetZ());\r
 813                                 \r
 814         Double_t resdR = 0;\r
 815         if(fV0Reader->GetNegativeMCParticle()->R() != 0){\r
 816           resdR = ((fV0Reader->GetXYRadius() - fV0Reader->GetNegativeMCParticle()->R())/fV0Reader->GetNegativeMCParticle()->R())*100;\r
 817         }\r
 818 \r
 819         fHistograms->FillHistogram("Resolution_dR", fV0Reader->GetNegativeMCParticle()->R(), resdR);\r
 820         fHistograms->FillHistogram("Resolution_MC_R", fV0Reader->GetNegativeMCParticle()->R());\r
 821         fHistograms->FillHistogram("Resolution_ESD_R", fV0Reader->GetXYRadius());\r
 822         fHistograms->FillHistogram("Resolution_dR_dPt", resdR, resdPt);\r
 823       }//if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22)\r
 824       else{\r
 825         fIsTrueReconstructedGammas.push_back(kFALSE);\r
 826       }\r
 827     }//if(fDoMCTruth)\r
 828   }//while(fV0Reader->NextV0)\r
 829   fHistograms->FillHistogram("ESD_NumberOfSurvivingV0s", nSurvivingV0s);\r
 830   fHistograms->FillHistogram("ESD_NumberOfV0s", fV0Reader->GetNumberOfV0s());\r
 831         \r
 832   //cout << "nearly at the end of doMCTruth" << endl;\r
 833         \r
 834 }\r
 835 \r
 836 void AliAnalysisTaskGammaConversion::ProcessGammasForNeutralMesonAnalysis(){\r
 837   // see header file for documentation\r
 838         \r
 839   for(UInt_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammas.size();firstGammaIndex++){\r
 840     for(UInt_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fKFReconstructedGammas.size();secondGammaIndex++){\r
 841                         \r
 842       AliKFParticle * twoGammaDecayCandidateDaughter0 = &fKFReconstructedGammas[firstGammaIndex];\r
 843       AliKFParticle * twoGammaDecayCandidateDaughter1 = &fKFReconstructedGammas[secondGammaIndex];\r
 844       \r
 845       if(electronv1[firstGammaIndex]==electronv1[secondGammaIndex] || electronv1[firstGammaIndex]==electronv2[secondGammaIndex]){\r
 846         continue;\r
 847       }\r
 848       if(electronv2[firstGammaIndex]==electronv1[secondGammaIndex] || electronv2[firstGammaIndex]==electronv2[secondGammaIndex]){\r
 849         continue;\r
 850       }\r
 851 \r
 852       /*\r
 853         if(fIsTrueReconstructedGammas[firstGammaIndex] == kFALSE || fIsTrueReconstructedGammas[secondGammaIndex] == kFALSE){\r
 854         continue;\r
 855         }\r
 856       */\r
 857 \r
 858       AliKFParticle *twoGammaCandidate = new AliKFParticle(*twoGammaDecayCandidateDaughter0,*twoGammaDecayCandidateDaughter1);\r
 859                         \r
 860       Double_t massTwoGammaCandidate = 0.;\r
 861       Double_t widthTwoGammaCandidate = 0.;\r
 862       Double_t chi2TwoGammaCandidate =10000.;   \r
 863       twoGammaCandidate->GetMass(massTwoGammaCandidate,widthTwoGammaCandidate);\r
 864       if(twoGammaCandidate->GetNDF()>0){\r
 865         chi2TwoGammaCandidate = twoGammaCandidate->GetChi2()/twoGammaCandidate->GetNDF();\r
 866                                 \r
 867         if(chi2TwoGammaCandidate>0 && chi2TwoGammaCandidate<fV0Reader->GetChi2CutMeson()){                                      \r
 868                                         \r
 869           TVector3 MomentumVectorTwoGammaCandidate(twoGammaCandidate->GetPx(),twoGammaCandidate->GetPy(),twoGammaCandidate->GetPz());\r
 870           TVector3 SpaceVectorTwoGammaCandidate(twoGammaCandidate->GetX(),twoGammaCandidate->GetY(),twoGammaCandidate->GetZ());\r
 871                                         \r
 872           Double_t openingAngleTwoGammaCandidate = twoGammaDecayCandidateDaughter0->GetAngle(*twoGammaDecayCandidateDaughter1);                                 \r
 873           Double_t rapidity;\r
 874           if(twoGammaCandidate->GetE() - twoGammaCandidate->GetPz() == 0 || twoGammaCandidate->GetE() + twoGammaCandidate->GetPz() == 0){\r
 875             rapidity=0;\r
 876           }\r
 877           else{\r
 878             rapidity = 0.5*(TMath::Log((twoGammaCandidate->GetE() +twoGammaCandidate->GetPz()) / (twoGammaCandidate->GetE()-twoGammaCandidate->GetPz())));\r
 879           }\r
 880                                         \r
 881           if(openingAngleTwoGammaCandidate < fMinOpeningAngleGhostCut) continue;   // minimum opening angle to avoid using ghosttracks\r
 882 \r
 883           fHistograms->FillHistogram("ESD_Mother_GammaDaughter_OpeningAngle", openingAngleTwoGammaCandidate);\r
 884           fHistograms->FillHistogram("ESD_Mother_Energy", twoGammaCandidate->GetE());\r
 885           fHistograms->FillHistogram("ESD_Mother_Pt", MomentumVectorTwoGammaCandidate.Pt());\r
 886           fHistograms->FillHistogram("ESD_Mother_Eta", MomentumVectorTwoGammaCandidate.Eta());\r
 887           fHistograms->FillHistogram("ESD_Mother_Rapid", rapidity);                                     \r
 888           fHistograms->FillHistogram("ESD_Mother_Phi", SpaceVectorTwoGammaCandidate.Phi());\r
 889           fHistograms->FillHistogram("ESD_Mother_Mass", massTwoGammaCandidate);\r
 890           fHistograms->FillHistogram("ESD_Mother_R", SpaceVectorTwoGammaCandidate.Pt());    // Pt in Space == R!!!\r
 891           fHistograms->FillHistogram("ESD_Mother_ZR", twoGammaCandidate->GetZ(), SpaceVectorTwoGammaCandidate.Pt());\r
 892           fHistograms->FillHistogram("ESD_Mother_XY", twoGammaCandidate->GetX(), twoGammaCandidate->GetY());\r
 893           fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt",massTwoGammaCandidate ,MomentumVectorTwoGammaCandidate.Pt());\r
 894         }\r
 895       }\r
 896       delete twoGammaCandidate;\r
 897                         \r
 898       //cout << "nearly at the end of processgamma for neutral meson ..." << endl;\r
 899                         \r
 900                         \r
 901     }\r
 902   }\r
 903 }\r
 904 \r
 905 void AliAnalysisTaskGammaConversion::CalculateBackground(){\r
 906   // see header file for documentation\r
 907         \r
 908   vector<AliKFParticle> vectorCurrentEventGoodV0s = fV0Reader->GetCurrentEventGoodV0s();\r
 909   vector<AliKFParticle> vectorPreviousEventGoodV0s = fV0Reader->GetPreviousEventGoodV0s();\r
 910   for(UInt_t iCurrent=0;iCurrent<vectorCurrentEventGoodV0s.size();iCurrent++){\r
 911     AliKFParticle * currentEventGoodV0 = &vectorCurrentEventGoodV0s.at(iCurrent);\r
 912     for(UInt_t iPrevious=0;iPrevious<vectorPreviousEventGoodV0s.size();iPrevious++){\r
 913       AliKFParticle * previousGoodV0 = &vectorPreviousEventGoodV0s.at(iPrevious);\r
 914 \r
 915       AliKFParticle *backgroundCandidate = new AliKFParticle(*currentEventGoodV0,*previousGoodV0);\r
 916                         \r
 917       Double_t massBG =0.;\r
 918       Double_t widthBG = 0.;\r
 919       Double_t chi2BG =10000.;  \r
 920       backgroundCandidate->GetMass(massBG,widthBG);\r
 921       if(backgroundCandidate->GetNDF()>0){\r
 922         chi2BG = backgroundCandidate->GetChi2()/backgroundCandidate->GetNDF();\r
 923         if(chi2BG>0 && chi2BG<fV0Reader->GetChi2CutMeson()){\r
 924                                         \r
 925           TVector3 MomentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());\r
 926           TVector3 SpaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());\r
 927                                         \r
 928           Double_t openingAngleBG = currentEventGoodV0->GetAngle(*previousGoodV0);\r
 929 \r
 930           Double_t rapidity;\r
 931           if(backgroundCandidate->GetE() - backgroundCandidate->GetPz() == 0 || backgroundCandidate->GetE() + backgroundCandidate->GetPz() == 0) rapidity=0;\r
 932           else rapidity = 0.5*(TMath::Log((backgroundCandidate->GetE() +backgroundCandidate->GetPz()) / (backgroundCandidate->GetE()-backgroundCandidate->GetPz())));\r
 933 \r
 934                                         \r
 935                                         \r
 936                                         \r
 937           if(openingAngleBG < fMinOpeningAngleGhostCut ) continue;   // minimum opening angle to avoid using ghosttracks\r
 938                         \r
 939                                         \r
 940           fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle", openingAngleBG);\r
 941           fHistograms->FillHistogram("ESD_Background_Energy", backgroundCandidate->GetE());\r
 942           fHistograms->FillHistogram("ESD_Background_Pt",  MomentumVectorbackgroundCandidate.Pt());\r
 943           fHistograms->FillHistogram("ESD_Background_Eta", MomentumVectorbackgroundCandidate.Eta());\r
 944           fHistograms->FillHistogram("ESD_Background_Rapidity", rapidity);\r
 945           fHistograms->FillHistogram("ESD_Background_Phi", SpaceVectorbackgroundCandidate.Phi());\r
 946           fHistograms->FillHistogram("ESD_Background_Mass", massBG);\r
 947           fHistograms->FillHistogram("ESD_Background_R", SpaceVectorbackgroundCandidate.Pt());  // Pt in Space == R!!!!\r
 948           fHistograms->FillHistogram("ESD_Background_ZR", backgroundCandidate->GetZ(), SpaceVectorbackgroundCandidate.Pt());\r
 949           fHistograms->FillHistogram("ESD_Background_XY", backgroundCandidate->GetX(), backgroundCandidate->GetY());\r
 950           fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt",massBG,MomentumVectorbackgroundCandidate.Pt());\r
 951         }\r
 952       }\r
 953       delete backgroundCandidate;   \r
 954       //cout << "nearly at the end of background" << endl;\r
 955                         \r
 956     }\r
 957   }\r
 958 }\r
 959 \r
 960 void AliAnalysisTaskGammaConversion::Terminate(Option_t */*option*/)\r
 961 {\r
 962   // Terminate analysis\r
 963   //\r
 964   AliDebug(1,"Do nothing in Terminate");\r
 965 }\r
 966 \r
 967 void AliAnalysisTaskGammaConversion::UserCreateOutputObjects()\r
 968 {\r
 969   // Create the output container\r
 970   if(fOutputContainer != NULL){\r
 971     delete fOutputContainer;\r
 972     fOutputContainer = NULL;\r
 973   }\r
 974   if(fOutputContainer == NULL){\r
 975     fOutputContainer = new TList();\r
 976   }\r
 977         \r
 978   //Adding the histograms to the output container\r
 979   fHistograms->GetOutputContainer(fOutputContainer);\r
 980         \r
 981         \r
 982   if(fWriteNtuple){\r
 983     if(fGammaNtuple == NULL){\r
 984       fGammaNtuple = new TNtuple("V0ntuple","V0ntuple","OnTheFly:HasVertex:NegPIDProb:PosPIDProb:X:Y:Z:R:MotherCandidateNDF:MotherCandidateChi2:MotherCandidateEnergy:MotherCandidateEta:MotherCandidatePt:MotherCandidateMass:MotherCandidateWidth:MCMotherCandidatePT:EPOpeningAngle:ElectronEnergy:ElectronPt:ElectronEta:ElectronPhi:PositronEnergy:PositronPt:PositronEta:PositronPhi:HasSameMCMother:MotherMCParticlePIDCode",50000);\r
 985     }\r
 986     if(fNeutralMesonNtuple == NULL){\r
 987       fNeutralMesonNtuple = new TNtuple("NeutralMesonNtuple","NeutralMesonNtuple","test");\r
 988     }\r
 989     TList * ntupleTList = new TList();\r
 990     ntupleTList->SetName("Ntuple");\r
 991     ntupleTList->Add((TNtuple*)fGammaNtuple);\r
 992     fOutputContainer->Add(ntupleTList);\r
 993   }\r
 994         \r
 995   fOutputContainer->SetName(GetName());\r
 996 }\r
 997 \r
 998 Double_t AliAnalysisTaskGammaConversion::GetMCOpeningAngle(TParticle* daughter0, TParticle* daughter1) const{\r
 999   //helper function\r
1000   TVector3 v3D0(daughter0->Px(),daughter0->Py(),daughter0->Pz());\r
1001   TVector3 v3D1(daughter1->Px(),daughter1->Py(),daughter1->Pz());\r
1002   return v3D0.Angle(v3D1);\r
1003 }\r