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 "TNtuple.h"\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   fESDEvent(NULL),      \r
  53   fOutputContainer(NULL),\r
  54   fHistograms(NULL),\r
  55   fDoMCTruth(kFALSE),\r
  56   fMCAllGammas(),\r
  57   fMCPi0s(),\r
  58   fMCEtas(),\r
  59   fMCGammaChic(),\r
  60   fKFReconstructedGammas(),\r
  61   fIsTrueReconstructedGammas(),\r
  62   fElectronv1(),\r
  63   fElectronv2(),\r
  64   fCurrentEventPosElectron(),\r
  65   fPreviousEventPosElectron(),\r
  66   fCurrentEventNegElectron(),\r
  67   fPreviousEventNegElectron(),\r
  68   fKFReconstructedGammasCut(),                  \r
  69   fPreviousEventTLVNegElectron(),\r
  70   fPreviousEventTLVPosElectron(),       \r
  71   fElectronMass(-1),\r
  72   fGammaMass(-1),\r
  73   fPi0Mass(-1),\r
  74   fEtaMass(-1),\r
  75   fGammaWidth(-1),\r
  76   fPi0Width(-1),\r
  77   fEtaWidth(-1),\r
  78   fMinOpeningAngleGhostCut(0.),\r
  79   fCalculateBackground(kFALSE),\r
  80   fWriteNtuple(kFALSE),\r
  81   fGammaNtuple(NULL),\r
  82   fNeutralMesonNtuple(NULL),\r
  83   fTotalNumberOfAddedNtupleEntries(0)\r
  84 {\r
  85   // Default constructor\r
  86   // Common I/O in slot 0\r
  87   DefineInput (0, TChain::Class());\r
  88   DefineOutput(0, TTree::Class());\r
  89         \r
  90   // Your private output\r
  91   DefineOutput(1, TList::Class());\r
  92 }\r
  93 \r
  94 AliAnalysisTaskGammaConversion::AliAnalysisTaskGammaConversion(const char* name):\r
  95   AliAnalysisTaskSE(name),\r
  96   fV0Reader(NULL),\r
  97   fStack(NULL),\r
  98   fESDEvent(NULL),      \r
  99   fOutputContainer(0x0),\r
 100   fHistograms(NULL),\r
 101   fDoMCTruth(kFALSE),\r
 102   fMCAllGammas(),\r
 103   fMCPi0s(),\r
 104   fMCEtas(),\r
 105   fMCGammaChic(),\r
 106   fKFReconstructedGammas(),\r
 107   fIsTrueReconstructedGammas(),\r
 108   fElectronv1(),\r
 109   fElectronv2(),\r
 110   fCurrentEventPosElectron(),\r
 111   fPreviousEventPosElectron(),\r
 112   fCurrentEventNegElectron(),\r
 113   fPreviousEventNegElectron(),\r
 114   fKFReconstructedGammasCut(),  \r
 115   fPreviousEventTLVNegElectron(),\r
 116   fPreviousEventTLVPosElectron(),\r
 117   fElectronMass(-1),\r
 118   fGammaMass(-1),\r
 119   fPi0Mass(-1),\r
 120   fEtaMass(-1),\r
 121   fGammaWidth(-1),\r
 122   fPi0Width(-1),\r
 123   fEtaWidth(-1),\r
 124   fMinOpeningAngleGhostCut(0.),\r
 125   fCalculateBackground(kFALSE),\r
 126   fWriteNtuple(kFALSE),\r
 127   fGammaNtuple(NULL),\r
 128   fNeutralMesonNtuple(NULL),\r
 129   fTotalNumberOfAddedNtupleEntries(0)\r
 130 {\r
 131   // Common I/O in slot 0\r
 132   DefineInput (0, TChain::Class());\r
 133   DefineOutput(0, TTree::Class());\r
 134         \r
 135   // Your private output\r
 136   DefineOutput(1, TList::Class());\r
 137 }\r
 138 \r
 139 AliAnalysisTaskGammaConversion::~AliAnalysisTaskGammaConversion() \r
 140 {\r
 141   // Remove all pointers\r
 142         \r
 143   if(fOutputContainer){\r
 144     fOutputContainer->Clear() ; \r
 145     delete fOutputContainer ;\r
 146   }\r
 147   if(fHistograms){\r
 148     delete fHistograms;\r
 149   }\r
 150   if(fV0Reader){\r
 151     delete fV0Reader;\r
 152   }\r
 153 }\r
 154 \r
 155 \r
 156 void AliAnalysisTaskGammaConversion::Init()\r
 157 {\r
 158   // Initialization\r
 159   AliLog::SetGlobalLogLevel(AliLog::kError);\r
 160 }\r
 161 \r
 162 \r
 163 void AliAnalysisTaskGammaConversion::Exec(Option_t */*option*/)\r
 164 {\r
 165   // Execute analysis for current event\r
 166         \r
 167   ConnectInputData("");\r
 168         \r
 169   //clear vectors\r
 170   fMCAllGammas.clear();\r
 171   fMCPi0s.clear();\r
 172   fMCEtas.clear();\r
 173   fMCGammaChic.clear();\r
 174         \r
 175   fKFReconstructedGammas.clear();\r
 176   fIsTrueReconstructedGammas.clear();\r
 177   fElectronv1.clear();\r
 178   fElectronv2.clear();\r
 179   fCurrentEventPosElectron.clear();\r
 180   fCurrentEventNegElectron.clear();     \r
 181   fKFReconstructedGammasCut.clear(); \r
 182         \r
 183   //Clear the data in the v0Reader\r
 184   fV0Reader->UpdateEventByEventData();\r
 185 \r
 186   \r
 187   // Process the MC information\r
 188   if(fDoMCTruth){\r
 189     ProcessMCData();\r
 190   }\r
 191   \r
 192   //Process the v0 information with no cuts\r
 193   ProcessV0sNoCut();\r
 194   \r
 195   // Process the v0 information\r
 196   ProcessV0s();\r
 197   \r
 198   //calculate background if flag is set\r
 199   if(fCalculateBackground){\r
 200     CalculateBackground();\r
 201   }\r
 202   \r
 203   // Process reconstructed gammas\r
 204   ProcessGammasForNeutralMesonAnalysis();\r
 205 \r
 206   CheckV0Efficiency();\r
 207   \r
 208   //Process reconstructed gammas electrons for Chi_c Analysis\r
 209   ProcessGammaElectronsForChicAnalysis();\r
 210 \r
 211   PostData(1, fOutputContainer);\r
 212         \r
 213 }\r
 214 \r
 215 void AliAnalysisTaskGammaConversion::ConnectInputData(Option_t */*option*/){\r
 216   // see header file for documentation\r
 217         \r
 218   if(fV0Reader == NULL){\r
 219     // Write warning here cuts and so on are default if this ever happens\r
 220   }\r
 221   fV0Reader->Initialize();\r
 222 }\r
 223 \r
 224 \r
 225 \r
 226 void AliAnalysisTaskGammaConversion::ProcessMCData(){\r
 227   // see header file for documentation\r
 228         \r
 229   fStack = fV0Reader->GetMCStack();\r
 230 \r
 231   if(fV0Reader->CheckForPrimaryVertex() == kFALSE){\r
 232     return; // aborts if the primary vertex does not have contributors.\r
 233   }\r
 234 \r
 235   for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++) {\r
 236     TParticle* particle = (TParticle *)fStack->Particle(iTracks);\r
 237                 \r
 238     if (!particle) {\r
 239       //print warning here\r
 240       continue;\r
 241     }\r
 242 \r
 243     ///////////////////////Begin Chic Analysis/////////////////////////////\r
 244 \r
 245 \r
 246     if(particle->GetPdgCode() == 443){//Is JPsi\r
 247 \r
 248       if(particle->GetNDaughters()==2){\r
 249         if(TMath::Abs(fStack->Particle(particle->GetFirstDaughter())->GetPdgCode()) == 11 &&\r
 250            TMath::Abs(fStack->Particle(particle->GetLastDaughter())->GetPdgCode()) == 11){\r
 251           TParticle* daug0 = fStack->Particle(particle->GetFirstDaughter());\r
 252           TParticle* daug1 = fStack->Particle(particle->GetLastDaughter());\r
 253           if(TMath::Abs(daug0->Eta()) < 0.9 && TMath::Abs(daug1->Eta()) < 0.9)\r
 254             fHistograms->FillTable("Table_Electrons",3);//e+ e-  from J/Psi inside acceptance\r
 255 \r
 256           if( TMath::Abs(daug0->Eta()) < 0.9){\r
 257             if(daug0->GetPdgCode() == -11)\r
 258               fHistograms->FillTable("Table_Electrons",1);//e+  from J/Psi inside acceptance\r
 259             else\r
 260               fHistograms->FillTable("Table_Electrons",2);//e-   from J/Psi inside acceptance\r
 261 \r
 262           }\r
 263           if(TMath::Abs(daug1->Eta()) < 0.9){\r
 264             if(daug1->GetPdgCode() == -11)\r
 265               fHistograms->FillTable("Table_Electrons",1);//e+  from J/Psi inside acceptance\r
 266             else\r
 267               fHistograms->FillTable("Table_Electrons",2);//e-   from J/Psi inside acceptance\r
 268           }\r
 269         }\r
 270       }\r
 271     }\r
 272     //              const int CHI_C0   = 10441;\r
 273     //              const int CHI_C1   = 20443;\r
 274     //              const int CHI_C2   = 445\r
 275     if(particle->GetPdgCode() == 22){//gamma from JPsi\r
 276       if(particle->GetMother(0) > -1){\r
 277         if(fStack->Particle(particle->GetMother(0))->GetPdgCode() == 10441 ||\r
 278            fStack->Particle(particle->GetMother(0))->GetPdgCode() == 20443 ||\r
 279            fStack->Particle(particle->GetMother(0))->GetPdgCode() == 445){\r
 280           if(TMath::Abs(particle->Eta()) < 1.2)\r
 281             fHistograms->FillTable("Table_Electrons",17);// gamma from chic inside accptance\r
 282         }\r
 283       }\r
 284     }\r
 285     if(particle->GetPdgCode() == 10441 || particle->GetPdgCode() == 20443 || particle->GetPdgCode() == 445){\r
 286       if( particle->GetNDaughters() == 2){\r
 287         TParticle* daug0 = fStack->Particle(particle->GetFirstDaughter());\r
 288         TParticle* daug1 = fStack->Particle(particle->GetLastDaughter());\r
 289 \r
 290         if( (daug0->GetPdgCode() == 443 || daug0->GetPdgCode() == 22) && (daug1->GetPdgCode() == 443 || daug1->GetPdgCode() == 22) ){\r
 291           if( daug0->GetPdgCode() == 443){\r
 292             TParticle* daugE0 = fStack->Particle(daug0->GetFirstDaughter());\r
 293             TParticle* daugE1 = fStack->Particle(daug0->GetLastDaughter());\r
 294             if( TMath::Abs(daug1->Eta()) < 1.2 && TMath::Abs(daugE0->Eta()) < 0.9 && TMath::Abs(daugE1->Eta()) < 0.9 )\r
 295               fHistograms->FillTable("Table_Electrons",18);\r
 296 \r
 297           }//if\r
 298           else if (daug1->GetPdgCode() == 443){\r
 299             TParticle* daugE0 = fStack->Particle(daug1->GetFirstDaughter());\r
 300             TParticle* daugE1 = fStack->Particle(daug1->GetLastDaughter());\r
 301             if( TMath::Abs(daug0->Eta()) < 1.2 && TMath::Abs(daugE0->Eta()) < 0.9 && TMath::Abs(daugE1->Eta()) < 0.9 )\r
 302               fHistograms->FillTable("Table_Electrons",18);\r
 303           }//else if\r
 304         }//gamma o Jpsi\r
 305       }//GetNDaughters\r
 306     }\r
 307 \r
 308 \r
 309     /////////////////////End Chic Analysis////////////////////////////\r
 310 \r
 311 \r
 312     if(TMath::Abs(particle->Eta())> fV0Reader->GetEtaCut() )    continue;\r
 313                                         \r
 314     if(particle->R()>fV0Reader->GetMaxRCut())   continue; // cuts on distance from collision point\r
 315                 \r
 316     Double_t tmpPhi=particle->Phi();\r
 317                 \r
 318     if(particle->Phi()> TMath::Pi()){\r
 319       tmpPhi = particle->Phi()-(2*TMath::Pi());\r
 320     }\r
 321                 \r
 322     Double_t rapidity;\r
 323     if(particle->Energy() - particle->Pz() == 0 || particle->Energy() + particle->Pz() == 0){\r
 324       rapidity=0;\r
 325     }\r
 326     else{\r
 327       rapidity = 0.5*(TMath::Log((particle->Energy()+particle->Pz()) / (particle->Energy()-particle->Pz())));\r
 328     }   \r
 329                 \r
 330     //process the gammas\r
 331     if (particle->GetPdgCode() == 22){\r
 332                         \r
 333       if(particle->GetMother(0) >-1 && fStack->Particle(particle->GetMother(0))->GetPdgCode() == 22){\r
 334         continue; // no photon as mothers!\r
 335       }\r
 336 \r
 337       if(particle->GetMother(0) >= fStack->GetNprimary()){\r
 338         continue; // the gamma has a mother, and it is not a primary particle\r
 339       }\r
 340 \r
 341       fMCAllGammas.push_back(particle);\r
 342                         \r
 343       fHistograms->FillHistogram("MC_allGamma_Energy", particle->Energy());\r
 344       fHistograms->FillHistogram("MC_allGamma_Pt", particle->Pt());\r
 345       fHistograms->FillHistogram("MC_allGamma_Eta", particle->Eta());\r
 346       fHistograms->FillHistogram("MC_allGamma_Phi", tmpPhi);\r
 347       fHistograms->FillHistogram("MC_allGamma_Rapid", rapidity);\r
 348                         \r
 349                         \r
 350       if(particle->GetMother(0) < 0){   // direct gamma\r
 351         fHistograms->FillHistogram("MC_allDirectGamma_Energy",particle->Energy());\r
 352         fHistograms->FillHistogram("MC_allDirectGamma_Pt", particle->Pt());\r
 353         fHistograms->FillHistogram("MC_allDirectGamma_Eta", particle->Eta());\r
 354         fHistograms->FillHistogram("MC_allDirectGamma_Phi", tmpPhi);\r
 355         fHistograms->FillHistogram("MC_allDirectGamma_Rapid", rapidity);                                \r
 356       }\r
 357                         \r
 358                         \r
 359       // looking for conversion (electron + positron from pairbuilding (= 5) )\r
 360       TParticle* ePos = NULL;\r
 361       TParticle* eNeg = NULL;\r
 362                         \r
 363       if(particle->GetNDaughters() >= 2){\r
 364         for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){\r
 365           TParticle *tmpDaughter = fStack->Particle(daughterIndex);\r
 366           if(tmpDaughter->GetUniqueID() == 5){\r
 367             if(tmpDaughter->GetPdgCode() == 11){\r
 368               eNeg = tmpDaughter;\r
 369             }\r
 370             else if(tmpDaughter->GetPdgCode() == -11){\r
 371               ePos = tmpDaughter;\r
 372             }\r
 373           }\r
 374         }\r
 375       }\r
 376                         \r
 377                         \r
 378       if(ePos == NULL || eNeg == NULL){ // means we do not have two daughters from pair production\r
 379         continue;\r
 380       }\r
 381                         \r
 382                         \r
 383       Double_t ePosPhi = ePos->Phi();\r
 384       if(ePos->Phi()> TMath::Pi()) ePosPhi = ePos->Phi()-(2*TMath::Pi());\r
 385                         \r
 386       Double_t eNegPhi = eNeg->Phi();\r
 387       if(eNeg->Phi()> TMath::Pi()) eNegPhi = eNeg->Phi()-(2*TMath::Pi());\r
 388                         \r
 389                         \r
 390       if(ePos->Pt()<fV0Reader->GetPtCut() || eNeg->Pt()<fV0Reader->GetPtCut()){\r
 391         continue; // no reconstruction below the Pt cut\r
 392       }\r
 393                                         \r
 394       if(TMath::Abs(ePos->Eta())> fV0Reader->GetEtaCut() || TMath::Abs(eNeg->Eta())> fV0Reader->GetEtaCut()){\r
 395         continue;\r
 396       } \r
 397                                 \r
 398       if(ePos->R()>fV0Reader->GetMaxRCut()){\r
 399         continue; // cuts on distance from collision point\r
 400       }\r
 401       \r
 402       \r
 403       if((TMath::Abs(ePos->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()  > ePos->R()){\r
 404         continue;               // line cut to exclude regions where we do not reconstruct\r
 405       }         \r
 406                 \r
 407       fHistograms->FillHistogram("MC_ConvGamma_Energy", particle->Energy());\r
 408       fHistograms->FillHistogram("MC_ConvGamma_Pt", particle->Pt());\r
 409       fHistograms->FillHistogram("MC_ConvGamma_Eta", particle->Eta());\r
 410       fHistograms->FillHistogram("MC_ConvGamma_Phi", tmpPhi);\r
 411       fHistograms->FillHistogram("MC_ConvGamma_Rapid", rapidity);\r
 412       fHistograms->FillHistogram("MC_ConvGamma_Pt_Eta", particle->Pt(),particle->Eta());\r
 413                         \r
 414       fHistograms->FillHistogram("MC_E_Energy", eNeg->Energy());\r
 415       fHistograms->FillHistogram("MC_E_Pt", eNeg->Pt());\r
 416       fHistograms->FillHistogram("MC_E_Eta", eNeg->Eta());\r
 417       fHistograms->FillHistogram("MC_E_Phi", eNegPhi);\r
 418                         \r
 419       fHistograms->FillHistogram("MC_P_Energy", ePos->Energy());\r
 420       fHistograms->FillHistogram("MC_P_Pt", ePos->Pt());\r
 421       fHistograms->FillHistogram("MC_P_Eta", ePos->Eta());\r
 422       fHistograms->FillHistogram("MC_P_Phi", ePosPhi);\r
 423                         \r
 424                         \r
 425                         \r
 426       //cout << "filled histos for converted gamma, ePos, eNeg" << endl;\r
 427                         \r
 428       // begin Mapping \r
 429       Int_t rBin    = fHistograms->GetRBin(ePos->R());\r
 430       Int_t phiBin  = fHistograms->GetPhiBin(particle->Phi());\r
 431                         \r
 432       TString nameMCMappingPhiR="";\r
 433       nameMCMappingPhiR.Form("MC_Conversion_Mapping-Phi%02d-R%02d",phiBin,rBin);\r
 434       fHistograms->FillHistogram(nameMCMappingPhiR, ePos->Vz(), particle->Eta());\r
 435                         \r
 436       TString nameMCMappingPhi="";\r
 437       nameMCMappingPhi.Form("MC_Conversion_Mapping-Phi%02d",phiBin);\r
 438       fHistograms->FillHistogram(nameMCMappingPhi, particle->Eta());\r
 439                         \r
 440       TString nameMCMappingR="";\r
 441       nameMCMappingR.Form("MC_Conversion_Mapping-R%02d",rBin);\r
 442       fHistograms->FillHistogram(nameMCMappingR, particle->Eta());\r
 443                         \r
 444       TString nameMCMappingPhiInR="";\r
 445       nameMCMappingPhiInR.Form("MC_Conversion_Mapping_Phi_R-%02d",rBin);\r
 446       fHistograms->FillHistogram(nameMCMappingPhiInR, tmpPhi);\r
 447       //end mapping\r
 448                         \r
 449       fHistograms->FillHistogram("MC_Conversion_R",ePos->R());\r
 450       fHistograms->FillHistogram("MC_Conversion_ZR",ePos->Vz(),ePos->R());\r
 451       fHistograms->FillHistogram("MC_Conversion_XY",ePos->Vx(),ePos->Vy());\r
 452       fHistograms->FillHistogram("MC_Conversion_OpeningAngle",GetMCOpeningAngle(ePos, eNeg));\r
 453                         \r
 454       //cout << "mapping is done" << endl;\r
 455                         \r
 456                         \r
 457       if(particle->GetMother(0) < 0){ // no mother = direct gamma, still inside converted\r
 458         fHistograms->FillHistogram("MC_ConvDirectGamma_Energy",particle->Energy());\r
 459         fHistograms->FillHistogram("MC_ConvDirectGamma_Pt", particle->Pt());\r
 460         fHistograms->FillHistogram("MC_ConvDirectGamma_Eta", particle->Eta());\r
 461         fHistograms->FillHistogram("MC_ConvDirectGamma_Phi", tmpPhi);\r
 462         fHistograms->FillHistogram("MC_ConvDirectGamma_Rapid", rapidity);\r
 463                                 \r
 464       } // end direct gamma\r
 465       else{   // mother exits \r
 466         if( fStack->Particle(particle->GetMother(0))->GetPdgCode()==10441 ||//chic0 \r
 467             fStack->Particle(particle->GetMother(0))->GetPdgCode()==20443 ||//psi2S\r
 468             fStack->Particle(particle->GetMother(0))->GetPdgCode()==445  //chic2\r
 469             ){ \r
 470           fMCGammaChic.push_back(particle);\r
 471         }\r
 472       }  // end if mother exits\r
 473     } // end if particle is a photon\r
 474                 \r
 475     if(particle->GetNDaughters() == 2){\r
 476                         \r
 477       TParticle* daughter0 = (TParticle*)fStack->Particle(particle->GetFirstDaughter());\r
 478       TParticle* daughter1 = (TParticle*)fStack->Particle(particle->GetLastDaughter());\r
 479                         \r
 480       if(daughter0->GetPdgCode() != 22 || daughter1->GetPdgCode() != 22) continue; //check for gamma gamma daughters\r
 481                         \r
 482                         \r
 483                         \r
 484       // check for conversions now -> have to pass eta and line cut!\r
 485       Bool_t daughter0Electron = kFALSE;\r
 486       Bool_t daughter0Positron = kFALSE;\r
 487       Bool_t daughter1Electron = kFALSE;\r
 488       Bool_t daughter1Positron = kFALSE;\r
 489                         \r
 490                         \r
 491                         \r
 492       if(daughter0->GetNDaughters() >= 2){\r
 493         for(Int_t TrackIndex=daughter0->GetFirstDaughter();TrackIndex<=daughter0->GetLastDaughter();TrackIndex++){\r
 494           TParticle *tmpDaughter = fStack->Particle(TrackIndex);\r
 495           if(tmpDaughter->GetUniqueID() == 5){\r
 496             if(tmpDaughter->GetPdgCode() == 11){\r
 497               if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){\r
 498                 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()  < tmpDaughter->R() ){\r
 499                   daughter0Electron = kTRUE;\r
 500                 }\r
 501                                                                 \r
 502               }\r
 503             }\r
 504             else if(tmpDaughter->GetPdgCode() == -11){\r
 505               if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){\r
 506                 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()  < tmpDaughter->R() ){\r
 507                   daughter0Positron = kTRUE;\r
 508                 }\r
 509                                                                 \r
 510               }\r
 511                                                         \r
 512             }\r
 513           }\r
 514         }\r
 515       }\r
 516                         \r
 517                         \r
 518         \r
 519       if(daughter1->GetNDaughters() >= 2){\r
 520         for(Int_t TrackIndex=daughter1->GetFirstDaughter();TrackIndex<=daughter1->GetLastDaughter();TrackIndex++){\r
 521           TParticle *tmpDaughter = fStack->Particle(TrackIndex);\r
 522           if(tmpDaughter->GetUniqueID() == 5){\r
 523             if(tmpDaughter->GetPdgCode() == 11){\r
 524               if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){\r
 525                 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()  < tmpDaughter->R() ){\r
 526                   daughter1Electron = kTRUE;\r
 527                 }\r
 528                                                                 \r
 529               }\r
 530             }\r
 531             else if(tmpDaughter->GetPdgCode() == -11){\r
 532               if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){\r
 533                 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()  < tmpDaughter->R() ){\r
 534                   daughter1Positron = kTRUE;\r
 535                 }\r
 536                                                                 \r
 537               }\r
 538                                                         \r
 539             }\r
 540           }\r
 541         }\r
 542       }\r
 543                         \r
 544                                                                                                 \r
 545                         \r
 546                         \r
 547       if(particle->GetPdgCode()==111){     //Pi0\r
 548         if( iTracks >= fStack->GetNprimary()){\r
 549           fHistograms->FillHistogram("MC_Pi0_Secondaries_Eta", particle->Eta());\r
 550           fHistograms->FillHistogram("MC_Pi0_Secondaries_Rapid", rapidity);\r
 551           fHistograms->FillHistogram("MC_Pi0_Secondaries_Phi", tmpPhi);\r
 552           fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt", particle->Pt());\r
 553           fHistograms->FillHistogram("MC_Pi0_Secondaries_Energy", particle->Energy());\r
 554           fHistograms->FillHistogram("MC_Pi0_Secondaries_R", particle->R());\r
 555           fHistograms->FillHistogram("MC_Pi0_Secondaries_ZR", particle->Vz(),particle->R());\r
 556           fHistograms->FillHistogram("MC_Pi0_Secondaries_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
 557           fHistograms->FillHistogram("MC_Pi0_Secondaries_XY", particle->Vx(),particle->Vy());//only fill from one daughter to avoid multiple filling\r
 558                                         \r
 559           if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){\r
 560             fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());\r
 561             fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);\r
 562             if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){\r
 563               fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());\r
 564               fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);\r
 565             }\r
 566           }\r
 567         }\r
 568         else{\r
 569           fHistograms->FillHistogram("MC_Pi0_Eta", particle->Eta());    \r
 570           fHistograms->FillHistogram("MC_Pi0_Rapid", rapidity);\r
 571           fHistograms->FillHistogram("MC_Pi0_Phi", tmpPhi);\r
 572           fHistograms->FillHistogram("MC_Pi0_Pt", particle->Pt());\r
 573           fHistograms->FillHistogram("MC_Pi0_Energy", particle->Energy());\r
 574           fHistograms->FillHistogram("MC_Pi0_R", particle->R());\r
 575           fHistograms->FillHistogram("MC_Pi0_ZR", particle->Vz(),particle->R());\r
 576           fHistograms->FillHistogram("MC_Pi0_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
 577           fHistograms->FillHistogram("MC_Pi0_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling\r
 578                                         \r
 579           if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){\r
 580             fHistograms->FillHistogram("MC_Pi0_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());\r
 581             fHistograms->FillHistogram("MC_Pi0_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);\r
 582             if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){\r
 583               fHistograms->FillHistogram("MC_Pi0_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());\r
 584               fHistograms->FillHistogram("MC_Pi0_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);\r
 585             }\r
 586           }\r
 587         }\r
 588       }\r
 589                         \r
 590       if(particle->GetPdgCode()==221){   //Eta\r
 591         fHistograms->FillHistogram("MC_Eta_Eta", particle->Eta());\r
 592         fHistograms->FillHistogram("MC_Eta_Rapid", rapidity);\r
 593         fHistograms->FillHistogram("MC_Eta_Phi",tmpPhi);\r
 594         fHistograms->FillHistogram("MC_Eta_Pt", particle->Pt());\r
 595         fHistograms->FillHistogram("MC_Eta_Energy", particle->Energy());\r
 596         fHistograms->FillHistogram("MC_Eta_R", particle->R());\r
 597         fHistograms->FillHistogram("MC_Eta_ZR", particle->Vz(),particle->R());\r
 598         fHistograms->FillHistogram("MC_Eta_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
 599         fHistograms->FillHistogram("MC_Eta_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling\r
 600                                 \r
 601         if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){\r
 602           fHistograms->FillHistogram("MC_Eta_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());\r
 603           fHistograms->FillHistogram("MC_Eta_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);\r
 604           if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){\r
 605             fHistograms->FillHistogram("MC_Eta_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());\r
 606             fHistograms->FillHistogram("MC_Eta_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);\r
 607           }\r
 608                                         \r
 609         }\r
 610                                 \r
 611       }\r
 612                         \r
 613       // all motherparticles with 2 gammas as daughters\r
 614       fHistograms->FillHistogram("MC_Mother_R", particle->R());\r
 615       fHistograms->FillHistogram("MC_Mother_ZR", particle->Vz(),particle->R());\r
 616       fHistograms->FillHistogram("MC_Mother_XY", particle->Vx(),particle->Vy());\r
 617       fHistograms->FillHistogram("MC_Mother_Mass", particle->GetCalcMass());\r
 618       fHistograms->FillHistogram("MC_Mother_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
 619       fHistograms->FillHistogram("MC_Mother_Energy", particle->Energy());\r
 620       fHistograms->FillHistogram("MC_Mother_Pt", particle->Pt());\r
 621       fHistograms->FillHistogram("MC_Mother_Eta", particle->Eta());\r
 622       fHistograms->FillHistogram("MC_Mother_Rapid", rapidity);\r
 623       fHistograms->FillHistogram("MC_Mother_Phi",tmpPhi);\r
 624       fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt",particle->GetMass(),particle->Pt());                 \r
 625       if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){\r
 626         fHistograms->FillHistogram("MC_Mother_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());\r
 627         fHistograms->FillHistogram("MC_Mother_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);\r
 628         fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt_withinAcceptance",particle->GetMass(),particle->Pt());                      \r
 629         if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){\r
 630           fHistograms->FillHistogram("MC_Mother_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());\r
 631           fHistograms->FillHistogram("MC_Mother_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);\r
 632           fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt_ConvGamma_withinAcceptance",particle->GetMass(),particle->Pt());                  \r
 633 \r
 634         }\r
 635                                 \r
 636                                 \r
 637       }\r
 638                         \r
 639       //cout << "mother histos are filled" << endl;\r
 640                         \r
 641     } // end if(particle->GetNDaughters() == 2)\r
 642                 \r
 643   }// end for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++)\r
 644         \r
 645   //cout << "right before the end of processMCdata" << endl;\r
 646         \r
 647 } // end ProcessMCData\r
 648 \r
 649 \r
 650 \r
 651 void AliAnalysisTaskGammaConversion::FillNtuple(){\r
 652   //Fills the ntuple with the different values\r
 653 \r
 654   if(fGammaNtuple == NULL){\r
 655     return;\r
 656   }\r
 657   Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();\r
 658   for(Int_t i=0;i<numberOfV0s;i++){\r
 659     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
 660     AliESDv0 * cV0 = fV0Reader->GetV0(i);\r
 661     Double_t negPID=0;\r
 662     Double_t posPID=0;\r
 663     fV0Reader->GetPIDProbability(negPID,posPID);\r
 664     values[0]=cV0->GetOnFlyStatus();\r
 665     values[1]=fV0Reader->CheckForPrimaryVertex();\r
 666     values[2]=negPID;\r
 667     values[3]=posPID;\r
 668     values[4]=fV0Reader->GetX();\r
 669     values[5]=fV0Reader->GetY();\r
 670     values[6]=fV0Reader->GetZ();\r
 671     values[7]=fV0Reader->GetXYRadius();\r
 672     values[8]=fV0Reader->GetMotherCandidateNDF();\r
 673     values[9]=fV0Reader->GetMotherCandidateChi2();\r
 674     values[10]=fV0Reader->GetMotherCandidateEnergy();\r
 675     values[11]=fV0Reader->GetMotherCandidateEta();\r
 676     values[12]=fV0Reader->GetMotherCandidatePt();\r
 677     values[13]=fV0Reader->GetMotherCandidateMass();\r
 678     values[14]=fV0Reader->GetMotherCandidateWidth();\r
 679     //    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
 680     values[16]=fV0Reader->GetOpeningAngle();\r
 681     values[17]=fV0Reader->GetNegativeTrackEnergy();\r
 682     values[18]=fV0Reader->GetNegativeTrackPt();\r
 683     values[19]=fV0Reader->GetNegativeTrackEta();\r
 684     values[20]=fV0Reader->GetNegativeTrackPhi();\r
 685     values[21]=fV0Reader->GetPositiveTrackEnergy();\r
 686     values[22]=fV0Reader->GetPositiveTrackPt();\r
 687     values[23]=fV0Reader->GetPositiveTrackEta();\r
 688     values[24]=fV0Reader->GetPositiveTrackPhi();\r
 689     values[25]=fV0Reader->HasSameMCMother();\r
 690     if(values[25] != 0){\r
 691       values[26]=fV0Reader->GetMotherMCParticlePDGCode();\r
 692       values[15]=fV0Reader->GetMotherMCParticle()->Pt();\r
 693     }\r
 694     fTotalNumberOfAddedNtupleEntries++;\r
 695     fGammaNtuple->Fill(values);\r
 696   }\r
 697   fV0Reader->ResetV0IndexNumber();\r
 698         \r
 699 }\r
 700 \r
 701 void AliAnalysisTaskGammaConversion::ProcessV0sNoCut(){\r
 702   // Process all the V0's without applying any cuts to it\r
 703 \r
 704   Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();\r
 705   for(Int_t i=0;i<numberOfV0s;i++){\r
 706     /*AliESDv0 * cV0 = */fV0Reader->GetV0(i);\r
 707 \r
 708     if(fV0Reader->CheckForPrimaryVertex() == kFALSE){\r
 709       return;\r
 710     }\r
 711     \r
 712     if(fDoMCTruth){\r
 713       \r
 714       if(fV0Reader->HasSameMCMother() == kFALSE){\r
 715         continue;\r
 716       }\r
 717                 \r
 718       TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();\r
 719       TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();\r
 720 \r
 721       if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){\r
 722         continue;\r
 723       }\r
 724       if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){\r
 725         continue;\r
 726       }\r
 727         \r
 728       if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){\r
 729       \r
 730         fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt", fV0Reader->GetMotherCandidatePt());\r
 731         fHistograms->FillHistogram("ESD_NoCutConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());\r
 732         fHistograms->FillHistogram("ESD_NoCutConvGamma_Eta", fV0Reader->GetMotherCandidateEta());                               \r
 733         fHistograms->FillHistogram("ESD_NoCutConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());\r
 734         fHistograms->FillHistogram("ESD_NoCutConvGamma_Mass", fV0Reader->GetMotherCandidateMass());\r
 735         fHistograms->FillHistogram("ESD_NoCutConvGamma_Width", fV0Reader->GetMotherCandidateWidth());\r
 736         fHistograms->FillHistogram("ESD_NoCutConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());\r
 737         fHistograms->FillHistogram("ESD_NoCutConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());\r
 738         fHistograms->FillHistogram("ESD_NoCutConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());\r
 739         fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());\r
 740         \r
 741         fHistograms->FillHistogram("ESD_NoCutConversion_XY", fV0Reader->GetX(),fV0Reader->GetY());\r
 742         fHistograms->FillHistogram("ESD_NoCutConversion_R", fV0Reader->GetXYRadius());\r
 743         fHistograms->FillHistogram("ESD_NoCutConversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());\r
 744         fHistograms->FillHistogram("ESD_NoCutConversion_OpeningAngle", fV0Reader->GetOpeningAngle());\r
 745         \r
 746         /*\r
 747           ESD_NoCutConvGamma_Pt\r
 748           ESD_NoCutConvGamma_Energy\r
 749           ESD_NoCutConvGamma_Eta\r
 750           ESD_NoCutConvGamma_Phi\r
 751           ESD_NoCutConvGamma_Mass\r
 752           ESD_NoCutConvGamma_Width\r
 753           ESD_NoCutConvGamma_Chi2\r
 754           ESD_NoCutConvGamma_NDF\r
 755           ESD_NoCutConvGamma_PtvsEta\r
 756           ESD_NoCutConversion_XY\r
 757           ESD_NoCutConversion_R\r
 758           ESD_NoCutConversion_ZR\r
 759           ESD_NoCutConversion_OpeningAngle\r
 760         */\r
 761       }\r
 762     }\r
 763   }\r
 764   fV0Reader->ResetV0IndexNumber();\r
 765 }\r
 766 \r
 767 void AliAnalysisTaskGammaConversion::ProcessV0s(){\r
 768   // see header file for documentation\r
 769         \r
 770   if(fWriteNtuple == kTRUE){\r
 771     FillNtuple();\r
 772   }\r
 773         \r
 774   Int_t nSurvivingV0s=0;\r
 775   while(fV0Reader->NextV0()){\r
 776     nSurvivingV0s++;\r
 777                 \r
 778                 \r
 779     //-------------------------- filling v0 information -------------------------------------\r
 780     fHistograms->FillHistogram("ESD_Conversion_R", fV0Reader->GetXYRadius());\r
 781     fHistograms->FillHistogram("ESD_Conversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());\r
 782     fHistograms->FillHistogram("ESD_Conversion_XY", fV0Reader->GetX(),fV0Reader->GetY());\r
 783     fHistograms->FillHistogram("ESD_Conversion_OpeningAngle", fV0Reader->GetOpeningAngle());    \r
 784                 \r
 785     fHistograms->FillHistogram("ESD_E_Energy", fV0Reader->GetNegativeTrackEnergy());\r
 786     fHistograms->FillHistogram("ESD_E_Pt", fV0Reader->GetNegativeTrackPt());\r
 787     fHistograms->FillHistogram("ESD_E_Eta", fV0Reader->GetNegativeTrackEta());\r
 788     fHistograms->FillHistogram("ESD_E_Phi", fV0Reader->GetNegativeTrackPhi());\r
 789                 \r
 790     fHistograms->FillHistogram("ESD_P_Energy", fV0Reader->GetPositiveTrackEnergy());\r
 791     fHistograms->FillHistogram("ESD_P_Pt", fV0Reader->GetPositiveTrackPt());\r
 792     fHistograms->FillHistogram("ESD_P_Eta", fV0Reader->GetPositiveTrackEta());\r
 793     fHistograms->FillHistogram("ESD_P_Phi", fV0Reader->GetPositiveTrackPhi());\r
 794                 \r
 795     fHistograms->FillHistogram("ESD_ConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());\r
 796     fHistograms->FillHistogram("ESD_ConvGamma_Pt", fV0Reader->GetMotherCandidatePt());\r
 797     fHistograms->FillHistogram("ESD_ConvGamma_Eta", fV0Reader->GetMotherCandidateEta());\r
 798     fHistograms->FillHistogram("ESD_ConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());\r
 799     fHistograms->FillHistogram("ESD_ConvGamma_Mass", fV0Reader->GetMotherCandidateMass());\r
 800     fHistograms->FillHistogram("ESD_ConvGamma_Width", fV0Reader->GetMotherCandidateWidth());\r
 801     fHistograms->FillHistogram("ESD_ConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());\r
 802     fHistograms->FillHistogram("ESD_ConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());\r
 803     fHistograms->FillHistogram("ESD_ConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());\r
 804     fHistograms->FillHistogram("ESD_ConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());\r
 805                 \r
 806                 \r
 807     // begin mapping\r
 808     Int_t rBin    = fHistograms->GetRBin(fV0Reader->GetXYRadius());\r
 809     Int_t phiBin  = fHistograms->GetPhiBin(fV0Reader->GetNegativeTrackPhi());\r
 810     Double_t motherCandidateEta= fV0Reader->GetMotherCandidateEta();\r
 811                 \r
 812     TString nameESDMappingPhiR="";\r
 813     nameESDMappingPhiR.Form("ESD_Conversion_Mapping-Phi%02d-R%02d",phiBin,rBin);\r
 814     fHistograms->FillHistogram(nameESDMappingPhiR, fV0Reader->GetZ(), motherCandidateEta);\r
 815                 \r
 816     TString nameESDMappingPhi="";\r
 817     nameESDMappingPhi.Form("ESD_Conversion_Mapping-Phi%02d",phiBin);\r
 818     fHistograms->FillHistogram(nameESDMappingPhi, fV0Reader->GetZ(), motherCandidateEta);\r
 819                 \r
 820     TString nameESDMappingR="";\r
 821     nameESDMappingR.Form("ESD_Conversion_Mapping-R%02d",rBin);\r
 822     fHistograms->FillHistogram(nameESDMappingR, fV0Reader->GetZ(), motherCandidateEta);  \r
 823                 \r
 824     TString nameESDMappingPhiInR="";\r
 825     nameESDMappingPhiInR.Form("ESD_Conversion_Mapping_Phi_R-%02d",rBin);\r
 826     fHistograms->FillHistogram(nameESDMappingPhiInR, fV0Reader->GetMotherCandidatePhi());\r
 827     // end mapping\r
 828                 \r
 829     fKFReconstructedGammas.push_back(*fV0Reader->GetMotherCandidateKFCombination());\r
 830     fElectronv1.push_back(fV0Reader->GetCurrentV0()->GetPindex());\r
 831     fElectronv2.push_back(fV0Reader->GetCurrentV0()->GetNindex());\r
 832 \r
 833                 \r
 834     //----------------------------------- checking for "real" conversions (MC match) --------------------------------------\r
 835     if(fDoMCTruth){\r
 836                         \r
 837       if(fV0Reader->HasSameMCMother() == kFALSE){\r
 838         fIsTrueReconstructedGammas.push_back(kFALSE);\r
 839         continue;\r
 840       }\r
 841       \r
 842                 \r
 843       TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();\r
 844       TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();\r
 845 \r
 846       if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){\r
 847         fIsTrueReconstructedGammas.push_back(kFALSE);\r
 848         continue;\r
 849       }\r
 850       if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){\r
 851         fIsTrueReconstructedGammas.push_back(kFALSE);\r
 852         continue;\r
 853       }\r
 854         \r
 855       if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){\r
 856         fIsTrueReconstructedGammas.push_back(kTRUE);\r
 857                                 \r
 858         fHistograms->FillHistogram("ESD_TrueConvGamma_Pt", fV0Reader->GetMotherCandidatePt());\r
 859         fHistograms->FillHistogram("ESD_TrueConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());\r
 860         fHistograms->FillHistogram("ESD_TrueConvGamma_Eta", fV0Reader->GetMotherCandidateEta());                                \r
 861         fHistograms->FillHistogram("ESD_TrueConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());\r
 862         fHistograms->FillHistogram("ESD_TrueConvGamma_Mass", fV0Reader->GetMotherCandidateMass());\r
 863         fHistograms->FillHistogram("ESD_TrueConvGamma_Width", fV0Reader->GetMotherCandidateWidth());\r
 864         fHistograms->FillHistogram("ESD_TrueConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());\r
 865         fHistograms->FillHistogram("ESD_TrueConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());\r
 866         fHistograms->FillHistogram("ESD_TrueConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());\r
 867         fHistograms->FillHistogram("ESD_TrueConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());\r
 868         fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLength", /*fV0Reader->GetNegativeTrackLength()*/fV0Reader->GetNegativeNTPCClusters());\r
 869         fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLength", /*fV0Reader->GetPositiveTrackLength()*/fV0Reader->GetPositiveNTPCClusters());\r
 870         fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLengthVSInvMass",/*fV0Reader->GetNegativeTrackLength()*/fV0Reader->GetNegativeNTPCClusters(),fV0Reader->GetMotherCandidateMass());\r
 871         fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLengthVSInvMass",/*fV0Reader->GetPositiveTrackLength()*/fV0Reader->GetPositiveNTPCClusters(),fV0Reader->GetMotherCandidateMass());\r
 872         \r
 873         fHistograms->FillHistogram("ESD_TrueConversion_XY", fV0Reader->GetX(),fV0Reader->GetY());\r
 874         fHistograms->FillHistogram("ESD_TrueConversion_R", fV0Reader->GetXYRadius());\r
 875         fHistograms->FillHistogram("ESD_TrueConversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());\r
 876         fHistograms->FillHistogram("ESD_TrueConversion_OpeningAngle", fV0Reader->GetOpeningAngle());\r
 877 \r
 878         \r
 879         /*\r
 880           fHistograms->FillHistogram("ESD_TrueConversion_XY", fV0Reader->GetX(),fV0Reader->GetY());\r
 881           fHistograms->FillHistogram("ESD_TrueConversion_R", fV0Reader->GetXYRadius());\r
 882           fHistograms->FillHistogram("ESD_TrueConversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());\r
 883           fHistograms->FillHistogram("ESD_TrueConversion_OpeningAngle", fV0Reader->GetOpeningAngle());\r
 884         */\r
 885 \r
 886 \r
 887 \r
 888         //resolution\r
 889         Double_t mcpt   = fV0Reader->GetMotherMCParticle()->Pt();\r
 890         Double_t esdpt  = fV0Reader->GetMotherCandidatePt();\r
 891         Double_t resdPt = 0;\r
 892         if(mcpt != 0){\r
 893           resdPt = ((esdpt - mcpt)/mcpt)*100;\r
 894         }\r
 895                                 \r
 896         fHistograms->FillHistogram("Resolution_dPt", mcpt, resdPt);\r
 897         fHistograms->FillHistogram("Resolution_MC_Pt", mcpt);\r
 898         fHistograms->FillHistogram("Resolution_ESD_Pt", esdpt);\r
 899                                 \r
 900         Double_t resdZ = 0;\r
 901         if(fV0Reader->GetNegativeMCParticle()->Vz() != 0){\r
 902           resdZ = ((fV0Reader->GetZ() -fV0Reader->GetNegativeMCParticle()->Vz())/fV0Reader->GetNegativeMCParticle()->Vz())*100;\r
 903         }\r
 904                                 \r
 905         fHistograms->FillHistogram("Resolution_dZ", fV0Reader->GetNegativeMCParticle()->Vz(), resdZ);\r
 906         fHistograms->FillHistogram("Resolution_MC_Z", fV0Reader->GetNegativeMCParticle()->Vz());\r
 907         fHistograms->FillHistogram("Resolution_ESD_Z", fV0Reader->GetZ());\r
 908                                 \r
 909         Double_t resdR = 0;\r
 910         if(fV0Reader->GetNegativeMCParticle()->R() != 0){\r
 911           resdR = ((fV0Reader->GetXYRadius() - fV0Reader->GetNegativeMCParticle()->R())/fV0Reader->GetNegativeMCParticle()->R())*100;\r
 912         }\r
 913 \r
 914         fHistograms->FillHistogram("Resolution_dR", fV0Reader->GetNegativeMCParticle()->R(), resdR);\r
 915         fHistograms->FillHistogram("Resolution_MC_R", fV0Reader->GetNegativeMCParticle()->R());\r
 916         fHistograms->FillHistogram("Resolution_ESD_R", fV0Reader->GetXYRadius());\r
 917         fHistograms->FillHistogram("Resolution_dR_dPt", resdR, resdPt);\r
 918       }//if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22)\r
 919       else{\r
 920         fIsTrueReconstructedGammas.push_back(kFALSE);\r
 921       }\r
 922     }//if(fDoMCTruth)\r
 923   }//while(fV0Reader->NextV0)\r
 924   fHistograms->FillHistogram("ESD_NumberOfSurvivingV0s", nSurvivingV0s);\r
 925   fHistograms->FillHistogram("ESD_NumberOfV0s", fV0Reader->GetNumberOfV0s());\r
 926         \r
 927   //cout << "nearly at the end of doMCTruth" << endl;\r
 928         \r
 929 }\r
 930 \r
 931 void AliAnalysisTaskGammaConversion::ProcessGammasForNeutralMesonAnalysis(){\r
 932   // see header file for documentation\r
 933         \r
 934   for(UInt_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammas.size();firstGammaIndex++){\r
 935     for(UInt_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fKFReconstructedGammas.size();secondGammaIndex++){\r
 936                         \r
 937       AliKFParticle * twoGammaDecayCandidateDaughter0 = &fKFReconstructedGammas[firstGammaIndex];\r
 938       AliKFParticle * twoGammaDecayCandidateDaughter1 = &fKFReconstructedGammas[secondGammaIndex];\r
 939       \r
 940       if(fElectronv1[firstGammaIndex]==fElectronv1[secondGammaIndex] || fElectronv1[firstGammaIndex]==fElectronv2[secondGammaIndex]){\r
 941         continue;\r
 942       }\r
 943       if(fElectronv2[firstGammaIndex]==fElectronv1[secondGammaIndex] || fElectronv2[firstGammaIndex]==fElectronv2[secondGammaIndex]){\r
 944         continue;\r
 945       }\r
 946 \r
 947       /*\r
 948         if(fIsTrueReconstructedGammas[firstGammaIndex] == kFALSE || fIsTrueReconstructedGammas[secondGammaIndex] == kFALSE){\r
 949         continue;\r
 950         }\r
 951       */\r
 952 \r
 953       AliKFParticle *twoGammaCandidate = new AliKFParticle(*twoGammaDecayCandidateDaughter0,*twoGammaDecayCandidateDaughter1);\r
 954                         \r
 955       Double_t massTwoGammaCandidate = 0.;\r
 956       Double_t widthTwoGammaCandidate = 0.;\r
 957       Double_t chi2TwoGammaCandidate =10000.;   \r
 958       twoGammaCandidate->GetMass(massTwoGammaCandidate,widthTwoGammaCandidate);\r
 959       if(twoGammaCandidate->GetNDF()>0){\r
 960         chi2TwoGammaCandidate = twoGammaCandidate->GetChi2()/twoGammaCandidate->GetNDF();\r
 961                                 \r
 962         if(chi2TwoGammaCandidate>0 && chi2TwoGammaCandidate<fV0Reader->GetChi2CutMeson()){                                      \r
 963                                         \r
 964           TVector3 momentumVectorTwoGammaCandidate(twoGammaCandidate->GetPx(),twoGammaCandidate->GetPy(),twoGammaCandidate->GetPz());\r
 965           TVector3 spaceVectorTwoGammaCandidate(twoGammaCandidate->GetX(),twoGammaCandidate->GetY(),twoGammaCandidate->GetZ());\r
 966                                         \r
 967           Double_t openingAngleTwoGammaCandidate = twoGammaDecayCandidateDaughter0->GetAngle(*twoGammaDecayCandidateDaughter1);                                 \r
 968           Double_t rapidity;\r
 969           if(twoGammaCandidate->GetE() - twoGammaCandidate->GetPz() == 0 || twoGammaCandidate->GetE() + twoGammaCandidate->GetPz() == 0){\r
 970             rapidity=0;\r
 971           }\r
 972           else{\r
 973             rapidity = 0.5*(TMath::Log((twoGammaCandidate->GetE() +twoGammaCandidate->GetPz()) / (twoGammaCandidate->GetE()-twoGammaCandidate->GetPz())));\r
 974           }\r
 975                                         \r
 976           if(openingAngleTwoGammaCandidate < fMinOpeningAngleGhostCut) continue;   // minimum opening angle to avoid using ghosttracks\r
 977 \r
 978           fHistograms->FillHistogram("ESD_Mother_GammaDaughter_OpeningAngle", openingAngleTwoGammaCandidate);\r
 979           fHistograms->FillHistogram("ESD_Mother_Energy", twoGammaCandidate->GetE());\r
 980           fHistograms->FillHistogram("ESD_Mother_Pt", momentumVectorTwoGammaCandidate.Pt());\r
 981           fHistograms->FillHistogram("ESD_Mother_Eta", momentumVectorTwoGammaCandidate.Eta());\r
 982           fHistograms->FillHistogram("ESD_Mother_Rapid", rapidity);                                     \r
 983           fHistograms->FillHistogram("ESD_Mother_Phi", spaceVectorTwoGammaCandidate.Phi());\r
 984           fHistograms->FillHistogram("ESD_Mother_Mass", massTwoGammaCandidate);\r
 985           fHistograms->FillHistogram("ESD_Mother_R", spaceVectorTwoGammaCandidate.Pt());    // Pt in Space == R!!!\r
 986           fHistograms->FillHistogram("ESD_Mother_ZR", twoGammaCandidate->GetZ(), spaceVectorTwoGammaCandidate.Pt());\r
 987           fHistograms->FillHistogram("ESD_Mother_XY", twoGammaCandidate->GetX(), twoGammaCandidate->GetY());\r
 988           fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());\r
 989           fHistograms->FillHistogram("ESD_Mother_InvMass",massTwoGammaCandidate);\r
 990         }\r
 991       }\r
 992       delete twoGammaCandidate;\r
 993                         \r
 994       //cout << "nearly at the end of processgamma for neutral meson ..." << endl;\r
 995                         \r
 996                         \r
 997     }\r
 998   }\r
 999 }\r
1000 \r
1001 void AliAnalysisTaskGammaConversion::CalculateBackground(){\r
1002   // see header file for documentation\r
1003         \r
1004   vector<AliKFParticle> vectorCurrentEventGoodV0s = fV0Reader->GetCurrentEventGoodV0s();\r
1005   vector<AliKFParticle> vectorPreviousEventGoodV0s = fV0Reader->GetPreviousEventGoodV0s();\r
1006   for(UInt_t iCurrent=0;iCurrent<vectorCurrentEventGoodV0s.size();iCurrent++){\r
1007     AliKFParticle * currentEventGoodV0 = &vectorCurrentEventGoodV0s.at(iCurrent);\r
1008     for(UInt_t iPrevious=0;iPrevious<vectorPreviousEventGoodV0s.size();iPrevious++){\r
1009       AliKFParticle * previousGoodV0 = &vectorPreviousEventGoodV0s.at(iPrevious);\r
1010 \r
1011       AliKFParticle *backgroundCandidate = new AliKFParticle(*currentEventGoodV0,*previousGoodV0);\r
1012                         \r
1013       Double_t massBG =0.;\r
1014       Double_t widthBG = 0.;\r
1015       Double_t chi2BG =10000.;  \r
1016       backgroundCandidate->GetMass(massBG,widthBG);\r
1017       if(backgroundCandidate->GetNDF()>0){\r
1018         chi2BG = backgroundCandidate->GetChi2()/backgroundCandidate->GetNDF();\r
1019         if(chi2BG>0 && chi2BG<fV0Reader->GetChi2CutMeson()){\r
1020                                         \r
1021           TVector3 MomentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());\r
1022           TVector3 SpaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());\r
1023                                         \r
1024           Double_t openingAngleBG = currentEventGoodV0->GetAngle(*previousGoodV0);\r
1025 \r
1026           Double_t rapidity;\r
1027           if(backgroundCandidate->GetE() - backgroundCandidate->GetPz() == 0 || backgroundCandidate->GetE() + backgroundCandidate->GetPz() == 0) rapidity=0;\r
1028           else rapidity = 0.5*(TMath::Log((backgroundCandidate->GetE() +backgroundCandidate->GetPz()) / (backgroundCandidate->GetE()-backgroundCandidate->GetPz())));\r
1029 \r
1030                                         \r
1031                                         \r
1032                                         \r
1033           if(openingAngleBG < fMinOpeningAngleGhostCut ) continue;   // minimum opening angle to avoid using ghosttracks\r
1034                         \r
1035                                         \r
1036           fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle", openingAngleBG);\r
1037           fHistograms->FillHistogram("ESD_Background_Energy", backgroundCandidate->GetE());\r
1038           fHistograms->FillHistogram("ESD_Background_Pt",  MomentumVectorbackgroundCandidate.Pt());\r
1039           fHistograms->FillHistogram("ESD_Background_Eta", MomentumVectorbackgroundCandidate.Eta());\r
1040           fHistograms->FillHistogram("ESD_Background_Rapidity", rapidity);\r
1041           fHistograms->FillHistogram("ESD_Background_Phi", SpaceVectorbackgroundCandidate.Phi());\r
1042           fHistograms->FillHistogram("ESD_Background_Mass", massBG);\r
1043           fHistograms->FillHistogram("ESD_Background_R", SpaceVectorbackgroundCandidate.Pt());  // Pt in Space == R!!!!\r
1044           fHistograms->FillHistogram("ESD_Background_ZR", backgroundCandidate->GetZ(), SpaceVectorbackgroundCandidate.Pt());\r
1045           fHistograms->FillHistogram("ESD_Background_XY", backgroundCandidate->GetX(), backgroundCandidate->GetY());\r
1046           fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt",massBG,MomentumVectorbackgroundCandidate.Pt());\r
1047           fHistograms->FillHistogram("ESD_Background_InvMass",massBG);\r
1048         }\r
1049       }\r
1050       delete backgroundCandidate;   \r
1051       //cout << "nearly at the end of background" << endl;\r
1052                         \r
1053     }\r
1054   }\r
1055 }\r
1056 \r
1057 void AliAnalysisTaskGammaConversion::Terminate(Option_t */*option*/)\r
1058 {\r
1059   // Terminate analysis\r
1060   //\r
1061   AliDebug(1,"Do nothing in Terminate");\r
1062 }\r
1063 \r
1064 void AliAnalysisTaskGammaConversion::UserCreateOutputObjects()\r
1065 {\r
1066   // Create the output container\r
1067   if(fOutputContainer != NULL){\r
1068     delete fOutputContainer;\r
1069     fOutputContainer = NULL;\r
1070   }\r
1071   if(fOutputContainer == NULL){\r
1072     fOutputContainer = new TList();\r
1073   }\r
1074         \r
1075   //Adding the histograms to the output container\r
1076   fHistograms->GetOutputContainer(fOutputContainer);\r
1077         \r
1078         \r
1079   if(fWriteNtuple){\r
1080     if(fGammaNtuple == NULL){\r
1081       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
1082     }\r
1083     if(fNeutralMesonNtuple == NULL){\r
1084       fNeutralMesonNtuple = new TNtuple("NeutralMesonNtuple","NeutralMesonNtuple","test");\r
1085     }\r
1086     TList * ntupleTList = new TList();\r
1087     ntupleTList->SetName("Ntuple");\r
1088     ntupleTList->Add((TNtuple*)fGammaNtuple);\r
1089     fOutputContainer->Add(ntupleTList);\r
1090   }\r
1091         \r
1092   fOutputContainer->SetName(GetName());\r
1093 }\r
1094 \r
1095 Double_t AliAnalysisTaskGammaConversion::GetMCOpeningAngle(TParticle* const daughter0, TParticle* const daughter1) const{\r
1096   //helper function\r
1097   TVector3 v3D0(daughter0->Px(),daughter0->Py(),daughter0->Pz());\r
1098   TVector3 v3D1(daughter1->Px(),daughter1->Py(),daughter1->Pz());\r
1099   return v3D0.Angle(v3D1);\r
1100 }\r
1101 \r
1102 void AliAnalysisTaskGammaConversion::CheckV0Efficiency(){\r
1103 \r
1104   vector<Int_t> indexOfGammaParticle;\r
1105 \r
1106   fStack = fV0Reader->GetMCStack();\r
1107 \r
1108   if(fV0Reader->CheckForPrimaryVertex() == kFALSE){\r
1109     return; // aborts if the primary vertex does not have contributors.\r
1110   }\r
1111 \r
1112   for (Int_t iTracks = 0; iTracks < fStack->GetNprimary(); iTracks++) {\r
1113     TParticle* particle = (TParticle *)fStack->Particle(iTracks);\r
1114     if(particle->GetPdgCode()==22){     //Gamma\r
1115       if(particle->GetNDaughters() >= 2){\r
1116         TParticle* electron=NULL;\r
1117         TParticle* positron=NULL; \r
1118         for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){\r
1119           TParticle *tmpDaughter = fStack->Particle(daughterIndex);\r
1120           if(tmpDaughter->GetUniqueID() == 5){\r
1121             if(tmpDaughter->GetPdgCode() == 11){\r
1122               electron = tmpDaughter;\r
1123             }\r
1124             else if(tmpDaughter->GetPdgCode() == -11){\r
1125               positron = tmpDaughter;\r
1126             }\r
1127           }\r
1128         }\r
1129         if(electron!=NULL && positron!=0){\r
1130           if(electron->R()<160){\r
1131             indexOfGammaParticle.push_back(iTracks);\r
1132           }\r
1133         }\r
1134       }\r
1135     }\r
1136   }\r
1137 \r
1138   Int_t nFoundGammas=0;\r
1139   Int_t nNotFoundGammas=0;\r
1140 \r
1141   Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();\r
1142   for(Int_t i=0;i<numberOfV0s;i++){\r
1143     fV0Reader->GetV0(i);\r
1144     \r
1145     if(fV0Reader->HasSameMCMother() == kFALSE){\r
1146       continue;\r
1147     }\r
1148     \r
1149     TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();\r
1150     TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();\r
1151     \r
1152     if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){\r
1153       continue;\r
1154     }\r
1155     if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){\r
1156       continue;\r
1157     }\r
1158     \r
1159     if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){\r
1160       //TParticle * v0Gamma = fV0Reader->GetMotherMCParticle();\r
1161       for(UInt_t mcIndex=0;mcIndex<indexOfGammaParticle.size();mcIndex++){\r
1162         if(negativeMC->GetFirstMother()==indexOfGammaParticle[mcIndex]){\r
1163           nFoundGammas++;\r
1164         }\r
1165         else{\r
1166           nNotFoundGammas++;\r
1167         }\r
1168       }\r
1169     }\r
1170   }\r
1171   //  cout<<"Found: "<<nFoundGammas<<"  of: "<<indexOfGammaParticle.size()<<endl;\r
1172 }\r
1173 \r
1174 \r
1175 void AliAnalysisTaskGammaConversion::ProcessGammaElectronsForChicAnalysis(){\r
1176 \r
1177 \r
1178   fESDEvent = fV0Reader->GetESDEvent();\r
1179 \r
1180 \r
1181   vector <AliESDtrack*> ESDeNegTemp(0);\r
1182   vector <AliESDtrack*> ESDePosTemp(0);\r
1183   vector <AliESDtrack*> ESDxNegTemp(0);\r
1184   vector <AliESDtrack*> ESDxPosTemp(0);\r
1185   vector <AliESDtrack*> ESDeNegNoJPsi(0);\r
1186   vector <AliESDtrack*> ESDePosNoJPsi(0); \r
1187 \r
1188 \r
1189 \r
1190   fHistograms->FillTable("Table_Electrons",0);//Count number of Events\r
1191 \r
1192   for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){\r
1193     AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);\r
1194 \r
1195     if(!curTrack){\r
1196       //print warning here\r
1197       continue;\r
1198     }\r
1199 \r
1200     double p[3];if(!curTrack->GetConstrainedPxPyPz(p))continue;\r
1201     double r[3];curTrack->GetConstrainedXYZ(r);\r
1202 \r
1203     TVector3 rXYZ(r);\r
1204 \r
1205     fHistograms->FillTable("Table_Electrons",4);//Count number of ESD tracks\r
1206 \r
1207     Bool_t flagKink       =  kTRUE;\r
1208     Bool_t flagTPCrefit   =  kTRUE;\r
1209     Bool_t flagTRDrefit   =  kTRUE;\r
1210     Bool_t flagITSrefit   =  kTRUE;\r
1211     Bool_t flagTRDout     =  kTRUE;\r
1212     Bool_t flagVertex     =  kTRUE;\r
1213 \r
1214 \r
1215     //Cuts ---------------------------------------------------------------\r
1216 \r
1217     if(curTrack->GetKinkIndex(0) > 0){\r
1218       fHistograms->FillHistogram("Table_Electrons",5);//Count kink\r
1219       flagKink = kFALSE;\r
1220     }\r
1221 \r
1222     ULong_t trkStatus = curTrack->GetStatus();\r
1223 \r
1224     ULong_t tpcRefit = (trkStatus & AliESDtrack::kTPCrefit);\r
1225 \r
1226     if(!tpcRefit){\r
1227       fHistograms->FillHistogram("Table_Electrons",9);//Count not TPCrefit\r
1228       flagTPCrefit = kFALSE;\r
1229     }\r
1230 \r
1231     ULong_t itsRefit = (trkStatus & AliESDtrack::kITSrefit);\r
1232     if(!itsRefit){\r
1233       fHistograms->FillHistogram("Table_Electrons",10);//Count not ITSrefit\r
1234       flagITSrefit = kFALSE;\r
1235     }\r
1236 \r
1237     ULong_t trdRefit = (trkStatus & AliESDtrack::kTRDrefit);\r
1238 \r
1239     if(!trdRefit){\r
1240       fHistograms->FillHistogram("Table_Electrons",8); //Count not TRDrefit\r
1241       flagTRDrefit = kFALSE;\r
1242     }\r
1243 \r
1244     ULong_t trdOut = (trkStatus & AliESDtrack::kTRDout);\r
1245 \r
1246     if(!trdOut) {\r
1247       fHistograms->FillHistogram("Table_Electrons",7); //Count not TRDout\r
1248       flagTRDout = kFALSE;\r
1249     }\r
1250 \r
1251     double nsigmaToVxt = GetSigmaToVertex(curTrack);\r
1252 \r
1253     if(nsigmaToVxt > 3){\r
1254       fHistograms->FillHistogram("Table_Electrons",6); //Count Tracks with number of sigmas > 3\r
1255       flagVertex = kFALSE;\r
1256     }\r
1257 \r
1258     if(! (flagKink && flagTPCrefit && flagITSrefit && flagTRDrefit && flagTRDout && flagVertex ) ) continue;\r
1259     fHistograms->FillHistogram("Table_Electrons",11);//Count Tracks passed Cuts\r
1260 \r
1261 \r
1262     Stat_t pid, weight;\r
1263     GetPID(curTrack, pid, weight);\r
1264 \r
1265     if(pid!=0){\r
1266       fHistograms->FillHistogram("Table_Electrons",12); //Count Tracks with pid != 0\r
1267     }\r
1268 \r
1269     if(pid == 0){\r
1270       fHistograms->FillHistogram("Table_Electrons",13); //Count Tracks with pid != 0\r
1271     }\r
1272 \r
1273 \r
1274 \r
1275 \r
1276     Int_t labelMC = TMath::Abs(curTrack->GetLabel());\r
1277     TParticle* curParticle = fStack->Particle(labelMC);\r
1278 \r
1279 \r
1280 \r
1281 \r
1282     TLorentzVector curElec;\r
1283     curElec.SetXYZM(p[0],p[1],p[2],fElectronMass);\r
1284 \r
1285 \r
1286 \r
1287 \r
1288     if(curTrack->GetSign() > 0){\r
1289 \r
1290       ESDxPosTemp.push_back(curTrack);\r
1291 \r
1292       if( pid == 0){\r
1293 \r
1294         fHistograms->FillHistogram("ESD_ElectronPosNegPt",curElec.Pt());\r
1295         fHistograms->FillHistogram("ESD_ElectronPosPt",curElec.Pt());\r
1296         fHistograms->FillHistogram("MC_ElectronPosNegPt",curParticle->Pt());\r
1297         fHistograms->FillHistogram("ESD_ElectronPosNegEta",curElec.Eta());\r
1298         fHistograms->FillHistogram("MC_ElectronPosNegEta",curParticle->Eta());\r
1299         ESDePosTemp.push_back(curTrack);\r
1300 \r
1301 \r
1302 \r
1303       }\r
1304 \r
1305     }\r
1306     else {\r
1307       ESDxNegTemp.push_back(curTrack);\r
1308 \r
1309       if( pid == 0){\r
1310 \r
1311         fHistograms->FillHistogram("ESD_ElectronPosNegPt",curElec.Pt());\r
1312         fHistograms->FillHistogram("ESD_ElectronNegPt",curElec.Pt());\r
1313         fHistograms->FillHistogram("MC_ElectronPosNegPt",curParticle->Pt());\r
1314         fHistograms->FillHistogram("ESD_ElectronPosNegEta",curElec.Eta());\r
1315         fHistograms->FillHistogram("MC_ElectronPosNegEta",curParticle->Eta());\r
1316         ESDeNegTemp.push_back(curTrack);\r
1317 \r
1318 \r
1319 \r
1320 \r
1321       }\r
1322     }\r
1323 \r
1324   }\r
1325 \r
1326 \r
1327   Bool_t ePosJPsi = kFALSE;\r
1328   Bool_t eNegJPsi = kFALSE;             \r
1329   Bool_t ePosPi0  = kFALSE;\r
1330   Bool_t eNegPi0  = kFALSE;\r
1331         \r
1332   UInt_t iePosJPsi=0,ieNegJPsi=0,iePosPi0=0,ieNegPi0=0;\r
1333  \r
1334   for(UInt_t iNeg=0; iNeg < ESDeNegTemp.size(); iNeg++){\r
1335     if(fStack->Particle(TMath::Abs(ESDeNegTemp[iNeg]->GetLabel()))->GetPdgCode() == 11)\r
1336       if(fStack->Particle(TMath::Abs(ESDeNegTemp[iNeg]->GetLabel()))->GetMother(0) > -1){\r
1337         Int_t labelMother = fStack->Particle(TMath::Abs(ESDeNegTemp[iNeg]->GetLabel()))->GetMother(0);\r
1338         TParticle* partMother = fStack ->Particle(labelMother);\r
1339         if (partMother->GetPdgCode() == 111){\r
1340           ieNegPi0 = iNeg;\r
1341           eNegPi0 = kTRUE;\r
1342         }\r
1343         if(partMother->GetPdgCode() == 443){ //Mother JPsi\r
1344           fHistograms->FillTable("Table_Electrons",14);\r
1345           ieNegJPsi = iNeg;\r
1346           eNegJPsi = kTRUE;\r
1347         }\r
1348         else{   \r
1349           ESDeNegNoJPsi.push_back(ESDeNegTemp[iNeg]);\r
1350           //            cout<<"ESD No Positivo JPsi "<<endl;\r
1351         }\r
1352 \r
1353       }\r
1354   }     \r
1355 \r
1356   for(UInt_t iPos=0; iPos < ESDePosTemp.size(); iPos++){\r
1357     if(fStack->Particle(TMath::Abs(ESDePosTemp[iPos]->GetLabel()))->GetPdgCode() == -11)\r
1358       if(fStack->Particle(TMath::Abs(ESDePosTemp[iPos]->GetLabel()))->GetMother(0) > -1){\r
1359         Int_t labelMother = fStack->Particle(TMath::Abs(ESDePosTemp[iPos]->GetLabel()))->GetMother(0);\r
1360         TParticle* partMother = fStack ->Particle(labelMother);\r
1361         if (partMother->GetPdgCode() == 111){\r
1362           iePosPi0 = iPos;\r
1363           ePosPi0 = kTRUE;\r
1364         }\r
1365         if(partMother->GetPdgCode() == 443){ //Mother JPsi\r
1366           fHistograms->FillTable("Table_Electrons",15);\r
1367           iePosJPsi = iPos;\r
1368           ePosJPsi = kTRUE;\r
1369         }\r
1370         else{\r
1371           ESDePosNoJPsi.push_back(ESDePosTemp[iPos]);\r
1372           //            cout<<"ESD No Negativo JPsi "<<endl;\r
1373         }\r
1374 \r
1375       }\r
1376   }\r
1377         \r
1378   if( eNegJPsi && ePosJPsi ){\r
1379     TVector3 tempeNegV,tempePosV;\r
1380     tempeNegV.SetXYZ(ESDeNegTemp[ieNegJPsi]->Px(),ESDeNegTemp[ieNegJPsi]->Py(),ESDeNegTemp[ieNegJPsi]->Pz());                   \r
1381     tempePosV.SetXYZ(ESDePosTemp[iePosJPsi]->Px(),ESDePosTemp[iePosJPsi]->Py(),ESDePosTemp[iePosJPsi]->Pz());\r
1382     fHistograms->FillTable("Table_Electrons",16);\r
1383     fHistograms->FillHistogram("ESD_ElectronPosNegJPsiAngle",tempeNegV.Angle(tempePosV));       \r
1384     fHistograms->FillHistogram("MC_ElectronPosNegJPsiAngle",GetMCOpeningAngle(fStack->Particle(TMath::Abs(ESDeNegTemp[ieNegJPsi]->GetLabel())),\r
1385                                                                               fStack->Particle(TMath::Abs(ESDePosTemp[iePosJPsi]->GetLabel()))));       \r
1386   }\r
1387         \r
1388   if( eNegPi0 && ePosPi0 ){\r
1389     TVector3 tempeNegV,tempePosV;\r
1390     tempeNegV.SetXYZ(ESDeNegTemp[ieNegPi0]->Px(),ESDeNegTemp[ieNegPi0]->Py(),ESDeNegTemp[ieNegPi0]->Pz());\r
1391     tempePosV.SetXYZ(ESDePosTemp[iePosPi0]->Px(),ESDePosTemp[iePosPi0]->Py(),ESDePosTemp[iePosPi0]->Pz());\r
1392     fHistograms->FillHistogram("ESD_ElectronPosNegPi0Angle",tempeNegV.Angle(tempePosV));\r
1393     fHistograms->FillHistogram("MC_ElectronPosNegPi0Angle",GetMCOpeningAngle(fStack->Particle(TMath::Abs(ESDeNegTemp[ieNegPi0]->GetLabel())),\r
1394                                                                              fStack->Particle(TMath::Abs(ESDePosTemp[iePosPi0]->GetLabel()))));   \r
1395   }\r
1396          \r
1397 \r
1398   FillAngle("ESD_eNegePosAngleBeforeCut",GetTLorentzVector(ESDeNegTemp),GetTLorentzVector(ESDePosTemp));\r
1399 \r
1400   CleanWithAngleCuts(ESDeNegTemp,ESDePosTemp,fKFReconstructedGammas);\r
1401         \r
1402   vector <TLorentzVector> CurrentTLVeNeg = GetTLorentzVector(fCurrentEventNegElectron);\r
1403   vector <TLorentzVector> CurrentTLVePos = GetTLorentzVector(fCurrentEventPosElectron);\r
1404 \r
1405 \r
1406   FillAngle("ESD_eNegePosAngleAfterCut",CurrentTLVeNeg,CurrentTLVePos);\r
1407 \r
1408  \r
1409 \r
1410 \r
1411   //FillAngle("ESD_eNegePosAngleAfterCut",CurrentTLVeNeg,CurrentTLVePos);\r
1412 \r
1413 \r
1414   FillElectronInvMass("ESD_InvMass_ePluseMinus",CurrentTLVeNeg,CurrentTLVePos);\r
1415   FillElectronInvMass("ESD_InvMass_xPlusxMinus",GetTLorentzVector(ESDxNegTemp),GetTLorentzVector(ESDxPosTemp));\r
1416 \r
1417        \r
1418 \r
1419   FillGammaElectronInvMass("ESD_InvMass_GammaePluseMinusChiC","ESD_InvMass_GammaePluseMinusChiCDiff",\r
1420                            fKFReconstructedGammasCut,CurrentTLVeNeg,CurrentTLVePos);\r
1421 \r
1422   FillGammaElectronInvMass("ESD_InvMass_GammaePluseMinusPi0","ESD_InvMass_GammaePluseMinusPi0Diff",\r
1423                            fKFReconstructedGammasCut,CurrentTLVeNeg,CurrentTLVePos);\r
1424 \r
1425   //BackGround\r
1426 \r
1427   //Like Sign e+e-\r
1428   ElectronBackground("ESD_ENegBackground",CurrentTLVeNeg);\r
1429   ElectronBackground("ESD_EPosBackground",CurrentTLVePos);\r
1430   ElectronBackground("ESD_EPosENegBackground",CurrentTLVeNeg);\r
1431   ElectronBackground("ESD_EPosENegBackground",CurrentTLVePos);\r
1432 \r
1433   //        Like Sign e+e- no JPsi\r
1434   ElectronBackground("ESD_EPosENegNoJPsiBG",GetTLorentzVector(ESDeNegNoJPsi));\r
1435   ElectronBackground("ESD_EPosENegNoJPsiBG",GetTLorentzVector(ESDePosNoJPsi));\r
1436 \r
1437   //Mixed Event\r
1438 \r
1439   if( fCurrentEventPosElectron.size() > 0 && fCurrentEventNegElectron.size() > 0 && fKFReconstructedGammasCut.size() > 0 ){\r
1440     FillGammaElectronInvMass("ESD_EPosENegGammaBackgroundMX","ESD_EPosENegGammaBackgroundMXDiff",\r
1441                              fKFReconstructedGammasCut,fPreviousEventTLVNegElectron,fPreviousEventTLVPosElectron);\r
1442     fPreviousEventTLVNegElectron = CurrentTLVeNeg;\r
1443     fPreviousEventTLVPosElectron = CurrentTLVePos;\r
1444 \r
1445   }\r
1446 \r
1447   /*\r
1448   //Photons P\r
1449   Double_t vtx[3];\r
1450   vtx[0]=0;vtx[1]=0;vtx[2]=0;\r
1451   for(UInt_t i=0;i<fKFReconstructedGammasChic.size();i++){\r
1452 \r
1453   //      if(fMCGammaChicTempCut[i]->GetMother(0) < 0) continue;\r
1454 \r
1455 \r
1456 \r
1457   Int_t tempLabel = fStack->Particle(fMCGammaChicTempCut[i]->GetMother(0))->GetPdgCode();\r
1458   //      cout<<" Label Pedro Gonzalez " <<tempLabel <<endl;\r
1459 \r
1460   //      cout<<" Label Distance"<<fKFReconstructedGammasChic[i].GetDistanceFromVertex(vtx)<<endl;\r
1461 \r
1462   if( tempLabel == 10441 || tempLabel == 20443 || tempLabel == 445 )\r
1463 \r
1464   fHistograms->FillHistogram("ESD_PhotonsMomentum",fKFReconstructedGammasChic[i].GetMomentum());\r
1465 \r
1466 \r
1467   }\r
1468 \r
1469 \r
1470   */\r
1471 \r
1472 \r
1473 }\r
1474 \r
1475 void AliAnalysisTaskGammaConversion::FillAngle(TString histoName,vector <TLorentzVector> TLVeNeg, vector <TLorentzVector> TLVePos){\r
1476   for( UInt_t iNeg=0; iNeg < TLVeNeg.size(); iNeg++){\r
1477     for (UInt_t iPos=0; iPos < TLVePos.size(); iPos++){\r
1478       fHistograms->FillHistogram(histoName.Data(),TLVeNeg[iNeg].Vect().Angle(TLVePos[iPos].Vect()));\r
1479     }\r
1480   }\r
1481 }\r
1482 void AliAnalysisTaskGammaConversion::FillElectronInvMass(TString histoName, vector <TLorentzVector> eNeg, vector <TLorentzVector> ePos){\r
1483 \r
1484   for( UInt_t n=0; n < eNeg.size(); n++){\r
1485 \r
1486     TLorentzVector en = eNeg.at(n);\r
1487     for (UInt_t p=0; p < ePos.size(); p++){\r
1488       TLorentzVector ep = ePos.at(p);\r
1489       TLorentzVector np = ep + en;\r
1490       fHistograms->FillHistogram(histoName.Data(),np.M());\r
1491     }\r
1492   }\r
1493 \r
1494 }\r
1495 \r
1496 void AliAnalysisTaskGammaConversion::FillGammaElectronInvMass(TString histoMass,TString histoDiff,vector <AliKFParticle> fKFGammas,\r
1497                                                               vector <TLorentzVector> TLVeNeg,vector<TLorentzVector> TLVePos)\r
1498 {\r
1499 \r
1500 \r
1501   for( UInt_t iNeg=0; iNeg < TLVeNeg.size(); iNeg++ ){\r
1502 \r
1503     for (UInt_t iPos=0; iPos < TLVePos.size(); iPos++){\r
1504 \r
1505       TLorentzVector xy = TLVePos[iPos] + TLVeNeg[iNeg];\r
1506 \r
1507       for (UInt_t iGam=0; iGam < fKFGammas.size(); iGam++){\r
1508 \r
1509         AliKFParticle * GammaCandidate = &fKFGammas[iGam];\r
1510         TLorentzVector g;\r
1511 \r
1512         g.SetXYZM(GammaCandidate->GetPx(),GammaCandidate->GetPy(),GammaCandidate->GetPz(),fGammaMass);\r
1513         TLorentzVector xyg = xy + g;\r
1514         fHistograms->FillHistogram(histoMass.Data(),xyg.M());\r
1515         fHistograms->FillHistogram(histoDiff.Data(),(xyg.M()-xy.M()));\r
1516       }\r
1517     }\r
1518   }\r
1519 \r
1520 }\r
1521 void AliAnalysisTaskGammaConversion::ElectronBackground(TString hBg, vector <TLorentzVector> e)\r
1522 {\r
1523   for(UInt_t i=0; i < e.size(); i++)\r
1524     {\r
1525       for (UInt_t j=i+1; j < e.size(); j++)\r
1526         {\r
1527           TLorentzVector ee = e[i] + e[j];\r
1528 \r
1529           fHistograms->FillHistogram(hBg.Data(),ee.M());\r
1530         }\r
1531     }\r
1532 }\r
1533 \r
1534 \r
1535 void AliAnalysisTaskGammaConversion::CleanWithAngleCuts(vector <AliESDtrack*> NegativeElectrons,\r
1536                                                         vector <AliESDtrack*> PositiveElectrons, vector <AliKFParticle> Gammas){\r
1537 \r
1538   UInt_t  Nsize = NegativeElectrons.size();\r
1539   UInt_t  Psize = PositiveElectrons.size();\r
1540   UInt_t  Gsize = Gammas.size();\r
1541 \r
1542 \r
1543 \r
1544   vector <Bool_t> xNegBand(Nsize);\r
1545   vector <Bool_t> xPosBand(Psize);\r
1546   vector <Bool_t> gammaBand(Gsize);\r
1547 \r
1548 \r
1549   for(UInt_t iNeg=0; iNeg < Nsize; iNeg++) xNegBand[iNeg]=kTRUE;\r
1550   for(UInt_t iPos=0; iPos < Psize; iPos++) xPosBand[iPos]=kTRUE;\r
1551   for(UInt_t iGam=0; iGam < Gsize; iGam++) gammaBand[iGam]=kTRUE;\r
1552 \r
1553 \r
1554   for(UInt_t iPos=0; iPos < Psize; iPos++){\r
1555         \r
1556     Double_t P[3]; PositiveElectrons[iPos]->GetConstrainedPxPyPz(P); \r
1557 \r
1558     TVector3 ePosV(P[0],P[1],P[2]);\r
1559 \r
1560     for(UInt_t iNeg=0; iNeg < Nsize; iNeg++){\r
1561         \r
1562       Double_t N[3]; NegativeElectrons[iNeg]->GetConstrainedPxPyPz(N); \r
1563       TVector3 eNegV(N[0],N[1],N[2]);\r
1564 \r
1565       if(ePosV.Angle(eNegV) < 0.05){ //e+e- from gamma\r
1566         xPosBand[iPos]=kFALSE;\r
1567         xNegBand[iNeg]=kFALSE;\r
1568       }\r
1569 \r
1570       for(UInt_t iGam=0; iGam < Gsize; iGam++){\r
1571         AliKFParticle* GammaCandidate = &Gammas[iGam];\r
1572         TVector3 GammaCandidateVector(GammaCandidate->Px(),GammaCandidate->Py(),GammaCandidate->Pz());\r
1573         if(ePosV.Angle(GammaCandidateVector) < 0.05 || eNegV.Angle(GammaCandidateVector) < 0.05)\r
1574           gammaBand[iGam]=kFALSE;\r
1575       }\r
1576     }\r
1577   }\r
1578 \r
1579 \r
1580 \r
1581 \r
1582   for(UInt_t iPos=0; iPos < Psize; iPos++){\r
1583     if(xPosBand[iPos]){\r
1584       fCurrentEventPosElectron.push_back(PositiveElectrons[iPos]);\r
1585     }\r
1586   }\r
1587   for(UInt_t iNeg=0;iNeg < Nsize; iNeg++){\r
1588     if(xNegBand[iNeg]){\r
1589       fCurrentEventNegElectron.push_back(NegativeElectrons[iNeg]);\r
1590     }\r
1591   }\r
1592   for(UInt_t iGam=0; iGam < Gsize; iGam++){\r
1593     if(gammaBand[iGam]){\r
1594       fKFReconstructedGammasCut.push_back(Gammas[iGam]);\r
1595     }\r
1596   }\r
1597 }\r
1598 \r
1599 \r
1600 void  AliAnalysisTaskGammaConversion::GetPID(AliESDtrack *track, Stat_t &pid, Stat_t &weight)\r
1601 {\r
1602   pid = -1;\r
1603   weight = -1;\r
1604 \r
1605   double wpart[5];\r
1606   double wpartbayes[5];\r
1607 \r
1608   //get probability of the diffenrent particle types\r
1609   track->GetESDpid(wpart);\r
1610 \r
1611   // Tentative particle type "concentrations"\r
1612   double c[5]={0.01, 0.01, 0.85, 0.10, 0.05};\r
1613 \r
1614   //Bayes' formula\r
1615   double rcc = 0.;\r
1616   for (int i = 0; i < 5; i++)\r
1617     {\r
1618       rcc+=(c[i] * wpart[i]);\r
1619     }\r
1620 \r
1621 \r
1622 \r
1623   for (int i=0; i<5; i++) {\r
1624     if( rcc!=0){\r
1625       wpartbayes[i] = c[i] * wpart[i] / rcc;\r
1626     }\r
1627   }\r
1628 \r
1629 \r
1630 \r
1631   Float_t max=0.;\r
1632   int ipid=-1;\r
1633   //find most probable particle in ESD pid\r
1634   //0:Electron - 1:Muon - 2:Pion - 3:Kaon - 4:Proton\r
1635   for (int i = 0; i < 5; i++)\r
1636     {\r
1637       if (wpartbayes[i] > max)\r
1638         {\r
1639           ipid = i;\r
1640           max = wpartbayes[i];\r
1641         }\r
1642     }\r
1643 \r
1644   pid = ipid;\r
1645   weight = max;\r
1646 }\r
1647 double AliAnalysisTaskGammaConversion::GetSigmaToVertex(AliESDtrack* t)\r
1648 {\r
1649   // Calculates the number of sigma to the vertex.\r
1650 \r
1651   Float_t b[2];\r
1652   Float_t bRes[2];\r
1653   Float_t bCov[3];\r
1654   t->GetImpactParameters(b,bCov);\r
1655   if (bCov[0]<=0 || bCov[2]<=0) {\r
1656     AliDebug(1, "Estimated b resolution lower or equal zero!");\r
1657     bCov[0]=0; bCov[2]=0;\r
1658   }\r
1659   bRes[0] = TMath::Sqrt(bCov[0]);\r
1660   bRes[1] = TMath::Sqrt(bCov[2]);\r
1661 \r
1662   // -----------------------------------\r
1663   // How to get to a n-sigma cut?\r
1664   //\r
1665   // The accumulated statistics from 0 to d is\r
1666   //\r
1667   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)\r
1668   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)\r
1669   //\r
1670   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)\r
1671   // Can this be expressed in a different way?\r
1672 \r
1673   if (bRes[0] == 0 || bRes[1] ==0)\r
1674     return -1;\r
1675 \r
1676   double d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));\r
1677 \r
1678   // stupid rounding problem screws up everything:\r
1679   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(\r
1680   if (TMath::Exp(-d * d / 2) < 1e-10)\r
1681     return 1000;\r
1682 \r
1683 \r
1684   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);\r
1685   return d;\r
1686 }\r
1687 vector <TLorentzVector> AliAnalysisTaskGammaConversion::GetTLorentzVector(vector <AliESDtrack*> ESDtrack){\r
1688 \r
1689   vector <TLorentzVector> TLVtrack(0);\r
1690 \r
1691   for(UInt_t itrack=0; itrack < ESDtrack.size(); itrack++){\r
1692     double P[3]; ESDtrack[itrack]->GetConstrainedPxPyPz(P);\r
1693     TLorentzVector CurrentTrack;\r
1694     CurrentTrack.SetXYZM(P[0],P[1],P[2],fElectronMass);\r
1695     TLVtrack.push_back(CurrentTrack);\r
1696   }\r
1697 \r
1698   return TLVtrack;\r
1699 }\r
1700 \r