--- /dev/null
+#ifndef ALIANALYSISTASKDIMUONCFCONTAINERBUILDER_CXX
+#define ALIANALYSISTASKDIMUONCFCONTAINERBUILDER_CXX
+
+#include "AliAnalysisTaskDimuonCFContainerBuilder.h"
+#include "AliHeader.h"
+#include "AliESDHeader.h"
+#include "AliStack.h"
+#include "TParticle.h"
+#include "TString.h"
+#include "TLorentzVector.h"
+#include "TFile.h"
+#include "TH1I.h"
+#include "TH2D.h"
+#include "AliMCEvent.h"
+#include "AliAnalysisManager.h"
+#include "AliESDEvent.h"
+#include "AliAODEvent.h"
+#include "AliCFManager.h"
+#include "AliCFCutBase.h"
+#include "AliCFContainer.h"
+#include "AliCFEffGrid.h"
+#include "AliCFDataGrid.h"
+#include "TChain.h"
+#include "AliESDtrack.h"
+#include "AliLog.h"
+#include "AliESDMuonTrack.h"
+#include "AliESDtrack.h"
+#include "AliESDInputHandler.h"
+#include "TCanvas.h"
+#include "AliAODInputHandler.h"
+#include "AliAODMCParticle.h"
+
+// Analysis task for the building of a dimuon CF container
+// Also some single-muon variables are stored
+// L. Bianchi - Universita' & INFN Torino
+
+
+ClassImp(AliAnalysisTaskDimuonCFContainerBuilder)
+
+//__________________________________________________________________________
+AliAnalysisTaskDimuonCFContainerBuilder::AliAnalysisTaskDimuonCFContainerBuilder() :
+ fReadAODData(0),
+ fReadMCInfo(kTRUE),
+ fIsAccProduction(kTRUE),
+ fCFManager(0x0),
+ fQAHistList(0x0),
+ fNevt(0),
+ fBeamEnergy(3500.),
+ fOutput(0x0),
+ fCutOnzVtxSPD(kFALSE),
+ fCutOnNContributors(kFALSE),
+ fTrigClassMuon(""),
+ fTrigClassInteraction(""),
+ fDistinguishTrigClass(kFALSE)
+{
+
+ //Default constructor
+
+ Double_t chilims[2]={0.,10000.};
+ Double_t ptlims[2]={0.,100.};
+ Double_t thetalims[2]={0.,180.};
+ Double_t vtxlims[2]={-1000.,1000.};
+ SetChi2Limits(chilims);
+ SetChi2MatchLimits(chilims);
+ SetPtSingMuLimits(ptlims);
+ SetThetaSingMuLimits(thetalims);
+ SetZprimVertLimits(vtxlims);
+ SetTrigClassMuonName();
+ SetTrigClassInteracName();
+ TString namemuonside[4]={"CMUS1A-","CMUS1B-","CMUS1C-","CMUS1-E-"};
+ TString nameinteractionside[4]={"CINT1A-","CINT1B-","CINT1C-","CINT1-E"};
+ SetTrigClassMuonSideName(namemuonside);
+ SetTrigClassInteracSideName(nameinteractionside);
+}
+//___________________________________________________________________________
+AliAnalysisTaskDimuonCFContainerBuilder::AliAnalysisTaskDimuonCFContainerBuilder(const Char_t* name, Bool_t readaod, Bool_t readMC, Bool_t isaccept, Double_t beamEn) :
+ AliAnalysisTaskSE(name),
+ fReadAODData(0),
+ fReadMCInfo(kTRUE),
+ fIsAccProduction(kTRUE),
+ fCFManager(0x0),
+ fQAHistList(0x0),
+ fNevt(0),
+ fBeamEnergy(3500.),
+ fOutput(0x0),
+ fCutOnzVtxSPD(kFALSE),
+ fCutOnNContributors(kFALSE),
+ fTrigClassMuon(""),
+ fTrigClassInteraction(""),
+ fDistinguishTrigClass(kFALSE)
+{
+ //
+ // Constructor. Initialization of Inputs and Outputs
+ //
+ Info("AliAnalysisTaskDimuonCFContainerBuilder","Calling Constructor");
+
+ SetReadAODData(readaod);
+ SetReadMCinfo(readMC);
+ SetIsAccProd(isaccept);
+ SetBeamEnergy(beamEn);
+
+ Double_t chilims[2]={0.,10000.};
+ Double_t ptlims[2]={0.,100.};
+ Double_t thetalims[2]={0.,180.};
+ Double_t vtxlims[2]={-1000.,1000.};
+ SetChi2Limits(chilims);
+ SetChi2MatchLimits(chilims);
+ SetPtSingMuLimits(ptlims);
+ SetThetaSingMuLimits(thetalims);
+ SetZprimVertLimits(vtxlims);
+ SetTrigClassMuonName();
+ SetTrigClassInteracName();
+ TString namemuonside[4]={"CMUS1A-","CMUS1B-","CMUS1C-","CMUS1-E-"};
+ TString nameinteractionside[4]={"CINT1A-","CINT1B-","CINT1C-","CINT1-E"};
+ SetTrigClassMuonSideName(namemuonside);
+ SetTrigClassInteracSideName(nameinteractionside);
+
+ DefineOutput(1,TList::Class());
+ DefineOutput(2,AliCFContainer::Class());
+
+}
+
+//___________________________________________________________________________
+AliAnalysisTaskDimuonCFContainerBuilder& AliAnalysisTaskDimuonCFContainerBuilder::operator=(const AliAnalysisTaskDimuonCFContainerBuilder& c)
+{
+ //
+ // Assignment operator
+ //
+ if (this!=&c) {
+ AliAnalysisTaskSE::operator=(c) ;
+ fReadAODData = c.fReadAODData ;
+ fCFManager = c.fCFManager;
+ fQAHistList = c.fQAHistList ;
+ fNevt = c.fNevt ;
+ }
+ return *this;
+}
+
+//___________________________________________________________________________
+AliAnalysisTaskDimuonCFContainerBuilder::AliAnalysisTaskDimuonCFContainerBuilder(const AliAnalysisTaskDimuonCFContainerBuilder& c) :
+ AliAnalysisTaskSE(c),
+ fReadAODData(c.fReadAODData),
+ fReadMCInfo(c.fReadMCInfo),
+ fIsAccProduction(c.fIsAccProduction),
+ fCFManager(c.fCFManager),
+ fQAHistList(c.fQAHistList),
+ fNevt(c.fNevt),
+ fBeamEnergy(c.fBeamEnergy),
+ fOutput(c.fOutput),
+ fCutOnzVtxSPD(c.fCutOnzVtxSPD),
+ fCutOnNContributors(c.fCutOnNContributors),
+ fTrigClassMuon(c.fTrigClassMuon),
+ fTrigClassInteraction(c.fTrigClassInteraction),
+ fDistinguishTrigClass(c.fDistinguishTrigClass)
+{
+
+ // Copy Constructor
+
+}
+
+//___________________________________________________________________________
+AliAnalysisTaskDimuonCFContainerBuilder::~AliAnalysisTaskDimuonCFContainerBuilder() {
+ //
+ //destructor
+ //
+ Info("~AliAnalysisTaskDimuonCFContainerBuilder","Calling Destructor");
+ if (fCFManager) delete fCFManager ;
+ if (fQAHistList) {fQAHistList->Clear(); delete fQAHistList;}
+}
+
+//___________________________________________________________________________
+void AliAnalysisTaskDimuonCFContainerBuilder::UserCreateOutputObjects(){
+
+ fOutput = new TList();
+ fOutput->SetOwner();
+
+ TH1D *hnevts = new TH1D("hnevts","hnevts",1,0,1); // Stat check histos
+ TH1D *JpsiMult = new TH1D("JpsiMult","JpsiMult",20,0,20);
+
+ TH1D *ptdimuREC = new TH1D("ptdimuREC","ptdimuREC",200,0,20); // Dimu check histos
+ TH1D *ydimuREC = new TH1D("ydimuREC","ydimuREC",200,-10.,10.);
+ TH1D *ydimuFail = new TH1D("ydimuFail","ydimuFail",10,660,670);
+ TH1D *costHEdimuREC= new TH1D("costHEdimuREC","costHEdimuREC",200,-1.,1.);
+ TH1D *costCSdimuREC= new TH1D("costCSdimuREC","costCSdimuREC",200,-1.,1.);
+ TH1D *costHEdimuFail= new TH1D("costHEdimuFail","costHEdimuFail",10,660.,670.);
+ TH1D *costCSdimuFail= new TH1D("costCSdimuFail","costCSdimuFail",10,660.,670.);
+ TH1D *phiHEdimuREC = new TH1D("phiHEdimuREC","phiHEdimuREC",100,0.,TMath::Pi());
+ TH1D *phiCSdimuREC = new TH1D("phiCSdimuREC","phiCSdimuREC",100,0.,TMath::Pi());
+ TH1D *phiHEdimuFail = new TH1D("phiHEdimuFail","phiHEdimuFail",10,660.,670.);
+ TH1D *phiCSdimuFail = new TH1D("phiCSdimuFail","phiCSdimuFail",10,660.,670.);
+ TH1D *imassTot = new TH1D("imassTot","imassTot",100,0,8);
+ TH1D *trigCond = new TH1D("trigCond","trigCond",40,0,4);
+
+ TH1D *EmuonREC = new TH1D("EmuonREC","EmuonREC",200,0.,20.); // Mu check histos
+ TH1D *ptmuonREC= new TH1D("ptmuonREC","ptmuonREC",200,0.,20.);
+ TH1D *ymuonREC = new TH1D("ymuonREC","ymuonREC",200,-10.,10.);
+ TH1D *hdca = new TH1D("hdca","hdca",20,0,200);
+ TH2D *hdcay = new TH2D("hdcay","hdcay",200,0,200,20,-5,0);
+
+ TH1D *zvSPD = new TH1D("zvSPD","zvSPD",100,-50,50.); // Event check histos
+ TH1D *zvSPDcut = new TH1D("zvSPDcut","zvSPDcut",100,-50,50.);
+ TH1D *NContrib = new TH1D("NContrib","NContrib",100,0,100);
+
+
+ fOutput->Add(hnevts);
+ fOutput->Add(JpsiMult);
+
+ fOutput->Add(ptdimuREC);
+ fOutput->Add(ydimuREC);
+ fOutput->Add(ydimuFail);
+ fOutput->Add(costHEdimuREC);
+ fOutput->Add(costCSdimuREC);
+ fOutput->Add(costHEdimuFail);
+ fOutput->Add(costCSdimuFail);
+ fOutput->Add(phiHEdimuREC);
+ fOutput->Add(phiCSdimuREC);
+ fOutput->Add(phiHEdimuFail);
+ fOutput->Add(phiCSdimuFail);
+ fOutput->Add(imassTot);
+ fOutput->Add(trigCond);
+
+ fOutput->Add(EmuonREC);
+ fOutput->Add(ptmuonREC);
+ fOutput->Add(ymuonREC);
+ fOutput->Add(hdca);
+ fOutput->Add(hdcay);
+
+ fOutput->Add(zvSPD);
+ fOutput->Add(zvSPDcut);
+ fOutput->Add(NContrib);
+
+
+ fOutput->ls();
+
+
+}
+
+
+//_________________________________________________
+void AliAnalysisTaskDimuonCFContainerBuilder::UserExec(Option_t *)
+{
+
+
+ //Info("UserExec","");
+ fNevt++;
+ ((TH1D*)(fOutput->FindObject("hnevts")))->Fill(0.5);
+
+ Double_t containerInput[15]={666,666,666,666,666,666,666,666,666,666,666,666,666,666,666}; //y, pT, costHE, phiHE, costCS, phiCS, mass, TrigCond
+ Int_t cutAccept =1;
+ Int_t NumbJpsis =0;
+
+ if (!fReadAODData){ // ESD-based ANALYSIS
+
+ if (fReadMCInfo){
+
+ if (!fMCEvent) {
+ Error("UserExec","NO MC EVENT FOUND!");
+ //return;
+ }
+
+ // MC part ---------------------------------------------------------------------------------------------------
+
+ //fCFManager->SetEventInfo(fMCEvent); CHANGES IN NEW VERSIONS - MANUAL CUT ON PDG
+ AliStack* stack = fMCEvent->Stack();
+
+ // loop on the MC event
+ for (Int_t ipart=0; ipart<fMCEvent->GetNumberOfTracks(); ipart++) {
+
+ AliMCParticle *mcPart = (AliMCParticle*) fMCEvent->GetTrack(ipart);
+
+ TParticle *part = mcPart->Particle();
+ TParticle *part0 = mcPart->Particle();
+ TParticle *part1 = mcPart->Particle();
+
+ // Mother kinematics
+ Double_t e = part->Energy();
+ Double_t pz = part->Pz();
+ Double_t rapmc = Rap(e,pz);
+
+ // Selection of the resonance
+ //if (!fCFManager->CheckParticleCuts(0,mcPart)) continue;
+ if (part->GetPdgCode()!=443) continue; // MANUAL CUT ON PDG
+ NumbJpsis++;
+
+ // Decays kinematics
+ Int_t p0 = part->GetDaughter(0);
+ part0 = stack->Particle(p0);
+ Int_t pdg0 = part0->GetPdgCode();
+
+ Int_t p1 = part->GetDaughter(1);
+ part1 = stack->Particle(p1);
+ Int_t pdg1 = part1->GetPdgCode();
+
+ Double_t e0 = part0->Energy();
+ Double_t pz0 = part0->Pz();
+ Double_t py0 = part0->Py();
+ Double_t px0 = part0->Px();
+ Double_t charge0 = part0->GetPDG()->Charge()/3;
+ Double_t theta0 = (180./TMath::Pi())*TMath::ATan2(TMath::Sqrt(px0*px0+py0*py0),pz0);
+ Double_t pt0 = TMath::Sqrt(px0*px0+py0*py0);
+ Double_t Mom0 = part0->P();
+
+ Double_t e1 = part1->Energy();
+ Double_t pz1 = part1->Pz();
+ Double_t py1 = part1->Py();
+ Double_t px1 = part1->Px();
+ //Double_t charge1 = part1->GetPDG()->Charge()/3;
+ Double_t theta1 = (180./TMath::Pi())*TMath::ATan2(TMath::Sqrt(px1*px1+py1*py1),pz1);
+ Double_t pt1 = TMath::Sqrt(px1*px1+py1*py1);
+ Double_t Mom1 = part1->P();
+
+
+ if(pdg0==13 || pdg1==13) {
+
+ Double_t ptmc = TMath::Sqrt((px0+px1)*(px0+px1)+(py0+py1)*(py0+py1));
+ Double_t imassmc = Imass(e0,px0,py0,pz0,e1,px1,py1,pz1);
+
+ Double_t costCS = CostCS(px0,py0,pz0,e0,charge0,px1,py1,pz1,e1);
+ Double_t costHE = CostHE(px0,py0,pz0,e0,charge0,px1,py1,pz1,e1);
+ Double_t phiCS = PhiCS(px0,py0,pz0,e0,charge0,px1,py1,pz1,e1);
+ Double_t phiHE = PhiHE(px0,py0,pz0,e0,charge0,px1,py1,pz1,e1);
+
+ containerInput[0] = rapmc ;
+ containerInput[1] = ptmc;
+ containerInput[2] = costHE;
+ containerInput[3] = TMath::Abs(phiHE);
+ containerInput[4] = costCS;
+ containerInput[5] = TMath::Abs(phiCS);
+ containerInput[6] = imassmc;
+ containerInput[7] = 1.; //for generated no trigger condition
+ if (pt0<pt1) {
+ containerInput[8]=pt0;
+ containerInput[9]=pt1;
+ } else {
+ containerInput[8]=pt1;
+ containerInput[9]=pt0;
+ }
+ if (theta0<theta1) {
+ containerInput[10]=theta0;
+ containerInput[11]=theta1;
+ } else {
+ containerInput[10]=theta1;
+ containerInput[11]=theta0;
+ }
+ if (Mom0<Mom1) {
+ containerInput[12]=Mom0;
+ containerInput[13]=Mom1;
+ } else {
+ containerInput[12]=Mom1;
+ containerInput[13]=Mom0;
+ }
+ containerInput[14]=1.;
+
+ // fill the container at the first step
+ fCFManager->GetParticleContainer()->Fill(containerInput,0);
+
+ if(fIsAccProduction){
+ // acceptance cuts on single mu
+ if(theta0<fThetaSingMuCut[0]||theta0>fThetaSingMuCut[1]||theta1<fThetaSingMuCut[0]||theta1>fThetaSingMuCut[1]) cutAccept=0;
+ if(pt0<fPtSingMuCut[0] || pt0>fPtSingMuCut[1] || pt1<fPtSingMuCut[0] || pt1>fPtSingMuCut[1]) cutAccept=0;
+ if (cutAccept ==1) fCFManager->GetParticleContainer()->Fill(containerInput,2);
+ }
+ }
+ }
+ } //fReadMCInfo
+
+
+ // ESD part ---------------------------------------------------------------------------------------------------
+
+ AliESDEvent *fESD;
+ AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
+ fESD = esdH->GetEvent();
+ Int_t mult1 = fESD->GetNumberOfMuonTracks() ;
+
+ Int_t trigfired=-1;
+ Int_t trigside=-1;
+ if(fDistinguishTrigClass){
+ //aod->GetHeader()->SetFiredTriggerClasses("CMU");
+ TString trigclass = fESD->GetFiredTriggerClasses();
+ printf("TrigClass: %s\n",trigclass.Data());
+ if(trigclass.Contains(fTrigClassMuon)) trigfired = 1;
+ else if (trigclass.Contains(fTrigClassInteraction)) trigfired = 0;
+ printf("Muon container: %d\n",trigfired);
+ for(Int_t i=0;i<4;i++){
+ if(trigfired==1 && trigclass.Contains(fTrigClassMuonSide[i])) {trigside=i;printf("Muon Side: %d\n",trigside);}
+ if(trigfired==0 && trigclass.Contains(fTrigClassInteractionSide[i])) {trigside=i;printf("Interaction Side: %d\n",trigside);}
+ }
+ }
+
+ ((TH1D*)(fOutput->FindObject("zvSPD")))->Fill(fESD->GetPrimaryVertexSPD()->GetZ());
+ ((TH1D*)(fOutput->FindObject("NContrib")))->Fill(fESD->GetPrimaryVertexSPD()->GetNContributors());
+ if (!fCutOnzVtxSPD || (fCutOnzVtxSPD && (fESD->GetPrimaryVertexSPD()->GetZ()>fzPrimVertexSPD[0] && fESD->GetPrimaryVertexSPD()->GetZ()<fzPrimVertexSPD[1]))){
+ ((TH1D*)(fOutput->FindObject("zvSPDcut")))->Fill(fESD->GetPrimaryVertexSPD()->GetZ());
+ if (!fCutOnNContributors || (fCutOnNContributors && (fESD->GetPrimaryVertexSPD()->GetNContributors()>fNContributors[0] && fESD->GetPrimaryVertexSPD()->GetNContributors()<fNContributors[1]))){
+ for (Int_t j = 0; j<mult1; j++) {
+
+ AliESDMuonTrack* mu1 = new AliESDMuonTrack(*(fESD->GetMuonTrack(j)));
+ if (!mu1->ContainTrackerData()) continue;
+ if (mu1->GetChi2()<fChi2Track[0] || mu1->GetChi2()>fChi2Track[1]) continue;
+ if (mu1->GetChi2MatchTrigger()<fChi2MatchTrig[0] || mu1->GetChi2MatchTrigger()>fChi2MatchTrig[1]) continue;
+ Double_t zr1 = mu1->Charge();
+ Double_t pxr1 = mu1->Px();
+ Double_t pyr1 = mu1->Py();
+ Double_t pzr1 = mu1->Pz();
+ Double_t ptmu1 = TMath::Sqrt(pxr1*pxr1+pyr1*pyr1);
+ Double_t er1 = mu1->E();
+ Double_t Mom1 = mu1->P();
+ Double_t theta1 = (180./TMath::Pi())*mu1->Theta();
+ Double_t rapiditymu1 = Rap(er1,pzr1);
+ ((TH1D*)(fOutput->FindObject("EmuonREC")))->Fill(er1);
+ ((TH1D*)(fOutput->FindObject("ptmuonREC")))->Fill(ptmu1);
+ ((TH1D*)(fOutput->FindObject("ymuonREC")))->Fill(rapiditymu1);
+
+ if(zr1<0){
+
+ for (Int_t jj = 0; jj<mult1; jj++) {
+
+ AliESDMuonTrack* mu2 = new AliESDMuonTrack(*(fESD->GetMuonTrack(jj)));
+ if (!mu2->ContainTrackerData()) continue;
+ if (mu2->GetChi2()<fChi2Track[0] || mu2->GetChi2()>fChi2Track[1]) continue;
+ if (mu2->GetChi2MatchTrigger()<fChi2MatchTrig[0] || mu2->GetChi2MatchTrigger()>fChi2MatchTrig[1]) continue;
+ Double_t zr2 = mu2->Charge();
+
+ if(zr2>0){
+ Double_t TrigCondition=0;
+ if (mu1->GetMatchTrigger()>=mu2->GetMatchTrigger()) TrigCondition = mu1->GetMatchTrigger()+0.1*mu2->GetMatchTrigger();
+ else TrigCondition = mu2->GetMatchTrigger()+0.1*mu1->GetMatchTrigger();
+ ((TH1D*)(fOutput->FindObject("trigCond")))->Fill(TrigCondition);
+ Double_t pxr2 = mu2->Px();
+ Double_t pyr2 = mu2->Py();
+ Double_t pzr2 = mu2->Pz();
+ Double_t ptmu2 = TMath::Sqrt(pxr2*pxr2+pyr2*pyr2);
+ Double_t er2 = mu2->E();
+ Double_t Mom2 = mu2->P();
+ Double_t theta2 = (180./TMath::Pi())*mu2->Theta();
+
+ Double_t ptrec = TMath::Sqrt((pxr1+pxr2)*(pxr1+pxr2)+(pyr1+pyr2)*(pyr1+pyr2));
+ ((TH1D*)(fOutput->FindObject("ptdimuREC")))->Fill(ptrec);
+ Double_t raprec= Rap((er1+er2),(pzr1+pzr2));
+ ((TH1D*)(fOutput->FindObject("ydimuREC")))->Fill(raprec);
+ Double_t imassrec = Imass(er1,pxr1,pyr1,pzr1,er2,pxr2,pyr2,pzr2);
+ ((TH1D*)(fOutput->FindObject("imassTot")))->Fill(imassrec);
+
+ if (imassrec>12.) continue;
+
+ Double_t costCSrec = CostCS(pxr1,pyr1,pzr1,er1,zr1,pxr2,pyr2,pzr2,er2);
+ ((TH1D*)(fOutput->FindObject("costCSdimuREC")))->Fill(costCSrec);
+ if(costCSrec==666) ((TH1D*)(fOutput->FindObject("costCSdimuFail")))->Fill(costCSrec);
+ Double_t costHErec = CostHE(pxr1,pyr1,pzr1,er1,zr1,pxr2,pyr2,pzr2,er2);
+ ((TH1D*)(fOutput->FindObject("costHEdimuREC")))->Fill(costHErec);
+ if(costHErec==666) ((TH1D*)(fOutput->FindObject("costHEdimuFail")))->Fill(costHErec);
+ Double_t phiCSrec = PhiCS(pxr1,pyr1,pzr1,er1,zr1,pxr2,pyr2,pzr2,er2);
+ ((TH1D*)(fOutput->FindObject("phiCSdimuREC")))->Fill(phiCSrec);
+ if(phiCSrec==666) ((TH1D*)(fOutput->FindObject("phiCSdimuFail")))->Fill(phiCSrec);
+ Double_t phiHErec = PhiHE(pxr1,pyr1,pzr1,er1,zr1,pxr2,pyr2,pzr2,er2);
+ ((TH1D*)(fOutput->FindObject("phiHEdimuREC")))->Fill(phiHErec);
+ if(phiHErec==666) ((TH1D*)(fOutput->FindObject("phiHEdimuFail")))->Fill(phiHErec);
+
+ containerInput[0] = raprec ;
+ containerInput[1] = ptrec;
+ containerInput[2] = costHErec;
+ containerInput[3] = TMath::Abs(phiHErec);
+ containerInput[4] = costCSrec;
+ containerInput[5] = TMath::Abs(phiCSrec);
+ containerInput[6] = imassrec;
+ containerInput[7] = TrigCondition+0.05;
+ if (ptmu1<ptmu2) {
+ containerInput[8]=ptmu1;
+ containerInput[9]=ptmu2;
+ } else {
+ containerInput[8]=ptmu2;
+ containerInput[9]=ptmu1;
+ }
+ if (theta1<theta2) {
+ containerInput[10]=theta1;
+ containerInput[11]=theta2;
+ } else {
+ containerInput[10]=theta2;
+ containerInput[11]=theta1;
+ }
+ if (Mom1<Mom2) {
+ containerInput[12]=Mom1;
+ containerInput[13]=Mom2;
+ } else {
+ containerInput[12]=Mom2;
+ containerInput[13]=Mom1;
+ }
+ if (fDistinguishTrigClass && trigside==0) containerInput[14]=0.;
+ else if (fDistinguishTrigClass && trigside==1) containerInput[14]=1.;
+ else if (fDistinguishTrigClass && trigside==2) containerInput[14]=2.;
+ else if (fDistinguishTrigClass && trigside==3) containerInput[14]=3.;
+ else containerInput[14]=0.;
+
+ if(fDistinguishTrigClass){
+ if (trigfired==1) fCFManager->GetParticleContainer()->Fill(containerInput,5);
+ else if (trigfired==0) fCFManager->GetParticleContainer()->Fill(containerInput,1);
+ } else {
+ fCFManager->GetParticleContainer()->Fill(containerInput,1);
+ if(fIsAccProduction){
+ if(cutAccept==1) fCFManager->GetParticleContainer()->Fill(containerInput,3);
+ }
+ }
+ } // mu+ Selection
+
+ } // second mu Loop
+ } // mu- Selection
+ }
+ }
+ }
+ } else { // AOD-based ANALYSIS
+
+ AliAODEvent *aod;
+ AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
+ aod = aodH->GetEvent();
+ Int_t ntracks=aod->GetNumberOfTracks();
+
+ // MC part ---------------------------------------------------------------------------------------------------
+
+ if (fReadMCInfo){
+ TClonesArray *mcarray = dynamic_cast<TClonesArray*> (aod->FindListObject(AliAODMCParticle::StdBranchName())); //array of MC particles in this event
+ for(int ii=0;ii<mcarray->GetEntries();ii++){
+ AliAODMCParticle *mctrack = (AliAODMCParticle*) mcarray->At(ii);
+ if(mctrack->GetPdgCode()!=13) continue;
+ Int_t NumbMother = mctrack->GetMother();
+ if (NumbMother==-1) continue;
+ AliAODMCParticle *mother = (AliAODMCParticle*) mcarray->At(NumbMother);
+ if (mother->GetPdgCode()!=443) continue;
+ NumbJpsis++;
+ Int_t daught0 = mother->GetDaughter(0);
+ Int_t daught1 = mother->GetDaughter(1);
+ AliAODMCParticle *mcDaughter0 = (AliAODMCParticle*) mcarray->At(daught0);
+ Double_t pxmc0 = mcDaughter0->Px();
+ Double_t pymc0 = mcDaughter0->Py();
+ Double_t pzmc0 = mcDaughter0->Pz();
+ Double_t ptmc0 = mcDaughter0->Pt();
+ Double_t Mommc0 = mcDaughter0->P();
+ Double_t Emc0 = mcDaughter0->E();
+ Double_t thetamc0 = (180./TMath::Pi())*mcDaughter0->Theta();
+ Int_t charge0 = (Int_t) mcDaughter0->Charge()/3;
+ AliAODMCParticle *mcDaughter1 = (AliAODMCParticle*) mcarray->At(daught1);
+ Double_t pxmc1 = mcDaughter1->Px();
+ Double_t pymc1 = mcDaughter1->Py();
+ Double_t pzmc1 = mcDaughter1->Pz();
+ Double_t ptmc1 = mcDaughter1->Pt();
+ Double_t Mommc1 = mcDaughter1->P();
+ Double_t Emc1 = mcDaughter1->E();
+ Double_t thetamc1 = (180./TMath::Pi())*mcDaughter1->Theta();
+ Int_t charge1 = (Int_t) mcDaughter1->Charge()/3;
+
+ if (charge0==charge1 || TMath::Abs(mcDaughter0->GetPdgCode())!=13 || TMath::Abs(mcDaughter1->GetPdgCode())!=13) continue;
+ //Double_t rapJpsi = mother->Y(); //PROBLEM!!!
+ Double_t rapJpsi = Rap(mother->E(),mother->Pz());
+ Double_t ptJpsi = mother->Pt();
+ Double_t costHE = CostHE(pxmc0,pymc0,pzmc0,Emc0,(Double_t)charge0,pxmc1,pymc1,pzmc1,Emc1);
+ Double_t phiHE = PhiHE(pxmc0,pymc0,pzmc0,Emc0,(Double_t)charge0,pxmc1,pymc1,pzmc1,Emc1);
+ Double_t costCS = CostCS(pxmc0,pymc0,pzmc0,Emc0,(Double_t)charge0,pxmc1,pymc1,pzmc1,Emc1);
+ Double_t phiCS = PhiCS(pxmc0,pymc0,pzmc0,Emc0,(Double_t)charge0,pxmc1,pymc1,pzmc1,Emc1);
+ Double_t massJpsi = mother->M();
+ containerInput[0] = rapJpsi;
+ containerInput[1] = ptJpsi;
+ containerInput[2] = costHE;
+ containerInput[3] = TMath::Abs(phiHE);
+ containerInput[4] = costCS;
+ containerInput[5] = TMath::Abs(phiCS);
+ containerInput[6] = massJpsi;
+ containerInput[7] = 1.; //for generated no trigger condition
+ if (ptmc0<ptmc1) {
+ containerInput[8]=ptmc0;
+ containerInput[9]=ptmc1;
+ } else {
+ containerInput[8]=ptmc1;
+ containerInput[9]=ptmc0;
+ }
+ if (thetamc0<thetamc1) {
+ containerInput[10]=thetamc0;
+ containerInput[11]=thetamc1;
+ } else {
+ containerInput[10]=thetamc1;
+ containerInput[11]=thetamc0;
+ }
+ if (Mommc0<Mommc1) {
+ containerInput[12]=Mommc0;
+ containerInput[13]=Mommc1;
+ } else {
+ containerInput[12]=Mommc1;
+ containerInput[13]=Mommc0;
+ }
+ containerInput[14]=1.;
+
+ if(rapJpsi!=666 && costHE!=666 && phiHE!=666 && costCS!=666 && phiCS!=666){
+ fCFManager->GetParticleContainer()->Fill(containerInput,0);
+ if(fIsAccProduction){
+ // acceptance cuts on single mu
+ if(thetamc0<fThetaSingMuCut[0]||thetamc0>fThetaSingMuCut[1]||thetamc1<fThetaSingMuCut[0]||thetamc1>fThetaSingMuCut[1]) cutAccept=0;
+ if(ptmc0<fPtSingMuCut[0] || ptmc0>fPtSingMuCut[1] || ptmc1<fPtSingMuCut[0] || ptmc1>fPtSingMuCut[1]) cutAccept=0;
+ if (cutAccept ==1) fCFManager->GetParticleContainer()->Fill(containerInput,2);
+ }
+
+ }
+
+ }
+ }
+
+
+ // AOD part ---------------------------------------------------------------------------------------------------
+
+ Int_t trigfired=-1;
+ Int_t trigside=-1;
+ if(fDistinguishTrigClass){
+ TString trigclass = aod->GetFiredTriggerClasses();
+ printf("TrigClass: %s\n",trigclass.Data());
+ if(trigclass.Contains(fTrigClassMuon)) trigfired = 1;
+ else if (trigclass.Contains(fTrigClassInteraction)) trigfired = 0;
+ printf("Muon container: %d\n",trigfired);
+ for(Int_t i=0;i<4;i++){
+ if(trigfired==1 && trigclass.Contains(fTrigClassMuonSide[i])) {trigside=i;printf("Muon Side: %d\n",trigside);}
+ if(trigfired==0 && trigclass.Contains(fTrigClassInteractionSide[i])) {trigside=i;printf("Interaction Side: %d\n",trigside);}
+ }
+ }
+
+ ((TH1D*)(fOutput->FindObject("zvSPD")))->Fill(aod->GetPrimaryVertex()->GetZ());
+ Int_t ncontr = aod->GetPrimaryVertex()->GetNContributors();
+ ((TH1D*)(fOutput->FindObject("NContrib")))->Fill(ncontr);
+ if (!fCutOnzVtxSPD || (fCutOnzVtxSPD && (aod->GetPrimaryVertex()->GetZ()>fzPrimVertexSPD[0] && aod->GetPrimaryVertex()->GetZ()<fzPrimVertexSPD[1]))){ //NOT REALLY SPD VERTEX
+ ((TH1D*)(fOutput->FindObject("zvSPDcut")))->Fill(aod->GetPrimaryVertex()->GetZ());
+ if (!fCutOnNContributors || (fCutOnNContributors && (aod->GetPrimaryVertex()->GetNContributors()>fNContributors[0] && aod->GetPrimaryVertex()->GetNContributors()<fNContributors[1]))){
+ for (Int_t j = 0; j<ntracks; j++) {
+ AliAODTrack *mu1 = aod->GetTrack(j);
+ if(!mu1->IsMuonTrack()) continue;
+ if (mu1->Chi2perNDF()<fChi2Track[0] || mu1->Chi2perNDF()>fChi2Track[1]) continue;
+ if (mu1->GetChi2MatchTrigger()<fChi2MatchTrig[0] || mu1->GetChi2MatchTrigger()>fChi2MatchTrig[1]) continue;
+ Double_t chargemu1 = mu1->Charge();
+ Double_t pxmu1 = mu1->Px();
+ Double_t pymu1 = mu1->Py();
+ Double_t pzmu1 = mu1->Pz();
+ Double_t emu1 = mu1->E();
+ Double_t pmu1 = mu1->P();
+ Double_t ptmu1 = mu1->Pt();
+ Double_t rapiditymu1 = Rap(emu1,pzmu1);
+ Double_t thetamu1 = (180./TMath::Pi())*mu1->Theta();
+ Double_t Rdcamu1 = mu1->DCA();
+ ((TH1D*)(fOutput->FindObject("EmuonREC")))->Fill(emu1);
+ ((TH1D*)(fOutput->FindObject("ptmuonREC")))->Fill(ptmu1);
+ ((TH1D*)(fOutput->FindObject("ymuonREC")))->Fill(rapiditymu1);
+ if(chargemu1<0){
+ for (Int_t jj = 0; jj<ntracks; jj++) {
+ AliAODTrack *mu2 = aod->GetTrack(jj);
+ if(!mu2->IsMuonTrack()) continue;
+ if (mu2->Chi2perNDF()<fChi2Track[0] || mu2->Chi2perNDF()>fChi2Track[1]) continue;
+ if (mu2->GetChi2MatchTrigger()<fChi2MatchTrig[0] || mu2->GetChi2MatchTrigger()>fChi2MatchTrig[1]) continue;
+ Double_t chargemu2 = mu2->Charge();
+ Double_t pxmu2 = mu2->Px();
+ Double_t pymu2 = mu2->Py();
+ Double_t pzmu2 = mu2->Pz();
+ Double_t emu2 = mu2->E();
+ Double_t pmu2 = mu2->P();
+ Double_t ptmu2 = mu2->Pt();
+ //Double_t rapiditymu2 = Rap(emu2,pzmu2);
+ Double_t thetamu2 = (180./TMath::Pi())*mu2->Theta();
+ Double_t Rdcamu2 = mu2->DCA();
+ if(chargemu2>0){
+ if (mu1->GetMatchTrigger()>0) printf("Mu1: charge: %f, match: %d\n",chargemu1,1);
+ else printf("Mu1: charge: %f, match: %d\n",chargemu1,0);
+ if (mu2->GetMatchTrigger()>0) printf("Mu2: charge: %f, match: %d\n",chargemu2,1);
+ else printf("Mu2: charge: %f, match: %d\n",chargemu2,0);
+ ((TH1D*)(fOutput->FindObject("hdca")))->Fill(Rdcamu2);
+ ((TH1D*)(fOutput->FindObject("hdca")))->Fill(Rdcamu1);
+ Double_t TrigCondition=0;
+ if (mu1->GetMatchTrigger()>=mu2->GetMatchTrigger()) TrigCondition = mu1->GetMatchTrigger()+0.1*mu2->GetMatchTrigger();
+ else TrigCondition = mu2->GetMatchTrigger()+0.1*mu1->GetMatchTrigger();
+ containerInput[0] = (Double_t) Rap((emu1+emu2),(pzmu1+pzmu2));
+ ((TH1D*)(fOutput->FindObject("ydimuREC")))->Fill(containerInput[0]);
+ ((TH2D*)(fOutput->FindObject("hdcay")))->Fill(TMath::Max(Rdcamu1,Rdcamu2),containerInput[0]);
+ //((TH2D*)(fOutput->FindObject("hdcay")))->Fill(Rdcamu2,containerInput[0]);
+ containerInput[1] = (Double_t) TMath::Sqrt((pxmu1+pxmu2)*(pxmu1+pxmu2)+(pymu1+pymu2)*(pymu1+pymu2));
+ ((TH1D*)(fOutput->FindObject("ptdimuREC")))->Fill(containerInput[1]);
+ containerInput[2] = CostHE(pxmu1,pymu1,pzmu1,emu1,chargemu1,pxmu2,pymu2,pzmu2,emu2);
+ ((TH1D*)(fOutput->FindObject("costHEdimuREC")))->Fill(containerInput[2]);
+ if(containerInput[2]==666) ((TH1D*)(fOutput->FindObject("costHEdimuFail")))->Fill(containerInput[2]);
+ containerInput[3] = TMath::Abs(PhiHE(pxmu1,pymu1,pzmu1,emu1,chargemu1,pxmu2,pymu2,pzmu2,emu2));
+ ((TH1D*)(fOutput->FindObject("phiHEdimuREC")))->Fill(containerInput[3]);
+ if(containerInput[3]==666) ((TH1D*)(fOutput->FindObject("phiHEdimuFail")))->Fill(containerInput[3]);
+ containerInput[4] = CostCS(pxmu1,pymu1,pzmu1,emu1,chargemu1,pxmu2,pymu2,pzmu2,emu2);
+ ((TH1D*)(fOutput->FindObject("costCSdimuREC")))->Fill(containerInput[4]);
+ if(containerInput[4]==666) ((TH1D*)(fOutput->FindObject("costCSdimuFail")))->Fill(containerInput[4]);
+ containerInput[5] = TMath::Abs(PhiCS(pxmu1,pymu1,pzmu1,emu1,chargemu1,pxmu2,pymu2,pzmu2,emu2));
+ ((TH1D*)(fOutput->FindObject("phiCSdimuREC")))->Fill(containerInput[5]);
+ if(containerInput[5]==666) ((TH1D*)(fOutput->FindObject("phiCSdimuFail")))->Fill(containerInput[5]);
+ containerInput[6] = Imass(emu1,pxmu1,pymu1,pzmu1,emu2,pxmu2,pymu2,pzmu2);
+ ((TH1D*)(fOutput->FindObject("imassTot")))->Fill(containerInput[6]);
+ containerInput[7] = TrigCondition+0.05;
+ ((TH1D*)(fOutput->FindObject("trigCond")))->Fill(containerInput[7]);
+ if (ptmu1<ptmu2) {
+ containerInput[8]=ptmu1;
+ containerInput[9]=ptmu2;
+ } else {
+ containerInput[8]=ptmu2;
+ containerInput[9]=ptmu1;
+ }
+ if (thetamu1<thetamu2) {
+ containerInput[10]=thetamu1;
+ containerInput[11]=thetamu2;
+ } else {
+ containerInput[10]=thetamu2;
+ containerInput[11]=thetamu1;
+ }
+ if (pmu1<pmu2) {
+ containerInput[12]=pmu1;
+ containerInput[13]=pmu2;
+ } else {
+ containerInput[12]=pmu2;
+ containerInput[13]=pmu1;
+ }
+ if (fDistinguishTrigClass && trigside==0) containerInput[14]=0.;
+ else if (fDistinguishTrigClass && trigside==1) containerInput[14]=1.;
+ else if (fDistinguishTrigClass && trigside==2) containerInput[14]=2.;
+ else if (fDistinguishTrigClass && trigside==3) containerInput[14]=3.;
+ else containerInput[14]=0.;
+
+ if(containerInput[2]!=666 && containerInput[3]!=666 && containerInput[4]!=666 && containerInput[5]!=666){
+ if(fDistinguishTrigClass){
+ if (trigfired==1) fCFManager->GetParticleContainer()->Fill(containerInput,5);
+ else if (trigfired==0) fCFManager->GetParticleContainer()->Fill(containerInput,1);
+ } else {
+ fCFManager->GetParticleContainer()->Fill(containerInput,1);
+ if(fIsAccProduction){
+ if (cutAccept ==1) fCFManager->GetParticleContainer()->Fill(containerInput,3);
+ }
+ }
+ } else if (containerInput[0]==666) ((TH1D*)(fOutput->FindObject("ydimuFail")))->Fill(containerInput[0]);
+
+ }
+
+ }
+
+ }
+
+ }
+ }
+ }
+
+
+ }
+
+
+// ----------
+ if (fReadMCInfo) ((TH1D*)(fOutput->FindObject("JpsiMult")))->Fill(NumbJpsis);
+
+ PostData(1,fOutput) ;
+ PostData(2,fCFManager->GetParticleContainer()) ;
+
+
+}
+
+//________________________________________________________________________
+Double_t AliAnalysisTaskDimuonCFContainerBuilder::Imass(Double_t e1, Double_t px1, Double_t py1, Double_t pz1,
+ Double_t e2, Double_t px2, Double_t py2, Double_t pz2)
+{
+// invariant mass calculation
+ Double_t imassrec = TMath::Sqrt((e1+e2)*(e1+e2)-((px1+px2)*(px1+px2)+
+ (py1+py2)*(py1+py2)+(pz1+pz2)*(pz1+pz2)));
+ return imassrec;
+}
+
+//________________________________________________________________________
+Double_t AliAnalysisTaskDimuonCFContainerBuilder::Rap(Double_t e, Double_t pz)
+{
+// calculate rapidity
+ Double_t rap;
+ if(e>TMath::Abs(pz)){ // in order to avoid problems with AODs
+ rap = 0.5*TMath::Log((e+pz)/(e-pz));
+ return rap;
+ }
+ else{
+ rap = 666.;
+ return rap;
+ }
+}
+
+//________________________________________________________________________
+Double_t AliAnalysisTaskDimuonCFContainerBuilder::CostCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
+Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
+{
+// Cosine of the theta decay angle (mu+) in the Collins-Soper frame
+
+ TLorentzVector PMu1CM, PMu2CM, PProjCM, PTargCM, PDimuCM; // In the CM. frame
+ TLorentzVector PMu1Dimu, PMu2Dimu, PProjDimu, PTargDimu; // In the dimuon rest frame
+ TVector3 beta,zaxisCS;
+ Double_t mp=0.93827231;
+ //
+ // --- Fill the Lorentz vector for projectile and target in the CM frame
+ //
+ PProjCM.SetPxPyPzE(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp));
+ PTargCM.SetPxPyPzE(0.,0.,fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp));
+ //
+ // --- Get the muons parameters in the CM frame
+ //
+ PMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
+ PMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
+ //
+ // --- Obtain the dimuon parameters in the CM frame
+ //
+ PDimuCM=PMu1CM+PMu2CM;
+ //
+ // --- Translate the dimuon parameters in the dimuon rest frame
+ //
+ beta=(-1./PDimuCM.E())*PDimuCM.Vect();
+ if(beta.Mag()>=1) return 666.;
+ PMu1Dimu=PMu1CM;
+ PMu2Dimu=PMu2CM;
+ PProjDimu=PProjCM;
+ PTargDimu=PTargCM;
+ PMu1Dimu.Boost(beta);
+ PMu2Dimu.Boost(beta);
+ PProjDimu.Boost(beta);
+ PTargDimu.Boost(beta);
+
+ //Debugging part -------------------------------------
+ Double_t debugProj[4]={0.,0.,0.,0.};
+ Double_t debugTarg[4]={0.,0.,0.,0.};
+ Double_t debugMu1[4]={0.,0.,0.,0.};
+ Double_t debugMu2[4]={0.,0.,0.,0.};
+ PMu1Dimu.GetXYZT(debugMu1);
+ PMu2Dimu.GetXYZT(debugMu2);
+ PProjDimu.GetXYZT(debugProj);
+ PTargDimu.GetXYZT(debugTarg);
+ if (debugProj[0]!=debugProj[0] ||debugProj[1]!=debugProj[1] || debugProj[2]!=debugProj[2] ||debugProj[3]!=debugProj[3]) return 666;
+ if (debugTarg[0]!=debugTarg[0] ||debugTarg[1]!=debugTarg[1] || debugTarg[2]!=debugTarg[2] ||debugTarg[3]!=debugTarg[3]) return 666;
+ if (debugMu1[0]!=debugMu1[0] ||debugMu1[1]!=debugMu1[1] || debugMu1[2]!=debugMu1[2] ||debugMu1[3]!=debugMu1[3]) return 666;
+ if (debugMu2[0]!=debugMu2[0] ||debugMu2[1]!=debugMu2[1] || debugMu2[2]!=debugMu2[2] ||debugMu2[3]!=debugMu2[3]) return 666;
+ //----------------------------------------------------
+
+ // --- Determine the z axis for the CS angle
+ zaxisCS=(((PProjDimu.Vect()).Unit())-((PTargDimu.Vect()).Unit())).Unit();
+
+ // --- Determine the CS angle (angle between mu+ and the z axis defined above)
+ Double_t cost;
+
+ if(charge1>0) {cost = zaxisCS.Dot((PMu1Dimu.Vect()).Unit());}
+ else {cost = zaxisCS.Dot((PMu2Dimu.Vect()).Unit());}
+
+ return cost;
+}
+
+//________________________________________________________________________
+Double_t AliAnalysisTaskDimuonCFContainerBuilder::CostHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
+Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
+{
+// Cosine of the theta decay angle (mu+) in the Helicity frame
+
+ TLorentzVector PMu1CM, PMu2CM, PDimuCM; // In the CM frame
+ TLorentzVector PMu1Dimu, PMu2Dimu; // In the dimuon rest frame
+ TVector3 beta,zaxisCS;
+ //
+ // --- Get the muons parameters in the CM frame
+ //
+ PMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
+ PMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
+ //
+ // --- Obtain the dimuon parameters in the CM frame
+ //
+ PDimuCM=PMu1CM+PMu2CM;
+ //
+ // --- Translate the muon parameters in the dimuon rest frame
+ //
+ beta=(-1./PDimuCM.E())*PDimuCM.Vect();
+ if(beta.Mag()>=1) return 666.;
+ PMu1Dimu=PMu1CM;
+ PMu2Dimu=PMu2CM;
+ PMu1Dimu.Boost(beta);
+ PMu2Dimu.Boost(beta);
+
+ //Debugging part -------------------------------------
+ Double_t debugMu1[4]={0.,0.,0.,0.};
+ Double_t debugMu2[4]={0.,0.,0.,0.};
+ PMu1Dimu.GetXYZT(debugMu1);
+ PMu2Dimu.GetXYZT(debugMu2);
+ if (debugMu1[0]!=debugMu1[0] ||debugMu1[1]!=debugMu1[1] || debugMu1[2]!=debugMu1[2] ||debugMu1[3]!=debugMu1[3]) return 666;
+ if (debugMu2[0]!=debugMu2[0] ||debugMu2[1]!=debugMu2[1] || debugMu2[2]!=debugMu2[2] ||debugMu2[3]!=debugMu2[3]) return 666;
+ //----------------------------------------------------
+
+ // --- Determine the z axis for the calculation of the polarization angle (i.e. the direction of the dimuon in the CM system)
+ TVector3 zaxis;
+ zaxis=(PDimuCM.Vect()).Unit();
+
+ // --- Calculation of the polarization angle (angle between mu+ and the z axis defined above)
+ Double_t cost;
+ if(charge1>0) {cost = zaxis.Dot((PMu1Dimu.Vect()).Unit());}
+ else {cost = zaxis.Dot((PMu2Dimu.Vect()).Unit());}
+ return cost;
+}
+
+//________________________________________________________________________
+Double_t AliAnalysisTaskDimuonCFContainerBuilder::PhiCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
+Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
+{
+// Phi decay angle (mu+) in the Collins-Soper frame
+
+ TLorentzVector PMu1CM, PMu2CM, PProjCM, PTargCM, PDimuCM; // In the CM frame
+ TLorentzVector PMu1Dimu, PMu2Dimu, PProjDimu, PTargDimu; // In the dimuon rest frame
+ TVector3 beta,yaxisCS, xaxisCS, zaxisCS;
+ Double_t mp=0.93827231;
+
+ // --- Fill the Lorentz vector for projectile and target in the CM frame
+ PProjCM.SetPxPyPzE(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp));
+ PTargCM.SetPxPyPzE(0.,0.,fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp));
+
+ // --- Get the muons parameters in the CM frame
+ PMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
+ PMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
+
+ // --- Obtain the dimuon parameters in the CM frame
+ PDimuCM=PMu1CM+PMu2CM;
+
+ // --- Translate the dimuon parameters in the dimuon rest frame
+ beta=(-1./PDimuCM.E())*PDimuCM.Vect();
+ if(beta.Mag()>=1) return 666.;
+ PMu1Dimu=PMu1CM;
+ PMu2Dimu=PMu2CM;
+ PProjDimu=PProjCM;
+ PTargDimu=PTargCM;
+ PMu1Dimu.Boost(beta);
+ PMu2Dimu.Boost(beta);
+ PProjDimu.Boost(beta);
+ PTargDimu.Boost(beta);
+
+ //Debugging part -------------------------------------
+ Double_t debugProj[4]={0.,0.,0.,0.};
+ Double_t debugTarg[4]={0.,0.,0.,0.};
+ Double_t debugMu1[4]={0.,0.,0.,0.};
+ Double_t debugMu2[4]={0.,0.,0.,0.};
+ PMu1Dimu.GetXYZT(debugMu1);
+ PMu2Dimu.GetXYZT(debugMu2);
+ PProjDimu.GetXYZT(debugProj);
+ PTargDimu.GetXYZT(debugTarg);
+ if (debugProj[0]!=debugProj[0] ||debugProj[1]!=debugProj[1] || debugProj[2]!=debugProj[2] ||debugProj[3]!=debugProj[3]) return 666;
+ if (debugTarg[0]!=debugTarg[0] ||debugTarg[1]!=debugTarg[1] || debugTarg[2]!=debugTarg[2] ||debugTarg[3]!=debugTarg[3]) return 666;
+ if (debugMu1[0]!=debugMu1[0] ||debugMu1[1]!=debugMu1[1] || debugMu1[2]!=debugMu1[2] ||debugMu1[3]!=debugMu1[3]) return 666;
+ if (debugMu2[0]!=debugMu2[0] ||debugMu2[1]!=debugMu2[1] || debugMu2[2]!=debugMu2[2] ||debugMu2[3]!=debugMu2[3]) return 666;
+ //----------------------------------------------------
+
+ // --- Determine the z axis for the CS angle
+ zaxisCS=(((PProjDimu.Vect()).Unit())-((PTargDimu.Vect()).Unit())).Unit();
+ yaxisCS=(((PProjDimu.Vect()).Unit()).Cross((PTargDimu.Vect()).Unit())).Unit();
+ xaxisCS=(yaxisCS.Cross(zaxisCS)).Unit();
+
+ Double_t phi=0.;
+ if(charge1>0) {
+ phi = TMath::ATan2((PMu1Dimu.Vect()).Dot(yaxisCS),((PMu1Dimu.Vect()).Dot(xaxisCS)));
+ } else {
+ phi = TMath::ATan2((PMu2Dimu.Vect()).Dot(yaxisCS),((PMu2Dimu.Vect()).Dot(xaxisCS)));
+ }
+ if (phi>TMath::Pi()) phi=phi-TMath::Pi();
+
+ return phi;
+}
+
+//________________________________________________________________________
+Double_t AliAnalysisTaskDimuonCFContainerBuilder::PhiHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
+Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
+{
+// Phi decay angle (mu+) in the Helicity frame
+ TLorentzVector PMu1Lab, PMu2Lab, PProjLab, PTargLab, PDimuLab; // In the lab. frame
+ TLorentzVector PMu1Dimu, PMu2Dimu, PProjDimu, PTargDimu; // In the dimuon rest frame
+ TVector3 beta,xaxis, yaxis,zaxis;
+ Double_t mp=0.93827231;
+
+ // --- Get the muons parameters in the LAB frame
+ PMu1Lab.SetPxPyPzE(px1,py1,pz1,e1);
+ PMu2Lab.SetPxPyPzE(px2,py2,pz2,e2);
+
+ // --- Obtain the dimuon parameters in the LAB frame
+ PDimuLab=PMu1Lab+PMu2Lab;
+ zaxis=(PDimuLab.Vect()).Unit();
+
+ // --- Translate the muon parameters in the dimuon rest frame
+ beta=(-1./PDimuLab.E())*PDimuLab.Vect();
+ if(beta.Mag()>=1.) return 666.;
+
+ PProjLab.SetPxPyPzE(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp)); // proiettile
+ PTargLab.SetPxPyPzE(0.,0.,fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp)); // bersaglio
+
+ PProjDimu=PProjLab;
+ PTargDimu=PTargLab;
+
+ PProjDimu.Boost(beta);
+ PTargDimu.Boost(beta);
+
+ yaxis=((PProjDimu.Vect()).Cross(PTargDimu.Vect())).Unit();
+ xaxis=(yaxis.Cross(zaxis)).Unit();
+
+ PMu1Dimu=PMu1Lab;
+ PMu2Dimu=PMu2Lab;
+ PMu1Dimu.Boost(beta);
+ PMu2Dimu.Boost(beta);
+
+ //Debugging part -------------------------------------
+ Double_t debugProj[4]={0.,0.,0.,0.};
+ Double_t debugTarg[4]={0.,0.,0.,0.};
+ Double_t debugMu1[4]={0.,0.,0.,0.};
+ Double_t debugMu2[4]={0.,0.,0.,0.};
+ PMu1Dimu.GetXYZT(debugMu1);
+ PMu2Dimu.GetXYZT(debugMu2);
+ PProjDimu.GetXYZT(debugProj);
+ PTargDimu.GetXYZT(debugTarg);
+ if (debugProj[0]!=debugProj[0] ||debugProj[1]!=debugProj[1] || debugProj[2]!=debugProj[2] ||debugProj[3]!=debugProj[3]) return 666;
+ if (debugTarg[0]!=debugTarg[0] ||debugTarg[1]!=debugTarg[1] || debugTarg[2]!=debugTarg[2] ||debugTarg[3]!=debugTarg[3]) return 666;
+ if (debugMu1[0]!=debugMu1[0] ||debugMu1[1]!=debugMu1[1] || debugMu1[2]!=debugMu1[2] ||debugMu1[3]!=debugMu1[3]) return 666;
+ if (debugMu2[0]!=debugMu2[0] ||debugMu2[1]!=debugMu2[1] || debugMu2[2]!=debugMu2[2] ||debugMu2[3]!=debugMu2[3]) return 666;
+ //----------------------------------------------------
+
+ Double_t phi=0.;
+ if(charge1 > 0) {
+ phi = TMath::ATan2((PMu1Dimu.Vect()).Dot(yaxis),(PMu1Dimu.Vect()).Dot(xaxis));
+ } else {
+ phi = TMath::ATan2((PMu2Dimu.Vect()).Dot(yaxis),(PMu2Dimu.Vect()).Dot(xaxis));
+ }
+ return phi;
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskDimuonCFContainerBuilder::Terminate(Option_t *)
+{
+
+}
+
+#endif
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff