AliAnalysisTaskSE(),\r
fV0Reader(NULL),\r
fStack(NULL),\r
+ fESDEvent(NULL), \r
fOutputContainer(NULL),\r
fHistograms(NULL),\r
fDoMCTruth(kFALSE),\r
fIsTrueReconstructedGammas(),\r
fElectronv1(),\r
fElectronv2(),\r
+ fCurrentEventPosElectron(),\r
+ fPreviousEventPosElectron(),\r
+ fCurrentEventNegElectron(),\r
+ fPreviousEventNegElectron(),\r
+ fKFReconstructedGammasCut(), \r
+ fPreviousEventTLVNegElectron(),\r
+ fPreviousEventTLVPosElectron(), \r
fElectronMass(-1),\r
fGammaMass(-1),\r
fPi0Mass(-1),\r
AliAnalysisTaskSE(name),\r
fV0Reader(NULL),\r
fStack(NULL),\r
+ fESDEvent(NULL), \r
fOutputContainer(0x0),\r
fHistograms(NULL),\r
fDoMCTruth(kFALSE),\r
fIsTrueReconstructedGammas(),\r
fElectronv1(),\r
fElectronv2(),\r
+ fCurrentEventPosElectron(),\r
+ fPreviousEventPosElectron(),\r
+ fCurrentEventNegElectron(),\r
+ fPreviousEventNegElectron(),\r
+ fKFReconstructedGammasCut(), \r
+ fPreviousEventTLVNegElectron(),\r
+ fPreviousEventTLVPosElectron(),\r
fElectronMass(-1),\r
fGammaMass(-1),\r
fPi0Mass(-1),\r
fIsTrueReconstructedGammas.clear();\r
fElectronv1.clear();\r
fElectronv2.clear();\r
+ fCurrentEventPosElectron.clear();\r
+ fCurrentEventNegElectron.clear(); \r
+ fKFReconstructedGammasCut.clear(); \r
\r
//Clear the data in the v0Reader\r
fV0Reader->UpdateEventByEventData();\r
\r
// Process reconstructed gammas\r
ProcessGammasForNeutralMesonAnalysis();\r
+\r
+ CheckV0Efficiency();\r
\r
+ //Process reconstructed gammas electrons for Chi_c Analysis\r
+ ProcessGammaElectronsForChicAnalysis();\r
+\r
PostData(1, fOutputContainer);\r
\r
}\r
continue;\r
}\r
\r
+ ///////////////////////Begin Chic Analysis/////////////////////////////\r
+\r
+\r
+ if(particle->GetPdgCode() == 443){//Is JPsi\r
+\r
+ if(particle->GetNDaughters()==2){\r
+ if(TMath::Abs(fStack->Particle(particle->GetFirstDaughter())->GetPdgCode()) == 11 &&\r
+ TMath::Abs(fStack->Particle(particle->GetLastDaughter())->GetPdgCode()) == 11){\r
+ TParticle* daug0 = fStack->Particle(particle->GetFirstDaughter());\r
+ TParticle* daug1 = fStack->Particle(particle->GetLastDaughter());\r
+ if(TMath::Abs(daug0->Eta()) < 0.9 && TMath::Abs(daug1->Eta()) < 0.9)\r
+ fHistograms->FillTable("Table_Electrons",3);//e+ e- from J/Psi inside acceptance\r
+\r
+ if( TMath::Abs(daug0->Eta()) < 0.9){\r
+ if(daug0->GetPdgCode() == -11)\r
+ fHistograms->FillTable("Table_Electrons",1);//e+ from J/Psi inside acceptance\r
+ else\r
+ fHistograms->FillTable("Table_Electrons",2);//e- from J/Psi inside acceptance\r
+\r
+ }\r
+ if(TMath::Abs(daug1->Eta()) < 0.9){\r
+ if(daug1->GetPdgCode() == -11)\r
+ fHistograms->FillTable("Table_Electrons",1);//e+ from J/Psi inside acceptance\r
+ else\r
+ fHistograms->FillTable("Table_Electrons",2);//e- from J/Psi inside acceptance\r
+ }\r
+ }\r
+ }\r
+ }\r
+ // const int CHI_C0 = 10441;\r
+ // const int CHI_C1 = 20443;\r
+ // const int CHI_C2 = 445\r
+ if(particle->GetPdgCode() == 22){//gamma from JPsi\r
+ if(particle->GetMother(0) > -1){\r
+ if(fStack->Particle(particle->GetMother(0))->GetPdgCode() == 10441 ||\r
+ fStack->Particle(particle->GetMother(0))->GetPdgCode() == 20443 ||\r
+ fStack->Particle(particle->GetMother(0))->GetPdgCode() == 445){\r
+ if(TMath::Abs(particle->Eta()) < 1.2)\r
+ fHistograms->FillTable("Table_Electrons",17);// gamma from chic inside accptance\r
+ }\r
+ }\r
+ }\r
+ if(particle->GetPdgCode() == 10441 || particle->GetPdgCode() == 20443 || particle->GetPdgCode() == 445){\r
+ if( particle->GetNDaughters() == 2){\r
+ TParticle* daug0 = fStack->Particle(particle->GetFirstDaughter());\r
+ TParticle* daug1 = fStack->Particle(particle->GetLastDaughter());\r
+\r
+ if( (daug0->GetPdgCode() == 443 || daug0->GetPdgCode() == 22) && (daug1->GetPdgCode() == 443 || daug1->GetPdgCode() == 22) ){\r
+ if( daug0->GetPdgCode() == 443){\r
+ TParticle* daugE0 = fStack->Particle(daug0->GetFirstDaughter());\r
+ TParticle* daugE1 = fStack->Particle(daug0->GetLastDaughter());\r
+ if( TMath::Abs(daug1->Eta()) < 1.2 && TMath::Abs(daugE0->Eta()) < 0.9 && TMath::Abs(daugE1->Eta()) < 0.9 )\r
+ fHistograms->FillTable("Table_Electrons",18);\r
+\r
+ }//if\r
+ else if (daug1->GetPdgCode() == 443){\r
+ TParticle* daugE0 = fStack->Particle(daug1->GetFirstDaughter());\r
+ TParticle* daugE1 = fStack->Particle(daug1->GetLastDaughter());\r
+ if( TMath::Abs(daug0->Eta()) < 1.2 && TMath::Abs(daugE0->Eta()) < 0.9 && TMath::Abs(daugE1->Eta()) < 0.9 )\r
+ fHistograms->FillTable("Table_Electrons",18);\r
+ }//else if\r
+ }//gamma o Jpsi\r
+ }//GetNDaughters\r
+ }\r
+\r
+\r
+ /////////////////////End Chic Analysis////////////////////////////\r
+\r
+\r
if(TMath::Abs(particle->Eta())> fV0Reader->GetEtaCut() ) continue;\r
\r
if(particle->R()>fV0Reader->GetMaxRCut()) continue; // cuts on distance from collision point\r
fHistograms->FillHistogram("ESD_Mother_ZR", twoGammaCandidate->GetZ(), spaceVectorTwoGammaCandidate.Pt());\r
fHistograms->FillHistogram("ESD_Mother_XY", twoGammaCandidate->GetX(), twoGammaCandidate->GetY());\r
fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());\r
+ fHistograms->FillHistogram("ESD_Mother_InvMass",massTwoGammaCandidate);\r
}\r
}\r
delete twoGammaCandidate;\r
fHistograms->FillHistogram("ESD_Background_ZR", backgroundCandidate->GetZ(), SpaceVectorbackgroundCandidate.Pt());\r
fHistograms->FillHistogram("ESD_Background_XY", backgroundCandidate->GetX(), backgroundCandidate->GetY());\r
fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt",massBG,MomentumVectorbackgroundCandidate.Pt());\r
+ fHistograms->FillHistogram("ESD_Background_InvMass",massBG);\r
}\r
}\r
delete backgroundCandidate; \r
TVector3 v3D1(daughter1->Px(),daughter1->Py(),daughter1->Pz());\r
return v3D0.Angle(v3D1);\r
}\r
+\r
+void AliAnalysisTaskGammaConversion::CheckV0Efficiency(){\r
+\r
+ vector<Int_t> indexOfGammaParticle;\r
+\r
+ fStack = fV0Reader->GetMCStack();\r
+\r
+ if(fV0Reader->CheckForPrimaryVertex() == kFALSE){\r
+ return; // aborts if the primary vertex does not have contributors.\r
+ }\r
+\r
+ for (Int_t iTracks = 0; iTracks < fStack->GetNprimary(); iTracks++) {\r
+ TParticle* particle = (TParticle *)fStack->Particle(iTracks);\r
+ if(particle->GetPdgCode()==22){ //Gamma\r
+ if(particle->GetNDaughters() >= 2){\r
+ TParticle* electron=NULL;\r
+ TParticle* positron=NULL; \r
+ for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){\r
+ TParticle *tmpDaughter = fStack->Particle(daughterIndex);\r
+ if(tmpDaughter->GetUniqueID() == 5){\r
+ if(tmpDaughter->GetPdgCode() == 11){\r
+ electron = tmpDaughter;\r
+ }\r
+ else if(tmpDaughter->GetPdgCode() == -11){\r
+ positron = tmpDaughter;\r
+ }\r
+ }\r
+ }\r
+ if(electron!=NULL && positron!=0){\r
+ if(electron->R()<160){\r
+ indexOfGammaParticle.push_back(iTracks);\r
+ }\r
+ }\r
+ }\r
+ }\r
+ }\r
+\r
+ Int_t nFoundGammas=0;\r
+ Int_t nNotFoundGammas=0;\r
+\r
+ Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();\r
+ for(Int_t i=0;i<numberOfV0s;i++){\r
+ fV0Reader->GetV0(i);\r
+ \r
+ if(fV0Reader->HasSameMCMother() == kFALSE){\r
+ continue;\r
+ }\r
+ \r
+ TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();\r
+ TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();\r
+ \r
+ if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){\r
+ continue;\r
+ }\r
+ if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){\r
+ continue;\r
+ }\r
+ \r
+ if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){\r
+ //TParticle * v0Gamma = fV0Reader->GetMotherMCParticle();\r
+ for(UInt_t mcIndex=0;mcIndex<indexOfGammaParticle.size();mcIndex++){\r
+ if(negativeMC->GetFirstMother()==indexOfGammaParticle[mcIndex]){\r
+ nFoundGammas++;\r
+ }\r
+ else{\r
+ nNotFoundGammas++;\r
+ }\r
+ }\r
+ }\r
+ }\r
+ // cout<<"Found: "<<nFoundGammas<<" of: "<<indexOfGammaParticle.size()<<endl;\r
+}\r
+\r
+\r
+void AliAnalysisTaskGammaConversion::ProcessGammaElectronsForChicAnalysis(){\r
+\r
+\r
+ fESDEvent = fV0Reader->GetESDEvent();\r
+\r
+\r
+ vector <AliESDtrack*> ESDeNegTemp(0);\r
+ vector <AliESDtrack*> ESDePosTemp(0);\r
+ vector <AliESDtrack*> ESDxNegTemp(0);\r
+ vector <AliESDtrack*> ESDxPosTemp(0);\r
+ vector <AliESDtrack*> ESDeNegNoJPsi(0);\r
+ vector <AliESDtrack*> ESDePosNoJPsi(0); \r
+\r
+\r
+\r
+ fHistograms->FillTable("Table_Electrons",0);//Count number of Events\r
+\r
+ for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){\r
+ AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);\r
+\r
+ if(!curTrack){\r
+ //print warning here\r
+ continue;\r
+ }\r
+\r
+ double p[3];if(!curTrack->GetConstrainedPxPyPz(p))continue;\r
+ double r[3];curTrack->GetConstrainedXYZ(r);\r
+\r
+ TVector3 rXYZ(r);\r
+\r
+ fHistograms->FillTable("Table_Electrons",4);//Count number of ESD tracks\r
+\r
+ Bool_t flagKink = kTRUE;\r
+ Bool_t flagTPCrefit = kTRUE;\r
+ Bool_t flagTRDrefit = kTRUE;\r
+ Bool_t flagITSrefit = kTRUE;\r
+ Bool_t flagTRDout = kTRUE;\r
+ Bool_t flagVertex = kTRUE;\r
+\r
+\r
+ //Cuts ---------------------------------------------------------------\r
+\r
+ if(curTrack->GetKinkIndex(0) > 0){\r
+ fHistograms->FillHistogram("Table_Electrons",5);//Count kink\r
+ flagKink = kFALSE;\r
+ }\r
+\r
+ ULong_t trkStatus = curTrack->GetStatus();\r
+\r
+ ULong_t tpcRefit = (trkStatus & AliESDtrack::kTPCrefit);\r
+\r
+ if(!tpcRefit){\r
+ fHistograms->FillHistogram("Table_Electrons",9);//Count not TPCrefit\r
+ flagTPCrefit = kFALSE;\r
+ }\r
+\r
+ ULong_t itsRefit = (trkStatus & AliESDtrack::kITSrefit);\r
+ if(!itsRefit){\r
+ fHistograms->FillHistogram("Table_Electrons",10);//Count not ITSrefit\r
+ flagITSrefit = kFALSE;\r
+ }\r
+\r
+ ULong_t trdRefit = (trkStatus & AliESDtrack::kTRDrefit);\r
+\r
+ if(!trdRefit){\r
+ fHistograms->FillHistogram("Table_Electrons",8); //Count not TRDrefit\r
+ flagTRDrefit = kFALSE;\r
+ }\r
+\r
+ ULong_t trdOut = (trkStatus & AliESDtrack::kTRDout);\r
+\r
+ if(!trdOut) {\r
+ fHistograms->FillHistogram("Table_Electrons",7); //Count not TRDout\r
+ flagTRDout = kFALSE;\r
+ }\r
+\r
+ double nsigmaToVxt = GetSigmaToVertex(curTrack);\r
+\r
+ if(nsigmaToVxt > 3){\r
+ fHistograms->FillHistogram("Table_Electrons",6); //Count Tracks with number of sigmas > 3\r
+ flagVertex = kFALSE;\r
+ }\r
+\r
+ if(! (flagKink && flagTPCrefit && flagITSrefit && flagTRDrefit && flagTRDout && flagVertex ) ) continue;\r
+ fHistograms->FillHistogram("Table_Electrons",11);//Count Tracks passed Cuts\r
+\r
+\r
+ Stat_t pid, weight;\r
+ GetPID(curTrack, pid, weight);\r
+\r
+ if(pid!=0){\r
+ fHistograms->FillHistogram("Table_Electrons",12); //Count Tracks with pid != 0\r
+ }\r
+\r
+ if(pid == 0){\r
+ fHistograms->FillHistogram("Table_Electrons",13); //Count Tracks with pid != 0\r
+ }\r
+\r
+\r
+\r
+\r
+ Int_t labelMC = TMath::Abs(curTrack->GetLabel());\r
+ TParticle* curParticle = fStack->Particle(labelMC);\r
+\r
+\r
+\r
+\r
+ TLorentzVector curElec;\r
+ curElec.SetXYZM(p[0],p[1],p[2],fElectronMass);\r
+\r
+\r
+\r
+\r
+ if(curTrack->GetSign() > 0){\r
+\r
+ ESDxPosTemp.push_back(curTrack);\r
+\r
+ if( pid == 0){\r
+\r
+ fHistograms->FillHistogram("ESD_ElectronPosNegPt",curElec.Pt());\r
+ fHistograms->FillHistogram("ESD_ElectronPosPt",curElec.Pt());\r
+ fHistograms->FillHistogram("MC_ElectronPosNegPt",curParticle->Pt());\r
+ fHistograms->FillHistogram("ESD_ElectronPosNegEta",curElec.Eta());\r
+ fHistograms->FillHistogram("MC_ElectronPosNegEta",curParticle->Eta());\r
+ ESDePosTemp.push_back(curTrack);\r
+\r
+\r
+\r
+ }\r
+\r
+ }\r
+ else {\r
+ ESDxNegTemp.push_back(curTrack);\r
+\r
+ if( pid == 0){\r
+\r
+ fHistograms->FillHistogram("ESD_ElectronPosNegPt",curElec.Pt());\r
+ fHistograms->FillHistogram("ESD_ElectronNegPt",curElec.Pt());\r
+ fHistograms->FillHistogram("MC_ElectronPosNegPt",curParticle->Pt());\r
+ fHistograms->FillHistogram("ESD_ElectronPosNegEta",curElec.Eta());\r
+ fHistograms->FillHistogram("MC_ElectronPosNegEta",curParticle->Eta());\r
+ ESDeNegTemp.push_back(curTrack);\r
+\r
+\r
+\r
+\r
+ }\r
+ }\r
+\r
+ }\r
+\r
+\r
+ Bool_t ePosJPsi = kFALSE;\r
+ Bool_t eNegJPsi = kFALSE; \r
+ Bool_t ePosPi0 = kFALSE;\r
+ Bool_t eNegPi0 = kFALSE;\r
+ \r
+ UInt_t iePosJPsi=0,ieNegJPsi=0,iePosPi0=0,ieNegPi0=0;\r
+ \r
+ for(UInt_t iNeg=0; iNeg < ESDeNegTemp.size(); iNeg++){\r
+ if(fStack->Particle(TMath::Abs(ESDeNegTemp[iNeg]->GetLabel()))->GetPdgCode() == 11)\r
+ if(fStack->Particle(TMath::Abs(ESDeNegTemp[iNeg]->GetLabel()))->GetMother(0) > -1){\r
+ Int_t labelMother = fStack->Particle(TMath::Abs(ESDeNegTemp[iNeg]->GetLabel()))->GetMother(0);\r
+ TParticle* partMother = fStack ->Particle(labelMother);\r
+ if (partMother->GetPdgCode() == 111){\r
+ ieNegPi0 = iNeg;\r
+ eNegPi0 = kTRUE;\r
+ }\r
+ if(partMother->GetPdgCode() == 443){ //Mother JPsi\r
+ fHistograms->FillTable("Table_Electrons",14);\r
+ ieNegJPsi = iNeg;\r
+ eNegJPsi = kTRUE;\r
+ }\r
+ else{ \r
+ ESDeNegNoJPsi.push_back(ESDeNegTemp[iNeg]);\r
+ // cout<<"ESD No Positivo JPsi "<<endl;\r
+ }\r
+\r
+ }\r
+ } \r
+\r
+ for(UInt_t iPos=0; iPos < ESDePosTemp.size(); iPos++){\r
+ if(fStack->Particle(TMath::Abs(ESDePosTemp[iPos]->GetLabel()))->GetPdgCode() == -11)\r
+ if(fStack->Particle(TMath::Abs(ESDePosTemp[iPos]->GetLabel()))->GetMother(0) > -1){\r
+ Int_t labelMother = fStack->Particle(TMath::Abs(ESDePosTemp[iPos]->GetLabel()))->GetMother(0);\r
+ TParticle* partMother = fStack ->Particle(labelMother);\r
+ if (partMother->GetPdgCode() == 111){\r
+ iePosPi0 = iPos;\r
+ ePosPi0 = kTRUE;\r
+ }\r
+ if(partMother->GetPdgCode() == 443){ //Mother JPsi\r
+ fHistograms->FillTable("Table_Electrons",15);\r
+ iePosJPsi = iPos;\r
+ ePosJPsi = kTRUE;\r
+ }\r
+ else{\r
+ ESDePosNoJPsi.push_back(ESDePosTemp[iPos]);\r
+ // cout<<"ESD No Negativo JPsi "<<endl;\r
+ }\r
+\r
+ }\r
+ }\r
+ \r
+ if( eNegJPsi && ePosJPsi ){\r
+ TVector3 tempeNegV,tempePosV;\r
+ tempeNegV.SetXYZ(ESDeNegTemp[ieNegJPsi]->Px(),ESDeNegTemp[ieNegJPsi]->Py(),ESDeNegTemp[ieNegJPsi]->Pz()); \r
+ tempePosV.SetXYZ(ESDePosTemp[iePosJPsi]->Px(),ESDePosTemp[iePosJPsi]->Py(),ESDePosTemp[iePosJPsi]->Pz());\r
+ fHistograms->FillTable("Table_Electrons",16);\r
+ fHistograms->FillHistogram("ESD_ElectronPosNegJPsiAngle",tempeNegV.Angle(tempePosV)); \r
+ fHistograms->FillHistogram("MC_ElectronPosNegJPsiAngle",GetMCOpeningAngle(fStack->Particle(TMath::Abs(ESDeNegTemp[ieNegJPsi]->GetLabel())),\r
+ fStack->Particle(TMath::Abs(ESDePosTemp[iePosJPsi]->GetLabel())))); \r
+ }\r
+ \r
+ if( eNegPi0 && ePosPi0 ){\r
+ TVector3 tempeNegV,tempePosV;\r
+ tempeNegV.SetXYZ(ESDeNegTemp[ieNegPi0]->Px(),ESDeNegTemp[ieNegPi0]->Py(),ESDeNegTemp[ieNegPi0]->Pz());\r
+ tempePosV.SetXYZ(ESDePosTemp[iePosPi0]->Px(),ESDePosTemp[iePosPi0]->Py(),ESDePosTemp[iePosPi0]->Pz());\r
+ fHistograms->FillHistogram("ESD_ElectronPosNegPi0Angle",tempeNegV.Angle(tempePosV));\r
+ fHistograms->FillHistogram("MC_ElectronPosNegPi0Angle",GetMCOpeningAngle(fStack->Particle(TMath::Abs(ESDeNegTemp[ieNegPi0]->GetLabel())),\r
+ fStack->Particle(TMath::Abs(ESDePosTemp[iePosPi0]->GetLabel())))); \r
+ }\r
+ \r
+\r
+ FillAngle("ESD_eNegePosAngleBeforeCut",GetTLorentzVector(ESDeNegTemp),GetTLorentzVector(ESDePosTemp));\r
+\r
+ CleanWithAngleCuts(ESDeNegTemp,ESDePosTemp,fKFReconstructedGammas);\r
+ \r
+ vector <TLorentzVector> CurrentTLVeNeg = GetTLorentzVector(fCurrentEventNegElectron);\r
+ vector <TLorentzVector> CurrentTLVePos = GetTLorentzVector(fCurrentEventPosElectron);\r
+\r
+\r
+ FillAngle("ESD_eNegePosAngleAfterCut",CurrentTLVeNeg,CurrentTLVePos);\r
+\r
+ \r
+\r
+\r
+ //FillAngle("ESD_eNegePosAngleAfterCut",CurrentTLVeNeg,CurrentTLVePos);\r
+\r
+\r
+ FillElectronInvMass("ESD_InvMass_ePluseMinus",CurrentTLVeNeg,CurrentTLVePos);\r
+ FillElectronInvMass("ESD_InvMass_xPlusxMinus",GetTLorentzVector(ESDxNegTemp),GetTLorentzVector(ESDxPosTemp));\r
+\r
+ \r
+\r
+ FillGammaElectronInvMass("ESD_InvMass_GammaePluseMinusChiC","ESD_InvMass_GammaePluseMinusChiCDiff",\r
+ fKFReconstructedGammasCut,CurrentTLVeNeg,CurrentTLVePos);\r
+\r
+ FillGammaElectronInvMass("ESD_InvMass_GammaePluseMinusPi0","ESD_InvMass_GammaePluseMinusPi0Diff",\r
+ fKFReconstructedGammasCut,CurrentTLVeNeg,CurrentTLVePos);\r
+\r
+ //BackGround\r
+\r
+ //Like Sign e+e-\r
+ ElectronBackground("ESD_ENegBackground",CurrentTLVeNeg);\r
+ ElectronBackground("ESD_EPosBackground",CurrentTLVePos);\r
+ ElectronBackground("ESD_EPosENegBackground",CurrentTLVeNeg);\r
+ ElectronBackground("ESD_EPosENegBackground",CurrentTLVePos);\r
+\r
+ // Like Sign e+e- no JPsi\r
+ ElectronBackground("ESD_EPosENegNoJPsiBG",GetTLorentzVector(ESDeNegNoJPsi));\r
+ ElectronBackground("ESD_EPosENegNoJPsiBG",GetTLorentzVector(ESDePosNoJPsi));\r
+\r
+ //Mixed Event\r
+\r
+ if( fCurrentEventPosElectron.size() > 0 && fCurrentEventNegElectron.size() > 0 && fKFReconstructedGammasCut.size() > 0 ){\r
+ FillGammaElectronInvMass("ESD_EPosENegGammaBackgroundMX","ESD_EPosENegGammaBackgroundMXDiff",\r
+ fKFReconstructedGammasCut,fPreviousEventTLVNegElectron,fPreviousEventTLVPosElectron);\r
+ fPreviousEventTLVNegElectron = CurrentTLVeNeg;\r
+ fPreviousEventTLVPosElectron = CurrentTLVePos;\r
+\r
+ }\r
+\r
+ /*\r
+ //Photons P\r
+ Double_t vtx[3];\r
+ vtx[0]=0;vtx[1]=0;vtx[2]=0;\r
+ for(UInt_t i=0;i<fKFReconstructedGammasChic.size();i++){\r
+\r
+ // if(fMCGammaChicTempCut[i]->GetMother(0) < 0) continue;\r
+\r
+\r
+\r
+ Int_t tempLabel = fStack->Particle(fMCGammaChicTempCut[i]->GetMother(0))->GetPdgCode();\r
+ // cout<<" Label Pedro Gonzalez " <<tempLabel <<endl;\r
+\r
+ // cout<<" Label Distance"<<fKFReconstructedGammasChic[i].GetDistanceFromVertex(vtx)<<endl;\r
+\r
+ if( tempLabel == 10441 || tempLabel == 20443 || tempLabel == 445 )\r
+\r
+ fHistograms->FillHistogram("ESD_PhotonsMomentum",fKFReconstructedGammasChic[i].GetMomentum());\r
+\r
+\r
+ }\r
+\r
+\r
+ */\r
+\r
+\r
+}\r
+\r
+void AliAnalysisTaskGammaConversion::FillAngle(TString histoName,vector <TLorentzVector> TLVeNeg, vector <TLorentzVector> TLVePos){\r
+ for( UInt_t iNeg=0; iNeg < TLVeNeg.size(); iNeg++){\r
+ for (UInt_t iPos=0; iPos < TLVePos.size(); iPos++){\r
+ fHistograms->FillHistogram(histoName.Data(),TLVeNeg[iNeg].Vect().Angle(TLVePos[iPos].Vect()));\r
+ }\r
+ }\r
+}\r
+void AliAnalysisTaskGammaConversion::FillElectronInvMass(TString histoName, vector <TLorentzVector> eNeg, vector <TLorentzVector> ePos){\r
+\r
+ for( UInt_t n=0; n < eNeg.size(); n++){\r
+\r
+ TLorentzVector en = eNeg.at(n);\r
+ for (UInt_t p=0; p < ePos.size(); p++){\r
+ TLorentzVector ep = ePos.at(p);\r
+ TLorentzVector np = ep + en;\r
+ fHistograms->FillHistogram(histoName.Data(),np.M());\r
+ }\r
+ }\r
+\r
+}\r
+\r
+void AliAnalysisTaskGammaConversion::FillGammaElectronInvMass(TString histoMass,TString histoDiff,vector <AliKFParticle> fKFGammas,\r
+ vector <TLorentzVector> TLVeNeg,vector<TLorentzVector> TLVePos)\r
+{\r
+\r
+\r
+ for( UInt_t iNeg=0; iNeg < TLVeNeg.size(); iNeg++ ){\r
+\r
+ for (UInt_t iPos=0; iPos < TLVePos.size(); iPos++){\r
+\r
+ TLorentzVector xy = TLVePos[iPos] + TLVeNeg[iNeg];\r
+\r
+ for (UInt_t iGam=0; iGam < fKFGammas.size(); iGam++){\r
+\r
+ AliKFParticle * GammaCandidate = &fKFGammas[iGam];\r
+ TLorentzVector g;\r
+\r
+ g.SetXYZM(GammaCandidate->GetPx(),GammaCandidate->GetPy(),GammaCandidate->GetPz(),fGammaMass);\r
+ TLorentzVector xyg = xy + g;\r
+ fHistograms->FillHistogram(histoMass.Data(),xyg.M());\r
+ fHistograms->FillHistogram(histoDiff.Data(),(xyg.M()-xy.M()));\r
+ }\r
+ }\r
+ }\r
+\r
+}\r
+void AliAnalysisTaskGammaConversion::ElectronBackground(TString hBg, vector <TLorentzVector> e)\r
+{\r
+ for(UInt_t i=0; i < e.size(); i++)\r
+ {\r
+ for (UInt_t j=i+1; j < e.size(); j++)\r
+ {\r
+ TLorentzVector ee = e[i] + e[j];\r
+\r
+ fHistograms->FillHistogram(hBg.Data(),ee.M());\r
+ }\r
+ }\r
+}\r
+\r
+\r
+void AliAnalysisTaskGammaConversion::CleanWithAngleCuts(vector <AliESDtrack*> NegativeElectrons,\r
+ vector <AliESDtrack*> PositiveElectrons, vector <AliKFParticle> Gammas){\r
+\r
+ UInt_t N = NegativeElectrons.size();\r
+ UInt_t P = PositiveElectrons.size();\r
+ UInt_t G = Gammas.size();\r
+\r
+\r
+\r
+ vector <Bool_t> xNegBand(N);\r
+ vector <Bool_t> xPosBand(P);\r
+ vector <Bool_t> gammaBand(G);\r
+\r
+\r
+ for(UInt_t iNeg=0; iNeg < N; iNeg++) xNegBand[iNeg]=kTRUE;\r
+ for(UInt_t iPos=0; iPos < P; iPos++) xPosBand[iPos]=kTRUE;\r
+ for(UInt_t iGam=0; iGam < G; iGam++) gammaBand[iGam]=kTRUE;\r
+\r
+\r
+ for(UInt_t iPos=0; iPos < P; iPos++){\r
+ \r
+ Double_t P[3]; PositiveElectrons[iPos]->GetConstrainedPxPyPz(P); \r
+\r
+ TVector3 ePosV(P[0],P[1],P[2]);\r
+\r
+ for(UInt_t iNeg=0; iNeg < N; iNeg++){\r
+ \r
+ Double_t N[3]; NegativeElectrons[iNeg]->GetConstrainedPxPyPz(N); \r
+ TVector3 eNegV(N[0],N[1],N[2]);\r
+\r
+ if(ePosV.Angle(eNegV) < 0.05){ //e+e- from gamma\r
+ xPosBand[iPos]=kFALSE;\r
+ xNegBand[iNeg]=kFALSE;\r
+ }\r
+\r
+ for(UInt_t iGam=0; iGam < G; iGam++){\r
+ AliKFParticle* GammaCandidate = &Gammas[iGam];\r
+ TVector3 GammaCandidateVector(GammaCandidate->Px(),GammaCandidate->Py(),GammaCandidate->Pz());\r
+ if(ePosV.Angle(GammaCandidateVector) < 0.05 || eNegV.Angle(GammaCandidateVector) < 0.05)\r
+ gammaBand[iGam]=kFALSE;\r
+ }\r
+ }\r
+ }\r
+\r
+\r
+\r
+\r
+ for(UInt_t iPos=0; iPos < P; iPos++){\r
+ if(xPosBand[iPos]){\r
+ fCurrentEventPosElectron.push_back(PositiveElectrons[iPos]);\r
+ }\r
+ }\r
+ for(UInt_t iNeg=0;iNeg < N; iNeg++){\r
+ if(xNegBand[iNeg]){\r
+ fCurrentEventNegElectron.push_back(NegativeElectrons[iNeg]);\r
+ }\r
+ }\r
+ for(UInt_t iGam=0; iGam < G; iGam++){\r
+ if(gammaBand[iGam]){\r
+ fKFReconstructedGammasCut.push_back(Gammas[iGam]);\r
+ }\r
+ }\r
+}\r
+\r
+\r
+void AliAnalysisTaskGammaConversion::GetPID(AliESDtrack *track, Stat_t &pid, Stat_t &weight)\r
+{\r
+ pid = -1;\r
+ weight = -1;\r
+\r
+ double wpart[5];\r
+ double wpartbayes[5];\r
+\r
+ //get probability of the diffenrent particle types\r
+ track->GetESDpid(wpart);\r
+\r
+ // Tentative particle type "concentrations"\r
+ double c[5]={0.01, 0.01, 0.85, 0.10, 0.05};\r
+\r
+ //Bayes' formula\r
+ double rcc = 0.;\r
+ for (int i = 0; i < 5; i++)\r
+ {\r
+ rcc+=(c[i] * wpart[i]);\r
+ }\r
+\r
+\r
+\r
+ for (int i=0; i<5; i++) {\r
+ if( rcc!=0){\r
+ wpartbayes[i] = c[i] * wpart[i] / rcc;\r
+ }\r
+ }\r
+\r
+\r
+\r
+ Float_t max=0.;\r
+ int ipid=-1;\r
+ //find most probable particle in ESD pid\r
+ //0:Electron - 1:Muon - 2:Pion - 3:Kaon - 4:Proton\r
+ for (int i = 0; i < 5; i++)\r
+ {\r
+ if (wpartbayes[i] > max)\r
+ {\r
+ ipid = i;\r
+ max = wpartbayes[i];\r
+ }\r
+ }\r
+\r
+ pid = ipid;\r
+ weight = max;\r
+}\r
+double AliAnalysisTaskGammaConversion::GetSigmaToVertex(AliESDtrack* t)\r
+{\r
+ // Calculates the number of sigma to the vertex.\r
+\r
+ Float_t b[2];\r
+ Float_t bRes[2];\r
+ Float_t bCov[3];\r
+ t->GetImpactParameters(b,bCov);\r
+ if (bCov[0]<=0 || bCov[2]<=0) {\r
+ AliDebug(1, "Estimated b resolution lower or equal zero!");\r
+ bCov[0]=0; bCov[2]=0;\r
+ }\r
+ bRes[0] = TMath::Sqrt(bCov[0]);\r
+ bRes[1] = TMath::Sqrt(bCov[2]);\r
+\r
+ // -----------------------------------\r
+ // How to get to a n-sigma cut?\r
+ //\r
+ // The accumulated statistics from 0 to d is\r
+ //\r
+ // -> Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)\r
+ // -> 1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)\r
+ //\r
+ // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)\r
+ // Can this be expressed in a different way?\r
+\r
+ if (bRes[0] == 0 || bRes[1] ==0)\r
+ return -1;\r
+\r
+ double d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));\r
+\r
+ // stupid rounding problem screws up everything:\r
+ // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(\r
+ if (TMath::Exp(-d * d / 2) < 1e-10)\r
+ return 1000;\r
+\r
+\r
+ d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);\r
+ return d;\r
+}\r
+vector <TLorentzVector> AliAnalysisTaskGammaConversion::GetTLorentzVector(vector <AliESDtrack*> ESDtrack){\r
+\r
+ vector <TLorentzVector> TLVtrack(0);\r
+\r
+ for(UInt_t itrack=0; itrack < ESDtrack.size(); itrack++){\r
+ double P[3]; ESDtrack[itrack]->GetConstrainedPxPyPz(P);\r
+ TLorentzVector CurrentTrack;\r
+ CurrentTrack.SetXYZM(P[0],P[1],P[2],fElectronMass);\r
+ TLVtrack.push_back(CurrentTrack);\r
+ }\r
+\r
+ return TLVtrack;\r
+}\r
+\r
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff