PWGGA/GammaConv/AliV0ReaderV1.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Authors: Svein Lindal, Daniel Lohner                                   *
   5  * Version 1.0                                                            *
   6  *                                                                        *
   7  *                                                                        *
   8  * based on: on older version (see aliroot up to v5-04-42-AN)             *
   9  *           AliV0Reader.cxx                                              *
  10  *           Authors: Kathrin Koch, Kenneth Aamodt, Ana Marin             *
  11  *                                                                        *
  12  * Permission to use, copy, modify and distribute this software and its   *
  13  * documentation strictly for non-commercial purposes is hereby granted   *
  14  * without fee, provided that the above copyright notice appears in all   *
  15  * copies and that both the copyright notice and this permission notice   *
  16  * appear in the supporting documentation. The authors make no claims     *
  17  * about the suitability of this software for any purpose. It is          *
  18  * provided "as is" without express or implied warranty.                  *
  19  **************************************************************************/
  20
  21 ////////////////////////////////////////////////
  22 //---------------------------------------------
  23 // Class reconstructing conversion photons from V0s
  24 //---------------------------------------------
  25 ////////////////////////////////////////////////
  26
  27 // The AliAODConversionPhotons will return the Position (== ESDTrack->GetID())of the positive (negative) ESDTrack
  28 // in the Array of ESDTracks stored in the ESDEvent via GetTrackLabelPositive() (GetTrackLabelNagative()).
  29 // Since it is the Position in the Array it is always positive.
  30 // After the conversion ESD --> AOD each AOD track will give you the former Position in the ESDArray via the
  31 // Function AODTrack->GetID(). AODTracks are stored for different TrackParameter (e.g. HybridTracks, TPC only ...). No Standard
  32 // AODTracks (not copies of AliExternalTrackParam) will be stored with negative values of AODTrack->GetID() (Formula: (ESDTrack->GetID()+1)*-1)
  33 // This leads to the possibility of more than one AOD track with the same negative ID. For that reason all AOD tracks are additionally flaged.
  34 // In this analysis we should therefore only use AODTracks with positive values of AODTrack->GetID().
  35
  36 // If you want to find the AODTrack corresponding to the daugher track of a AliAODConversionPhoton you have to find the AODTrack with
  37 // AODTrack->GetID() == GetTrackLabelPositive() (GetTrackLabelNagative()).
  38
  39
  40 #include "AliV0ReaderV1.h"
  41 #include "AliKFParticle.h"
  42 #include "AliAODv0.h"
  43 #include "AliESDv0.h"
  44 #include "AliAODEvent.h"
  45 #include "AliESDEvent.h"
  46 #include "AliPID.h"
  47 #include "AliMCEvent.h"
  48 #include "AliStack.h"
  49 #include "AliMCEventHandler.h"
  50 #include "AliESDpid.h"
  51 #include "AliESDtrackCuts.h"
  52 #include "TRandom3.h"
  53 #include "AliGenCocktailEventHeader.h"
  54 #include "TList.h"
  55 #include "AliKFConversionPhoton.h"
  56 #include "AliAODConversionPhoton.h"
  57 #include "AliConversionPhotonBase.h"
  58 #include "TVector.h"
  59 #include "AliKFVertex.h"
  60 #include "AliAODTrack.h"
  61 #include "AliESDtrack.h"
  62 #include "AliAnalysisManager.h"
  63 #include "AliInputEventHandler.h"
  64 #include "AliAODHandler.h"
  65 #include "AliPIDResponse.h"
  66 #include "TChain.h"
  67 #include "TFile.h"
  68 #include "TString.h"
  69 #include "TObjArray.h"
  70
  71 class iostream;
  72
  73
  74 using namespace std;
  75
  76 ClassImp(AliV0ReaderV1)
  77
  78 //________________________________________________________________________
  79 AliV0ReaderV1::AliV0ReaderV1(const char *name) : AliAnalysisTaskSE(name),
  80    fConversionCuts(NULL),
  81    fConversionGammas(NULL),
  82    fUseImprovedVertex(kTRUE),
  83    fUseOwnXYZCalculation(kTRUE),
  84    fUseConstructGamma(kFALSE),
  85    kUseAODConversionPhoton(kTRUE),
  86    fCreateAOD(kFALSE),
  87    fDeltaAODBranchName("GammaConv"),
  88    fDeltaAODFilename("AliAODGammaConversion.root"),
  89    fRelabelAODs(kFALSE),
  90    fEventIsSelected(kFALSE),
  91    fNumberOfPrimaryTracks(0),
  92    fPeriodName("")
  93 {
  94    // Default constructor
  95
  96    DefineInput(0, TChain::Class());
  97 }
  98
  99 //________________________________________________________________________
 100 AliV0ReaderV1::~AliV0ReaderV1()
 101 {
 102    // default deconstructor
 103
 104    if(fConversionGammas){
 105       fConversionGammas->Delete();// Clear Objects
 106       delete fConversionGammas;
 107       fConversionGammas=0x0;
 108    }
 109 }
 110 /*
 111 //________________________________________________________________________
 112 AliV0ReaderV1::AliV0ReaderV1(AliV0ReaderV1 &original) : AliAnalysisTaskSE(original),
 113 fConversionCuts(NULL),
 114 fConversionGammas(NULL),
 115 fUseImprovedVertex(original.fUseImprovedVertex),
 116 fUseOwnXYZCalculation(original.fUseOwnXYZCalculation),
 117 fUseConstructGamma(original.fUseConstructGamma),
 118 kUseAODConversionPhoton(original.kUseAODConversionPhoton),
 119 fCreateAOD(original.fCreateAOD),
 120 fDeltaAODBranchName(original.fDeltaAODBranchName),
 121 fDeltaAODFilename(original.fDeltaAODFilename),
 122 fEventIsSelected(original.fEventIsSelected)
 123 {
 124 // Default constructor
 125
 126 DefineInput(0, TChain::Class());
 127 }
 128
 129 //____________________________________________________________
 130 AliV0ReaderV1 &AliV0ReaderV1::operator=(const AliV0ReaderV1 &ref){
 131 //
 132 // Assignment operator
 133 // Only copies pointers, object is not the owner of the references
 134 //
 135 if(this != &ref){
 136 AliAnalysisTaskSE::operator=(ref);
 137 fUseImprovedVertex=ref.fUseImprovedVertex;
 138 fUseOwnXYZCalculation=ref.fUseOwnXYZCalculation;
 139 fUseConstructGamma=ref.fUseConstructGamma;
 140 kUseAODConversionPhoton=ref.kUseAODConversionPhoton;
 141 fCreateAOD=ref.fCreateAOD;
 142 fDeltaAODBranchName=ref.fDeltaAODBranchName;
 143 fDeltaAODFilename=ref.fDeltaAODFilename;
 144 fEventIsSelected=ref.fEventIsSelected;
 145 }
 146 return *this;
 147 }
 148 */
 149 //________________________________________________________________________
 150 void AliV0ReaderV1::Init()
 151 {
 152    // Initialize function to be called once before analysis
 153    if(fConversionCuts==NULL){
 154       if(fConversionCuts==NULL)AliError("No Cut Selection initialized");
 155    }
 156    if(fCreateAOD){kUseAODConversionPhoton=kTRUE;}
 157
 158    if(fConversionGammas != NULL){
 159       delete fConversionGammas;
 160       fConversionGammas=NULL;
 161    }
 162
 163    if(fConversionGammas == NULL){
 164       if(kUseAODConversionPhoton){
 165          fConversionGammas = new TClonesArray("AliAODConversionPhoton",100);}
 166       else{
 167          fConversionGammas = new TClonesArray("AliKFConversionPhoton",100);}
 168    }
 169    fConversionGammas->Delete();//Reset the TClonesArray
 170 }
 171
 172 //________________________________________________________________________
 173 void AliV0ReaderV1::UserCreateOutputObjects()
 174 {
 175    // Create AODs
 176
 177    if(fCreateAOD){
 178       if(fConversionCuts){
 179          fDeltaAODBranchName.Append("_");
 180          fDeltaAODBranchName.Append(fConversionCuts->GetCutNumber());
 181          fDeltaAODBranchName.Append("_gamma");
 182       }
 183       fConversionGammas->SetName(fDeltaAODBranchName.Data());
 184
 185       AddAODBranch("TClonesArray", &fConversionGammas, fDeltaAODFilename.Data());
 186       AliAnalysisManager::GetAnalysisManager()->RegisterExtraFile(fDeltaAODFilename.Data());
 187    }
 188
 189 }
 190 //________________________________________________________________________
 191 Bool_t AliV0ReaderV1::Notify()
 192 {
 193    if (fPeriodName.CompareTo("") == 0){
 194       AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 195       if(man) {
 196          AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 197          if (inputHandler){
 198             TTree* tree = (TTree*) inputHandler->GetTree();
 199             TFile* file = (TFile*) tree->GetCurrentFile();
 200             TString fileName(file->GetName());
 201             TObjArray *arr = fileName.Tokenize("/");
 202             for (Int_t i = 0; i < arr->GetEntriesFast();i++ ){
 203                TObjString* testObjString = (TObjString*)arr->At(i);
 204                if (testObjString->GetString().Contains("LHC")){
 205                   fPeriodName = testObjString->GetString();
 206                   i = arr->GetEntriesFast();
 207                }
 208             }
 209          }
 210       }
 211       if(!fConversionCuts->GetDoEtaShift()) return kTRUE; // No Eta Shift requested, continue
 212       if(fConversionCuts->GetEtaShift() == 0.0){ // Eta Shift requested but not set, get shift automatically
 213          fConversionCuts->GetCorrectEtaShiftFromPeriod(fPeriodName);
 214          fConversionCuts->DoEtaShift(kFALSE); // Eta Shift Set, make sure that it is called only once
 215          return kTRUE;
 216       }
 217       else{
 218          printf(" Gamma Conversion Reader %s :: Eta Shift Manually Set to %f \n\n",
 219                 (fConversionCuts->GetCutNumber()).Data(),fConversionCuts->GetEtaShift());
 220          fConversionCuts->DoEtaShift(kFALSE); // Eta Shift Set, make sure that it is called only once
 221       }
 222    }
 223
 224    return kTRUE;
 225 }
 226 //________________________________________________________________________
 227 void AliV0ReaderV1::UserExec(Option_t *option){
 228
 229    // Check if correctly initialized
 230    if(!fConversionGammas)Init();
 231
 232    // User Exec
 233    fEventIsSelected=ProcessEvent(fInputEvent,fMCEvent);
 234
 235    // AOD Output
 236    FillAODOutput();
 237
 238 }
 239
 240 //________________________________________________________________________
 241 Bool_t AliV0ReaderV1::ProcessEvent(AliVEvent *inputEvent,AliMCEvent *mcEvent)
 242 {
 243
 244    //Reset the TClonesArray
 245    fConversionGammas->Delete();
 246
 247    fInputEvent=inputEvent;
 248    fMCEvent=mcEvent;
 249
 250    if(!fInputEvent){
 251       AliError("No Input event");
 252       return kFALSE;
 253    }
 254    if(!fConversionCuts){AliError("No ConversionCuts");return kFALSE;}
 255
 256    // Count Primary Tracks Event
 257    CountTracks();
 258
 259    // Event Cuts
 260    if(!fConversionCuts->EventIsSelected(fInputEvent,fMCEvent))return kFALSE;
 261
 262    // Set Magnetic Field
 263    AliKFParticle::SetField(fInputEvent->GetMagneticField());
 264
 265    if(fInputEvent->IsA()==AliESDEvent::Class()){
 266       ProcessESDV0s();
 267    }
 268    if(fInputEvent->IsA()==AliAODEvent::Class()){
 269       GetAODConversionGammas();
 270    }
 271
 272
 273
 274    return kTRUE;
 275 }
 276 ///________________________________________________________________________
 277 void AliV0ReaderV1::FillAODOutput()
 278 {
 279    // Fill AOD Output with reconstructed Photons
 280
 281    if(fInputEvent->IsA()==AliESDEvent::Class()){
 282       ///Make sure delta aod is filled if standard aod is filled (for synchronization when reading aod with standard aod)
 283       if(fCreateAOD) {
 284          AliAODHandler * aodhandler = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 285          if (aodhandler && aodhandler->GetFillAOD()) {
 286             AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillExtension(kTRUE);
 287             //PostData(0, fConversionGammas);
 288
 289          }
 290       }
 291    }
 292 }
 293
 294 ///________________________________________________________________________
 295 const AliExternalTrackParam *AliV0ReaderV1::GetExternalTrackParam(AliESDv0 *fCurrentV0,Int_t &tracklabel,Int_t charge){
 296
 297    // Get External Track Parameter with given charge
 298
 299    if(!(charge==1||charge==-1)){AliError("Charge not defined");return 0x0;}
 300
 301    // Check for sign flip
 302    if(fCurrentV0){
 303       if(!fCurrentV0->GetParamN()||!fCurrentV0->GetParamP())return 0x0;
 304       if(!fConversionCuts->GetTrack(fInputEvent,fCurrentV0->GetNindex())||!fConversionCuts->GetTrack(fInputEvent,fCurrentV0->GetPindex()))return 0x0;
 305       if((fConversionCuts->GetTrack(fInputEvent,fCurrentV0->GetPindex()))->Charge()==charge){
 306          tracklabel=fCurrentV0->GetPindex();
 307          return fCurrentV0->GetParamP();}
 308       if((fConversionCuts->GetTrack(fInputEvent,fCurrentV0->GetNindex()))->Charge()==charge){
 309          tracklabel=fCurrentV0->GetNindex();
 310          return fCurrentV0->GetParamN();}
 311    }
 312    return 0x0;
 313 }
 314
 315 ///________________________________________________________________________
 316 Bool_t AliV0ReaderV1::ProcessESDV0s()
 317 {
 318    // Process ESD V0s for conversion photon reconstruction
 319    AliESDEvent *fESDEvent=dynamic_cast<AliESDEvent*>(fInputEvent);
 320
 321    AliKFConversionPhoton *fCurrentMotherKFCandidate=NULL;
 322
 323    if(fESDEvent){
 324
 325       for(Int_t currentV0Index=0;currentV0Index<fESDEvent->GetNumberOfV0s();currentV0Index++){
 326          AliESDv0 *fCurrentV0=(AliESDv0*)(fESDEvent->GetV0(currentV0Index));
 327          if(!fCurrentV0){
 328             printf("Requested V0 does not exist");
 329             continue;
 330          }
 331
 332          fCurrentMotherKFCandidate=ReconstructV0(fCurrentV0,currentV0Index);
 333
 334          if(fCurrentMotherKFCandidate){
 335             // Add Gamma to the TClonesArray
 336
 337             if(kUseAODConversionPhoton){
 338                new((*fConversionGammas)[fConversionGammas->GetEntriesFast()]) AliAODConversionPhoton(fCurrentMotherKFCandidate);
 339                AliAODConversionPhoton * currentConversionPhoton = (AliAODConversionPhoton*)(fConversionGammas->At(fConversionGammas->GetEntriesFast()-1));
 340                currentConversionPhoton->SetMass(fCurrentMotherKFCandidate->M());
 341                currentConversionPhoton->SetMassToZero();
 342
 343             }
 344             else{
 345                new((*fConversionGammas)[fConversionGammas->GetEntriesFast()]) AliKFConversionPhoton(*fCurrentMotherKFCandidate);
 346             }
 347
 348             delete fCurrentMotherKFCandidate;
 349             fCurrentMotherKFCandidate=NULL;
 350          }
 351       }
 352    }
 353    return kTRUE;
 354 }
 355
 356 ///________________________________________________________________________
 357 AliKFConversionPhoton *AliV0ReaderV1::ReconstructV0(AliESDv0 *fCurrentV0,Int_t currentV0Index)
 358 {
 359    // Reconstruct conversion photon from ESD v0
 360    fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kPhotonIn);
 361
 362    //checks if on the fly mode is set
 363    if(!fConversionCuts->SelectV0Finder(fCurrentV0->GetOnFlyStatus())){
 364       fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kOnFly);
 365       return 0x0;
 366    }
 367
 368    // TrackLabels
 369    Int_t currentTrackLabels[2]={-1,-1};
 370
 371    // Get Daughter KF Particles
 372
 373    const AliExternalTrackParam *fCurrentExternalTrackParamPositive=GetExternalTrackParamP(fCurrentV0,currentTrackLabels[0]);
 374    const AliExternalTrackParam *fCurrentExternalTrackParamNegative=GetExternalTrackParamN(fCurrentV0,currentTrackLabels[1]);
 375
 376    if(!fCurrentExternalTrackParamPositive||!fCurrentExternalTrackParamNegative)return 0x0;
 377
 378    // Apply some Cuts before Reconstruction
 379
 380    AliVTrack * posTrack = fConversionCuts->GetTrack(fInputEvent,currentTrackLabels[0]);
 381    AliVTrack * negTrack = fConversionCuts->GetTrack(fInputEvent,currentTrackLabels[1]);
 382
 383    if(!negTrack || !posTrack) {
 384       fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kNoTracks);
 385       return 0x0;
 386    }
 387
 388    // Track Cuts
 389    if(!fConversionCuts->TracksAreSelected(negTrack, posTrack)){
 390       fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kTrackCuts);
 391       return 0x0;
 392    }
 393
 394    // PID Cuts
 395    if(!fConversionCuts->dEdxCuts(negTrack) || !fConversionCuts->dEdxCuts(posTrack)) {
 396       fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kdEdxCuts);
 397       return 0x0;
 398    }
 399
 400    // Reconstruct Photon
 401
 402    AliKFConversionPhoton *fCurrentMotherKF=NULL;
 403
 404    AliKFParticle fCurrentNegativeKFParticle(*(fCurrentExternalTrackParamNegative),11);
 405    AliKFParticle fCurrentPositiveKFParticle(*(fCurrentExternalTrackParamPositive),-11);
 406
 407    // Reconstruct Gamma
 408
 409    if(fUseConstructGamma){
 410       fCurrentMotherKF = new AliKFConversionPhoton();
 411       fCurrentMotherKF->ConstructGamma(fCurrentNegativeKFParticle,fCurrentPositiveKFParticle);
 412    }else{
 413       fCurrentMotherKF = new AliKFConversionPhoton(fCurrentNegativeKFParticle,fCurrentPositiveKFParticle);
 414       fCurrentMotherKF->SetMassConstraint(0,0.0001);
 415    }
 416
 417    // Set Track Labels
 418
 419    fCurrentMotherKF->SetTrackLabels(currentTrackLabels[0],currentTrackLabels[1]);
 420
 421    // Set V0 index
 422
 423    fCurrentMotherKF->SetV0Index(currentV0Index);
 424
 425    //Set MC Label
 426
 427    if(fMCEvent){
 428
 429       AliStack *fMCStack= fMCEvent->Stack();
 430
 431       Int_t labelp=TMath::Abs(fConversionCuts->GetTrack(fInputEvent,fCurrentMotherKF->GetTrackLabelPositive())->GetLabel());
 432       Int_t labeln=TMath::Abs(fConversionCuts->GetTrack(fInputEvent,fCurrentMotherKF->GetTrackLabelNegative())->GetLabel());
 433
 434       TParticle *fNegativeMCParticle = fMCStack->Particle(labeln);
 435       TParticle *fPositiveMCParticle = fMCStack->Particle(labelp);
 436
 437       if(fPositiveMCParticle&&fNegativeMCParticle){
 438          fCurrentMotherKF->SetMCLabelPositive(labelp);
 439          fCurrentMotherKF->SetMCLabelNegative(labeln);
 440       }
 441    }
 442
 443    // Update Vertex (moved for same eta compared to old)
 444    //      cout << currentV0Index <<" \t before: \t" << fCurrentMotherKF->GetPx() << "\t" << fCurrentMotherKF->GetPy() << "\t" << fCurrentMotherKF->GetPz()  << endl;
 445    if(fUseImprovedVertex == kTRUE){
 446       AliKFVertex primaryVertexImproved(*fInputEvent->GetPrimaryVertex());
 447       //        cout << "Prim Vtx: " << primaryVertexImproved.GetX() << "\t" << primaryVertexImproved.GetY() << "\t" << primaryVertexImproved.GetZ() << endl;
 448       primaryVertexImproved+=*fCurrentMotherKF;
 449       fCurrentMotherKF->SetProductionVertex(primaryVertexImproved);
 450    }
 451    // SetPsiPair
 452
 453    Double_t PsiPair=GetPsiPair(fCurrentV0,fCurrentExternalTrackParamPositive,fCurrentExternalTrackParamNegative);
 454    fCurrentMotherKF->SetPsiPair(PsiPair);
 455
 456
 457    // Recalculate ConversionPoint
 458    Double_t dca[2]={0,0};
 459    if(fUseOwnXYZCalculation){
 460       Double_t convpos[3]={0,0,0};
 461       if(!GetConversionPoint(fCurrentExternalTrackParamPositive,fCurrentExternalTrackParamNegative,convpos,dca)){
 462          fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kConvPointFail);
 463          delete fCurrentMotherKF;
 464          fCurrentMotherKF=NULL;
 465          return 0x0;
 466       }
 467
 468       fCurrentMotherKF->SetConversionPoint(convpos);
 469    }
 470
 471    if(fCurrentMotherKF->GetNDF() > 0.)
 472       fCurrentMotherKF->SetChi2perNDF(fCurrentMotherKF->GetChi2()/fCurrentMotherKF->GetNDF());   //->Photon is created before all chi2 relevant changes are performed, set it "by hand"
 473
 474
 475    // Set Dilepton Mass (moved down for same eta compared to old)
 476    fCurrentMotherKF->SetMass(fCurrentMotherKF->M());
 477
 478    // Apply Photon Cuts
 479
 480    if(!fConversionCuts->PhotonCuts(fCurrentMotherKF,fInputEvent)){
 481       fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kPhotonCuts);
 482       delete fCurrentMotherKF;
 483       fCurrentMotherKF=NULL;
 484       return 0x0;
 485    }
 486
 487    //     cout << currentV0Index <<" \t after: \t" <<fCurrentMotherKF->GetPx() << "\t" << fCurrentMotherKF->GetPy() << "\t" << fCurrentMotherKF->GetPz()  << endl;
 488
 489    fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kPhotonOut);
 490    return fCurrentMotherKF;
 491 }
 492
 493 ///________________________________________________________________________
 494 Double_t AliV0ReaderV1::GetPsiPair(const AliESDv0* v0, const AliExternalTrackParam *positiveparam,const AliExternalTrackParam *negativeparam) const {
 495    //
 496    // Angle between daughter momentum plane and plane
 497    //
 498
 499    AliExternalTrackParam nt(*negativeparam);
 500    AliExternalTrackParam pt(*positiveparam);
 501
 502    Float_t magField = fInputEvent->GetMagneticField();
 503
 504    Double_t xyz[3] = {0.,0.,0.};
 505    v0->GetXYZ(xyz[0],xyz[1],xyz[2]);
 506
 507    // Double_t pPlus[3]  = {pt.Px(),pt.Py(),pt.Pz()};
 508    // Double_t pMinus[3] = {nt.Px(),nt.Py(),nt.Pz()};
 509
 510    // Double_t u[3] = {pPlus[0]+pMinus[0],pPlus[1]+pMinus[1],pPlus[2]+pMinus[2]};
 511    // Double_t normu = sqrt( (u[0]*u[0]) + (u[1]*u[1]) + (u[2]*u[2]) );
 512
 513    // u[0] = u[0] / normu;
 514    // u[1] = u[1] / normu;
 515    // u[2] = u[2] / normu;
 516
 517    // Double_t normpPlus = sqrt( (pPlus[0]*pPlus[0]) + (pPlus[1]*pPlus[1]) + (pPlus[2]*pPlus[2]) );
 518    // Double_t normpMinus = sqrt( (pMinus[0]*pMinus[0]) + (pMinus[1]*pMinus[1]) + (pMinus[2]*pMinus[2]) );
 519
 520    // pPlus[0] = pPlus[0] / normpPlus;
 521    // pPlus[1] = pPlus[1] / normpPlus;
 522    // pPlus[2] = pPlus[2] / normpPlus;
 523
 524    // pMinus[0] = pMinus[0] / normpMinus;
 525    // pMinus[1] = pMinus[1] / normpMinus;
 526    // pMinus[2] = pMinus[2] / normpMinus;
 527
 528    // Double_t v[3] = {0,0,0}; // pPlus X pMinus
 529    // v[0] = (pPlus[1]*pMinus[2]) - (pPlus[2]*pMinus[1]);
 530    // v[1] = (pPlus[2]*pMinus[0]) - (pPlus[0]*pMinus[2]);
 531    // v[2] = (pPlus[0]*pMinus[1]) - (pPlus[1]*pMinus[0]);
 532
 533    // Double_t w[3] = {0,0,0}; // u X v
 534    // w[0] = (u[1]*v[2]) - (u[2]*v[1]);
 535    // w[1] = (u[2]*v[0]) - (u[0]*v[2]);
 536    // w[2] = (u[0]*v[1]) - (u[1]*v[0]);
 537
 538    // Double_t z[3] = {0,0,1};
 539    // Double_t wc[3] = {0,0,0}; // u X v
 540    // wc[0] = (u[1]*z[2]) - (u[2]*z[1]);
 541    // wc[1] = (u[2]*z[0]) - (u[0]*z[2]);
 542    // wc[2] = (u[0]*z[1]) - (u[1]*z[0]);
 543
 544    // Double_t PhiV = TMath::ACos((w[0]*wc[0]) + (w[1]*wc[1]) + (w[2]*wc[2]));
 545    //return abs(PhiV);
 546
 547
 548    // TVector3 pPlus(pt.Px(),pt.Py(),pt.Pz());
 549    // TVector3 pMinus(nt.Px(),nt.Py(),nt.Pz());
 550
 551    // TVector3 u = pMinus + pPlus;
 552    // u = u*(1/u.Mag());
 553
 554    // TVector3 pHPlus = pPlus*(1/pPlus.Mag());
 555    // TVector3 pHMinus = pMinus*(1/pMinus.Mag());
 556
 557    // TVector3 v = pHPlus.Cross(pHMinus);
 558    // TVector3 w = u.Cross(v);
 559    // TVector3 z(0,0,1);
 560    // TVector3 wc = u.Cross(z);
 561
 562    // Double_t PhiV = w * wc;
 563
 564    Double_t mn[3] = {0,0,0};
 565    Double_t mp[3] = {0,0,0};
 566
 567    v0->GetNPxPyPz(mn[0],mn[1],mn[2]);//reconstructed cartesian momentum components of negative daughter;
 568    v0->GetPPxPyPz(mp[0],mp[1],mp[2]);//reconstructed cartesian momentum components of positive daughter;
 569
 570    Double_t deltat = 1.;
 571    deltat = TMath::ATan(mp[2]/(TMath::Sqrt(mp[0]*mp[0] + mp[1]*mp[1])+1.e-13)) - TMath::ATan(mn[2]/(TMath::Sqrt(mn[0]*mn[0] + mn[1]*mn[1])+1.e-13));//difference of angles of the two daughter tracks with z-axis
 572    Double_t radiussum = TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]) + 50;//radius to which tracks shall be propagated
 573
 574    Double_t momPosProp[3] = {0,0,0};
 575    Double_t momNegProp[3] = {0,0,0};
 576
 577    Double_t psiPair = 4.;
 578    if(nt.PropagateTo(radiussum,magField) == 0) return psiPair; //propagate tracks to the outside -> Better Purity and Efficiency
 579    if(pt.PropagateTo(radiussum,magField) == 0) return psiPair; //propagate tracks to the outside -> Better Purity and Efficiency
 580
 581    pt.GetPxPyPz(momPosProp);//Get momentum vectors of tracks after propagation
 582    nt.GetPxPyPz(momNegProp);
 583
 584    Double_t pEle =
 585       TMath::Sqrt(momNegProp[0]*momNegProp[0]+momNegProp[1]*momNegProp[1]+momNegProp[2]*momNegProp[2]);//absolute momentum value of negative daughter
 586
 587    Double_t pPos =
 588       TMath::Sqrt(momPosProp[0]*momPosProp[0]+momPosProp[1]*momPosProp[1]+momPosProp[2]*momPosProp[2]);//absolute momentum value of positive daughter
 589
 590    Double_t scalarproduct =
 591       momPosProp[0]*momNegProp[0]+momPosProp[1]*momNegProp[1]+momPosProp[2]*momNegProp[2];//scalar product of propagated positive and negative daughters' momenta
 592
 593    Double_t chipair = TMath::ACos(scalarproduct/(pEle*pPos));//Angle between propagated daughter tracks
 594
 595 //    psiPair =  TMath::Abs(TMath::ASin(deltat/chipair));
 596    psiPair = TMath::ASin(deltat/chipair);
 597    return psiPair;
 598 }
 599
 600 ///________________________________________________________________________
 601 Bool_t AliV0ReaderV1::GetHelixCenter(const AliExternalTrackParam *track,Double_t center[2]){
 602
 603    // Get Center of the helix track parametrization
 604
 605    Int_t charge=track->Charge();
 606    Double_t b=fInputEvent->GetMagneticField();
 607
 608    Double_t     helix[6];
 609    track->GetHelixParameters(helix,b);
 610
 611    Double_t xpos =      helix[5];
 612    Double_t ypos =      helix[0];
 613    Double_t radius = TMath::Abs(1./helix[4]);
 614    Double_t phi = helix[2];
 615
 616    if(phi < 0){
 617       phi = phi + 2*TMath::Pi();
 618    }
 619
 620    phi -= TMath::Pi()/2.;
 621    Double_t xpoint =    radius * TMath::Cos(phi);
 622    Double_t ypoint =    radius * TMath::Sin(phi);
 623
 624    if(b<0){
 625       if(charge > 0){
 626          xpoint = - xpoint;
 627          ypoint = - ypoint;
 628       }
 629
 630       if(charge < 0){
 631          xpoint =       xpoint;
 632          ypoint =       ypoint;
 633       }
 634    }
 635    if(b>0){
 636       if(charge > 0){
 637          xpoint =       xpoint;
 638          ypoint =       ypoint;
 639       }
 640
 641       if(charge < 0){
 642          xpoint = - xpoint;
 643          ypoint = - ypoint;
 644       }
 645    }
 646    center[0] =  xpos + xpoint;
 647    center[1] =  ypos + ypoint;
 648
 649    return 1;
 650 }
 651 ///________________________________________________________________________
 652 Bool_t AliV0ReaderV1::GetConversionPoint(const AliExternalTrackParam *pparam,const AliExternalTrackParam *nparam,Double_t convpos[3],Double_t dca[2]){
 653
 654    // Recalculate Conversion Point
 655
 656    if(!pparam||!nparam)return kFALSE;
 657
 658    Double_t helixcenterpos[2];
 659    GetHelixCenter(pparam,helixcenterpos);
 660
 661    Double_t helixcenterneg[2];
 662    GetHelixCenter(nparam,helixcenterneg);
 663
 664    Double_t helixpos[6];
 665    pparam->GetHelixParameters(helixpos,fInputEvent->GetMagneticField());
 666    Double_t posradius = TMath::Abs(1./helixpos[4]);
 667
 668    Double_t helixneg[6];
 669    nparam->GetHelixParameters(helixneg,fInputEvent->GetMagneticField());
 670    Double_t negradius = TMath::Abs(1./helixneg[4]);
 671
 672    // Calculate xy-position
 673
 674    Double_t xpos = helixcenterpos[0];
 675    Double_t ypos = helixcenterpos[1];
 676    Double_t xneg = helixcenterneg[0];
 677    Double_t yneg = helixcenterneg[1];
 678
 679    convpos[0] = (xpos*negradius + xneg*posradius)/(negradius+posradius);
 680    convpos[1] = (ypos*negradius+ yneg*posradius)/(negradius+posradius);
 681
 682
 683    // Calculate Track XY vertex position
 684
 685    Double_t deltaXPos = convpos[0] -    xpos;
 686    Double_t deltaYPos = convpos[1] -    ypos;
 687
 688    Double_t deltaXNeg = convpos[0] -    xneg;
 689    Double_t deltaYNeg = convpos[1] -    yneg;
 690
 691    Double_t alphaPos =  TMath::Pi() + TMath::ATan2(-deltaYPos,-deltaXPos);
 692    Double_t alphaNeg =  TMath::Pi() + TMath::ATan2(-deltaYNeg,-deltaXNeg);
 693
 694    Double_t vertexXNeg =        xneg +  TMath::Abs(negradius)*TMath::Cos(alphaNeg);
 695    Double_t vertexYNeg =        yneg +  TMath::Abs(negradius)*TMath::Sin(alphaNeg);
 696
 697    Double_t vertexXPos =        xpos +  TMath::Abs(posradius)*TMath::Cos(alphaPos);
 698    Double_t vertexYPos =        ypos +  TMath::Abs(posradius)*TMath::Sin(alphaPos);
 699
 700    AliExternalTrackParam p(*pparam);
 701    AliExternalTrackParam n(*nparam);
 702
 703    TVector2 vertexPos(vertexXPos,vertexYPos);
 704    TVector2 vertexNeg(vertexXNeg,vertexYNeg);
 705
 706    // Convert to local coordinate system
 707    TVector2 vertexPosRot=vertexPos.Rotate(-p.GetAlpha());
 708    TVector2 vertexNegRot=vertexNeg.Rotate(-n.GetAlpha());
 709
 710    // Propagate Track Params to Vertex
 711
 712    if(!p.PropagateTo(vertexPosRot.X(),fInputEvent->GetMagneticField()))return kFALSE;
 713    if(!n.PropagateTo(vertexNegRot.X(),fInputEvent->GetMagneticField()))return kFALSE;
 714
 715    // Check whether propagation was sucessful
 716
 717    if(TMath::Abs(vertexPos.Mod()-TMath::Sqrt(p.GetX()*p.GetX()+p.GetY()*p.GetY()))>0.01)return kFALSE;
 718    if(TMath::Abs(vertexNeg.Mod()-TMath::Sqrt(n.GetX()*n.GetX()+n.GetY()*n.GetY()))>0.01)return kFALSE;
 719
 720    // Calculate z position
 721
 722    convpos[2] = (p.GetZ()*negradius+n.GetZ()*posradius)/(negradius+posradius);
 723
 724    // Calculate DCA
 725    TVector2 vdca=vertexPos-vertexNeg;
 726    dca[0]=vdca.Mod();
 727    dca[1]=TMath::Abs(n.GetZ()-p.GetZ());
 728
 729    return kTRUE;
 730 }
 731 //________________________________________________________________________
 732 Bool_t AliV0ReaderV1::GetAODConversionGammas(){
 733
 734    // Get reconstructed conversion photons from satellite AOD file
 735
 736    AliAODEvent *fAODEvent=dynamic_cast<AliAODEvent*>(fInputEvent);
 737
 738    if(fAODEvent){
 739
 740       if(fConversionGammas == NULL){
 741          fConversionGammas = new TClonesArray("AliAODConversionPhoton",100);
 742       }
 743       fConversionGammas->Delete();//Reset the TClonesArray
 744
 745       //Get Gammas from satellite AOD gamma branch
 746
 747       AliAODConversionPhoton *gamma=0x0;
 748
 749       TClonesArray *fInputGammas=dynamic_cast<TClonesArray*>(fAODEvent->FindListObject(fDeltaAODBranchName.Data()));
 750       if(!fInputGammas){
 751          FindDeltaAODBranchName();
 752          fInputGammas=dynamic_cast<TClonesArray*>(fAODEvent->FindListObject(fDeltaAODBranchName.Data()));}
 753       if(!fInputGammas){AliError("No Gamma Satellites found");return kFALSE;}
 754       // Apply Selection Cuts to Gammas and create local working copy
 755       if(fInputGammas){
 756          for(Int_t i=0;i<fInputGammas->GetEntriesFast();i++){
 757             gamma=dynamic_cast<AliAODConversionPhoton*>(fInputGammas->At(i));
 758             if(gamma){
 759                if(fRelabelAODs)RelabelAODPhotonCandidates(gamma);
 760                if(fConversionCuts->PhotonIsSelected(gamma,fInputEvent)){
 761                   new((*fConversionGammas)[fConversionGammas->GetEntriesFast()]) AliAODConversionPhoton(*gamma);}
 762             }
 763          }
 764       }
 765    }
 766
 767    if(fConversionGammas->GetEntries()){return kTRUE;}
 768
 769    return kFALSE;
 770 }
 771
 772 //________________________________________________________________________
 773 void AliV0ReaderV1::FindDeltaAODBranchName(){
 774
 775    // Find delta AOD branchname containing reconstructed photons
 776
 777    TList *list=fInputEvent->GetList();
 778    for(Int_t ii=0;ii<list->GetEntries();ii++){
 779       TString name((list->At(ii))->GetName());
 780       if(name.BeginsWith(fDeltaAODBranchName)&&name.EndsWith("gamma")){
 781          fDeltaAODBranchName=name;
 782          AliInfo(Form("Set DeltaAOD BranchName to: %s",fDeltaAODBranchName.Data()));
 783          return;
 784       }
 785    }
 786 }
 787 //________________________________________________________________________
 788 void AliV0ReaderV1::RelabelAODPhotonCandidates(AliAODConversionPhoton *PhotonCandidate){
 789
 790    // Relabeling For AOD Event
 791    // ESDiD -> AODiD
 792    // MCLabel -> AODMCLabel
 793    Bool_t AODLabelPos = kFALSE;
 794    Bool_t AODLabelNeg = kFALSE;
 795
 796    for(Int_t i = 0; i<fInputEvent->GetNumberOfTracks();i++){
 797       AliAODTrack *tempDaughter = static_cast<AliAODTrack*>(fInputEvent->GetTrack(i));
 798       if(!AODLabelPos){
 799          if( tempDaughter->GetID() == PhotonCandidate->GetTrackLabelPositive() ){
 800             PhotonCandidate->SetMCLabelPositive(abs(tempDaughter->GetLabel()));
 801             PhotonCandidate->SetLabelPositive(i);
 802             AODLabelPos = kTRUE;
 803          }
 804       }
 805       if(!AODLabelNeg){
 806          if( tempDaughter->GetID() == PhotonCandidate->GetTrackLabelNegative()){
 807             PhotonCandidate->SetMCLabelNegative(abs(tempDaughter->GetLabel()));
 808             PhotonCandidate->SetLabelNegative(i);
 809             AODLabelNeg = kTRUE;
 810          }
 811       }
 812       if(AODLabelNeg && AODLabelPos){
 813          return;
 814       }
 815    }
 816    if(!AODLabelPos || !AODLabelNeg){
 817       AliError(Form("NO AOD Daughters Found Pos: %i %i Neg: %i %i",AODLabelPos,PhotonCandidate->GetTrackLabelPositive(),AODLabelNeg,PhotonCandidate->GetTrackLabelNegative()));
 818    }
 819
 820 }
 821
 822 //________________________________________________________________________
 823 void AliV0ReaderV1::CountTracks(){
 824
 825    // In Case of MC count only MB tracks
 826    // if(fMCEvent && fInputEvent->IsA()==AliESDEvent::Class()){
 827    //    fConversionCuts->GetNotRejectedParticles(1,NULL,fMCEvent);
 828    // }
 829    // if(fMCEvent && fInputEvent->IsA()==AliAODEvent::Class()){
 830    //    fConversionCuts->GetNotRejectedParticles(1,NULL,fInputEvent);
 831    // }
 832
 833    if(fInputEvent->IsA()==AliESDEvent::Class()){
 834       // Using standard function for setting Cuts
 835       Bool_t selectPrimaries=kTRUE;
 836       AliESDtrackCuts *EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(selectPrimaries);
 837       EsdTrackCuts->SetMaxDCAToVertexZ(2);
 838       EsdTrackCuts->SetEtaRange(-0.8, 0.8);
 839       EsdTrackCuts->SetPtRange(0.15);
 840       fNumberOfPrimaryTracks = 0;
 841       for(Int_t iTracks = 0; iTracks < fInputEvent->GetNumberOfTracks(); iTracks++){
 842          AliESDtrack* curTrack = (AliESDtrack*) fInputEvent->GetTrack(iTracks);
 843          if(!curTrack) continue;
 844          if(!EsdTrackCuts->AcceptTrack(curTrack)) continue;
 845          //if(fMCEvent && !(fConversionCuts->IsParticleFromBGEvent(abs(curTrack->GetLabel()),fMCEvent->Stack(),fInputEvent))) continue;
 846          fNumberOfPrimaryTracks++;
 847       }
 848       delete EsdTrackCuts;
 849       EsdTrackCuts=0x0;
 850    }
 851    else if(fInputEvent->IsA()==AliAODEvent::Class()){
 852       fNumberOfPrimaryTracks = 0;
 853       for(Int_t iTracks = 0; iTracks<fInputEvent->GetNumberOfTracks(); iTracks++){
 854          AliAODTrack* curTrack = (AliAODTrack*) fInputEvent->GetTrack(iTracks);
 855          if(curTrack->GetID()<0) continue; // Avoid double counting of tracks
 856          if(!curTrack->IsPrimaryCandidate()) continue;
 857          if(abs(curTrack->Eta())>0.8) continue;
 858          if(curTrack->Pt()<0.15) continue;
 859          //if(fMCEvent && !(fConversionCuts->IsParticleFromBGEvent(abs(curTrack->GetLabel()),NULL,fInputEvent))) continue;
 860          //if(abs(curTrack->ZAtDCA())>2) continue; // Only Set For TPCOnly tracks
 861          fNumberOfPrimaryTracks++;
 862       }
 863    }
 864
 865    return;
 866 }
 867
 868 //________________________________________________________________________
 869 void AliV0ReaderV1::Terminate(Option_t *)
 870 {
 871
 872 }