PWG3/muon/AliAnalysisTaskDimuonCFContainerBuilder.cxx

   1 #ifndef ALIANALYSISTASKDIMUONCFCONTAINERBUILDER_CXX
   2 #define ALIANALYSISTASKDIMUONCFCONTAINERBUILDER_CXX
   3
   4 #include "AliAnalysisTaskDimuonCFContainerBuilder.h"
   5 #include "AliStack.h"
   6 #include "TParticle.h"
   7 #include "TString.h"
   8 #include "TLorentzVector.h"
   9 #include "TH2D.h"
  10 #include "AliAnalysisManager.h"
  11 #include "AliESDEvent.h"
  12 #include "AliAODEvent.h"
  13 #include "AliCFManager.h"
  14 #include "AliCFContainer.h"
  15 #include "AliESDMuonTrack.h"
  16 #include "AliESDInputHandler.h"
  17 #include "AliAODInputHandler.h"
  18 #include "AliAODMCParticle.h"
  19
  20 //      Analysis task for the building of a dimuon CF container
  21 //      Also some single-muon variables are stored
  22 //      All the variable ranges and binning are set in the AddTask macro
  23 //
  24 //      author: L. Bianchi - Universita' & INFN Torino
  25
  26
  27 ClassImp(AliAnalysisTaskDimuonCFContainerBuilder)
  28
  29 //__________________________________________________________________________
  30 AliAnalysisTaskDimuonCFContainerBuilder::AliAnalysisTaskDimuonCFContainerBuilder() :
  31   fReadAODData(0),
  32   fReadMCInfo(kTRUE),
  33   fIsAccProduction(kTRUE),
  34   fCFManager(0x0),
  35   fQAHistList(0x0),
  36   fNevt(0),
  37   fBeamEnergy(3500.),
  38   fOutput(0x0),
  39   fCutOnzVtxSPD(kFALSE),
  40   fCutOnNContributors(kFALSE),
  41   fTrigClassMuon(""),
  42   fTrigClassInteraction(""),
  43   fDistinguishTrigClass(kFALSE)
  44 {
  45
  46   //Default constructor
  47
  48   Double_t chilims[2]={0.,10000.};
  49   Double_t ptlims[2]={0.,100.};
  50   Double_t thetalims[2]={0.,180.};
  51   Double_t vtxlims[2]={-1000.,1000.};
  52   SetChi2Limits(chilims);
  53   SetChi2MatchLimits(chilims);
  54   SetPtSingMuLimits(ptlims);
  55   SetThetaSingMuLimits(thetalims);
  56   SetZprimVertLimits(vtxlims);
  57   SetTrigClassMuonName();
  58   SetTrigClassInteracName();
  59   TString namemuonside[4]={"CMUS1A-","CMUS1B-","CMUS1C-","CMUS1-E-"};
  60   TString nameinteractionside[4]={"CINT1A-","CINT1B-","CINT1C-","CINT1-E"};
  61   SetTrigClassMuonSideName(namemuonside);
  62   SetTrigClassInteracSideName(nameinteractionside);
  63 }
  64 //___________________________________________________________________________
  65 AliAnalysisTaskDimuonCFContainerBuilder::AliAnalysisTaskDimuonCFContainerBuilder(const Char_t* name, Bool_t readaod, Bool_t readMC, Bool_t isaccept, Double_t beamEn) :
  66   AliAnalysisTaskSE(name),
  67   fReadAODData(0),
  68   fReadMCInfo(kTRUE),
  69   fIsAccProduction(kTRUE),
  70   fCFManager(0x0),
  71   fQAHistList(0x0),
  72   fNevt(0),
  73   fBeamEnergy(3500.),
  74   fOutput(0x0),
  75   fCutOnzVtxSPD(kFALSE),
  76   fCutOnNContributors(kFALSE),
  77   fTrigClassMuon(""),
  78   fTrigClassInteraction(""),
  79   fDistinguishTrigClass(kFALSE)
  80 {
  81   //
  82   // Constructor. Initialization of Inputs and Outputs
  83   //
  84   Info("AliAnalysisTaskDimuonCFContainerBuilder","Calling Constructor");
  85
  86   SetReadAODData(readaod);
  87   SetReadMCinfo(readMC);
  88   SetIsAccProd(isaccept);
  89   SetBeamEnergy(beamEn);
  90
  91   Double_t chilims[2]={0.,10000.};
  92   Double_t ptlims[2]={0.,100.};
  93   Double_t thetalims[2]={0.,180.};
  94   Double_t vtxlims[2]={-1000.,1000.};
  95   SetChi2Limits(chilims);
  96   SetChi2MatchLimits(chilims);
  97   SetPtSingMuLimits(ptlims);
  98   SetThetaSingMuLimits(thetalims);
  99   SetZprimVertLimits(vtxlims);
 100   SetTrigClassMuonName();
 101   SetTrigClassInteracName();
 102   TString namemuonside[4]={"CMUS1A-","CMUS1B-","CMUS1C-","CMUS1-E-"};
 103   TString nameinteractionside[4]={"CINT1A-","CINT1B-","CINT1C-","CINT1-E"};
 104   SetTrigClassMuonSideName(namemuonside);
 105   SetTrigClassInteracSideName(nameinteractionside);
 106
 107   DefineOutput(1,TList::Class());
 108   DefineOutput(2,AliCFContainer::Class());
 109
 110 }
 111
 112 //___________________________________________________________________________
 113 AliAnalysisTaskDimuonCFContainerBuilder& AliAnalysisTaskDimuonCFContainerBuilder::operator=(const AliAnalysisTaskDimuonCFContainerBuilder& c)
 114 {
 115   //
 116   // Assignment operator
 117   //
 118   if (this!=&c) {
 119     AliAnalysisTaskSE::operator=(c) ;
 120     fReadAODData = c.fReadAODData ;
 121     fCFManager  = c.fCFManager;
 122     fQAHistList = c.fQAHistList ;
 123     fNevt = c.fNevt ;
 124   }
 125   return *this;
 126 }
 127
 128 //___________________________________________________________________________
 129 AliAnalysisTaskDimuonCFContainerBuilder::AliAnalysisTaskDimuonCFContainerBuilder(const AliAnalysisTaskDimuonCFContainerBuilder& c) :
 130   AliAnalysisTaskSE(c),
 131   fReadAODData(c.fReadAODData),
 132   fReadMCInfo(c.fReadMCInfo),
 133   fIsAccProduction(c.fIsAccProduction),
 134   fCFManager(c.fCFManager),
 135   fQAHistList(c.fQAHistList),
 136   fNevt(c.fNevt),
 137   fBeamEnergy(c.fBeamEnergy),
 138   fOutput(c.fOutput),
 139   fCutOnzVtxSPD(c.fCutOnzVtxSPD),
 140   fCutOnNContributors(c.fCutOnNContributors),
 141   fTrigClassMuon(c.fTrigClassMuon),
 142   fTrigClassInteraction(c.fTrigClassInteraction),
 143   fDistinguishTrigClass(c.fDistinguishTrigClass)
 144 {
 145
 146   // Copy Constructor
 147
 148 }
 149
 150 //___________________________________________________________________________
 151 AliAnalysisTaskDimuonCFContainerBuilder::~AliAnalysisTaskDimuonCFContainerBuilder() {
 152   //
 153   //destructor
 154   //
 155   Info("~AliAnalysisTaskDimuonCFContainerBuilder","Calling Destructor");
 156   if (fCFManager)           delete fCFManager ;
 157   if (fQAHistList) {fQAHistList->Clear(); delete fQAHistList;}
 158 }
 159
 160 //___________________________________________________________________________
 161 void AliAnalysisTaskDimuonCFContainerBuilder::UserCreateOutputObjects(){
 162  //     UserCreateOutputObjects
 163  fOutput = new TList();
 164  fOutput->SetOwner();
 165
 166  TH1D *hnevts   = new TH1D("hnevts","hnevts",1,0,1);                                    // Stat check histos
 167  TH1D *jpsiMult = new TH1D("jpsiMult","jpsiMult",20,0,20);
 168
 169  TH1D *ptdimuREC    = new TH1D("ptdimuREC","ptdimuREC",200,0,20);                       // Dimu check histos
 170  TH1D *ydimuREC     = new TH1D("ydimuREC","ydimuREC",200,-10.,10.);
 171  TH1D *ydimuFail     = new TH1D("ydimuFail","ydimuFail",10,660,670);
 172  TH1D *costHEdimuREC= new TH1D("costHEdimuREC","costHEdimuREC",200,-1.,1.);
 173  TH1D *costCSdimuREC= new TH1D("costCSdimuREC","costCSdimuREC",200,-1.,1.);
 174  TH1D *costHEdimuFail= new TH1D("costHEdimuFail","costHEdimuFail",10,660.,670.);
 175  TH1D *costCSdimuFail= new TH1D("costCSdimuFail","costCSdimuFail",10,660.,670.);
 176  TH1D *phiHEdimuREC = new TH1D("phiHEdimuREC","phiHEdimuREC",100,0.,TMath::Pi());
 177  TH1D *phiCSdimuREC = new TH1D("phiCSdimuREC","phiCSdimuREC",100,0.,TMath::Pi());
 178  TH1D *phiHEdimuFail = new TH1D("phiHEdimuFail","phiHEdimuFail",10,660.,670.);
 179  TH1D *phiCSdimuFail = new TH1D("phiCSdimuFail","phiCSdimuFail",10,660.,670.);
 180  TH1D *imassTot     = new TH1D("imassTot","imassTot",100,0,8);
 181  TH1D *trigCond     = new TH1D("trigCond","trigCond",40,0,4);
 182
 183  TH1D *emuonREC = new TH1D("emuonREC","emuonREC",200,0.,20.);                           // Mu check histos
 184  TH1D *ptmuonREC= new TH1D("ptmuonREC","ptmuonREC",200,0.,20.);
 185  TH1D *ymuonREC = new TH1D("ymuonREC","ymuonREC",200,-10.,10.);
 186  TH1D *hdca     = new TH1D("hdca","hdca",20,0,200);
 187  TH2D *hdcay    = new TH2D("hdcay","hdcay",200,0,200,20,-5,0);
 188
 189  TH1D *zvSPD    = new TH1D("zvSPD","zvSPD",100,-50,50.);                                // Event check histos
 190  TH1D *zvSPDcut = new TH1D("zvSPDcut","zvSPDcut",100,-50,50.);
 191  TH1D *nContrib = new TH1D("nContrib","nContrib",100,0,100);
 192
 193
 194  fOutput->Add(hnevts);
 195  fOutput->Add(jpsiMult);
 196
 197  fOutput->Add(ptdimuREC);
 198  fOutput->Add(ydimuREC);
 199  fOutput->Add(ydimuFail);
 200  fOutput->Add(costHEdimuREC);
 201  fOutput->Add(costCSdimuREC);
 202  fOutput->Add(costHEdimuFail);
 203  fOutput->Add(costCSdimuFail);
 204  fOutput->Add(phiHEdimuREC);
 205  fOutput->Add(phiCSdimuREC);
 206  fOutput->Add(phiHEdimuFail);
 207  fOutput->Add(phiCSdimuFail);
 208  fOutput->Add(imassTot);
 209  fOutput->Add(trigCond);
 210
 211  fOutput->Add(emuonREC);
 212  fOutput->Add(ptmuonREC);
 213  fOutput->Add(ymuonREC);
 214  fOutput->Add(hdca);
 215  fOutput->Add(hdcay);
 216
 217  fOutput->Add(zvSPD);
 218  fOutput->Add(zvSPDcut);
 219  fOutput->Add(nContrib);
 220
 221
 222  fOutput->ls();
 223
 224
 225 }
 226
 227
 228 //_________________________________________________
 229 void AliAnalysisTaskDimuonCFContainerBuilder::UserExec(Option_t *)
 230 {
 231
 232
 233   //Info("UserExec","");
 234   fNevt++;
 235   ((TH1D*)(fOutput->FindObject("hnevts")))->Fill(0.5);
 236
 237   Double_t containerInput[15]={666,666,666,666,666,666,666,666,666,666,666,666,666,666,666};   //y, pT, costHE, phiHE, costCS, phiCS, mass, TrigCond
 238   Int_t cutAccept =1;
 239   Int_t numbJpsis =0;
 240
 241   if (!fReadAODData){     // ESD-based ANALYSIS
 242
 243     if (fReadMCInfo){
 244
 245       if (!fMCEvent) {
 246         Error("UserExec","NO MC EVENT FOUND!");
 247         //return;
 248       }
 249
 250       // MC part  ---------------------------------------------------------------------------------------------------
 251
 252       //fCFManager->SetEventInfo(fMCEvent);                                     CHANGES IN NEW VERSIONS - MANUAL CUT ON PDG
 253       AliStack* stack = fMCEvent->Stack();
 254
 255       // loop on the MC event
 256       for (Int_t ipart=0; ipart<fMCEvent->GetNumberOfTracks(); ipart++) {
 257
 258         AliMCParticle *mcPart  = (AliMCParticle*) fMCEvent->GetTrack(ipart);
 259
 260         TParticle *part = mcPart->Particle();
 261         TParticle *part0 = mcPart->Particle();
 262         TParticle *part1 = mcPart->Particle();
 263
 264         // Mother kinematics
 265         Double_t e = part->Energy();
 266         Double_t pz = part->Pz();
 267         Double_t rapmc = Rap(e,pz);
 268
 269         // Selection of the resonance
 270         //if (!fCFManager->CheckParticleCuts(0,mcPart)) continue;
 271         if (part->GetPdgCode()!=443) continue;                                  // MANUAL CUT ON PDG
 272         numbJpsis++;
 273
 274         // Decays kinematics
 275         Int_t p0 = part->GetDaughter(0);
 276         part0 = stack->Particle(p0);
 277         Int_t pdg0 = part0->GetPdgCode();
 278
 279         Int_t p1 = part->GetDaughter(1);
 280         part1 = stack->Particle(p1);
 281         Int_t pdg1 = part1->GetPdgCode();
 282
 283         Double_t e0 = part0->Energy();
 284         Double_t pz0 = part0->Pz();
 285         Double_t py0 = part0->Py();
 286         Double_t px0 = part0->Px();
 287         Double_t charge0 = part0->GetPDG()->Charge()/3;
 288         Double_t theta0 = (180./TMath::Pi())*TMath::ATan2(TMath::Sqrt(px0*px0+py0*py0),pz0);
 289         Double_t pt0 = TMath::Sqrt(px0*px0+py0*py0);
 290         Double_t mom0 = part0->P();
 291
 292         Double_t e1 = part1->Energy();
 293         Double_t pz1 = part1->Pz();
 294         Double_t py1 = part1->Py();
 295         Double_t px1 = part1->Px();
 296         //Double_t charge1 = part1->GetPDG()->Charge()/3;
 297         Double_t theta1 = (180./TMath::Pi())*TMath::ATan2(TMath::Sqrt(px1*px1+py1*py1),pz1);
 298         Double_t pt1 = TMath::Sqrt(px1*px1+py1*py1);
 299         Double_t mom1 = part1->P();
 300
 301
 302         if(pdg0==13 || pdg1==13) {
 303
 304           Double_t ptmc = TMath::Sqrt((px0+px1)*(px0+px1)+(py0+py1)*(py0+py1));
 305           Double_t imassmc = Imass(e0,px0,py0,pz0,e1,px1,py1,pz1);
 306
 307           Double_t costCS = CostCS(px0,py0,pz0,e0,charge0,px1,py1,pz1,e1);
 308           Double_t costHE = CostHE(px0,py0,pz0,e0,charge0,px1,py1,pz1,e1);
 309           Double_t phiCS  = PhiCS(px0,py0,pz0,e0,charge0,px1,py1,pz1,e1);
 310           Double_t phiHE  = PhiHE(px0,py0,pz0,e0,charge0,px1,py1,pz1,e1);
 311
 312           containerInput[0]  = rapmc ;
 313           containerInput[1]  = ptmc;
 314           containerInput[2]  = costHE;
 315           containerInput[3]  = TMath::Abs(phiHE);
 316           containerInput[4]  = costCS;
 317           containerInput[5]  = TMath::Abs(phiCS);
 318           containerInput[6]  = imassmc;
 319           containerInput[7]  = 1.;              //for generated no trigger condition
 320           if (pt0<pt1) {
 321             containerInput[8]=pt0;
 322             containerInput[9]=pt1;
 323           } else {
 324             containerInput[8]=pt1;
 325             containerInput[9]=pt0;
 326           }
 327           if (theta0<theta1) {
 328             containerInput[10]=theta0;
 329             containerInput[11]=theta1;
 330           } else {
 331             containerInput[10]=theta1;
 332             containerInput[11]=theta0;
 333           }
 334           if (mom0<mom1) {
 335             containerInput[12]=mom0;
 336             containerInput[13]=mom1;
 337           } else {
 338             containerInput[12]=mom1;
 339             containerInput[13]=mom0;
 340           }
 341           containerInput[14]=1.;
 342
 343           // fill the container at the first step
 344           fCFManager->GetParticleContainer()->Fill(containerInput,0);
 345
 346           if(fIsAccProduction){
 347             // acceptance cuts on single mu
 348             if(theta0<fThetaSingMuCut[0]||theta0>fThetaSingMuCut[1]||theta1<fThetaSingMuCut[0]||theta1>fThetaSingMuCut[1]) cutAccept=0;
 349             if(pt0<fPtSingMuCut[0] || pt0>fPtSingMuCut[1] || pt1<fPtSingMuCut[0] || pt1>fPtSingMuCut[1]) cutAccept=0;
 350             if (cutAccept ==1) fCFManager->GetParticleContainer()->Fill(containerInput,2);
 351           }
 352         }
 353       }
 354     } //fReadMCInfo
 355
 356
 357       // ESD part  ---------------------------------------------------------------------------------------------------
 358
 359       AliESDEvent *fESD;
 360       AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 361       fESD = esdH->GetEvent();
 362       Int_t mult1 = fESD->GetNumberOfMuonTracks() ;
 363
 364       Int_t trigfired=-1;
 365       Int_t trigside=-1;
 366       if(fDistinguishTrigClass){
 367         //aod->GetHeader()->SetFiredTriggerClasses("CMU");
 368         TString trigclass = fESD->GetFiredTriggerClasses();
 369         printf("TrigClass: %s\n",trigclass.Data());
 370         if(trigclass.Contains(fTrigClassMuon)) trigfired = 1;
 371         else if (trigclass.Contains(fTrigClassInteraction)) trigfired = 0;
 372         printf("Muon container: %d\n",trigfired);
 373         for(Int_t i=0;i<4;i++){
 374           if(trigfired==1 && trigclass.Contains(fTrigClassMuonSide[i])) {trigside=i;printf("Muon Side: %d\n",trigside);}
 375           if(trigfired==0 && trigclass.Contains(fTrigClassInteractionSide[i])) {trigside=i;printf("Interaction Side: %d\n",trigside);}
 376         }
 377       }
 378
 379       ((TH1D*)(fOutput->FindObject("zvSPD")))->Fill(fESD->GetPrimaryVertexSPD()->GetZ());
 380       ((TH1D*)(fOutput->FindObject("nContrib")))->Fill(fESD->GetPrimaryVertexSPD()->GetNContributors());
 381       if (!fCutOnzVtxSPD || (fCutOnzVtxSPD && (fESD->GetPrimaryVertexSPD()->GetZ()>fzPrimVertexSPD[0] && fESD->GetPrimaryVertexSPD()->GetZ()<fzPrimVertexSPD[1]))){
 382       ((TH1D*)(fOutput->FindObject("zvSPDcut")))->Fill(fESD->GetPrimaryVertexSPD()->GetZ());
 383       if (!fCutOnNContributors || (fCutOnNContributors && (fESD->GetPrimaryVertexSPD()->GetNContributors()>fNContributors[0] && fESD->GetPrimaryVertexSPD()->GetNContributors()<fNContributors[1]))){
 384       for (Int_t j = 0; j<mult1; j++) {
 385
 386         AliESDMuonTrack* mu1 = new AliESDMuonTrack(*(fESD->GetMuonTrack(j)));
 387         if (!mu1->ContainTrackerData()) continue;
 388         if (mu1->GetChi2()<fChi2Track[0] || mu1->GetChi2()>fChi2Track[1]) continue;
 389         if (mu1->GetChi2MatchTrigger()<fChi2MatchTrig[0] || mu1->GetChi2MatchTrigger()>fChi2MatchTrig[1]) continue;
 390         Double_t zr1 = mu1->Charge();
 391         Double_t pxr1  = mu1->Px();
 392         Double_t pyr1  = mu1->Py();
 393         Double_t pzr1  = mu1->Pz();
 394         Double_t ptmu1 = TMath::Sqrt(pxr1*pxr1+pyr1*pyr1);
 395         Double_t er1 = mu1->E();
 396         Double_t mom1 = mu1->P();
 397         Double_t theta1 = (180./TMath::Pi())*mu1->Theta();
 398         Double_t rapiditymu1 = Rap(er1,pzr1);
 399         ((TH1D*)(fOutput->FindObject("emuonREC")))->Fill(er1);
 400         ((TH1D*)(fOutput->FindObject("ptmuonREC")))->Fill(ptmu1);
 401         ((TH1D*)(fOutput->FindObject("ymuonREC")))->Fill(rapiditymu1);
 402
 403         if(zr1<0){
 404
 405             for (Int_t jj = 0; jj<mult1; jj++) {
 406
 407               AliESDMuonTrack* mu2 = new AliESDMuonTrack(*(fESD->GetMuonTrack(jj)));
 408               if (!mu2->ContainTrackerData()) continue;
 409               if (mu2->GetChi2()<fChi2Track[0] || mu2->GetChi2()>fChi2Track[1]) continue;
 410               if (mu2->GetChi2MatchTrigger()<fChi2MatchTrig[0] || mu2->GetChi2MatchTrigger()>fChi2MatchTrig[1]) continue;
 411               Double_t zr2 = mu2->Charge();
 412
 413               if(zr2>0){
 414                 Double_t trigCondition=0;
 415                 if (mu1->GetMatchTrigger()>=mu2->GetMatchTrigger()) trigCondition = mu1->GetMatchTrigger()+0.1*mu2->GetMatchTrigger();
 416                 else trigCondition = mu2->GetMatchTrigger()+0.1*mu1->GetMatchTrigger();
 417                 ((TH1D*)(fOutput->FindObject("trigCond")))->Fill(trigCondition);
 418                 Double_t pxr2 = mu2->Px();
 419                 Double_t pyr2 = mu2->Py();
 420                 Double_t pzr2 = mu2->Pz();
 421                 Double_t ptmu2 = TMath::Sqrt(pxr2*pxr2+pyr2*pyr2);
 422                 Double_t er2 = mu2->E();
 423                 Double_t mom2 = mu2->P();
 424                 Double_t theta2 = (180./TMath::Pi())*mu2->Theta();
 425
 426                 Double_t ptrec = TMath::Sqrt((pxr1+pxr2)*(pxr1+pxr2)+(pyr1+pyr2)*(pyr1+pyr2));
 427                 ((TH1D*)(fOutput->FindObject("ptdimuREC")))->Fill(ptrec);
 428                 Double_t raprec= Rap((er1+er2),(pzr1+pzr2));
 429                 ((TH1D*)(fOutput->FindObject("ydimuREC")))->Fill(raprec);
 430                 Double_t imassrec = Imass(er1,pxr1,pyr1,pzr1,er2,pxr2,pyr2,pzr2);
 431                 ((TH1D*)(fOutput->FindObject("imassTot")))->Fill(imassrec);
 432
 433                 if (imassrec>12.) continue;
 434
 435                 Double_t costCSrec = CostCS(pxr1,pyr1,pzr1,er1,zr1,pxr2,pyr2,pzr2,er2);
 436                 ((TH1D*)(fOutput->FindObject("costCSdimuREC")))->Fill(costCSrec);
 437                 if((Int_t)costCSrec==666) ((TH1D*)(fOutput->FindObject("costCSdimuFail")))->Fill(costCSrec);
 438                 Double_t costHErec = CostHE(pxr1,pyr1,pzr1,er1,zr1,pxr2,pyr2,pzr2,er2);
 439                 ((TH1D*)(fOutput->FindObject("costHEdimuREC")))->Fill(costHErec);
 440                 if((Int_t)costHErec==666) ((TH1D*)(fOutput->FindObject("costHEdimuFail")))->Fill(costHErec);
 441                 Double_t phiCSrec  = PhiCS(pxr1,pyr1,pzr1,er1,zr1,pxr2,pyr2,pzr2,er2);
 442                 ((TH1D*)(fOutput->FindObject("phiCSdimuREC")))->Fill(phiCSrec);
 443                 if((Int_t)phiCSrec==666) ((TH1D*)(fOutput->FindObject("phiCSdimuFail")))->Fill(phiCSrec);
 444                 Double_t phiHErec  = PhiHE(pxr1,pyr1,pzr1,er1,zr1,pxr2,pyr2,pzr2,er2);
 445                 ((TH1D*)(fOutput->FindObject("phiHEdimuREC")))->Fill(phiHErec);
 446                 if((Int_t)phiHErec==666) ((TH1D*)(fOutput->FindObject("phiHEdimuFail")))->Fill(phiHErec);
 447
 448                 containerInput[0] = raprec ;
 449                 containerInput[1] = ptrec;
 450                 containerInput[2] = costHErec;
 451                 containerInput[3] = TMath::Abs(phiHErec);
 452                 containerInput[4] = costCSrec;
 453                 containerInput[5] = TMath::Abs(phiCSrec);
 454                 containerInput[6] = imassrec;
 455                 containerInput[7] = trigCondition+0.05;
 456                 if (ptmu1<ptmu2) {
 457                   containerInput[8]=ptmu1;
 458                   containerInput[9]=ptmu2;
 459                 } else {
 460                   containerInput[8]=ptmu2;
 461                   containerInput[9]=ptmu1;
 462                 }
 463                 if (theta1<theta2) {
 464                   containerInput[10]=theta1;
 465                   containerInput[11]=theta2;
 466                 } else {
 467                   containerInput[10]=theta2;
 468                   containerInput[11]=theta1;
 469                 }
 470                 if (mom1<mom2) {
 471                   containerInput[12]=mom1;
 472                   containerInput[13]=mom2;
 473                 } else {
 474                   containerInput[12]=mom2;
 475                   containerInput[13]=mom1;
 476                 }
 477                 if (fDistinguishTrigClass && trigside==0) containerInput[14]=0.;
 478                 else if (fDistinguishTrigClass && trigside==1) containerInput[14]=1.;
 479                 else if (fDistinguishTrigClass && trigside==2) containerInput[14]=2.;
 480                 else if (fDistinguishTrigClass && trigside==3) containerInput[14]=3.;
 481                 else containerInput[14]=0.;
 482
 483                 if(fDistinguishTrigClass){
 484                   if (trigfired==1) fCFManager->GetParticleContainer()->Fill(containerInput,5);
 485                   else if (trigfired==0) fCFManager->GetParticleContainer()->Fill(containerInput,1);
 486                 } else {
 487                   fCFManager->GetParticleContainer()->Fill(containerInput,1);
 488                   if(fIsAccProduction){
 489                     if(cutAccept==1) fCFManager->GetParticleContainer()->Fill(containerInput,3);
 490                   }
 491                 }
 492               }  // mu+ Selection
 493
 494             }      // second mu Loop
 495          }          // mu- Selection
 496      }
 497         }
 498       }
 499   } else {     // AOD-based ANALYSIS
 500
 501     AliAODEvent *aod;
 502     AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 503     aod = aodH->GetEvent();
 504     Int_t ntracks=aod->GetNumberOfTracks();
 505
 506     // MC part  ---------------------------------------------------------------------------------------------------
 507
 508     if (fReadMCInfo){
 509       TClonesArray *mcarray = dynamic_cast<TClonesArray*> (aod->FindListObject(AliAODMCParticle::StdBranchName()));  //array of MC particles in this event
 510       for(int ii=0;ii<mcarray->GetEntries();ii++){
 511         AliAODMCParticle *mctrack = (AliAODMCParticle*) mcarray->At(ii);
 512         if(mctrack->GetPdgCode()!=13) continue;
 513         Int_t numbMother = mctrack->GetMother();
 514         if (numbMother==-1) continue;
 515         AliAODMCParticle *mother = (AliAODMCParticle*) mcarray->At(numbMother);
 516         if (mother->GetPdgCode()!=443) continue;
 517         numbJpsis++;
 518         Int_t daught0 = mother->GetDaughter(0);
 519         Int_t daught1 = mother->GetDaughter(1);
 520         AliAODMCParticle *mcDaughter0 = (AliAODMCParticle*) mcarray->At(daught0);
 521         Double_t pxmc0 = mcDaughter0->Px();
 522         Double_t pymc0 = mcDaughter0->Py();
 523         Double_t pzmc0 = mcDaughter0->Pz();
 524         Double_t ptmc0 = mcDaughter0->Pt();
 525         Double_t mommc0 = mcDaughter0->P();
 526         Double_t emc0  = mcDaughter0->E();
 527         Double_t thetamc0 = (180./TMath::Pi())*mcDaughter0->Theta();
 528         Int_t charge0 = (Int_t) mcDaughter0->Charge()/3;
 529         AliAODMCParticle *mcDaughter1 = (AliAODMCParticle*) mcarray->At(daught1);
 530         Double_t pxmc1 = mcDaughter1->Px();
 531         Double_t pymc1 = mcDaughter1->Py();
 532         Double_t pzmc1 = mcDaughter1->Pz();
 533         Double_t ptmc1 = mcDaughter1->Pt();
 534         Double_t mommc1 = mcDaughter1->P();
 535         Double_t emc1  = mcDaughter1->E();
 536         Double_t thetamc1 = (180./TMath::Pi())*mcDaughter1->Theta();
 537         Int_t charge1 = (Int_t) mcDaughter1->Charge()/3;
 538
 539         if (charge0==charge1 || TMath::Abs(mcDaughter0->GetPdgCode())!=13 || TMath::Abs(mcDaughter1->GetPdgCode())!=13) continue;
 540         Double_t rapJpsi = Rap(mother->E(),mother->Pz());
 541         Double_t ptJpsi = mother->Pt();
 542         Double_t costHE = CostHE(pxmc0,pymc0,pzmc0,emc0,(Double_t)charge0,pxmc1,pymc1,pzmc1,emc1);
 543         Double_t phiHE  = PhiHE(pxmc0,pymc0,pzmc0,emc0,(Double_t)charge0,pxmc1,pymc1,pzmc1,emc1);
 544         Double_t costCS = CostCS(pxmc0,pymc0,pzmc0,emc0,(Double_t)charge0,pxmc1,pymc1,pzmc1,emc1);
 545         Double_t phiCS  = PhiCS(pxmc0,pymc0,pzmc0,emc0,(Double_t)charge0,pxmc1,pymc1,pzmc1,emc1);
 546         Double_t massJpsi = mother->M();
 547         containerInput[0] = rapJpsi;
 548         containerInput[1] = ptJpsi;
 549         containerInput[2] = costHE;
 550         containerInput[3] = TMath::Abs(phiHE);
 551         containerInput[4] = costCS;
 552         containerInput[5] = TMath::Abs(phiCS);
 553         containerInput[6] = massJpsi;
 554         containerInput[7] = 1.;         //for generated no trigger condition
 555         if (ptmc0<ptmc1) {
 556           containerInput[8]=ptmc0;
 557           containerInput[9]=ptmc1;
 558         } else {
 559           containerInput[8]=ptmc1;
 560           containerInput[9]=ptmc0;
 561         }
 562         if (thetamc0<thetamc1) {
 563           containerInput[10]=thetamc0;
 564           containerInput[11]=thetamc1;
 565         } else {
 566           containerInput[10]=thetamc1;
 567           containerInput[11]=thetamc0;
 568         }
 569         if (mommc0<mommc1) {
 570           containerInput[12]=mommc0;
 571           containerInput[13]=mommc1;
 572         } else {
 573           containerInput[12]=mommc1;
 574           containerInput[13]=mommc0;
 575         }
 576         containerInput[14]=1.;
 577
 578         if((Int_t)rapJpsi!=666 && (Int_t)costHE!=666 && (Int_t)phiHE!=666 && (Int_t)costCS!=666 && (Int_t)phiCS!=666){
 579           fCFManager->GetParticleContainer()->Fill(containerInput,0);
 580           if(fIsAccProduction){
 581             // acceptance cuts on single mu
 582             if(thetamc0<fThetaSingMuCut[0]||thetamc0>fThetaSingMuCut[1]||thetamc1<fThetaSingMuCut[0]||thetamc1>fThetaSingMuCut[1]) cutAccept=0;
 583             if(ptmc0<fPtSingMuCut[0] || ptmc0>fPtSingMuCut[1] || ptmc1<fPtSingMuCut[0] || ptmc1>fPtSingMuCut[1]) cutAccept=0;
 584             if (cutAccept ==1) fCFManager->GetParticleContainer()->Fill(containerInput,2);
 585           }
 586
 587         }
 588
 589         }
 590       }
 591
 592
 593       // AOD part  ---------------------------------------------------------------------------------------------------
 594
 595       Int_t trigfired=-1;
 596       Int_t trigside=-1;
 597       if(fDistinguishTrigClass){
 598         TString trigclass = aod->GetFiredTriggerClasses();
 599         printf("TrigClass: %s\n",trigclass.Data());
 600         if(trigclass.Contains(fTrigClassMuon)) trigfired = 1;
 601         else if (trigclass.Contains(fTrigClassInteraction)) trigfired = 0;
 602         printf("Muon container: %d\n",trigfired);
 603         for(Int_t i=0;i<4;i++){
 604           if(trigfired==1 && trigclass.Contains(fTrigClassMuonSide[i])) {trigside=i;printf("Muon Side: %d\n",trigside);}
 605           if(trigfired==0 && trigclass.Contains(fTrigClassInteractionSide[i])) {trigside=i;printf("Interaction Side: %d\n",trigside);}
 606         }
 607       }
 608
 609       ((TH1D*)(fOutput->FindObject("zvSPD")))->Fill(aod->GetPrimaryVertex()->GetZ());
 610       Int_t ncontr = aod->GetPrimaryVertex()->GetNContributors();
 611       ((TH1D*)(fOutput->FindObject("nContrib")))->Fill(ncontr);
 612       if (!fCutOnzVtxSPD || (fCutOnzVtxSPD && (aod->GetPrimaryVertex()->GetZ()>fzPrimVertexSPD[0] && aod->GetPrimaryVertex()->GetZ()<fzPrimVertexSPD[1]))){     //NOT REALLY SPD VERTEX
 613       ((TH1D*)(fOutput->FindObject("zvSPDcut")))->Fill(aod->GetPrimaryVertex()->GetZ());
 614       if (!fCutOnNContributors || (fCutOnNContributors && (aod->GetPrimaryVertex()->GetNContributors()>fNContributors[0] && aod->GetPrimaryVertex()->GetNContributors()<fNContributors[1]))){
 615       for (Int_t j = 0; j<ntracks; j++) {
 616         AliAODTrack *mu1 = aod->GetTrack(j);
 617         if(!mu1->IsMuonTrack()) continue;
 618         if (mu1->Chi2perNDF()<fChi2Track[0] || mu1->Chi2perNDF()>fChi2Track[1]) continue;
 619         if (mu1->GetChi2MatchTrigger()<fChi2MatchTrig[0] || mu1->GetChi2MatchTrigger()>fChi2MatchTrig[1]) continue;
 620         Double_t chargemu1 = mu1->Charge();
 621         Double_t pxmu1 = mu1->Px();
 622         Double_t pymu1 = mu1->Py();
 623         Double_t pzmu1 = mu1->Pz();
 624         Double_t emu1 = mu1->E();
 625         Double_t pmu1 = mu1->P();
 626         Double_t ptmu1 = mu1->Pt();
 627         Double_t rapiditymu1 = Rap(emu1,pzmu1);
 628         Double_t thetamu1 = (180./TMath::Pi())*mu1->Theta();
 629         Double_t rdcamu1 = mu1->DCA();
 630         ((TH1D*)(fOutput->FindObject("emuonREC")))->Fill(emu1);
 631         ((TH1D*)(fOutput->FindObject("ptmuonREC")))->Fill(ptmu1);
 632         ((TH1D*)(fOutput->FindObject("ymuonREC")))->Fill(rapiditymu1);
 633         if(chargemu1<0){
 634           for (Int_t jj = 0; jj<ntracks; jj++) {
 635             AliAODTrack *mu2 = aod->GetTrack(jj);
 636             if(!mu2->IsMuonTrack()) continue;
 637             if (mu2->Chi2perNDF()<fChi2Track[0] || mu2->Chi2perNDF()>fChi2Track[1]) continue;
 638             if (mu2->GetChi2MatchTrigger()<fChi2MatchTrig[0] || mu2->GetChi2MatchTrigger()>fChi2MatchTrig[1]) continue;
 639             Double_t chargemu2 = mu2->Charge();
 640             Double_t pxmu2 = mu2->Px();
 641             Double_t pymu2 = mu2->Py();
 642             Double_t pzmu2 = mu2->Pz();
 643             Double_t emu2 = mu2->E();
 644             Double_t pmu2 = mu2->P();
 645             Double_t ptmu2 = mu2->Pt();
 646             //Double_t rapiditymu2 = Rap(emu2,pzmu2);
 647             Double_t thetamu2 = (180./TMath::Pi())*mu2->Theta();
 648             Double_t rdcamu2 = mu2->DCA();
 649             if(chargemu2>0){
 650             if (mu1->GetMatchTrigger()>0) printf("Mu1: charge: %f, match: %d\n",chargemu1,1);
 651             else  printf("Mu1: charge: %f, match: %d\n",chargemu1,0);
 652             if (mu2->GetMatchTrigger()>0) printf("Mu2: charge: %f, match: %d\n",chargemu2,1);
 653             else printf("Mu2: charge: %f, match: %d\n",chargemu2,0);
 654               ((TH1D*)(fOutput->FindObject("hdca")))->Fill(rdcamu2);
 655               ((TH1D*)(fOutput->FindObject("hdca")))->Fill(rdcamu1);
 656               Double_t trigCondition=0;
 657               if (mu1->GetMatchTrigger()>=mu2->GetMatchTrigger()) trigCondition = mu1->GetMatchTrigger()+0.1*mu2->GetMatchTrigger();
 658                 else trigCondition = mu2->GetMatchTrigger()+0.1*mu1->GetMatchTrigger();
 659               containerInput[0] = (Double_t) Rap((emu1+emu2),(pzmu1+pzmu2));
 660               ((TH1D*)(fOutput->FindObject("ydimuREC")))->Fill(containerInput[0]);
 661               ((TH2D*)(fOutput->FindObject("hdcay")))->Fill(TMath::Max(rdcamu1,rdcamu2),containerInput[0]);
 662               //((TH2D*)(fOutput->FindObject("hdcay")))->Fill(Rdcamu2,containerInput[0]);
 663               containerInput[1] = (Double_t) TMath::Sqrt((pxmu1+pxmu2)*(pxmu1+pxmu2)+(pymu1+pymu2)*(pymu1+pymu2));
 664               ((TH1D*)(fOutput->FindObject("ptdimuREC")))->Fill(containerInput[1]);
 665               containerInput[2] = CostHE(pxmu1,pymu1,pzmu1,emu1,chargemu1,pxmu2,pymu2,pzmu2,emu2);
 666               ((TH1D*)(fOutput->FindObject("costHEdimuREC")))->Fill(containerInput[2]);
 667               if(containerInput[2]==666) ((TH1D*)(fOutput->FindObject("costHEdimuFail")))->Fill(containerInput[2]);
 668               containerInput[3] = TMath::Abs(PhiHE(pxmu1,pymu1,pzmu1,emu1,chargemu1,pxmu2,pymu2,pzmu2,emu2));
 669               ((TH1D*)(fOutput->FindObject("phiHEdimuREC")))->Fill(containerInput[3]);
 670               if(containerInput[3]==666) ((TH1D*)(fOutput->FindObject("phiHEdimuFail")))->Fill(containerInput[3]);
 671               containerInput[4] = CostCS(pxmu1,pymu1,pzmu1,emu1,chargemu1,pxmu2,pymu2,pzmu2,emu2);
 672               ((TH1D*)(fOutput->FindObject("costCSdimuREC")))->Fill(containerInput[4]);
 673               if(containerInput[4]==666) ((TH1D*)(fOutput->FindObject("costCSdimuFail")))->Fill(containerInput[4]);
 674               containerInput[5] = TMath::Abs(PhiCS(pxmu1,pymu1,pzmu1,emu1,chargemu1,pxmu2,pymu2,pzmu2,emu2));
 675               ((TH1D*)(fOutput->FindObject("phiCSdimuREC")))->Fill(containerInput[5]);
 676               if(containerInput[5]==666) ((TH1D*)(fOutput->FindObject("phiCSdimuFail")))->Fill(containerInput[5]);
 677               containerInput[6] = Imass(emu1,pxmu1,pymu1,pzmu1,emu2,pxmu2,pymu2,pzmu2);
 678               ((TH1D*)(fOutput->FindObject("imassTot")))->Fill(containerInput[6]);
 679               containerInput[7] = trigCondition+0.05;
 680               ((TH1D*)(fOutput->FindObject("trigCond")))->Fill(containerInput[7]);
 681               if (ptmu1<ptmu2) {
 682                 containerInput[8]=ptmu1;
 683                 containerInput[9]=ptmu2;
 684               } else {
 685                 containerInput[8]=ptmu2;
 686                 containerInput[9]=ptmu1;
 687               }
 688               if (thetamu1<thetamu2) {
 689                 containerInput[10]=thetamu1;
 690                 containerInput[11]=thetamu2;
 691               } else {
 692                 containerInput[10]=thetamu2;
 693                 containerInput[11]=thetamu1;
 694               }
 695               if (pmu1<pmu2) {
 696                 containerInput[12]=pmu1;
 697                 containerInput[13]=pmu2;
 698               } else {
 699                 containerInput[12]=pmu2;
 700                 containerInput[13]=pmu1;
 701               }
 702               if (fDistinguishTrigClass && trigside==0) containerInput[14]=0.;
 703               else if (fDistinguishTrigClass && trigside==1) containerInput[14]=1.;
 704               else if (fDistinguishTrigClass && trigside==2) containerInput[14]=2.;
 705               else if (fDistinguishTrigClass && trigside==3) containerInput[14]=3.;
 706               else containerInput[14]=0.;
 707
 708               if(containerInput[2]!=666 && containerInput[3]!=666 && containerInput[4]!=666 && containerInput[5]!=666){
 709                 if(fDistinguishTrigClass){
 710                   if (trigfired==1) fCFManager->GetParticleContainer()->Fill(containerInput,5);
 711                   else if (trigfired==0) fCFManager->GetParticleContainer()->Fill(containerInput,1);
 712                 } else {
 713                   fCFManager->GetParticleContainer()->Fill(containerInput,1);
 714                   if(fIsAccProduction){
 715                     if (cutAccept ==1) fCFManager->GetParticleContainer()->Fill(containerInput,3);
 716                   }
 717                 }
 718               } else if (containerInput[0]==666) ((TH1D*)(fOutput->FindObject("ydimuFail")))->Fill(containerInput[0]);
 719
 720           }
 721
 722         }
 723
 724       }
 725
 726     }
 727     }
 728     }
 729
 730
 731   }
 732
 733
 734 //  ----------
 735   if (fReadMCInfo) ((TH1D*)(fOutput->FindObject("jpsiMult")))->Fill(numbJpsis);
 736
 737   PostData(1,fOutput) ;
 738   PostData(2,fCFManager->GetParticleContainer()) ;
 739
 740
 741 }
 742
 743 //________________________________________________________________________
 744 Double_t AliAnalysisTaskDimuonCFContainerBuilder::Imass(Double_t e1, Double_t px1, Double_t py1, Double_t pz1,
 745                                    Double_t e2, Double_t px2, Double_t py2, Double_t pz2) const
 746 {
 747 // invariant mass calculation
 748     Double_t imassrec = TMath::Sqrt((e1+e2)*(e1+e2)-((px1+px2)*(px1+px2)+
 749                                     (py1+py2)*(py1+py2)+(pz1+pz2)*(pz1+pz2)));
 750     return imassrec;
 751 }
 752
 753 //________________________________________________________________________
 754 Double_t AliAnalysisTaskDimuonCFContainerBuilder::Rap(Double_t e, Double_t pz) const
 755 {
 756 // calculate rapidity
 757     Double_t rap;
 758     if(e>TMath::Abs(pz)){                         // in order to avoid problems with AODs
 759         rap = 0.5*TMath::Log((e+pz)/(e-pz));
 760         return rap;
 761     }
 762     else{
 763         rap = 666.;
 764         return rap;
 765     }
 766 }
 767
 768 //________________________________________________________________________
 769 Double_t AliAnalysisTaskDimuonCFContainerBuilder::CostCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
 770 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
 771 {
 772 // Cosine of the theta decay angle (mu+) in the Collins-Soper frame
 773
 774   TLorentzVector pMu1CM, pMu2CM, pProjCM, pTargCM, pDimuCM; // In the CM. frame
 775   TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame
 776   TVector3 beta,zaxisCS;
 777   Double_t mp=0.93827231;
 778   //
 779   // --- Fill the Lorentz vector for projectile and target in the CM frame
 780   //
 781   pProjCM.SetPxPyPzE(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp));
 782   pTargCM.SetPxPyPzE(0.,0.,fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp));
 783   //
 784   // --- Get the muons parameters in the CM frame
 785   //
 786   pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
 787   pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
 788   //
 789   // --- Obtain the dimuon parameters in the CM frame
 790   //
 791   pDimuCM=pMu1CM+pMu2CM;
 792   //
 793   // --- Translate the dimuon parameters in the dimuon rest frame
 794   //
 795   beta=(-1./pDimuCM.E())*pDimuCM.Vect();
 796   if(beta.Mag()>=1) return 666.;
 797   pMu1Dimu=pMu1CM;
 798   pMu2Dimu=pMu2CM;
 799   pProjDimu=pProjCM;
 800   pTargDimu=pTargCM;
 801   pMu1Dimu.Boost(beta);
 802   pMu2Dimu.Boost(beta);
 803   pProjDimu.Boost(beta);
 804   pTargDimu.Boost(beta);
 805
 806   //Debugging part -------------------------------------
 807   Double_t debugProj[4]={0.,0.,0.,0.};
 808   Double_t debugTarg[4]={0.,0.,0.,0.};
 809   Double_t debugMu1[4]={0.,0.,0.,0.};
 810   Double_t debugMu2[4]={0.,0.,0.,0.};
 811   pMu1Dimu.GetXYZT(debugMu1);
 812   pMu2Dimu.GetXYZT(debugMu2);
 813   pProjDimu.GetXYZT(debugProj);
 814   pTargDimu.GetXYZT(debugTarg);
 815   if (debugProj[0]!=debugProj[0] ||debugProj[1]!=debugProj[1] || debugProj[2]!=debugProj[2] ||debugProj[3]!=debugProj[3]) return 666;
 816   if (debugTarg[0]!=debugTarg[0] ||debugTarg[1]!=debugTarg[1] || debugTarg[2]!=debugTarg[2] ||debugTarg[3]!=debugTarg[3]) return 666;
 817   if (debugMu1[0]!=debugMu1[0] ||debugMu1[1]!=debugMu1[1] || debugMu1[2]!=debugMu1[2] ||debugMu1[3]!=debugMu1[3]) return 666;
 818   if (debugMu2[0]!=debugMu2[0] ||debugMu2[1]!=debugMu2[1] || debugMu2[2]!=debugMu2[2] ||debugMu2[3]!=debugMu2[3]) return 666;
 819   //----------------------------------------------------
 820
 821   // --- Determine the z axis for the CS angle
 822   zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit();
 823
 824   // --- Determine the CS angle (angle between mu+ and the z axis defined above)
 825   Double_t cost;
 826
 827   if(charge1>0) {cost = zaxisCS.Dot((pMu1Dimu.Vect()).Unit());}
 828   else {cost = zaxisCS.Dot((pMu2Dimu.Vect()).Unit());}
 829
 830   return cost;
 831 }
 832
 833 //________________________________________________________________________
 834 Double_t AliAnalysisTaskDimuonCFContainerBuilder::CostHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
 835 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
 836 {
 837 // Cosine of the theta decay angle (mu+) in the Helicity frame
 838
 839   TLorentzVector pMu1CM, pMu2CM, pDimuCM; // In the CM frame
 840   TLorentzVector pMu1Dimu, pMu2Dimu; // In the dimuon rest frame
 841   TVector3 beta,zaxisCS;
 842   //
 843   // --- Get the muons parameters in the CM frame
 844   //
 845   pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
 846   pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
 847   //
 848   // --- Obtain the dimuon parameters in the CM frame
 849   //
 850   pDimuCM=pMu1CM+pMu2CM;
 851   //
 852   // --- Translate the muon parameters in the dimuon rest frame
 853   //
 854   beta=(-1./pDimuCM.E())*pDimuCM.Vect();
 855   if(beta.Mag()>=1) return 666.;
 856   pMu1Dimu=pMu1CM;
 857   pMu2Dimu=pMu2CM;
 858   pMu1Dimu.Boost(beta);
 859   pMu2Dimu.Boost(beta);
 860
 861   //Debugging part -------------------------------------
 862   Double_t debugMu1[4]={0.,0.,0.,0.};
 863   Double_t debugMu2[4]={0.,0.,0.,0.};
 864   pMu1Dimu.GetXYZT(debugMu1);
 865   pMu2Dimu.GetXYZT(debugMu2);
 866   if (debugMu1[0]!=debugMu1[0] ||debugMu1[1]!=debugMu1[1] || debugMu1[2]!=debugMu1[2] ||debugMu1[3]!=debugMu1[3]) return 666;
 867   if (debugMu2[0]!=debugMu2[0] ||debugMu2[1]!=debugMu2[1] || debugMu2[2]!=debugMu2[2] ||debugMu2[3]!=debugMu2[3]) return 666;
 868   //----------------------------------------------------
 869
 870   // --- Determine the z axis for the calculation of the polarization angle (i.e. the direction of the dimuon in the CM system)
 871   TVector3 zaxis;
 872   zaxis=(pDimuCM.Vect()).Unit();
 873
 874   // --- Calculation of the polarization angle (angle between mu+ and the z axis defined above)
 875   Double_t cost;
 876   if(charge1>0) {cost = zaxis.Dot((pMu1Dimu.Vect()).Unit());}
 877   else {cost = zaxis.Dot((pMu2Dimu.Vect()).Unit());}
 878   return cost;
 879 }
 880
 881 //________________________________________________________________________
 882 Double_t AliAnalysisTaskDimuonCFContainerBuilder::PhiCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
 883 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
 884 {
 885 // Phi decay angle (mu+) in the Collins-Soper frame
 886
 887    TLorentzVector pMu1CM, pMu2CM, pProjCM, pTargCM, pDimuCM; // In the CM frame
 888    TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame
 889    TVector3 beta,yaxisCS, xaxisCS, zaxisCS;
 890    Double_t mp=0.93827231;
 891
 892    // --- Fill the Lorentz vector for projectile and target in the CM frame
 893    pProjCM.SetPxPyPzE(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp));
 894    pTargCM.SetPxPyPzE(0.,0.,fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp));
 895
 896    // --- Get the muons parameters in the CM frame
 897    pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
 898    pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
 899
 900    // --- Obtain the dimuon parameters in the CM frame
 901    pDimuCM=pMu1CM+pMu2CM;
 902
 903    // --- Translate the dimuon parameters in the dimuon rest frame
 904    beta=(-1./pDimuCM.E())*pDimuCM.Vect();
 905    if(beta.Mag()>=1) return 666.;
 906    pMu1Dimu=pMu1CM;
 907    pMu2Dimu=pMu2CM;
 908    pProjDimu=pProjCM;
 909    pTargDimu=pTargCM;
 910    pMu1Dimu.Boost(beta);
 911    pMu2Dimu.Boost(beta);
 912    pProjDimu.Boost(beta);
 913    pTargDimu.Boost(beta);
 914
 915    //Debugging part -------------------------------------
 916    Double_t debugProj[4]={0.,0.,0.,0.};
 917    Double_t debugTarg[4]={0.,0.,0.,0.};
 918    Double_t debugMu1[4]={0.,0.,0.,0.};
 919    Double_t debugMu2[4]={0.,0.,0.,0.};
 920    pMu1Dimu.GetXYZT(debugMu1);
 921    pMu2Dimu.GetXYZT(debugMu2);
 922    pProjDimu.GetXYZT(debugProj);
 923    pTargDimu.GetXYZT(debugTarg);
 924    if (debugProj[0]!=debugProj[0] ||debugProj[1]!=debugProj[1] || debugProj[2]!=debugProj[2] ||debugProj[3]!=debugProj[3]) return 666;
 925    if (debugTarg[0]!=debugTarg[0] ||debugTarg[1]!=debugTarg[1] || debugTarg[2]!=debugTarg[2] ||debugTarg[3]!=debugTarg[3]) return 666;
 926    if (debugMu1[0]!=debugMu1[0] ||debugMu1[1]!=debugMu1[1] || debugMu1[2]!=debugMu1[2] ||debugMu1[3]!=debugMu1[3]) return 666;
 927    if (debugMu2[0]!=debugMu2[0] ||debugMu2[1]!=debugMu2[1] || debugMu2[2]!=debugMu2[2] ||debugMu2[3]!=debugMu2[3]) return 666;
 928    //----------------------------------------------------
 929
 930    // --- Determine the z axis for the CS angle
 931    zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit();
 932    yaxisCS=(((pProjDimu.Vect()).Unit()).Cross((pTargDimu.Vect()).Unit())).Unit();
 933    xaxisCS=(yaxisCS.Cross(zaxisCS)).Unit();
 934
 935    Double_t phi=0.;
 936    if(charge1>0) {
 937        phi = TMath::ATan2((pMu1Dimu.Vect()).Dot(yaxisCS),((pMu1Dimu.Vect()).Dot(xaxisCS)));
 938    } else {
 939        phi = TMath::ATan2((pMu2Dimu.Vect()).Dot(yaxisCS),((pMu2Dimu.Vect()).Dot(xaxisCS)));
 940    }
 941    if (phi>TMath::Pi()) phi=phi-TMath::Pi();
 942
 943    return phi;
 944 }
 945
 946 //________________________________________________________________________
 947 Double_t AliAnalysisTaskDimuonCFContainerBuilder::PhiHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
 948 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
 949 {
 950 // Phi decay angle (mu+) in the Helicity frame
 951   TLorentzVector pMu1Lab, pMu2Lab, pProjLab, pTargLab, pDimuLab; // In the lab. frame
 952   TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame
 953   TVector3 beta,xaxis, yaxis,zaxis;
 954   Double_t mp=0.93827231;
 955
 956   // --- Get the muons parameters in the LAB frame
 957   pMu1Lab.SetPxPyPzE(px1,py1,pz1,e1);
 958   pMu2Lab.SetPxPyPzE(px2,py2,pz2,e2);
 959
 960   // --- Obtain the dimuon parameters in the LAB frame
 961   pDimuLab=pMu1Lab+pMu2Lab;
 962   zaxis=(pDimuLab.Vect()).Unit();
 963
 964   // --- Translate the muon parameters in the dimuon rest frame
 965   beta=(-1./pDimuLab.E())*pDimuLab.Vect();
 966   if(beta.Mag()>=1.) return 666.;
 967
 968   pProjLab.SetPxPyPzE(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp)); // proiettile
 969   pTargLab.SetPxPyPzE(0.,0.,fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp)); // bersaglio
 970
 971   pProjDimu=pProjLab;
 972   pTargDimu=pTargLab;
 973
 974   pProjDimu.Boost(beta);
 975   pTargDimu.Boost(beta);
 976
 977   yaxis=((pProjDimu.Vect()).Cross(pTargDimu.Vect())).Unit();
 978   xaxis=(yaxis.Cross(zaxis)).Unit();
 979
 980   pMu1Dimu=pMu1Lab;
 981   pMu2Dimu=pMu2Lab;
 982   pMu1Dimu.Boost(beta);
 983   pMu2Dimu.Boost(beta);
 984
 985   //Debugging part -------------------------------------
 986   Double_t debugProj[4]={0.,0.,0.,0.};
 987   Double_t debugTarg[4]={0.,0.,0.,0.};
 988   Double_t debugMu1[4]={0.,0.,0.,0.};
 989   Double_t debugMu2[4]={0.,0.,0.,0.};
 990   pMu1Dimu.GetXYZT(debugMu1);
 991   pMu2Dimu.GetXYZT(debugMu2);
 992   pProjDimu.GetXYZT(debugProj);
 993   pTargDimu.GetXYZT(debugTarg);
 994   if (debugProj[0]!=debugProj[0] ||debugProj[1]!=debugProj[1] || debugProj[2]!=debugProj[2] ||debugProj[3]!=debugProj[3]) return 666;
 995   if (debugTarg[0]!=debugTarg[0] ||debugTarg[1]!=debugTarg[1] || debugTarg[2]!=debugTarg[2] ||debugTarg[3]!=debugTarg[3]) return 666;
 996   if (debugMu1[0]!=debugMu1[0] ||debugMu1[1]!=debugMu1[1] || debugMu1[2]!=debugMu1[2] ||debugMu1[3]!=debugMu1[3]) return 666;
 997   if (debugMu2[0]!=debugMu2[0] ||debugMu2[1]!=debugMu2[1] || debugMu2[2]!=debugMu2[2] ||debugMu2[3]!=debugMu2[3]) return 666;
 998   //----------------------------------------------------
 999
1000   Double_t phi=0.;
1001    if(charge1 > 0) {
1002       phi = TMath::ATan2((pMu1Dimu.Vect()).Dot(yaxis),(pMu1Dimu.Vect()).Dot(xaxis));
1003      } else {
1004       phi = TMath::ATan2((pMu2Dimu.Vect()).Dot(yaxis),(pMu2Dimu.Vect()).Dot(xaxis));
1005    }
1006    return phi;
1007 }
1008
1009 //________________________________________________________________________
1010 void AliAnalysisTaskDimuonCFContainerBuilder::Terminate(Option_t *)
1011 {
1012
1013 }
1014
1015 #endif