PWG/muon/AliAnalysisTaskMuonTreeBuilder.cxx

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