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