PWG3/muon/AliAnalysisTaskMuonTreeBuilder.cxx

   1 #ifndef ALIANALYSISTASKMUONTREEBUILDER_CXX
   2 #define ALIANALYSISTASKMUONTREEBUILDER_CXX
   3
   4 #include "AliAnalysisTaskMuonTreeBuilder.h"
   5 #include "AliMCEvent.h"
   6 #include "AliESDMuonTrack.h"
   7 #include "AliESDVertex.h"
   8 #include "AliStack.h"
   9 #include "AliHeader.h"
  10 #include "AliESDHeader.h"
  11 #include "TParticle.h"
  12 #include "TLorentzVector.h"
  13 #include "TFile.h"
  14 #include "TH1I.h"
  15 #include "AliAnalysisManager.h"
  16 #include "AliESDEvent.h"
  17 #include "AliAODEvent.h"
  18 #include "TChain.h"
  19 #include "AliESDtrack.h"
  20 #include "AliLog.h"
  21 #include "AliESDtrack.h"
  22 #include "AliESDInputHandler.h"
  23 #include "TCanvas.h"
  24 #include "AliPhysicsSelection.h"
  25 #include "AliMultiplicity.h"
  26
  27 //      Analysis task for muon-dimuon analysis
  28 //      Works for real and MC events
  29 //      Works with the corresponding AddTask macro
  30 //      Includes a flag for physics selection
  31 //
  32 //      author: L. Bianchi - Universita' & INFN Torino
  33
  34 ClassImp(AliAnalysisTaskMuonTreeBuilder)
  35
  36 //__________________________________________________________________________
  37 AliAnalysisTaskMuonTreeBuilder::AliAnalysisTaskMuonTreeBuilder() :
  38   fNevt(0),
  39   fBeamEnergy(5000.),
  40   fOutput(0x0),
  41   fOutputTree(0x0),
  42   fIsMC(kFALSE),
  43   fIsSelected(kFALSE),
  44   fNumMuonTracks(0),
  45   fNumSPDTracklets(0),
  46   fNumContributors(0),
  47   fNumDimuons(0)
  48 {
  49   //
  50   //Default ctor
  51   //
  52 }
  53 //___________________________________________________________________________
  54 AliAnalysisTaskMuonTreeBuilder::AliAnalysisTaskMuonTreeBuilder(const Char_t* name) :
  55   AliAnalysisTaskSE(name),
  56   fNevt(0),
  57   fBeamEnergy(5000.),
  58   fOutput(0x0),
  59   fOutputTree(0x0),
  60   fIsMC(kFALSE),
  61   fIsSelected(kFALSE),
  62   fNumMuonTracks(0),
  63   fNumSPDTracklets(0),
  64   fNumContributors(0),
  65   fNumDimuons(0)
  66 {
  67   //
  68   // Constructor. Initialization of Inputs and Outputs
  69   //
  70   Info("AliAnalysisTaskMuonTreeBuilder","Calling Constructor");
  71
  72   DefineOutput(1,TTree::Class());
  73 }
  74
  75 //___________________________________________________________________________
  76 AliAnalysisTaskMuonTreeBuilder& AliAnalysisTaskMuonTreeBuilder::operator=(const AliAnalysisTaskMuonTreeBuilder& c)
  77 {
  78   //
  79   // Assignment operator
  80   //
  81   if (this!=&c) {
  82     AliAnalysisTaskSE::operator=(c) ;
  83     fNevt = c.fNevt ;
  84   }
  85   return *this;
  86 }
  87
  88 //___________________________________________________________________________
  89 AliAnalysisTaskMuonTreeBuilder::AliAnalysisTaskMuonTreeBuilder(const AliAnalysisTaskMuonTreeBuilder& c) :
  90   AliAnalysisTaskSE(c),
  91   fNevt(c.fNevt),
  92   fBeamEnergy(c.fBeamEnergy),
  93   fOutput(c.fOutput),
  94   fOutputTree(c.fOutputTree),
  95   fIsMC(kFALSE),
  96   fIsSelected(kFALSE),
  97   fNumMuonTracks(c.fNumMuonTracks),
  98   fNumSPDTracklets(c.fNumSPDTracklets),
  99   fNumContributors(c.fNumContributors),
 100   fNumDimuons(c.fNumDimuons)
 101 {
 102   //
 103   // Copy Constructor                                                                           FIDUCIAL REGIONS?
 104   //
 105 }
 106
 107 //___________________________________________________________________________
 108 AliAnalysisTaskMuonTreeBuilder::~AliAnalysisTaskMuonTreeBuilder() {
 109   //
 110   //destructor
 111   //
 112   Info("~AliAnalysisTaskMuonTreeBuilder","Calling Destructor");
 113 }
 114
 115 //___________________________________________________________________________
 116 void AliAnalysisTaskMuonTreeBuilder::UserCreateOutputObjects(){
 117
 118 //
 119 //  Creating User-Defined Output Objects
 120 //
 121
 122  // TREE OUTPUT----------------------------------------------------------
 123  OpenFile(1);
 124  fOutputTree = new TTree("krec","Tree of reconstructed muons");
 125
 126  fOutputTree->Branch("IsSelected",&fIsSelected,"IsSelected/O");
 127  fOutputTree->Branch("FiredTriggerClasses",fTrigClass,"FiredTriggerClasses/C");
 128
 129  fOutputTree->Branch("NumMuonTracks",&fNumMuonTracks,"NumMuonTracks/I");
 130  fOutputTree->Branch("NumContributors",&fNumContributors,"NumContributors/I");
 131  fOutputTree->Branch("NumSPDTracklets",&fNumSPDTracklets,"NumSPDTraclets/I");
 132  fOutputTree->Branch("Vertex",fVertex,"Vertex[3]/D");
 133  fOutputTree->Branch("pT",fpT,"pT[10]/D");
 134  fOutputTree->Branch("E",fE,"E[10]/D");
 135  fOutputTree->Branch("px",fpx,"px[10]/D");
 136  fOutputTree->Branch("py",fpy,"py[10]/D");
 137  fOutputTree->Branch("pz",fpz,"pz[10]/D");
 138  fOutputTree->Branch("pxUncorr",fpxUncorr,"pxUncorr[10]/D");
 139  fOutputTree->Branch("pyUncorr",fpyUncorr,"pyUncorr[10]/D");
 140  fOutputTree->Branch("pzUncorr",fpzUncorr,"pzUncorr[10]/D");
 141  fOutputTree->Branch("y",fy,"y[10]/D");
 142  fOutputTree->Branch("eta",feta,"eta[10]/D");
 143  fOutputTree->Branch("phi",fphi,"phi[10]/D");
 144  fOutputTree->Branch("MatchTrig",fMatchTrig,"MatchTrig[10]/I");
 145  fOutputTree->Branch("TrackChi2",fTrackChi2,"TrackChi2[10]/D");
 146  fOutputTree->Branch("MatchTrigChi2",fMatchTrigChi2,"MatchTrigChi2[10]/D");
 147  fOutputTree->Branch("DCA",fDCA,"DCA[10]/D");
 148  fOutputTree->Branch("Charge",fCharge,"Charge[10]/S");
 149  fOutputTree->Branch("MuFamily",fMuFamily,"MuFamily[10]/I");
 150  fOutputTree->Branch("RAtAbsEnd",fRAtAbsEnd,"RAtAbsEnd[10]/D");
 151
 152  fOutputTree->Branch("NumDimuons",&fNumDimuons,"NumDimuons/I");
 153  fOutputTree->Branch("DimuConstituent",fDimuonConstituent,"DimuonConstituent[45][2]/I");
 154  fOutputTree->Branch("pTdimuon",fpTdimuon,"pTdimuon[45]/D");
 155  fOutputTree->Branch("pxdimuon",fpxdimuon,"pxdimuon[45]/D");
 156  fOutputTree->Branch("pydimuon",fpydimuon,"pydimuon[45]/D");
 157  fOutputTree->Branch("pzdimuon",fpzdimuon,"pzdimuon[45]/D");
 158  fOutputTree->Branch("ydimuon",fydimuon,"ydimuon[45]/D");
 159  fOutputTree->Branch("Imassdimuon",fiMassdimuon,"iMassdimuon[45]/D");
 160  fOutputTree->Branch("costCS",fcostCS,"costCS[45]/D");
 161  fOutputTree->Branch("costHE",fcostHE,"costHE[45]/D");
 162  fOutputTree->Branch("phiCS",fphiCS,"phiCS[45]/D");
 163  fOutputTree->Branch("phiHE",fphiHE,"phiHE[45]/D");
 164
 165  fOutputTree->Branch("PDG",fPDG,"PDG[10]/I");
 166  fOutputTree->Branch("PDGmother",fPDGmother,"PDGmother[10]/I");
 167  fOutputTree->Branch("PDGdimu",fPDGdimu,"PDGdimu[45]/I");
 168
 169 }
 170
 171
 172
 173 //_________________________________________________
 174 void AliAnalysisTaskMuonTreeBuilder::UserExec(Option_t *)
 175 {
 176 //
 177 //  User Exec
 178 //
 179
 180   fNevt++;
 181
 182
 183   fNumMuonTracks=0;
 184   fNumSPDTracklets=666.;
 185   fNumContributors=666.;
 186   fNumDimuons=0;
 187   fIsSelected=kFALSE;
 188   fVertex[0]=666.; fVertex[1]=666.; fVertex[2]=666.;
 189   for(Int_t i=0; i<10;i++){
 190     fpT[i]=666.;
 191     fE[i]=666.;
 192     fpx[i]=666;
 193     fpy[i]=666;
 194     fpz[i]=666;
 195     fpxUncorr[i]=666;
 196     fpyUncorr[i]=666;
 197     fpzUncorr[i]=666;
 198     fy[i]=666.;
 199     feta[i]=666.;
 200     fphi[i]=666.;
 201     fMatchTrig[i]=666.;
 202     fTrackChi2[i]=666.;
 203     fMatchTrigChi2[i]=666.;
 204     fDCA[i]=666.;
 205     fPDG[i]=666;
 206     fPDGmother[i]=666;
 207     fCharge[i]=666;
 208     fMuFamily[i]=666;
 209     fRAtAbsEnd[i]=666;
 210   }
 211   for(Int_t i=0; i<45;i++){
 212     fDimuonConstituent[i][0]=666;  fDimuonConstituent[i][1]=666;
 213     fpTdimuon[i]=666.;
 214     fpxdimuon[i]=666.;
 215     fpydimuon[i]=666.;
 216     fpzdimuon[i]=666.;
 217     fydimuon[i]=666.;
 218     fiMassdimuon[i]=666.;
 219     fcostCS[i]=666.;
 220     fcostHE[i]=666.;
 221     fphiCS[i]=666.;
 222     fphiHE[i]=666.;
 223     fPDGdimu[i]=666;
 224   }
 225
 226
 227 ////////
 228 //// ESD
 229 ////////
 230
 231   AliESDEvent *fESD = 0x0;
 232   AliMCEvent*  mcEvent  = 0x0;
 233
 234   if(fIsMC){
 235     if (!fMCEvent) {
 236       Error("UserExec","NO MC EVENT FOUND!");
 237       return;
 238     }
 239   }
 240
 241   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 242
 243   fIsSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 244
 245   fNumMuonTracks = fESD->GetNumberOfMuonTracks() ;
 246    Int_t loopEnd = fNumMuonTracks;
 247   if(!fIsMC) {
 248     TString cla = fESD->GetFiredTriggerClasses();
 249     sprintf(fTrigClass,"%s",cla.Data());
 250   }
 251
 252   if(fNumMuonTracks>0 && fIsMC){
 253     mcEvent = MCEvent();
 254   }
 255   fNumContributors = fESD->GetPrimaryVertexSPD()->GetNContributors();
 256   fNumSPDTracklets = fESD->GetMultiplicity()->GetNumberOfTracklets();
 257   fVertex[0]=fESD->GetPrimaryVertexSPD()->GetX();
 258   fVertex[1]=fESD->GetPrimaryVertexSPD()->GetY();
 259   fVertex[2]=fESD->GetPrimaryVertexSPD()->GetZ();
 260   printf("fVertex : %f - %f - %f\n",fVertex[0],fVertex[1],fVertex[2]);
 261
 262   Int_t numdimu = 0;
 263   for (Int_t j = 0; j<loopEnd; j++) {
 264     AliESDMuonTrack* mu1 = new AliESDMuonTrack(*(fESD->GetMuonTrack(j)));
 265     if (!mu1->ContainTrackerData()) {fNumMuonTracks=fNumMuonTracks-1; continue;}
 266     Double_t charge1 = mu1->Charge();
 267     fCharge[j] = mu1->Charge();
 268     fpT[j]   = mu1->Pt();
 269     fpx[j]  = mu1->Px();
 270     fpy[j]  = mu1->Py();
 271     fpz[j]  = mu1->Pz();
 272     fpxUncorr[j]  = mu1->PxUncorrected();
 273     fpyUncorr[j]  = mu1->PyUncorrected();
 274     fpzUncorr[j]  = mu1->PzUncorrected();
 275     fy[j]    = mu1->Y();
 276     feta[j]    = mu1->Eta();
 277     fphi[j]  = Phideg(mu1->Phi());
 278     Double_t emu1    = mu1->E();
 279     fE[j] = emu1;
 280     fDCA[j] = mu1->GetDCA();
 281     fTrackChi2[j] = mu1->GetChi2();
 282     fMatchTrig[j]   = mu1->GetMatchTrigger();
 283     fMatchTrigChi2[j]= mu1->GetChi2MatchTrigger();
 284     fRAtAbsEnd[j]=mu1->GetRAtAbsorberEnd();
 285
 286     AliMCParticle* mcTrack = 0x0;
 287     if(fIsMC){
 288     if(mu1->GetLabel()==-1) continue;
 289     mcTrack = (AliMCParticle*)mcEvent->GetTrack(mu1->GetLabel());
 290     fPDG[j] = mcTrack->PdgCode();
 291     if (mcTrack->GetMother()==-1) continue;
 292     fPDGmother[j] = ((AliMCParticle*)mcEvent->GetTrack(mcTrack->GetMother()))->PdgCode();
 293     if (TMath::Abs(fPDG[j])==13) fMuFamily[j] = FindMuFamily(mcTrack,mcEvent);
 294     }
 295     for (Int_t jj = j+1; jj<loopEnd; jj++) {
 296       AliESDMuonTrack* mu2 = new AliESDMuonTrack(*(fESD->GetMuonTrack(jj)));
 297       if (!mu2->ContainTrackerData()) continue;
 298
 299       Double_t pxmu2  = mu2->Px();
 300       Double_t pymu2  = mu2->Py();
 301       Double_t pzmu2  = mu2->Pz();
 302       Double_t emu2   = mu2->E();
 303       //Double_t charge2= mu2->Charge();
 304
 305       fpTdimuon[numdimu] = TMath::Sqrt((fpx[j]+pxmu2)*(fpx[j]+pxmu2)+(fpy[j]+pymu2)*(fpy[j]+pymu2));
 306       fpxdimuon[numdimu] = fpx[j]+pxmu2;
 307       fpydimuon[numdimu] = fpy[j]+pymu2;
 308       fpzdimuon[numdimu] = fpz[j]+pzmu2;
 309       fydimuon[numdimu] = Rap((emu1+emu2),(fpz[j]+pzmu2));
 310       fiMassdimuon[numdimu] = Imass(emu1,fpx[j],fpy[j],fpz[j],emu2,pxmu2,pymu2,pzmu2);
 311       fcostCS[numdimu]=CostCS(fpx[j],fpy[j],fpz[j],emu1,charge1,pxmu2,pymu2,pzmu2,emu2);
 312       fcostHE[numdimu]=CostHE(fpx[j],fpy[j],fpz[j],emu1,charge1,pxmu2,pymu2,pzmu2,emu2);
 313       fphiCS[numdimu] = PhiCS(fpx[j],fpy[j],fpz[j],emu1,charge1,pxmu2,pymu2,pzmu2,emu2);
 314       fphiHE[numdimu] = PhiHE(fpx[j],fpy[j],fpz[j],emu1,charge1,pxmu2,pymu2,pzmu2,emu2);
 315
 316       fDimuonConstituent[numdimu][0]=j;  fDimuonConstituent[numdimu][1]=jj;
 317
 318       numdimu++;
 319
 320       if(fIsMC){
 321       if(mu2->GetLabel()==-1) continue;
 322       if(TMath::Abs(fPDG[j])==13 && TMath::Abs(((AliMCParticle*)mcEvent->GetTrack(mu2->GetLabel()))->PdgCode())==13) fPDGdimu[numdimu-1]=FindDimuFamily(mcTrack,(AliMCParticle*)mcEvent->GetTrack(mu2->GetLabel()),mcEvent);
 323       else fPDGdimu[numdimu-1]=-3;
 324       }
 325
 326       delete mu2;
 327     }
 328     fNumDimuons=numdimu;
 329
 330
 331     delete mu1;
 332     }
 333   fOutputTree->Fill();
 334   PostData(1,fOutputTree);
 335 }
 336
 337 //________________________________________________________________________
 338 Int_t AliAnalysisTaskMuonTreeBuilder::FindDimuFamily(AliMCParticle* mcTrack1,AliMCParticle* mcTrack2, AliMCEvent* mcEvent) const
 339 {
 340   // finds the family of the dimuon (works only if the 2 muons are real muons (not hadrons))
 341   Int_t familynumber;
 342
 343   if(mcTrack1->GetMother()==mcTrack2->GetMother()) familynumber = TMath::Abs(((AliMCParticle*)mcEvent->GetTrack(mcTrack1->GetMother()))->PdgCode());
 344   else{
 345     Int_t familymu1 = FindMuFamily(mcTrack1,mcEvent);
 346     Int_t familymu2 = FindMuFamily(mcTrack2,mcEvent);
 347     if(familymu1==5 && familymu2==5) familynumber=5;                                                    //bb dimuon
 348     else if(familymu1==4 && familymu2==4) familynumber=4;                                               //cc dimuon
 349     else if((familymu1==4 && familymu2==5)||(familymu2==4 && familymu1==5)) familynumber=45;            //bc dimuon
 350     else if (familymu1==-2 || familymu2==-2 || familymu1==-3 || familymu2==-3) familynumber=-2;         //hadron dimuon (at least 1 hadron involved)
 351     else familynumber=-1;
 352   }
 353
 354   return familynumber;
 355 }
 356
 357 //________________________________________________________________________
 358 Int_t AliAnalysisTaskMuonTreeBuilder::FindMuFamily(AliMCParticle* mcTrack, AliMCEvent* mcEvent) const
 359 {
 360   // finds the family of the muon
 361   Int_t imother = mcTrack->GetMother();
 362   if ( imother<0 ) return -1; // Drell-Yan Muon
 363
 364   Int_t igrandma = imother;
 365
 366   AliMCParticle* motherPart = (AliMCParticle*)mcEvent->GetTrack(imother);
 367   Int_t motherPdg = motherPart->PdgCode();
 368
 369   // Track is an heavy flavor muon
 370   Int_t absPdg = TMath::Abs(motherPdg);
 371   if(absPdg/100==5 || absPdg/1000==5) return 5;
 372   if(absPdg/100==4 || absPdg/1000==4){
 373     Int_t newMother = -1;
 374     igrandma = imother;
 375     Int_t absGrandMotherPdg = TMath::Abs(motherPart->PdgCode());
 376     while ( absGrandMotherPdg > 10 ) {
 377       igrandma = ((AliMCParticle*)mcEvent->GetTrack(igrandma))->GetMother();
 378       if( igrandma < 0 ) break;
 379       absGrandMotherPdg = TMath::Abs(((AliMCParticle*)mcEvent->GetTrack(igrandma))->PdgCode());
 380     }
 381
 382     if (absGrandMotherPdg==5) newMother = 5; // Charm from beauty
 383     else if (absGrandMotherPdg==4) newMother = 4;
 384
 385     if(newMother<0) {
 386       //AliWarning("Mother not correctly found! Set to charm!\n");
 387       newMother = 4;
 388     }
 389
 390     return newMother;
 391   }
 392
 393   Int_t nPrimaries = mcEvent->Stack()->GetNprimary();
 394   // Track is a bkg. muon
 395   if (imother<nPrimaries) {
 396     return -2;                  //is a primary
 397   }
 398   else {
 399     return -3;                  //is a secondary
 400   }
 401
 402 }
 403
 404 //________________________________________________________________________
 405 Double_t AliAnalysisTaskMuonTreeBuilder::Imass(Double_t e1, Double_t px1, Double_t py1, Double_t pz1,
 406                                    Double_t e2, Double_t px2, Double_t py2, Double_t pz2) const
 407 {
 408 // invariant mass calculation
 409     Double_t imassrec = TMath::Sqrt((e1+e2)*(e1+e2)-((px1+px2)*(px1+px2)+
 410                                     (py1+py2)*(py1+py2)+(pz1+pz2)*(pz1+pz2)));
 411     return imassrec;
 412 }
 413
 414 //________________________________________________________________________
 415 Double_t AliAnalysisTaskMuonTreeBuilder::Rap(Double_t e, Double_t pz) const
 416 {
 417 // calculate rapidity
 418     Double_t rap;
 419     if(e>TMath::Abs(pz)){
 420         rap = 0.5*TMath::Log((e+pz)/(e-pz));
 421         return rap;
 422     }
 423     else{
 424         rap = 666.;
 425         return rap;
 426     }
 427 }
 428
 429 //________________________________________________________________________
 430 Double_t AliAnalysisTaskMuonTreeBuilder::Phideg(Double_t phi) const
 431 {
 432 // calculate Phi in range [-180,180]
 433     Double_t phideg;
 434
 435         phideg = phi-TMath::Pi();
 436         phideg = phideg*57.296;
 437         return phideg;
 438 }
 439
 440 //________________________________________________________________________
 441 Double_t AliAnalysisTaskMuonTreeBuilder::CostCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
 442 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
 443 {
 444 // Cosine of the theta decay angle (mu+) in the Collins-Soper frame
 445
 446   TLorentzVector pMu1CM, pMu2CM, pProjCM, pTargCM, pDimuCM; // In the CM. frame
 447   TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame
 448   TVector3 beta,zaxisCS;
 449   Double_t mp=0.93827231;
 450   //
 451   // --- Fill the Lorentz vector for projectile and target in the CM frame
 452   //
 453   pProjCM.SetPxPyPzE(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp));
 454   pTargCM.SetPxPyPzE(0.,0.,fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp));
 455   //
 456   // --- Get the muons parameters in the CM frame
 457   //
 458   pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
 459   pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
 460   //
 461   // --- Obtain the dimuon parameters in the CM frame
 462   //
 463   pDimuCM=pMu1CM+pMu2CM;
 464   //
 465   // --- Translate the dimuon parameters in the dimuon rest frame
 466   //
 467   beta=(-1./pDimuCM.E())*pDimuCM.Vect();
 468   if(beta.Mag()>=1) return 666.;
 469   pMu1Dimu=pMu1CM;
 470   pMu2Dimu=pMu2CM;
 471   pProjDimu=pProjCM;
 472   pTargDimu=pTargCM;
 473   pMu1Dimu.Boost(beta);
 474   pMu2Dimu.Boost(beta);
 475   pProjDimu.Boost(beta);
 476   pTargDimu.Boost(beta);
 477
 478   //Debugging part -------------------------------------
 479   Double_t debugProj[4]={0.,0.,0.,0.};
 480   Double_t debugTarg[4]={0.,0.,0.,0.};
 481   Double_t debugMu1[4]={0.,0.,0.,0.};
 482   Double_t debugMu2[4]={0.,0.,0.,0.};
 483   pMu1Dimu.GetXYZT(debugMu1);
 484   pMu2Dimu.GetXYZT(debugMu2);
 485   pProjDimu.GetXYZT(debugProj);
 486   pTargDimu.GetXYZT(debugTarg);
 487   if (debugProj[0]!=debugProj[0] ||debugProj[1]!=debugProj[1] || debugProj[2]!=debugProj[2] ||debugProj[3]!=debugProj[3]) return 666;
 488   if (debugTarg[0]!=debugTarg[0] ||debugTarg[1]!=debugTarg[1] || debugTarg[2]!=debugTarg[2] ||debugTarg[3]!=debugTarg[3]) return 666;
 489   if (debugMu1[0]!=debugMu1[0] ||debugMu1[1]!=debugMu1[1] || debugMu1[2]!=debugMu1[2] ||debugMu1[3]!=debugMu1[3]) return 666;
 490   if (debugMu2[0]!=debugMu2[0] ||debugMu2[1]!=debugMu2[1] || debugMu2[2]!=debugMu2[2] ||debugMu2[3]!=debugMu2[3]) return 666;
 491   //----------------------------------------------------
 492
 493   // --- Determine the z axis for the CS angle
 494   zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit();
 495
 496   // --- Determine the CS angle (angle between mu+ and the z axis defined above)
 497   Double_t cost;
 498
 499   if(charge1>0) {cost = zaxisCS.Dot((pMu1Dimu.Vect()).Unit());}
 500   else {cost = zaxisCS.Dot((pMu2Dimu.Vect()).Unit());}
 501
 502   return cost;
 503 }
 504
 505 //________________________________________________________________________
 506 Double_t AliAnalysisTaskMuonTreeBuilder::CostHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
 507 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
 508 {
 509 // Cosine of the theta decay angle (mu+) in the Helicity frame
 510
 511   TLorentzVector pMu1CM, pMu2CM, pDimuCM; // In the CM frame
 512   TLorentzVector pMu1Dimu, pMu2Dimu; // In the dimuon rest frame
 513   TVector3 beta,zaxisCS;
 514   //
 515   // --- Get the muons parameters in the CM frame
 516   //
 517   pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
 518   pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
 519   //
 520   // --- Obtain the dimuon parameters in the CM frame
 521   //
 522   pDimuCM=pMu1CM+pMu2CM;
 523   //
 524   // --- Translate the muon parameters in the dimuon rest frame
 525   //
 526   beta=(-1./pDimuCM.E())*pDimuCM.Vect();
 527   if(beta.Mag()>=1) return 666.;
 528   pMu1Dimu=pMu1CM;
 529   pMu2Dimu=pMu2CM;
 530   pMu1Dimu.Boost(beta);
 531   pMu2Dimu.Boost(beta);
 532
 533   //Debugging part -------------------------------------
 534   Double_t debugMu1[4]={0.,0.,0.,0.};
 535   Double_t debugMu2[4]={0.,0.,0.,0.};
 536   pMu1Dimu.GetXYZT(debugMu1);
 537   pMu2Dimu.GetXYZT(debugMu2);
 538   if (debugMu1[0]!=debugMu1[0] ||debugMu1[1]!=debugMu1[1] || debugMu1[2]!=debugMu1[2] ||debugMu1[3]!=debugMu1[3]) return 666;
 539   if (debugMu2[0]!=debugMu2[0] ||debugMu2[1]!=debugMu2[1] || debugMu2[2]!=debugMu2[2] ||debugMu2[3]!=debugMu2[3]) return 666;
 540   //----------------------------------------------------
 541
 542   // --- Determine the z axis for the calculation of the polarization angle (i.e. the direction of the dimuon in the CM system)
 543   TVector3 zaxis;
 544   zaxis=(pDimuCM.Vect()).Unit();
 545
 546   // --- Calculation of the polarization angle (angle between mu+ and the z axis defined above)
 547   Double_t cost;
 548   if(charge1>0) {cost = zaxis.Dot((pMu1Dimu.Vect()).Unit());}
 549   else {cost = zaxis.Dot((pMu2Dimu.Vect()).Unit());}
 550   return cost;
 551 }
 552
 553 //________________________________________________________________________
 554 Double_t AliAnalysisTaskMuonTreeBuilder::PhiCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
 555 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
 556 {
 557 // Phi decay angle (mu+) in the Collins-Soper frame
 558
 559    TLorentzVector pMu1CM, pMu2CM, pProjCM, pTargCM, pDimuCM; // In the CM frame
 560    TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame
 561    TVector3 beta,yaxisCS, xaxisCS, zaxisCS;
 562    Double_t mp=0.93827231;
 563
 564    // --- Fill the Lorentz vector for projectile and target in the CM frame
 565    pProjCM.SetPxPyPzE(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp));
 566    pTargCM.SetPxPyPzE(0.,0.,fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp));
 567
 568    // --- Get the muons parameters in the CM frame
 569    pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
 570    pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
 571
 572    // --- Obtain the dimuon parameters in the CM frame
 573    pDimuCM=pMu1CM+pMu2CM;
 574
 575    // --- Translate the dimuon parameters in the dimuon rest frame
 576    beta=(-1./pDimuCM.E())*pDimuCM.Vect();
 577    if(beta.Mag()>=1) return 666.;
 578    pMu1Dimu=pMu1CM;
 579    pMu2Dimu=pMu2CM;
 580    pProjDimu=pProjCM;
 581    pTargDimu=pTargCM;
 582    pMu1Dimu.Boost(beta);
 583    pMu2Dimu.Boost(beta);
 584    pProjDimu.Boost(beta);
 585    pTargDimu.Boost(beta);
 586
 587    //Debugging part -------------------------------------
 588    Double_t debugProj[4]={0.,0.,0.,0.};
 589    Double_t debugTarg[4]={0.,0.,0.,0.};
 590    Double_t debugMu1[4]={0.,0.,0.,0.};
 591    Double_t debugMu2[4]={0.,0.,0.,0.};
 592    pMu1Dimu.GetXYZT(debugMu1);
 593    pMu2Dimu.GetXYZT(debugMu2);
 594    pProjDimu.GetXYZT(debugProj);
 595    pTargDimu.GetXYZT(debugTarg);
 596    if (debugProj[0]!=debugProj[0] ||debugProj[1]!=debugProj[1] || debugProj[2]!=debugProj[2] ||debugProj[3]!=debugProj[3]) return 666;
 597    if (debugTarg[0]!=debugTarg[0] ||debugTarg[1]!=debugTarg[1] || debugTarg[2]!=debugTarg[2] ||debugTarg[3]!=debugTarg[3]) return 666;
 598    if (debugMu1[0]!=debugMu1[0] ||debugMu1[1]!=debugMu1[1] || debugMu1[2]!=debugMu1[2] ||debugMu1[3]!=debugMu1[3]) return 666;
 599    if (debugMu2[0]!=debugMu2[0] ||debugMu2[1]!=debugMu2[1] || debugMu2[2]!=debugMu2[2] ||debugMu2[3]!=debugMu2[3]) return 666;
 600    //----------------------------------------------------
 601
 602    // --- Determine the z axis for the CS angle
 603    zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit();
 604    yaxisCS=(((pProjDimu.Vect()).Unit()).Cross((pTargDimu.Vect()).Unit())).Unit();
 605    xaxisCS=(yaxisCS.Cross(zaxisCS)).Unit();
 606
 607    Double_t phi=0.;
 608    if(charge1>0) {
 609        phi = TMath::ATan2((pMu1Dimu.Vect()).Dot(yaxisCS),((pMu1Dimu.Vect()).Dot(xaxisCS)));
 610    } else {
 611        phi = TMath::ATan2((pMu2Dimu.Vect()).Dot(yaxisCS),((pMu2Dimu.Vect()).Dot(xaxisCS)));
 612    }
 613    if (phi>TMath::Pi()) phi=phi-TMath::Pi();
 614
 615    return phi;
 616 }
 617
 618 //________________________________________________________________________
 619 Double_t AliAnalysisTaskMuonTreeBuilder::PhiHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
 620 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
 621 {
 622 // Phi decay angle (mu+) in the Helicity frame
 623   TLorentzVector pMu1Lab, pMu2Lab, pProjLab, pTargLab, pDimuLab; // In the lab. frame
 624   TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame
 625   TVector3 beta,xaxis, yaxis,zaxis;
 626   Double_t mp=0.93827231;
 627
 628   // --- Get the muons parameters in the LAB frame
 629   pMu1Lab.SetPxPyPzE(px1,py1,pz1,e1);
 630   pMu2Lab.SetPxPyPzE(px2,py2,pz2,e2);
 631
 632   // --- Obtain the dimuon parameters in the LAB frame
 633   pDimuLab=pMu1Lab+pMu2Lab;
 634   zaxis=(pDimuLab.Vect()).Unit();
 635
 636   // --- Translate the muon parameters in the dimuon rest frame
 637   beta=(-1./pDimuLab.E())*pDimuLab.Vect();
 638   if(beta.Mag()>=1.) return 666.;
 639
 640   pProjLab.SetPxPyPzE(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp)); // proiettile
 641   pTargLab.SetPxPyPzE(0.,0.,fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp)); // bersaglio
 642
 643   pProjDimu=pProjLab;
 644   pTargDimu=pTargLab;
 645
 646   pProjDimu.Boost(beta);
 647   pTargDimu.Boost(beta);
 648
 649   yaxis=((pProjDimu.Vect()).Cross(pTargDimu.Vect())).Unit();
 650   xaxis=(yaxis.Cross(zaxis)).Unit();
 651
 652   pMu1Dimu=pMu1Lab;
 653   pMu2Dimu=pMu2Lab;
 654   pMu1Dimu.Boost(beta);
 655   pMu2Dimu.Boost(beta);
 656
 657   //Debugging part -------------------------------------
 658   Double_t debugProj[4]={0.,0.,0.,0.};
 659   Double_t debugTarg[4]={0.,0.,0.,0.};
 660   Double_t debugMu1[4]={0.,0.,0.,0.};
 661   Double_t debugMu2[4]={0.,0.,0.,0.};
 662   pMu1Dimu.GetXYZT(debugMu1);
 663   pMu2Dimu.GetXYZT(debugMu2);
 664   pProjDimu.GetXYZT(debugProj);
 665   pTargDimu.GetXYZT(debugTarg);
 666   if (debugProj[0]!=debugProj[0] ||debugProj[1]!=debugProj[1] || debugProj[2]!=debugProj[2] ||debugProj[3]!=debugProj[3]) return 666;
 667   if (debugTarg[0]!=debugTarg[0] ||debugTarg[1]!=debugTarg[1] || debugTarg[2]!=debugTarg[2] ||debugTarg[3]!=debugTarg[3]) return 666;
 668   if (debugMu1[0]!=debugMu1[0] ||debugMu1[1]!=debugMu1[1] || debugMu1[2]!=debugMu1[2] ||debugMu1[3]!=debugMu1[3]) return 666;
 669   if (debugMu2[0]!=debugMu2[0] ||debugMu2[1]!=debugMu2[1] || debugMu2[2]!=debugMu2[2] ||debugMu2[3]!=debugMu2[3]) return 666;
 670   //----------------------------------------------------
 671
 672   Double_t phi=0.;
 673    if(charge1 > 0) {
 674       phi = TMath::ATan2((pMu1Dimu.Vect()).Dot(yaxis),(pMu1Dimu.Vect()).Dot(xaxis));
 675      } else {
 676       phi = TMath::ATan2((pMu2Dimu.Vect()).Dot(yaxis),(pMu2Dimu.Vect()).Dot(xaxis));
 677    }
 678    return phi;
 679 }
 680
 681 //________________________________________________________________________
 682 void AliAnalysisTaskMuonTreeBuilder::Terminate(Option_t *)
 683 {
 684 // Terminate
 685 }
 686
 687 #endif