PWG3/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
 245   fIsSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 246
 247   fNumMuonTracks = fESD->GetNumberOfMuonTracks() ;
 248    Int_t loopEnd = fNumMuonTracks;
 249   if(!fIsMC) {
 250     TString cla = fESD->GetFiredTriggerClasses();
 251     sprintf(fTrigClass,"%s",cla.Data());
 252   }
 253
 254   if(fNumMuonTracks>0 && fIsMC){
 255     mcEvent = MCEvent();
 256   }
 257   fNumContributors = fESD->GetPrimaryVertexSPD()->GetNContributors();
 258   fNumSPDTracklets = fESD->GetMultiplicity()->GetNumberOfTracklets();
 259   fVertex[0]=fESD->GetPrimaryVertexSPD()->GetX();
 260   fVertex[1]=fESD->GetPrimaryVertexSPD()->GetY();
 261   fVertex[2]=fESD->GetPrimaryVertexSPD()->GetZ();
 262   printf("fVertex : %f - %f - %f\n",fVertex[0],fVertex[1],fVertex[2]);
 263
 264   Int_t numdimu = 0;
 265   for (Int_t j = 0; j<loopEnd; j++) {
 266     AliESDMuonTrack* mu1 = new AliESDMuonTrack(*(fESD->GetMuonTrack(j)));
 267     if (!mu1->ContainTrackerData()) {fNumMuonTracks=fNumMuonTracks-1; continue;}
 268     Double_t charge1 = mu1->Charge();
 269     fCharge[j] = mu1->Charge();
 270     fpT[j]   = mu1->Pt();
 271     fpx[j]  = mu1->Px();
 272     fpy[j]  = mu1->Py();
 273     fpz[j]  = mu1->Pz();
 274     fpxUncorr[j]  = mu1->PxUncorrected();
 275     fpyUncorr[j]  = mu1->PyUncorrected();
 276     fpzUncorr[j]  = mu1->PzUncorrected();
 277     fy[j]    = mu1->Y();
 278     feta[j]    = mu1->Eta();
 279     fphi[j]  = Phideg(mu1->Phi());
 280     Double_t emu1    = mu1->E();
 281     fE[j] = emu1;
 282     fDCA[j] = mu1->GetDCA();
 283     fTrackChi2[j] = mu1->GetChi2();
 284     fMatchTrig[j]   = mu1->GetMatchTrigger();
 285     fMatchTrigChi2[j]= mu1->GetChi2MatchTrigger();
 286     fRAtAbsEnd[j]=mu1->GetRAtAbsorberEnd();
 287
 288     AliMCParticle* mcTrack = 0x0;
 289     if(fIsMC){
 290     if(mu1->GetLabel()==-1) continue;
 291     mcTrack = (AliMCParticle*)mcEvent->GetTrack(mu1->GetLabel());
 292     fPDG[j] = mcTrack->PdgCode();
 293     if (mcTrack->GetMother()==-1) continue;
 294     fPDGmother[j] = ((AliMCParticle*)mcEvent->GetTrack(mcTrack->GetMother()))->PdgCode();
 295     if (TMath::Abs(fPDG[j])==13) fMuFamily[j] = FindMuFamily(mcTrack,mcEvent);
 296     }
 297     for (Int_t jj = j+1; jj<loopEnd; jj++) {
 298       AliESDMuonTrack* mu2 = new AliESDMuonTrack(*(fESD->GetMuonTrack(jj)));
 299       if (!mu2->ContainTrackerData()) continue;
 300
 301       Double_t pxmu2  = mu2->Px();
 302       Double_t pymu2  = mu2->Py();
 303       Double_t pzmu2  = mu2->Pz();
 304       Double_t emu2   = mu2->E();
 305       //Double_t charge2= mu2->Charge();
 306
 307       fpTdimuon[numdimu] = TMath::Sqrt((fpx[j]+pxmu2)*(fpx[j]+pxmu2)+(fpy[j]+pymu2)*(fpy[j]+pymu2));
 308       fpxdimuon[numdimu] = fpx[j]+pxmu2;
 309       fpydimuon[numdimu] = fpy[j]+pymu2;
 310       fpzdimuon[numdimu] = fpz[j]+pzmu2;
 311       fydimuon[numdimu] = Rap((emu1+emu2),(fpz[j]+pzmu2));
 312       fiMassdimuon[numdimu] = Imass(emu1,fpx[j],fpy[j],fpz[j],emu2,pxmu2,pymu2,pzmu2);
 313       fcostCS[numdimu]=CostCS(fpx[j],fpy[j],fpz[j],emu1,charge1,pxmu2,pymu2,pzmu2,emu2);
 314       fcostHE[numdimu]=CostHE(fpx[j],fpy[j],fpz[j],emu1,charge1,pxmu2,pymu2,pzmu2,emu2);
 315       fphiCS[numdimu] = PhiCS(fpx[j],fpy[j],fpz[j],emu1,charge1,pxmu2,pymu2,pzmu2,emu2);
 316       fphiHE[numdimu] = PhiHE(fpx[j],fpy[j],fpz[j],emu1,charge1,pxmu2,pymu2,pzmu2,emu2);
 317
 318       fDimuonConstituent[numdimu][0]=j;  fDimuonConstituent[numdimu][1]=jj;
 319
 320       numdimu++;
 321
 322       if(fIsMC){
 323       if(mu2->GetLabel()==-1) continue;
 324       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);
 325       else fPDGdimu[numdimu-1]=-3;
 326       }
 327
 328       delete mu2;
 329     }
 330     fNumDimuons=numdimu;
 331
 332
 333     delete mu1;
 334     }
 335   fOutputTree->Fill();
 336   PostData(1,fOutputTree);
 337 }
 338
 339 //________________________________________________________________________
 340 Int_t AliAnalysisTaskMuonTreeBuilder::FindDimuFamily(AliMCParticle* mcTrack1,AliMCParticle* mcTrack2, AliMCEvent* mcEvent) const
 341 {
 342   // finds the family of the dimuon (works only if the 2 muons are real muons (not hadrons))
 343   Int_t familynumber;
 344
 345   if(mcTrack1->GetMother()==mcTrack2->GetMother()) familynumber = TMath::Abs(((AliMCParticle*)mcEvent->GetTrack(mcTrack1->GetMother()))->PdgCode());
 346   else{
 347     Int_t familymu1 = FindMuFamily(mcTrack1,mcEvent);
 348     Int_t familymu2 = FindMuFamily(mcTrack2,mcEvent);
 349     if(familymu1==5 && familymu2==5) familynumber=5;                                                    //bb dimuon
 350     else if(familymu1==4 && familymu2==4) familynumber=4;                                               //cc dimuon
 351     else if((familymu1==4 && familymu2==5)||(familymu2==4 && familymu1==5)) familynumber=45;            //bc dimuon
 352     else if (familymu1==-2 || familymu2==-2 || familymu1==-3 || familymu2==-3) familynumber=-2;         //hadron dimuon (at least 1 hadron involved)
 353     else familynumber=-1;
 354   }
 355
 356   return familynumber;
 357 }
 358
 359 //________________________________________________________________________
 360 Int_t AliAnalysisTaskMuonTreeBuilder::FindMuFamily(AliMCParticle* mcTrack, AliMCEvent* mcEvent) const
 361 {
 362   // finds the family of the muon
 363   Int_t imother = mcTrack->GetMother();
 364   if ( imother<0 ) return -1; // Drell-Yan Muon
 365
 366   Int_t igrandma = imother;
 367
 368   AliMCParticle* motherPart = (AliMCParticle*)mcEvent->GetTrack(imother);
 369   Int_t motherPdg = motherPart->PdgCode();
 370
 371   // Track is an heavy flavor muon
 372   Int_t absPdg = TMath::Abs(motherPdg);
 373   if(absPdg/100==5 || absPdg/1000==5) return 5;
 374   if(absPdg/100==4 || absPdg/1000==4){
 375     Int_t newMother = -1;
 376     igrandma = imother;
 377     Int_t absGrandMotherPdg = TMath::Abs(motherPart->PdgCode());
 378     while ( absGrandMotherPdg > 10 ) {
 379       igrandma = ((AliMCParticle*)mcEvent->GetTrack(igrandma))->GetMother();
 380       if( igrandma < 0 ) break;
 381       absGrandMotherPdg = TMath::Abs(((AliMCParticle*)mcEvent->GetTrack(igrandma))->PdgCode());
 382     }
 383
 384     if (absGrandMotherPdg==5) newMother = 5; // Charm from beauty
 385     else if (absGrandMotherPdg==4) newMother = 4;
 386
 387     if(newMother<0) {
 388       //AliWarning("Mother not correctly found! Set to charm!\n");
 389       newMother = 4;
 390     }
 391
 392     return newMother;
 393   }
 394
 395   Int_t nPrimaries = mcEvent->Stack()->GetNprimary();
 396   // Track is a bkg. muon
 397   if (imother<nPrimaries) {
 398     return -2;                  //is a primary
 399   }
 400   else {
 401     return -3;                  //is a secondary
 402   }
 403
 404 }
 405
 406 //________________________________________________________________________
 407 Double_t AliAnalysisTaskMuonTreeBuilder::Imass(Double_t e1, Double_t px1, Double_t py1, Double_t pz1,
 408                                    Double_t e2, Double_t px2, Double_t py2, Double_t pz2) const
 409 {
 410 // invariant mass calculation
 411     Double_t imassrec = TMath::Sqrt((e1+e2)*(e1+e2)-((px1+px2)*(px1+px2)+
 412                                     (py1+py2)*(py1+py2)+(pz1+pz2)*(pz1+pz2)));
 413     return imassrec;
 414 }
 415
 416 //________________________________________________________________________
 417 Double_t AliAnalysisTaskMuonTreeBuilder::Rap(Double_t e, Double_t pz) const
 418 {
 419 // calculate rapidity
 420     Double_t rap;
 421     if(e>TMath::Abs(pz)){
 422         rap = 0.5*TMath::Log((e+pz)/(e-pz));
 423         return rap;
 424     }
 425     else{
 426         rap = 666.;
 427         return rap;
 428     }
 429 }
 430
 431 //________________________________________________________________________
 432 Double_t AliAnalysisTaskMuonTreeBuilder::Phideg(Double_t phi) const
 433 {
 434 // calculate Phi in range [-180,180]
 435     Double_t phideg;
 436
 437         phideg = phi-TMath::Pi();
 438         phideg = phideg*57.296;
 439         return phideg;
 440 }
 441
 442 //________________________________________________________________________
 443 Double_t AliAnalysisTaskMuonTreeBuilder::CostCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
 444 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
 445 {
 446 // Cosine of the theta decay angle (mu+) in the Collins-Soper frame
 447
 448   TLorentzVector pMu1CM, pMu2CM, pProjCM, pTargCM, pDimuCM; // In the CM. frame
 449   TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame
 450   TVector3 beta,zaxisCS;
 451   Double_t mp=0.93827231;
 452   //
 453   // --- Fill the Lorentz vector for projectile and target in the CM frame
 454   //
 455   pProjCM.SetPxPyPzE(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp));
 456   pTargCM.SetPxPyPzE(0.,0.,fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp));
 457   //
 458   // --- Get the muons parameters in the CM frame
 459   //
 460   pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
 461   pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
 462   //
 463   // --- Obtain the dimuon parameters in the CM frame
 464   //
 465   pDimuCM=pMu1CM+pMu2CM;
 466   //
 467   // --- Translate the dimuon parameters in the dimuon rest frame
 468   //
 469   beta=(-1./pDimuCM.E())*pDimuCM.Vect();
 470   if(beta.Mag()>=1) return 666.;
 471   pMu1Dimu=pMu1CM;
 472   pMu2Dimu=pMu2CM;
 473   pProjDimu=pProjCM;
 474   pTargDimu=pTargCM;
 475   pMu1Dimu.Boost(beta);
 476   pMu2Dimu.Boost(beta);
 477   pProjDimu.Boost(beta);
 478   pTargDimu.Boost(beta);
 479
 480   //Debugging part -------------------------------------
 481   Double_t debugProj[4]={0.,0.,0.,0.};
 482   Double_t debugTarg[4]={0.,0.,0.,0.};
 483   Double_t debugMu1[4]={0.,0.,0.,0.};
 484   Double_t debugMu2[4]={0.,0.,0.,0.};
 485   pMu1Dimu.GetXYZT(debugMu1);
 486   pMu2Dimu.GetXYZT(debugMu2);
 487   pProjDimu.GetXYZT(debugProj);
 488   pTargDimu.GetXYZT(debugTarg);
 489   if (debugProj[0]!=debugProj[0] ||debugProj[1]!=debugProj[1] || debugProj[2]!=debugProj[2] ||debugProj[3]!=debugProj[3]) return 666;
 490   if (debugTarg[0]!=debugTarg[0] ||debugTarg[1]!=debugTarg[1] || debugTarg[2]!=debugTarg[2] ||debugTarg[3]!=debugTarg[3]) return 666;
 491   if (debugMu1[0]!=debugMu1[0] ||debugMu1[1]!=debugMu1[1] || debugMu1[2]!=debugMu1[2] ||debugMu1[3]!=debugMu1[3]) return 666;
 492   if (debugMu2[0]!=debugMu2[0] ||debugMu2[1]!=debugMu2[1] || debugMu2[2]!=debugMu2[2] ||debugMu2[3]!=debugMu2[3]) return 666;
 493   //----------------------------------------------------
 494
 495   // --- Determine the z axis for the CS angle
 496   zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit();
 497
 498   // --- Determine the CS angle (angle between mu+ and the z axis defined above)
 499   Double_t cost;
 500
 501   if(charge1>0) {cost = zaxisCS.Dot((pMu1Dimu.Vect()).Unit());}
 502   else {cost = zaxisCS.Dot((pMu2Dimu.Vect()).Unit());}
 503
 504   return cost;
 505 }
 506
 507 //________________________________________________________________________
 508 Double_t AliAnalysisTaskMuonTreeBuilder::CostHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
 509 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
 510 {
 511 // Cosine of the theta decay angle (mu+) in the Helicity frame
 512
 513   TLorentzVector pMu1CM, pMu2CM, pDimuCM; // In the CM frame
 514   TLorentzVector pMu1Dimu, pMu2Dimu; // In the dimuon rest frame
 515   TVector3 beta,zaxisCS;
 516   //
 517   // --- Get the muons parameters in the CM frame
 518   //
 519   pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
 520   pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
 521   //
 522   // --- Obtain the dimuon parameters in the CM frame
 523   //
 524   pDimuCM=pMu1CM+pMu2CM;
 525   //
 526   // --- Translate the muon parameters in the dimuon rest frame
 527   //
 528   beta=(-1./pDimuCM.E())*pDimuCM.Vect();
 529   if(beta.Mag()>=1) return 666.;
 530   pMu1Dimu=pMu1CM;
 531   pMu2Dimu=pMu2CM;
 532   pMu1Dimu.Boost(beta);
 533   pMu2Dimu.Boost(beta);
 534
 535   //Debugging part -------------------------------------
 536   Double_t debugMu1[4]={0.,0.,0.,0.};
 537   Double_t debugMu2[4]={0.,0.,0.,0.};
 538   pMu1Dimu.GetXYZT(debugMu1);
 539   pMu2Dimu.GetXYZT(debugMu2);
 540   if (debugMu1[0]!=debugMu1[0] ||debugMu1[1]!=debugMu1[1] || debugMu1[2]!=debugMu1[2] ||debugMu1[3]!=debugMu1[3]) return 666;
 541   if (debugMu2[0]!=debugMu2[0] ||debugMu2[1]!=debugMu2[1] || debugMu2[2]!=debugMu2[2] ||debugMu2[3]!=debugMu2[3]) return 666;
 542   //----------------------------------------------------
 543
 544   // --- Determine the z axis for the calculation of the polarization angle (i.e. the direction of the dimuon in the CM system)
 545   TVector3 zaxis;
 546   zaxis=(pDimuCM.Vect()).Unit();
 547
 548   // --- Calculation of the polarization angle (angle between mu+ and the z axis defined above)
 549   Double_t cost;
 550   if(charge1>0) {cost = zaxis.Dot((pMu1Dimu.Vect()).Unit());}
 551   else {cost = zaxis.Dot((pMu2Dimu.Vect()).Unit());}
 552   return cost;
 553 }
 554
 555 //________________________________________________________________________
 556 Double_t AliAnalysisTaskMuonTreeBuilder::PhiCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
 557 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
 558 {
 559 // Phi decay angle (mu+) in the Collins-Soper frame
 560
 561    TLorentzVector pMu1CM, pMu2CM, pProjCM, pTargCM, pDimuCM; // In the CM frame
 562    TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame
 563    TVector3 beta,yaxisCS, xaxisCS, zaxisCS;
 564    Double_t mp=0.93827231;
 565
 566    // --- Fill the Lorentz vector for projectile and target in the CM frame
 567    pProjCM.SetPxPyPzE(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp));
 568    pTargCM.SetPxPyPzE(0.,0.,fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp));
 569
 570    // --- Get the muons parameters in the CM frame
 571    pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
 572    pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
 573
 574    // --- Obtain the dimuon parameters in the CM frame
 575    pDimuCM=pMu1CM+pMu2CM;
 576
 577    // --- Translate the dimuon parameters in the dimuon rest frame
 578    beta=(-1./pDimuCM.E())*pDimuCM.Vect();
 579    if(beta.Mag()>=1) return 666.;
 580    pMu1Dimu=pMu1CM;
 581    pMu2Dimu=pMu2CM;
 582    pProjDimu=pProjCM;
 583    pTargDimu=pTargCM;
 584    pMu1Dimu.Boost(beta);
 585    pMu2Dimu.Boost(beta);
 586    pProjDimu.Boost(beta);
 587    pTargDimu.Boost(beta);
 588
 589    //Debugging part -------------------------------------
 590    Double_t debugProj[4]={0.,0.,0.,0.};
 591    Double_t debugTarg[4]={0.,0.,0.,0.};
 592    Double_t debugMu1[4]={0.,0.,0.,0.};
 593    Double_t debugMu2[4]={0.,0.,0.,0.};
 594    pMu1Dimu.GetXYZT(debugMu1);
 595    pMu2Dimu.GetXYZT(debugMu2);
 596    pProjDimu.GetXYZT(debugProj);
 597    pTargDimu.GetXYZT(debugTarg);
 598    if (debugProj[0]!=debugProj[0] ||debugProj[1]!=debugProj[1] || debugProj[2]!=debugProj[2] ||debugProj[3]!=debugProj[3]) return 666;
 599    if (debugTarg[0]!=debugTarg[0] ||debugTarg[1]!=debugTarg[1] || debugTarg[2]!=debugTarg[2] ||debugTarg[3]!=debugTarg[3]) return 666;
 600    if (debugMu1[0]!=debugMu1[0] ||debugMu1[1]!=debugMu1[1] || debugMu1[2]!=debugMu1[2] ||debugMu1[3]!=debugMu1[3]) return 666;
 601    if (debugMu2[0]!=debugMu2[0] ||debugMu2[1]!=debugMu2[1] || debugMu2[2]!=debugMu2[2] ||debugMu2[3]!=debugMu2[3]) return 666;
 602    //----------------------------------------------------
 603
 604    // --- Determine the z axis for the CS angle
 605    zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit();
 606    yaxisCS=(((pProjDimu.Vect()).Unit()).Cross((pTargDimu.Vect()).Unit())).Unit();
 607    xaxisCS=(yaxisCS.Cross(zaxisCS)).Unit();
 608
 609    Double_t phi=0.;
 610    if(charge1>0) {
 611        phi = TMath::ATan2((pMu1Dimu.Vect()).Dot(yaxisCS),((pMu1Dimu.Vect()).Dot(xaxisCS)));
 612    } else {
 613        phi = TMath::ATan2((pMu2Dimu.Vect()).Dot(yaxisCS),((pMu2Dimu.Vect()).Dot(xaxisCS)));
 614    }
 615    if (phi>TMath::Pi()) phi=phi-TMath::Pi();
 616
 617    return phi;
 618 }
 619
 620 //________________________________________________________________________
 621 Double_t AliAnalysisTaskMuonTreeBuilder::PhiHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
 622 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
 623 {
 624 // Phi decay angle (mu+) in the Helicity frame
 625   TLorentzVector pMu1Lab, pMu2Lab, pProjLab, pTargLab, pDimuLab; // In the lab. frame
 626   TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame
 627   TVector3 beta,xaxis, yaxis,zaxis;
 628   Double_t mp=0.93827231;
 629
 630   // --- Get the muons parameters in the LAB frame
 631   pMu1Lab.SetPxPyPzE(px1,py1,pz1,e1);
 632   pMu2Lab.SetPxPyPzE(px2,py2,pz2,e2);
 633
 634   // --- Obtain the dimuon parameters in the LAB frame
 635   pDimuLab=pMu1Lab+pMu2Lab;
 636   zaxis=(pDimuLab.Vect()).Unit();
 637
 638   // --- Translate the muon parameters in the dimuon rest frame
 639   beta=(-1./pDimuLab.E())*pDimuLab.Vect();
 640   if(beta.Mag()>=1.) return 666.;
 641
 642   pProjLab.SetPxPyPzE(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp)); // proiettile
 643   pTargLab.SetPxPyPzE(0.,0.,fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp)); // bersaglio
 644
 645   pProjDimu=pProjLab;
 646   pTargDimu=pTargLab;
 647
 648   pProjDimu.Boost(beta);
 649   pTargDimu.Boost(beta);
 650
 651   yaxis=((pProjDimu.Vect()).Cross(pTargDimu.Vect())).Unit();
 652   xaxis=(yaxis.Cross(zaxis)).Unit();
 653
 654   pMu1Dimu=pMu1Lab;
 655   pMu2Dimu=pMu2Lab;
 656   pMu1Dimu.Boost(beta);
 657   pMu2Dimu.Boost(beta);
 658
 659   //Debugging part -------------------------------------
 660   Double_t debugProj[4]={0.,0.,0.,0.};
 661   Double_t debugTarg[4]={0.,0.,0.,0.};
 662   Double_t debugMu1[4]={0.,0.,0.,0.};
 663   Double_t debugMu2[4]={0.,0.,0.,0.};
 664   pMu1Dimu.GetXYZT(debugMu1);
 665   pMu2Dimu.GetXYZT(debugMu2);
 666   pProjDimu.GetXYZT(debugProj);
 667   pTargDimu.GetXYZT(debugTarg);
 668   if (debugProj[0]!=debugProj[0] ||debugProj[1]!=debugProj[1] || debugProj[2]!=debugProj[2] ||debugProj[3]!=debugProj[3]) return 666;
 669   if (debugTarg[0]!=debugTarg[0] ||debugTarg[1]!=debugTarg[1] || debugTarg[2]!=debugTarg[2] ||debugTarg[3]!=debugTarg[3]) return 666;
 670   if (debugMu1[0]!=debugMu1[0] ||debugMu1[1]!=debugMu1[1] || debugMu1[2]!=debugMu1[2] ||debugMu1[3]!=debugMu1[3]) return 666;
 671   if (debugMu2[0]!=debugMu2[0] ||debugMu2[1]!=debugMu2[1] || debugMu2[2]!=debugMu2[2] ||debugMu2[3]!=debugMu2[3]) return 666;
 672   //----------------------------------------------------
 673
 674   Double_t phi=0.;
 675    if(charge1 > 0) {
 676       phi = TMath::ATan2((pMu1Dimu.Vect()).Dot(yaxis),(pMu1Dimu.Vect()).Dot(xaxis));
 677      } else {
 678       phi = TMath::ATan2((pMu2Dimu.Vect()).Dot(yaxis),(pMu2Dimu.Vect()).Dot(xaxis));
 679    }
 680    return phi;
 681 }
 682
 683 //________________________________________________________________________
 684 void AliAnalysisTaskMuonTreeBuilder::Terminate(Option_t *)
 685 {
 686 // Terminate
 687 }
 688
 689 #endif