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