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