PWG3/muon/AliAnalysisTaskDimuonCFContainerBuilder.cxx

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