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 #include "AliV0ReaderV1.h"
  28 #include "AliKFParticle.h"
  29 #include "AliAODv0.h"
  30 #include "AliESDv0.h"
  31 #include "AliAODEvent.h"
  32 #include "AliESDEvent.h"
  33 #include "AliPID.h"
  34 #include "AliMCEvent.h"
  35 #include "AliStack.h"
  36 #include "AliMCEventHandler.h"
  37 #include "AliESDpid.h"
  38 #include "AliESDtrackCuts.h"
  39 #include "TRandom3.h"
  40 #include "AliGenCocktailEventHeader.h"
  41 #include "TList.h"
  42 #include "AliKFConversionPhoton.h"
  43 #include "AliAODConversionPhoton.h"
  44 #include "AliConversionPhotonBase.h"
  45 #include "TVector.h"
  46 #include "AliKFVertex.h"
  47 #include "AliAODTrack.h"
  48 #include "AliESDtrack.h"
  49 #include "AliAnalysisManager.h"
  50 #include "AliInputEventHandler.h"
  51 #include "AliAODHandler.h"
  52 #include "AliPIDResponse.h"
  53 #include "TChain.h"
  54 #include "TFile.h"
  55 #include "TString.h"
  56 #include "TObjArray.h"
  57
  58 class iostream;
  59
  60
  61 using namespace std;
  62
  63 ClassImp(AliV0ReaderV1)
  64
  65 //________________________________________________________________________
  66 AliV0ReaderV1::AliV0ReaderV1(const char *name) : AliAnalysisTaskSE(name),
  67     fConversionCuts(NULL),
  68     fConversionGammas(NULL),
  69     fUseImprovedVertex(kTRUE),
  70     fUseOwnXYZCalculation(kTRUE),
  71     fUseConstructGamma(kFALSE),
  72     kUseAODConversionPhoton(kTRUE),
  73     fCreateAOD(kFALSE),
  74     fDeltaAODBranchName("GammaConv"),
  75     fDeltaAODFilename("AliAODGammaConversion.root"),
  76     fCheckAODConsistenty(kFALSE),
  77     fEventIsSelected(kFALSE),
  78     fPeriodName("")
  79 {
  80     // Default constructor
  81
  82     DefineInput(0, TChain::Class());
  83 }
  84
  85 //________________________________________________________________________
  86 AliV0ReaderV1::~AliV0ReaderV1()
  87 {
  88     // default deconstructor
  89
  90     if(fConversionGammas){
  91         fConversionGammas->Delete();// Clear Objects
  92         delete fConversionGammas;
  93         fConversionGammas=0x0;
  94     }
  95 }
  96 /*
  97 //________________________________________________________________________
  98 AliV0ReaderV1::AliV0ReaderV1(AliV0ReaderV1 &original) : AliAnalysisTaskSE(original),
  99     fConversionCuts(NULL),
 100     fConversionGammas(NULL),
 101     fUseImprovedVertex(original.fUseImprovedVertex),
 102     fUseOwnXYZCalculation(original.fUseOwnXYZCalculation),
 103     fUseConstructGamma(original.fUseConstructGamma),
 104     kUseAODConversionPhoton(original.kUseAODConversionPhoton),
 105     fCreateAOD(original.fCreateAOD),
 106     fDeltaAODBranchName(original.fDeltaAODBranchName),
 107     fDeltaAODFilename(original.fDeltaAODFilename),
 108     fEventIsSelected(original.fEventIsSelected)
 109 {
 110     // Default constructor
 111
 112     DefineInput(0, TChain::Class());
 113 }
 114
 115 //____________________________________________________________
 116 AliV0ReaderV1 &AliV0ReaderV1::operator=(const AliV0ReaderV1 &ref){
 117         //
 118         // Assignment operator
 119         // Only copies pointers, object is not the owner of the references
 120         //
 121     if(this != &ref){
 122         AliAnalysisTaskSE::operator=(ref);
 123         fUseImprovedVertex=ref.fUseImprovedVertex;
 124         fUseOwnXYZCalculation=ref.fUseOwnXYZCalculation;
 125         fUseConstructGamma=ref.fUseConstructGamma;
 126         kUseAODConversionPhoton=ref.kUseAODConversionPhoton;
 127         fCreateAOD=ref.fCreateAOD;
 128         fDeltaAODBranchName=ref.fDeltaAODBranchName;
 129         fDeltaAODFilename=ref.fDeltaAODFilename;
 130         fEventIsSelected=ref.fEventIsSelected;
 131     }
 132     return *this;
 133 }
 134 */
 135 //________________________________________________________________________
 136 void AliV0ReaderV1::Init()
 137 {
 138     // Initialize function to be called once before analysis
 139     if(fConversionCuts==NULL){
 140         if(fConversionCuts==NULL)AliError("No Cut Selection initialized");
 141     }
 142     if(fCreateAOD){kUseAODConversionPhoton=kTRUE;}
 143
 144     if(fConversionGammas != NULL){
 145         delete fConversionGammas;
 146         fConversionGammas=NULL;
 147     }
 148
 149     if(fConversionGammas == NULL){
 150         if(kUseAODConversionPhoton){
 151             fConversionGammas = new TClonesArray("AliAODConversionPhoton",100);}
 152         else{
 153             fConversionGammas = new TClonesArray("AliKFConversionPhoton",100);}
 154     }
 155     fConversionGammas->Delete();//Reset the TClonesArray
 156 }
 157
 158 //________________________________________________________________________
 159 void AliV0ReaderV1::UserCreateOutputObjects()
 160 {
 161     // Create AODs
 162
 163     if(fCreateAOD){
 164         if(fConversionCuts){
 165             fDeltaAODBranchName.Append("_");
 166             fDeltaAODBranchName.Append(fConversionCuts->GetCutNumber());
 167             fDeltaAODBranchName.Append("_gamma");
 168         }
 169         fConversionGammas->SetName(fDeltaAODBranchName.Data());
 170
 171         AddAODBranch("TClonesArray", &fConversionGammas, fDeltaAODFilename.Data());
 172         AliAnalysisManager::GetAnalysisManager()->RegisterExtraFile(fDeltaAODFilename.Data());
 173     }
 174
 175 }
 176 //________________________________________________________________________
 177 Bool_t AliV0ReaderV1::Notify()
 178 {
 179    if (fPeriodName.CompareTo("") == 0){
 180       AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 181       if(man) {
 182          AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 183          if (inputHandler){
 184             TTree* tree = (TTree*) inputHandler->GetTree();
 185             TFile* file = (TFile*) tree->GetCurrentFile();
 186             TString fileName(file->GetName());
 187             TObjArray *arr = fileName.Tokenize("/");
 188             for (Int_t i = 0; i < arr->GetEntriesFast();i++ ){
 189                TObjString* testObjString = (TObjString*)arr->At(i);
 190                if (testObjString->GetString().Contains("LHC")){
 191                   fPeriodName = testObjString->GetString();
 192                   i = arr->GetEntriesFast();
 193                }
 194             }
 195          }
 196       }
 197       if(!fConversionCuts->GetDoEtaShift()) return kTRUE; // No Eta Shift requested, continue
 198       if(fConversionCuts->GetEtaShift() == 0.0){ // Eta Shift requested but not set, get shift automatically
 199          fConversionCuts->GetCorrectEtaShiftFromPeriod(fPeriodName);
 200          fConversionCuts->DoEtaShift(kFALSE); // Eta Shift Set, make sure that it is called only once
 201          return kTRUE;
 202       }
 203       else{
 204          printf(" Gamma Conversion Reader %s :: Eta Shift Manually Set to %f \n\n",
 205                 (fConversionCuts->GetCutNumber()).Data(),fConversionCuts->GetEtaShift());
 206          fConversionCuts->DoEtaShift(kFALSE); // Eta Shift Set, make sure that it is called only once
 207       }
 208    }
 209
 210    return kTRUE;
 211 }
 212 //________________________________________________________________________
 213 void AliV0ReaderV1::UserExec(Option_t *option){
 214
 215     // Check if correctly initialized
 216     if(!fConversionGammas)Init();
 217
 218     // User Exec
 219     fEventIsSelected=ProcessEvent(fInputEvent,fMCEvent);
 220 }
 221
 222 //________________________________________________________________________
 223 Bool_t AliV0ReaderV1::ProcessEvent(AliVEvent *inputEvent,AliMCEvent *mcEvent)
 224 {
 225
 226     //Reset the TClonesArray
 227     fConversionGammas->Delete();
 228
 229     fInputEvent=inputEvent;
 230     fMCEvent=mcEvent;
 231
 232     if(!fInputEvent){
 233         AliError("No Input event");
 234         return kFALSE;
 235     }
 236
 237     if(!fConversionCuts){AliError("No ConversionCuts");return kFALSE;}
 238
 239     // Event Cuts
 240     if(!fConversionCuts->EventIsSelected(fInputEvent,fMCEvent))return kFALSE;
 241
 242     // Set Magnetic Field
 243     AliKFParticle::SetField(fInputEvent->GetMagneticField());
 244
 245     if(fInputEvent->IsA()==AliESDEvent::Class()){
 246         ProcessESDV0s();
 247     }
 248     if(fInputEvent->IsA()==AliAODEvent::Class()){
 249         GetAODConversionGammas();
 250     }
 251
 252     // AOD Output
 253     FillAODOutput();
 254
 255     return kTRUE;
 256 }
 257 ///________________________________________________________________________
 258 void AliV0ReaderV1::FillAODOutput()
 259 {
 260     // Fill AOD Output with reconstructed Photons
 261
 262     if(fInputEvent->IsA()==AliESDEvent::Class()){
 263         ///Make sure delta aod is filled if standard aod is filled (for synchronization when reading aod with standard aod)
 264         if(fCreateAOD) {
 265             AliAODHandler * aodhandler = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 266             if (aodhandler && aodhandler->GetFillAOD()) {
 267                AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillExtension(kTRUE);
 268               //PostData(0, fConversionGammas);
 269
 270             }
 271         }
 272     }
 273 }
 274
 275 ///________________________________________________________________________
 276 const AliExternalTrackParam *AliV0ReaderV1::GetExternalTrackParam(AliESDv0 *fCurrentV0,Int_t &tracklabel,Int_t charge){
 277
 278     // Get External Track Parameter with given charge
 279
 280     if(!(charge==1||charge==-1)){AliError("Charge not defined");return 0x0;}
 281
 282     // Check for sign flip
 283     if(fCurrentV0){
 284         if(!fCurrentV0->GetParamN()||!fCurrentV0->GetParamP())return 0x0;
 285         if(!fConversionCuts->GetTrack(fInputEvent,fCurrentV0->GetNindex())||!fConversionCuts->GetTrack(fInputEvent,fCurrentV0->GetPindex()))return 0x0;
 286         if((fConversionCuts->GetTrack(fInputEvent,fCurrentV0->GetPindex()))->Charge()==charge){
 287             tracklabel=fCurrentV0->GetPindex();
 288             return fCurrentV0->GetParamP();}
 289         if((fConversionCuts->GetTrack(fInputEvent,fCurrentV0->GetNindex()))->Charge()==charge){
 290             tracklabel=fCurrentV0->GetNindex();
 291             return fCurrentV0->GetParamN();}
 292     }
 293     return 0x0;
 294 }
 295
 296 ///________________________________________________________________________
 297 Bool_t AliV0ReaderV1::ProcessESDV0s()
 298 {
 299     // Process ESD V0s for conversion photon reconstruction
 300     AliESDEvent *fESDEvent=dynamic_cast<AliESDEvent*>(fInputEvent);
 301
 302     AliKFConversionPhoton *fCurrentMotherKFCandidate=NULL;
 303
 304     if(fESDEvent){
 305
 306         for(Int_t currentV0Index=0;currentV0Index<fESDEvent->GetNumberOfV0s();currentV0Index++){
 307             AliESDv0 *fCurrentV0=(AliESDv0*)(fESDEvent->GetV0(currentV0Index));
 308             if(!fCurrentV0){
 309                 printf("Requested V0 does not exist");
 310                 continue;
 311             }
 312
 313             fCurrentMotherKFCandidate=ReconstructV0(fCurrentV0,currentV0Index);
 314
 315             if(fCurrentMotherKFCandidate){
 316
 317                //Bool_t aodV0Found = kFALSE;
 318                if(AODEvent() && fCheckAODConsistenty){
 319                   // for(Int_t i = 0; i<AODEvent()->GetNumberOfV0s();i++){
 320                   //    AliAODv0 *currebtAODV0 = AODEvent()->GetV0(i);
 321                   //    cout<<currebtAODV0->GetID()<<endl;
 322                   //    if(currebtAODV0->GetID() == currentV0Index)
 323                   //       aodV0Found = kTRUE;
 324                   // }
 325                   // if(!aodV0Found)
 326                   //    AliError(Form("AODV0 not Found belonging to ESDV0 %i",currentV0Index));
 327
 328
 329                   if(!(fConversionCuts->GetTrack(AODEvent(),fCurrentMotherKFCandidate->GetTrackLabelPositive())) || !(fConversionCuts->GetTrack(AODEvent(),fCurrentMotherKFCandidate->GetTrackLabelNegative()))){
 330                      fConversionGammas->Delete(); // Reset Gamma Array
 331                      AliError(Form("AOD Tracks not found for current ESD V0!!! V0 index %i, posESDtrack %i negESDtrack %i, Run Number: %i, Period Number: %i, NTracks: ESD %i AOD %i",
 332                                    currentV0Index,
 333                                    fCurrentMotherKFCandidate->GetTrackLabelPositive(),fCurrentMotherKFCandidate->GetTrackLabelNegative(),
 334                                    fInputEvent->GetRunNumber(),fInputEvent->GetPeriodNumber(),
 335                                    fInputEvent->GetNumberOfTracks(),AODEvent()->GetNumberOfTracks()));
 336                      return kTRUE;
 337                   }
 338                }
 339                 // Add Gamma to the TClonesArray
 340
 341                 if(kUseAODConversionPhoton){
 342                     new((*fConversionGammas)[fConversionGammas->GetEntriesFast()]) AliAODConversionPhoton(fCurrentMotherKFCandidate);
 343                                 AliAODConversionPhoton * currentConversionPhoton = (AliAODConversionPhoton*)(fConversionGammas->At(fConversionGammas->GetEntriesFast()-1));
 344                                 currentConversionPhoton->SetMass(fCurrentMotherKFCandidate->M());
 345                                 currentConversionPhoton->SetMassToZero();
 346
 347                 }
 348                 else{
 349                     new((*fConversionGammas)[fConversionGammas->GetEntriesFast()]) AliKFConversionPhoton(*fCurrentMotherKFCandidate);
 350                 }
 351
 352                 delete fCurrentMotherKFCandidate;
 353                 fCurrentMotherKFCandidate=NULL;
 354             }
 355         }
 356     }
 357     return kTRUE;
 358 }
 359
 360 ///________________________________________________________________________
 361 AliKFConversionPhoton *AliV0ReaderV1::ReconstructV0(AliESDv0 *fCurrentV0,Int_t currentV0Index)
 362 {
 363         // Reconstruct conversion photon from ESD v0
 364     fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kPhotonIn);
 365
 366     //checks if on the fly mode is set
 367     if(!fConversionCuts->SelectV0Finder(fCurrentV0->GetOnFlyStatus())){
 368        fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kOnFly);
 369        return 0x0;
 370     }
 371
 372     // TrackLabels
 373     Int_t currentTrackLabels[2]={-1,-1};
 374
 375     // Get Daughter KF Particles
 376
 377     const AliExternalTrackParam *fCurrentExternalTrackParamPositive=GetExternalTrackParamP(fCurrentV0,currentTrackLabels[0]);
 378     const AliExternalTrackParam *fCurrentExternalTrackParamNegative=GetExternalTrackParamN(fCurrentV0,currentTrackLabels[1]);
 379
 380     if(!fCurrentExternalTrackParamPositive||!fCurrentExternalTrackParamNegative)return 0x0;
 381
 382     // Apply some Cuts before Reconstruction
 383
 384     AliVTrack * posTrack = fConversionCuts->GetTrack(fInputEvent,currentTrackLabels[0]);
 385     AliVTrack * negTrack = fConversionCuts->GetTrack(fInputEvent,currentTrackLabels[1]);
 386
 387     if(!negTrack || !posTrack) {
 388        fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kNoTracks);
 389        return 0x0;
 390     }
 391
 392     // Track Cuts
 393     if(!fConversionCuts->TracksAreSelected(negTrack, posTrack)){
 394        fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kTrackCuts);
 395        return 0x0;
 396     }
 397
 398     // PID Cuts
 399     if(!fConversionCuts->dEdxCuts(negTrack) || !fConversionCuts->dEdxCuts(posTrack)) {
 400        fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kdEdxCuts);
 401        return 0x0;
 402     }
 403
 404     // Reconstruct Photon
 405
 406     AliKFConversionPhoton *fCurrentMotherKF=NULL;
 407
 408     AliKFParticle fCurrentNegativeKFParticle(*(fCurrentExternalTrackParamNegative),11);
 409     AliKFParticle fCurrentPositiveKFParticle(*(fCurrentExternalTrackParamPositive),-11);
 410
 411     // Reconstruct Gamma
 412
 413     if(fUseConstructGamma){
 414                         fCurrentMotherKF = new AliKFConversionPhoton();
 415                         fCurrentMotherKF->ConstructGamma(fCurrentNegativeKFParticle,fCurrentPositiveKFParticle);
 416     }else{
 417                         fCurrentMotherKF = new AliKFConversionPhoton(fCurrentNegativeKFParticle,fCurrentPositiveKFParticle);
 418                         fCurrentMotherKF->SetMassConstraint(0,0.0001);
 419     }
 420
 421     // Set Track Labels
 422
 423     fCurrentMotherKF->SetTrackLabels(currentTrackLabels[0],currentTrackLabels[1]);
 424
 425     // Set V0 index
 426
 427     fCurrentMotherKF->SetV0Index(currentV0Index);
 428
 429     //Set MC Label
 430
 431     if(fMCEvent){
 432
 433         AliStack *fMCStack= fMCEvent->Stack();
 434
 435         Int_t labelp=TMath::Abs(fConversionCuts->GetTrack(fInputEvent,fCurrentMotherKF->GetTrackLabelPositive())->GetLabel());
 436         Int_t labeln=TMath::Abs(fConversionCuts->GetTrack(fInputEvent,fCurrentMotherKF->GetTrackLabelNegative())->GetLabel());
 437
 438         TParticle *fNegativeMCParticle = fMCStack->Particle(labeln);
 439         TParticle *fPositiveMCParticle = fMCStack->Particle(labelp);
 440
 441         if(fPositiveMCParticle&&fNegativeMCParticle){
 442             fCurrentMotherKF->SetMCLabelPositive(labelp);
 443             fCurrentMotherKF->SetMCLabelNegative(labeln);
 444         }
 445     }
 446
 447     // Update Vertex (moved for same eta compared to old)
 448     if(fUseImprovedVertex == kTRUE){
 449        AliKFVertex primaryVertexImproved(*fInputEvent->GetPrimaryVertex());
 450        primaryVertexImproved+=*fCurrentMotherKF;
 451        fCurrentMotherKF->SetProductionVertex(primaryVertexImproved);
 452     }
 453
 454     // SetPsiPair
 455
 456     Double_t PsiPair=GetPsiPair(fCurrentV0,fCurrentExternalTrackParamPositive,fCurrentExternalTrackParamNegative);
 457     fCurrentMotherKF->SetPsiPair(PsiPair);
 458
 459
 460     // Recalculate ConversionPoint
 461     Double_t dca[2]={0,0};
 462     if(fUseOwnXYZCalculation){
 463         Double_t convpos[3]={0,0,0};
 464         if(!GetConversionPoint(fCurrentExternalTrackParamPositive,fCurrentExternalTrackParamNegative,convpos,dca)){
 465            fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kConvPointFail);
 466            delete fCurrentMotherKF;
 467            fCurrentMotherKF=NULL;
 468            return 0x0;
 469         }
 470
 471         fCurrentMotherKF->SetConversionPoint(convpos);
 472     }
 473
 474     if(fCurrentMotherKF->GetNDF() > 0.)
 475         fCurrentMotherKF->SetChi2perNDF(fCurrentMotherKF->GetChi2()/fCurrentMotherKF->GetNDF());   //->Photon is created before all chi2 relevant changes are performed, set it "by hand"
 476
 477
 478     // Set Dilepton Mass (moved down for same eta compared to old)
 479     fCurrentMotherKF->SetMass(fCurrentMotherKF->M());
 480
 481     // Apply Photon Cuts
 482
 483     if(!fConversionCuts->PhotonCuts(fCurrentMotherKF,fInputEvent)){
 484         fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kPhotonCuts);
 485         delete fCurrentMotherKF;
 486         fCurrentMotherKF=NULL;
 487         return 0x0;
 488     }
 489
 490     fConversionCuts->FillPhotonCutIndex(AliConversionCuts::kPhotonOut);
 491     return fCurrentMotherKF;
 492 }
 493
 494 ///________________________________________________________________________
 495 Double_t AliV0ReaderV1::GetPsiPair(const AliESDv0* v0, const AliExternalTrackParam *positiveparam,const AliExternalTrackParam *negativeparam) const {
 496     //
 497     // Angle between daughter momentum plane and plane
 498     //
 499
 500    AliExternalTrackParam nt(*negativeparam);
 501    AliExternalTrackParam pt(*positiveparam);
 502
 503    Float_t magField = fInputEvent->GetMagneticField();
 504
 505    Double_t xyz[3] = {0.,0.,0.};
 506    v0->GetXYZ(xyz[0],xyz[1],xyz[2]);
 507
 508    // Double_t pPlus[3]  = {pt.Px(),pt.Py(),pt.Pz()};
 509    // Double_t pMinus[3] = {nt.Px(),nt.Py(),nt.Pz()};
 510
 511    // Double_t u[3] = {pPlus[0]+pMinus[0],pPlus[1]+pMinus[1],pPlus[2]+pMinus[2]};
 512    // Double_t normu = sqrt( (u[0]*u[0]) + (u[1]*u[1]) + (u[2]*u[2]) );
 513
 514    // u[0] = u[0] / normu;
 515    // u[1] = u[1] / normu;
 516    // u[2] = u[2] / normu;
 517
 518    // Double_t normpPlus = sqrt( (pPlus[0]*pPlus[0]) + (pPlus[1]*pPlus[1]) + (pPlus[2]*pPlus[2]) );
 519    // Double_t normpMinus = sqrt( (pMinus[0]*pMinus[0]) + (pMinus[1]*pMinus[1]) + (pMinus[2]*pMinus[2]) );
 520
 521    // pPlus[0] = pPlus[0] / normpPlus;
 522    // pPlus[1] = pPlus[1] / normpPlus;
 523    // pPlus[2] = pPlus[2] / normpPlus;
 524
 525    // pMinus[0] = pMinus[0] / normpMinus;
 526    // pMinus[1] = pMinus[1] / normpMinus;
 527    // pMinus[2] = pMinus[2] / normpMinus;
 528
 529    // Double_t v[3] = {0,0,0}; // pPlus X pMinus
 530    // v[0] = (pPlus[1]*pMinus[2]) - (pPlus[2]*pMinus[1]);
 531    // v[1] = (pPlus[2]*pMinus[0]) - (pPlus[0]*pMinus[2]);
 532    // v[2] = (pPlus[0]*pMinus[1]) - (pPlus[1]*pMinus[0]);
 533
 534    // Double_t w[3] = {0,0,0}; // u X v
 535    // w[0] = (u[1]*v[2]) - (u[2]*v[1]);
 536    // w[1] = (u[2]*v[0]) - (u[0]*v[2]);
 537    // w[2] = (u[0]*v[1]) - (u[1]*v[0]);
 538
 539    // Double_t z[3] = {0,0,1};
 540    // Double_t wc[3] = {0,0,0}; // u X v
 541    // wc[0] = (u[1]*z[2]) - (u[2]*z[1]);
 542    // wc[1] = (u[2]*z[0]) - (u[0]*z[2]);
 543    // wc[2] = (u[0]*z[1]) - (u[1]*z[0]);
 544
 545    // Double_t PhiV = TMath::ACos((w[0]*wc[0]) + (w[1]*wc[1]) + (w[2]*wc[2]));
 546    //return abs(PhiV);
 547
 548
 549    // TVector3 pPlus(pt.Px(),pt.Py(),pt.Pz());
 550    // TVector3 pMinus(nt.Px(),nt.Py(),nt.Pz());
 551
 552    // TVector3 u = pMinus + pPlus;
 553    // u = u*(1/u.Mag());
 554
 555    // TVector3 pHPlus = pPlus*(1/pPlus.Mag());
 556    // TVector3 pHMinus = pMinus*(1/pMinus.Mag());
 557
 558    // TVector3 v = pHPlus.Cross(pHMinus);
 559    // TVector3 w = u.Cross(v);
 560    // TVector3 z(0,0,1);
 561    // TVector3 wc = u.Cross(z);
 562
 563    // Double_t PhiV = w * wc;
 564
 565    Double_t mn[3] = {0,0,0};
 566    Double_t mp[3] = {0,0,0};
 567
 568    v0->GetNPxPyPz(mn[0],mn[1],mn[2]);//reconstructed cartesian momentum components of negative daughter;
 569    v0->GetPPxPyPz(mp[0],mp[1],mp[2]);//reconstructed cartesian momentum components of positive daughter;
 570
 571    Double_t deltat = 1.;
 572    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
 573    Double_t radiussum = TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]) + 50;//radius to which tracks shall be propagated
 574
 575    Double_t momPosProp[3] = {0,0,0};
 576    Double_t momNegProp[3] = {0,0,0};
 577
 578    Double_t psiPair = 4.;
 579    if(nt.PropagateTo(radiussum,magField) == 0) return psiPair; //propagate tracks to the outside -> Better Purity and Efficiency
 580    if(pt.PropagateTo(radiussum,magField) == 0) return psiPair; //propagate tracks to the outside -> Better Purity and Efficiency
 581
 582    pt.GetPxPyPz(momPosProp);//Get momentum vectors of tracks after propagation
 583    nt.GetPxPyPz(momNegProp);
 584
 585    Double_t pEle =
 586        TMath::Sqrt(momNegProp[0]*momNegProp[0]+momNegProp[1]*momNegProp[1]+momNegProp[2]*momNegProp[2]);//absolute momentum value of negative daughter
 587
 588    Double_t pPos =
 589        TMath::Sqrt(momPosProp[0]*momPosProp[0]+momPosProp[1]*momPosProp[1]+momPosProp[2]*momPosProp[2]);//absolute momentum value of positive daughter
 590
 591    Double_t scalarproduct =
 592        momPosProp[0]*momNegProp[0]+momPosProp[1]*momNegProp[1]+momPosProp[2]*momNegProp[2];//scalar product of propagated positive and negative daughters' momenta
 593
 594    Double_t chipair = TMath::ACos(scalarproduct/(pEle*pPos));//Angle between propagated daughter tracks
 595
 596    psiPair =  TMath::Abs(TMath::ASin(deltat/chipair));
 597
 598    return psiPair;
 599 }
 600
 601 ///________________________________________________________________________
 602 Bool_t AliV0ReaderV1::GetHelixCenter(const AliExternalTrackParam *track,Double_t center[2]){
 603
 604     // Get Center of the helix track parametrization
 605
 606     Int_t charge=track->Charge();
 607     Double_t b=fInputEvent->GetMagneticField();
 608
 609     Double_t    helix[6];
 610     track->GetHelixParameters(helix,b);
 611
 612     Double_t xpos =     helix[5];
 613     Double_t ypos =     helix[0];
 614     Double_t radius = TMath::Abs(1./helix[4]);
 615     Double_t phi = helix[2];
 616
 617     if(phi < 0){
 618         phi = phi + 2*TMath::Pi();
 619     }
 620
 621     phi -= TMath::Pi()/2.;
 622     Double_t xpoint =   radius * TMath::Cos(phi);
 623     Double_t ypoint =   radius * TMath::Sin(phi);
 624
 625     if(b<0){
 626         if(charge > 0){
 627             xpoint = - xpoint;
 628             ypoint = - ypoint;
 629         }
 630
 631         if(charge < 0){
 632             xpoint =    xpoint;
 633             ypoint =    ypoint;
 634         }
 635     }
 636     if(b>0){
 637         if(charge > 0){
 638             xpoint =    xpoint;
 639             ypoint =    ypoint;
 640         }
 641
 642         if(charge < 0){
 643             xpoint = - xpoint;
 644             ypoint = - ypoint;
 645         }
 646     }
 647     center[0] = xpos + xpoint;
 648     center[1] = ypos + ypoint;
 649
 650     return 1;
 651 }
 652 ///________________________________________________________________________
 653 Bool_t AliV0ReaderV1::GetConversionPoint(const AliExternalTrackParam *pparam,const AliExternalTrackParam *nparam,Double_t convpos[3],Double_t dca[2]){
 654
 655     // Recalculate Conversion Point
 656
 657     if(!pparam||!nparam)return kFALSE;
 658
 659     Double_t helixcenterpos[2];
 660     GetHelixCenter(pparam,helixcenterpos);
 661
 662     Double_t helixcenterneg[2];
 663     GetHelixCenter(nparam,helixcenterneg);
 664
 665         Double_t helixpos[6];
 666         pparam->GetHelixParameters(helixpos,fInputEvent->GetMagneticField());
 667         Double_t posradius = TMath::Abs(1./helixpos[4]);
 668
 669         Double_t helixneg[6];
 670         nparam->GetHelixParameters(helixneg,fInputEvent->GetMagneticField());
 671         Double_t negradius = TMath::Abs(1./helixneg[4]);
 672
 673         // Calculate xy-position
 674
 675         Double_t xpos = helixcenterpos[0];
 676         Double_t ypos = helixcenterpos[1];
 677         Double_t xneg = helixcenterneg[0];
 678         Double_t yneg = helixcenterneg[1];
 679
 680         convpos[0] = (xpos*negradius + xneg*posradius)/(negradius+posradius);
 681         convpos[1] = (ypos*negradius+ yneg*posradius)/(negradius+posradius);
 682
 683
 684         // Calculate Track XY vertex position
 685
 686         Double_t deltaXPos = convpos[0] -       xpos;
 687         Double_t deltaYPos = convpos[1] -       ypos;
 688
 689         Double_t deltaXNeg = convpos[0] -       xneg;
 690         Double_t deltaYNeg = convpos[1] -       yneg;
 691
 692         Double_t alphaPos =     TMath::Pi() + TMath::ATan2(-deltaYPos,-deltaXPos);
 693         Double_t alphaNeg =     TMath::Pi() + TMath::ATan2(-deltaYNeg,-deltaXNeg);
 694
 695         Double_t vertexXNeg =   xneg +  TMath::Abs(negradius)*TMath::Cos(alphaNeg);
 696         Double_t vertexYNeg =   yneg +  TMath::Abs(negradius)*TMath::Sin(alphaNeg);
 697
 698         Double_t vertexXPos =   xpos +  TMath::Abs(posradius)*TMath::Cos(alphaPos);
 699         Double_t vertexYPos =   ypos +  TMath::Abs(posradius)*TMath::Sin(alphaPos);
 700
 701         AliExternalTrackParam p(*pparam);
 702         AliExternalTrackParam n(*nparam);
 703
 704         TVector2 vertexPos(vertexXPos,vertexYPos);
 705         TVector2 vertexNeg(vertexXNeg,vertexYNeg);
 706
 707         // Convert to local coordinate system
 708         TVector2 vertexPosRot=vertexPos.Rotate(-p.GetAlpha());
 709         TVector2 vertexNegRot=vertexNeg.Rotate(-n.GetAlpha());
 710
 711         // Propagate Track Params to Vertex
 712
 713         if(!p.PropagateTo(vertexPosRot.X(),fInputEvent->GetMagneticField()))return kFALSE;
 714         if(!n.PropagateTo(vertexNegRot.X(),fInputEvent->GetMagneticField()))return kFALSE;
 715
 716         // Check whether propagation was sucessful
 717
 718         if(TMath::Abs(vertexPos.Mod()-TMath::Sqrt(p.GetX()*p.GetX()+p.GetY()*p.GetY()))>0.01)return kFALSE;
 719         if(TMath::Abs(vertexNeg.Mod()-TMath::Sqrt(n.GetX()*n.GetX()+n.GetY()*n.GetY()))>0.01)return kFALSE;
 720
 721         // Calculate z position
 722
 723         convpos[2] = (p.GetZ()*negradius+n.GetZ()*posradius)/(negradius+posradius);
 724
 725         // Calculate DCA
 726         TVector2 vdca=vertexPos-vertexNeg;
 727         dca[0]=vdca.Mod();
 728         dca[1]=TMath::Abs(n.GetZ()-p.GetZ());
 729
 730         return kTRUE;
 731 }
 732 //________________________________________________________________________
 733 Bool_t AliV0ReaderV1::GetAODConversionGammas(){
 734
 735     // Get reconstructed conversion photons from satellite AOD file
 736
 737     AliAODEvent *fAODEvent=dynamic_cast<AliAODEvent*>(fInputEvent);
 738
 739     if(fAODEvent){
 740
 741        if(fConversionGammas == NULL){
 742           fConversionGammas = new TClonesArray("AliAODConversionPhoton",100);
 743        }
 744         fConversionGammas->Delete();//Reset the TClonesArray
 745
 746         //Get Gammas from satellite AOD gamma branch
 747
 748         AliAODConversionPhoton *gamma=0x0;
 749
 750         TClonesArray *fInputGammas=dynamic_cast<TClonesArray*>(fAODEvent->FindListObject(fDeltaAODBranchName.Data()));
 751         if(!fInputGammas){
 752            FindDeltaAODBranchName();
 753            fInputGammas=dynamic_cast<TClonesArray*>(fAODEvent->FindListObject(fDeltaAODBranchName.Data()));}
 754         if(!fInputGammas){AliError("No Gamma Satellites found");return kFALSE;}
 755         // Apply Selection Cuts to Gammas and create local working copy
 756         if(fInputGammas){
 757            for(Int_t i=0;i<fInputGammas->GetEntriesFast();i++){
 758               gamma=dynamic_cast<AliAODConversionPhoton*>(fInputGammas->At(i));
 759               if(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::Terminate(Option_t *)
 789 {
 790
 791 }