New Analysis task for upsilon anlaysis (Sang-un)
[u/mrichter/AliRoot.git] / PWG3 / muon / AliCFMuonResUpsilon.cxx
CommitLineData
014044ba 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16//-----------------------------------------------------------------------
17// Example of task (running locally, on AliEn and CAF),
18// which provides standard way of calculating acceptance and efficiency
19// between different steps of the procedure.
20// The ouptut of the task is a AliCFContainer from which the efficiencies
21// can be calculated
22//-----------------------------------------------------------------------
23// Author : R. Vernet, Consorzio Cometa - Catania (it)
24//-----------------------------------------------------------------------
25// Modification done by X. Lopez - LPC Clermont (fr)
26//-----------------------------------------------------------------------
27// Modification done by S. Ahn - LPC Clermont (fr), Konkuk University (kr)
28//-----------------------------------------------------------------------
29
30
31
32#ifndef ALICFMUONRESUPSILON_CXX
33#define ALICFMUONRESUPSILON_CXX
34
35#include "TH1.h"
36#include "TParticle.h"
37#include "TChain.h"
38#include "TLorentzVector.h"
39#include "TCanvas.h"
40
41#include "AliHeader.h"
42#include "AliStack.h"
43#include "AliMCEvent.h"
44#include "AliAnalysisManager.h"
45#include "AliLog.h"
46
47#include "AliVEvent.h"
48#include "AliPhysicsSelection.h"
49
50#include "AliCFMuonResUpsilon.h"
51#include "AliCFManager.h"
52#include "AliCFCutBase.h"
53#include "AliCFContainer.h"
54
55#include "AliESDEvent.h"
56#include "AliESDHeader.h"
57#include "AliESDtrack.h"
58#include "AliESDMuonTrack.h"
59#include "AliESDtrack.h"
60#include "AliESDInputHandler.h"
61
62#include "AliAODEvent.h"
63#include "AliAODInputHandler.h"
64#include "AliAODMCParticle.h"
65
66ClassImp(AliCFMuonResUpsilon)
67
68//__________________________________________________________________________
69AliCFMuonResUpsilon::AliCFMuonResUpsilon() :
70 AliAnalysisTaskSE(""),
71 fCFManager(0x0),
72 hnevts(0x0),
73 fIsPhysSelMB(kFALSE),
74 fIsPhysSelMUON(kFALSE),
75 fQAHistList(0x0),
76 fPDG(0),
77 fPtMin(0),
78 fPtMax(0),
79 fYMin(0),
80 fYMax(0),
81 fTrigClassMuon(""),
82 fTrigClassInteraction(""),
83 fDistinguishTrigClass(kTRUE)
84{
85
86 //Default ctor
87
88}
89//___________________________________________________________________________
90AliCFMuonResUpsilon::AliCFMuonResUpsilon(const Char_t* name) :
91 AliAnalysisTaskSE(name),
92 fCFManager(0x0),
93 hnevts(0x0),
94 fIsPhysSelMB(kFALSE),
95 fIsPhysSelMUON(kFALSE),
96 fQAHistList(0x0),
97 fPDG(0),
98 fPtMin(0),
99 fPtMax(0),
100 fYMin(0),
101 fYMax(0),
102 fTrigClassMuon(""),
103 fTrigClassInteraction(""),
104 fDistinguishTrigClass(kTRUE)
105{
106
107 // Constructor. Initialization of Inputs and Outputs
108
109 Info("AliCFMuonResUpsilon","Calling Constructor");
110
111 hnevts = new TH1D("hnevts","hnevts",5,0,5); // nevent CINT1B,CMUS1B + PhysSel
112
113 TString nameside[3]={"AC","B","E"};
114
115 SetTrigClassMuonName();
116 SetTrigClassInteracName();
117 SetTrigClassSideName(nameside);
118
119 DefineInput(0, TChain::Class());
120 DefineOutput(1,TH1D::Class());
121 DefineOutput(2,AliCFContainer::Class());
122}
123
124//___________________________________________________________________________
125AliCFMuonResUpsilon& AliCFMuonResUpsilon::operator=(const AliCFMuonResUpsilon& c)
126{
127
128 // Assignment operator
129
130 if (this!=&c) {
131 AliAnalysisTaskSE::operator=(c);
132 fReadAODData = c.fReadAODData;
133 fCFManager = c.fCFManager;
134 fQAHistList = c.fQAHistList;
135 hnevts = c.hnevts;
136 fPDG = c.fPDG;
137 fPtMin = c.fPtMin;
138 fPtMax = c.fPtMax;
139 fYMin = c.fYMin;
140 fYMax = c.fYMax;
141 }
142 return *this;
143}
144
145//___________________________________________________________________________
146AliCFMuonResUpsilon::AliCFMuonResUpsilon(const AliCFMuonResUpsilon& c) :
147 AliAnalysisTaskSE(c),
148 fReadAODData(c.fReadAODData),
149 fReadMCInfo(c.fReadMCInfo),
150 fCFManager(c.fCFManager),
151 hnevts(c.hnevts),
152 fIsPhysSelMB(kFALSE),
153 fIsPhysSelMUON(kFALSE),
154 fQAHistList(c.fQAHistList),
155 fPDG(c.fPDG),
156 fPtMin(c.fPtMin),
157 fPtMax(c.fPtMax),
158 fYMin(c.fYMin),
159 fYMax(c.fYMax),
160 fTrigClassMuon(c.fTrigClassMuon),
161 fTrigClassInteraction(c.fTrigClassInteraction),
162 fDistinguishTrigClass(c.fDistinguishTrigClass)
163{
164
165 // Copy Constructor
166
167}
168
169//___________________________________________________________________________
170AliCFMuonResUpsilon::~AliCFMuonResUpsilon() {
171
172 //destructor
173
174 Info("~AliCFMuonResUpsilon","Calling Destructor");
175 if (fCFManager) delete fCFManager ;
176 if (hnevts) delete hnevts ;
177 if (fQAHistList) {fQAHistList->Clear(); delete fQAHistList;}
178}
179//___________________________________________________________________________
180void AliCFMuonResUpsilon::UserCreateOutputObjects() {
181 // UserCreateOutputObjects
182
183 //TH1D *hnevts = new TH1D("hnevts","hnevts",5,0,5); // nevent CINT1B,CMUS1B + PhysSel
184
185 PostData(1,hnevts) ;
186 PostData(2,fCFManager->GetParticleContainer()) ;
187}
188
189//_________________________________________________
190void AliCFMuonResUpsilon::UserExec(Option_t *)
191{
192
193 // Main loop function
194
195 // variables
196 fIsPhysSelMB=kFALSE; // physics selection : MB
197 fIsPhysSelMUON=kFALSE; // physics selection : MUON
198
199 Double_t containerInput[14] = {0,} ;
200
201 hnevts->Fill(0.5);
202 containerInput[0] = 0.5;
203
204 if(!fReadAODData) { // ESD-based ANALYSIS
205
206 Bool_t flag=kFALSE;
207
208 if(fReadMCInfo) {
209 if (!fMCEvent) {
210 Error("UserExec","NO MC EVENT FOUND!");
211 return;
212 }
213
214 // MC part ----------------------------------------------------------------------------------
215 fCFManager->SetMCEventInfo(fMCEvent);
216 AliStack *stack = fMCEvent->Stack();
217 Int_t npart=fMCEvent->GetNumberOfTracks();
218
219 //printf("ESD: npart=%d\n", npart);
220
221 for (Int_t ipart=0; ipart<npart; ipart++) {
222
223 AliMCParticle *mcPart = (AliMCParticle*)fMCEvent->GetTrack(ipart);
224
225 // Mother kinematics
226 TParticle *part = mcPart->Particle();
227 Double_t e = part->Energy();
228 Double_t pz = part->Pz();
229 Double_t rapmc = Rap(e,pz);
230
231 // Selection of the resonance
232 if(!fCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,mcPart)) continue;
233 if(part->Pt()<0 || part->Pt()>1000) continue;
234 if(rapmc<-4 || rapmc>-2.4) continue;
235
236 // Decays kinematics
237 Int_t p0 = part->GetDaughter(0);
238 TParticle *part0 = stack->Particle(p0);
239 AliMCParticle *mcpart0 = new AliMCParticle(part0);
240 Int_t pdg0 = part0->GetPdgCode();
241
242 Int_t p1 = part->GetDaughter(1);
243 TParticle *part1 = stack->Particle(p1);
244 AliMCParticle *mcpart1 = new AliMCParticle(part1);
245 Int_t pdg1 = part1->GetPdgCode();
246
247 Double_t e0 = part0->Energy();
248 Double_t pz0 = part0->Pz();
249 Double_t py0 = part0->Py();
250 Double_t px0 = part0->Px();
251 Double_t charge0 = part0->GetPDG()->Charge()/3;
252 Double_t pt0 = TMath::Sqrt(px0*px0+py0*py0);
253 Double_t theta0 = (180./TMath::Pi())*TMath::ATan2(TMath::Sqrt(px0*px0+py0*py0),pz0);
254 Double_t eta0 = part0->Eta();
255 Double_t vz0 = part0->Vz();
256 Double_t mom0 = part0->P();
257
258 if(!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mcpart0)) continue;
259 if(pt0<0.5) continue;
260 if(theta0<171. || theta0>178.) continue;
261
262 Double_t e1 = part1->Energy();
263 Double_t pz1 = part1->Pz();
264 Double_t py1 = part1->Py();
265 Double_t px1 = part1->Px();
266 Double_t charge1 = part1->GetPDG()->Charge()/3;
267 Double_t pt1 = TMath::Sqrt(px1*px1+py1*py1);
268 Double_t theta1 = (180./TMath::Pi())*TMath::ATan2(TMath::Sqrt(px1*px1+py1*py1),pz1);
269 Double_t eta1 = part1->Eta();
270 Double_t vz1 = part1->Vz();
271 Double_t mom1 = part1->P();
272
273 if(!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mcpart1)) continue;
274 if(pt1<0.5) continue;
275 if(theta1<171. || theta1>178.) continue;
276
277 if(pdg0==-13 || pdg1==-13) {
278
279 Double_t ptmc = TMath::Sqrt((px0+px1)*(px0+px1)+(py0+py1)*(py0+py1));
280 Double_t imassmc = Imass(e0,px0,py0,pz0,e1,px1,py1,pz1);
281 if(!imassmc) continue;
282
283 containerInput[1] = rapmc ;
284 containerInput[2] = ptmc ;
285 containerInput[3] = imassmc ;
286 containerInput[4] = 1. ;
287 containerInput[5] = pt0 ;
288 containerInput[5] = pt1 ;
289 containerInput[6] = mom0 ;
290 containerInput[6] = mom1 ;
291 containerInput[7] = 1. ;
292 containerInput[8] = 0. ;
293 containerInput[8] = 0. ;
294 containerInput[9] = charge0+charge1 ;
295 containerInput[10] = eta0;
296 containerInput[10] = eta1;
297 containerInput[11] = theta0;
298 containerInput[11] = theta1;
299 containerInput[12] = vz0;
300 containerInput[12] = vz1;
301 containerInput[13] = 0;
302 containerInput[13] = 0;
303
304 // fill the container at the first step
305 fCFManager->GetParticleContainer()->Fill(containerInput,0); // MC container
306 flag=kTRUE;
307 } // second mu loop
308 } // first mu loop
309 } // end fReadMCInfo
310
311 if(fReadMCInfo && !flag) return;
312 // RECO : ESD part ----------------------------------------------------------------------------------
313
314 // ESD handler
315 AliESDEvent *fESD = 0x0;
316
317 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
318 if( ! esdH) {
319 AliError("Cannot get input event handler");
320 return;
321 }
322 fESD = esdH->GetEvent();
323
324
325 Int_t ntrk=fESD->GetNumberOfMuonTracks();
326
327 Int_t trigfired=-1;
328 Int_t trigside=-1;
329
330 if(!fReadMCInfo) {
331 fIsPhysSelMB=(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kMB));
332 fIsPhysSelMUON=(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kMUON));
333
334 if(fDistinguishTrigClass) {
335 TString trigclass = fESD->GetFiredTriggerClasses();
336 if(trigclass.Contains(fTrigClassMuon)) trigfired = 1;
337 else if(trigclass.Contains(fTrigClassInteraction)) trigfired = 0;
338 for(Int_t i=0;i<3;i++) {
339 if(trigfired==1 && trigclass.Contains(fTrigClassSide[i])) {trigside=i;}
340 if(trigfired==0 && trigclass.Contains(fTrigClassSide[i])) {trigside=i;}
341 }
342 }
343
344 if(trigside==1) { // Beam-Beam event
345 if(trigfired) { // CMUS1-B
346 containerInput[0]=3.5;
347 hnevts->Fill(3.5);
348 if(fIsPhysSelMUON) {
349 containerInput[0]=4.5;
350 hnevts->Fill(4.5); // PhysSel
351 }
352 } else { // CINT1-B
353 containerInput[0]=1.5;
354 hnevts->Fill(1.5);
355 if(fIsPhysSelMB) {
356 containerInput[0]=2.5;
357 hnevts->Fill(2.5); // PhysSel
358 }
359 }
360 }
361 }
362 else trigside = 0; // for MC
363
364 for(Int_t j=0; j<ntrk-1; j++) {
365 for(Int_t jj=j+1; jj<ntrk; jj++) {
366 AliESDMuonTrack *mu1 = new AliESDMuonTrack(*(fESD->GetMuonTrack(j)));
367 AliESDMuonTrack *mu2 = new AliESDMuonTrack(*(fESD->GetMuonTrack(jj)));
368
369 Double_t trigCondition = 0;
370 if (mu1->GetMatchTrigger()>=mu2->GetMatchTrigger()) trigCondition = mu1->GetMatchTrigger()+0.1*mu2->GetMatchTrigger();
371 else trigCondition = mu2->GetMatchTrigger()+0.1*mu1->GetMatchTrigger();
372
373 // mu1
374 Double_t zr1 = mu1->Charge();
375 Double_t pxr1 = mu1->Px();
376 Double_t pyr1 = mu1->Py();
377 Double_t pzr1 = mu1->Pz();
378 Double_t ptr1 = TMath::Sqrt(pxr1*pxr1+pyr1*pyr1);
379 Double_t er1 = mu1->E();
380 Double_t momr1 = mu1->P();
381 Double_t vzr1 = mu1->GetZ();
382 Double_t dcar1 = mu1->GetDCA();
383 Double_t rabsr1 = mu1->GetRAtAbsorberEnd();
384 Double_t theta_absr1 = 180.*(1.-TMath::ATan(rabsr1/505.)/TMath::Pi()); // [deg]
385 Double_t theta_absr1_2 = (theta_absr1*TMath::Pi()/180.)/2.;
386 Double_t tan_theta_absr1 = TMath::Tan(theta_absr1_2);
387 Double_t etar1 = -TMath::Log(tan_theta_absr1);
388
389 if(Rap(er1,pzr1)<-4 || Rap(er1,pzr1)>-2.5) continue;
390
391 // mu2
392 Double_t zr2 = mu2->Charge();
393 Double_t pxr2 = mu2->Px();
394 Double_t pyr2 = mu2->Py();
395 Double_t pzr2 = mu2->Pz();
396 Double_t ptr2 = TMath::Sqrt(pxr2*pxr2+pyr2*pyr2);
397 Double_t er2 = mu2->E();
398 Double_t momr2 = mu2->P();
399 Double_t vzr2 = mu2->GetZ();
400 Double_t dcar2 = mu2->GetDCA();
401 Double_t rabsr2 = mu2->GetRAtAbsorberEnd();
402 Double_t theta_absr2 = 180.*(1.-TMath::ATan(rabsr2/505.)/TMath::Pi()); // [deg]
403 Double_t theta_absr2_2 = (theta_absr2*TMath::Pi()/180.)/2.;
404 Double_t tan_theta_absr2 = TMath::Tan(theta_absr2_2);
405 Double_t etar2 = -TMath::Log(tan_theta_absr2);
406
407 if(Rap(er2,pzr2)<-4 || Rap(er2,pzr2)>-2.5) continue;
408
409 // for MC
410 if(fReadMCInfo) {
411 //mu1
412 rabsr1 = 0;
413 Double_t thetar1 = (180./TMath::Pi())*TMath::ATan2(TMath::Sqrt(pxr1*pxr1+pyr1*pyr1),pzr1);
414 theta_absr1 = thetar1;
415 etar1 = mu1->Eta();
416 //mu2
417 rabsr2 = 0;
418 Double_t thetar2 = (180./TMath::Pi())*TMath::ATan2(TMath::Sqrt(pxr2*pxr2+pyr2*pyr2),pzr2);
419 theta_absr2 = thetar2;
420 etar2 = mu2->Eta();
421 }
422
423 if(TMath::Abs(etar1) > 8 || TMath::Abs(etar2) > 8) continue;
424
425 // dimuon
426 Double_t ptrec = TMath::Sqrt((pxr1+pxr2)*(pxr1+pxr2)+(pyr1+pyr2)*(pyr1+pyr2));
427 Double_t raprec= Rap((er1+er2),(pzr1+pzr2));
428 Double_t imassrec = Imass(er1,pxr1,pyr1,pzr1,er2,pxr2,pyr2,pzr2);
429 if(!imassrec) continue;
430
431 containerInput[1] = raprec ;
432 containerInput[2] = ptrec ;
433 containerInput[3] = imassrec ;
434 containerInput[4] = trigCondition+0.05 ;
435 containerInput[5] = ptr1 ;
436 containerInput[5] = ptr2 ;
437 containerInput[6] = momr1;
438 containerInput[6] = momr2;
439 containerInput[7] = trigside;
440 containerInput[8] = rabsr1;
441 containerInput[8] = rabsr2;
442 containerInput[9] = zr1 + zr2;
443 containerInput[10] = etar1;
444 containerInput[10] = etar2;
445 containerInput[11] = theta_absr1;
446 containerInput[11] = theta_absr2;
447 containerInput[12] = vzr1;
448 containerInput[12] = vzr2;
449 containerInput[13] = dcar1;
450 containerInput[13] = dcar2;
451
452 if (trigfired==1 && trigside==1) fCFManager->GetParticleContainer()->Fill(containerInput,4);
453 if (trigfired==0 && trigside==1) fCFManager->GetParticleContainer()->Fill(containerInput,1);
454 if(fReadMCInfo) fCFManager->GetParticleContainer()->Fill(containerInput,1); // Rec container
455 } // second mu loop
456 } // first mu loop
457 } // end ESD-based ANALYSIS
458
459 else { // AOD-based ANALYSIS
460
461 Bool_t flag=kTRUE;
462
463 // AOD handler
464 AliAODEvent *fAOD = 0x0;
465
466 AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
467 fAOD = aodH->GetEvent();
468 if( ! aodH) {
469 AliError("Cannot get input event handler");
470 return;
471 }
472
473 // MC part -----------------------------------------------------------------------------------
474
475 if(fReadMCInfo) {
476 TClonesArray *mcArr = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
477
478 Int_t npart=mcArr->GetEntries();
479
480 for(Int_t i=0; i<npart; i++) {
481 AliAODMCParticle *mctrack = (AliAODMCParticle*) mcArr->At(i);
482 // select resonance
483 if(mctrack->GetPdgCode()!=fPDG) continue;
484 // cuts on resonance
485 if(!(mctrack->Pt()>0 && mctrack->Pt()<1000) || !(Rap(mctrack->E(),mctrack->Pz())>-4 && Rap(mctrack->E(),mctrack->Pz())<-2.4)) continue;
486 Int_t daug0 = mctrack->GetDaughter(0);
487 Int_t daug1 = mctrack->GetDaughter(1);
488 // daughter1
489 AliAODMCParticle *mcdaug0 = (AliAODMCParticle*) mcArr->At(daug0);
490 Double_t pt0 = mcdaug0->Pt();
491 Double_t mom0 = mcdaug0->P();
492 Double_t theta0 = (180./TMath::Pi())*mcdaug0->Theta();
493 Double_t eta0 = mcdaug0->Eta();
494 Int_t charge0 = (Int_t) mcdaug0->Charge()/3;
495 Double_t vz0 = mcdaug0->Zv();
496 if(!(pt0>0.5) || !(theta0>171. && theta0<178.)) continue;
497 // daughter2
498 AliAODMCParticle *mcdaug1 = (AliAODMCParticle*) mcArr->At(daug1);
499 Double_t pt1 = mcdaug1->Pt();
500 Double_t mom1 = mcdaug1->P();
501 Double_t theta1 = (180./TMath::Pi())*mcdaug1->Theta();
502 Double_t eta1 = mcdaug1->Eta();
503 Int_t charge1 = (Int_t) mcdaug1->Charge()/3;
504 Double_t vz1 = mcdaug1->Zv();
505 if(!(pt1>0.5) || !(theta1>171. && theta1<178.)) continue;
506
507 if (TMath::Abs(mcdaug0->GetPdgCode())!=13 || TMath::Abs(mcdaug1->GetPdgCode())!=13) continue;
508
509 Double_t rapmc = Rap(mctrack->E(), mctrack->Pz());
510 Double_t ptmc = mctrack->Pt();
511 Double_t imassmc = mctrack->M();
512 if(!imassmc) continue;
513
514 containerInput[1] = rapmc ;
515 containerInput[2] = ptmc ;
516 containerInput[3] = imassmc ;
517 containerInput[4] = 1. ;
518 containerInput[5] = pt0 ;
519 containerInput[5] = pt1 ;
520 containerInput[6] = mom0 ;
521 containerInput[6] = mom1 ;
522 containerInput[7] = 1. ;
523 containerInput[8] = 0. ;
524 containerInput[8] = 0. ;
525 containerInput[9] = charge0+charge1 ;
526 containerInput[10] = eta0;
527 containerInput[10] = eta1;
528 containerInput[11] = theta0;
529 containerInput[11] = theta1;
530 containerInput[12] = vz0;
531 containerInput[12] = vz1;
532 containerInput[13] = 0;
533 containerInput[13] = 0;
534
535 fCFManager->GetParticleContainer()->Fill(containerInput,0); // MC container
536 flag=kTRUE;
537 }
538 }
539
540if(fReadMCInfo && !flag) return;
541
542 // RECO : AOD part ----------------------------------------------------------------------------------
543
544
545 Int_t trigfired=-1;
546 Int_t trigside=-1;
547
548 Int_t ntrk = fAOD->GetNumberOfTracks();
549
550 if(!fReadMCInfo) {
551 fIsPhysSelMB=(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kMB));
552 fIsPhysSelMUON=(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kMUON));
553
554 if(fDistinguishTrigClass) {
555 TString trigclass = fAOD->GetFiredTriggerClasses();
556 if(trigclass.Contains(fTrigClassMuon)) trigfired = 1;
557 else if(trigclass.Contains(fTrigClassInteraction)) trigfired = 0;
558 for(Int_t i=0;i<3;i++) {
559 if(trigfired==1 && trigclass.Contains(fTrigClassSide[i].Data())) {trigside=i;}
560 if(trigfired==0 && trigclass.Contains(fTrigClassSide[i].Data())) {trigside=i;}
561 }
562 }
563
564 if(trigside==1) { // Beam-Beam event
565 if(trigfired) { // CMUS1-B
566 containerInput[0]=3.5;
567 hnevts->Fill(3.5);
568 if(fIsPhysSelMUON) {
569 containerInput[0]=4.5;
570 hnevts->Fill(4.5); // PhysSel
571 }
572 } else { // CINT1-B
573 containerInput[0]=1.5;
574 hnevts->Fill(1.5);
575 if(fIsPhysSelMB) {
576 containerInput[0]=2.5;
577 hnevts->Fill(2.5); // PhysSel
578 }
579 }
580 }
581 }
582 else trigside = 0; // for MC
583
584
585 for(Int_t j=0; j<ntrk; j++) {
586 for(Int_t jj=j+1; jj<ntrk; jj++) {
587 AliAODTrack *mu1 = fAOD->GetTrack(j);
588 if(!mu1->IsMuonTrack() || !(mu1->Y()>-4 && mu1->Y()<-2.5)) continue;
589 AliAODTrack *mu2 = fAOD->GetTrack(jj);
590 if(!mu2->IsMuonTrack() || !(mu2->Y()>-4 && mu2->Y()<-2.5)) continue;
591
592 Double_t trigCondition=0;
593 if (mu1->GetMatchTrigger()>=mu2->GetMatchTrigger()) trigCondition = mu1->GetMatchTrigger()+0.1*mu2->GetMatchTrigger();
594 else trigCondition = mu2->GetMatchTrigger()+0.1*mu1->GetMatchTrigger();
595
596 // mu1
597 Double_t zr1 = mu1->Charge();
598 Double_t pxr1 = mu1->Px();
599 Double_t pyr1 = mu1->Py();
600 Double_t pzr1 = mu1->Pz();
601 Double_t ptr1 = TMath::Sqrt(pxr1*pxr1+pyr1*pyr1);
602 Double_t er1 = mu1->E();
603 Double_t momr1 = mu1->P();
604 Double_t vzr1 = mu1->Zv();
605 Double_t dcar1 = mu1->DCA();
606 Double_t rabsr1 = mu1->GetRAtAbsorberEnd();
607 Double_t theta_absr1 = 180.*(1.-TMath::ATan(rabsr1/505.)/TMath::Pi()); // [deg]
608 Double_t theta_absr1_2 = (theta_absr1*TMath::Pi()/180.)/2.;
609 Double_t tan_theta_absr1 = TMath::Tan(theta_absr1_2);
610 Double_t etar1 = -TMath::Log(tan_theta_absr1);
611
612 // mu2
613 Double_t zr2 = mu2->Charge();
614 Double_t pxr2 = mu2->Px();
615 Double_t pyr2 = mu2->Py();
616 Double_t pzr2 = mu2->Pz();
617 Double_t ptr2 = TMath::Sqrt(pxr2*pxr2+pyr2*pyr2);
618 Double_t er2 = mu2->E();
619 Double_t momr2 = mu2->P();
620 Double_t vzr2 = mu2->Zv();
621 Double_t dcar2 = mu2->DCA();
622 Double_t rabsr2 = mu1->GetRAtAbsorberEnd();
623 Double_t theta_absr2 = 180.*(1.-TMath::ATan(rabsr2/505.)/TMath::Pi()); // [deg]
624 Double_t theta_absr2_2 = (theta_absr2*TMath::Pi()/180.)/2.;
625 Double_t tan_theta_absr2 = TMath::Tan(theta_absr2_2);
626 Double_t etar2 = -TMath::Log(tan_theta_absr2);
627
628 // for MC
629 if(fReadMCInfo) {
630 //mu1
631 rabsr1 = 0;
632 Double_t thetar1 = (180./TMath::Pi())*TMath::ATan2(TMath::Sqrt(pxr1*pxr1+pyr1*pyr1),pzr1);
633 theta_absr1 = thetar1;
634 etar1 = mu1->Eta();
635 //mu2
636 rabsr2 = 0;
637 Double_t thetar2 = (180./TMath::Pi())*TMath::ATan2(TMath::Sqrt(pxr2*pxr2+pyr2*pyr2),pzr2);
638 theta_absr2 = thetar2;
639 etar2 = mu2->Eta();
640 }
641
642 if(TMath::Abs(etar1) > 8 || TMath::Abs(etar2) > 8) continue;
643
644 // dimuon
645 Double_t ptrec = TMath::Sqrt((pxr1+pxr2)*(pxr1+pxr2)+(pyr1+pyr2)*(pyr1+pyr2));
646 Double_t raprec= Rap((er1+er2),(pzr1+pzr2));
647 Double_t imassrec = Imass(er1,pxr1,pyr1,pzr1,er2,pxr2,pyr2,pzr2);
648 if(!imassrec) continue;
649
650 containerInput[1] = raprec ;
651 containerInput[2] = ptrec ;
652 containerInput[3] = imassrec ;
653 containerInput[4] = trigCondition+0.05 ;
654 containerInput[5] = ptr1 ;
655 containerInput[5] = ptr2 ;
656 containerInput[6] = momr1;
657 containerInput[6] = momr2;
658 containerInput[7] = trigside;
659 containerInput[8] = rabsr1;
660 containerInput[8] = rabsr2;
661 containerInput[9] = zr1 + zr2;
662 containerInput[10] = etar1;
663 containerInput[10] = etar2;
664 containerInput[11] = theta_absr1;
665 containerInput[11] = theta_absr2;
666 containerInput[12] = vzr1;
667 containerInput[12] = vzr2;
668 containerInput[13] = dcar1;
669 containerInput[13] = dcar2;
670
671 if (trigfired==1 && trigside==1) fCFManager->GetParticleContainer()->Fill(containerInput,4);
672 if (trigfired==0 && trigside==1) fCFManager->GetParticleContainer()->Fill(containerInput,1);
673 if(fReadMCInfo) fCFManager->GetParticleContainer()->Fill(containerInput,1); // Rec container
674
675 } // second mu loop
676 } // first mu loop
677 } // end AOD-based ANALYSIS
678
679 // end User analysis loop
680}
681//________________________________________________________________________
682Double_t AliCFMuonResUpsilon::Imass(Float_t e1, Float_t px1, Float_t py1, Float_t pz1,
683 Float_t e2, Float_t px2, Float_t py2, Float_t pz2) const
684{
685// invariant mass calculation
686 Float_t imass = (e1+e2)*(e1+e2)-((px1+px2)*(px1+px2)+(py1+py2)*(py1+py2)+(pz1+pz2)*(pz1+pz2));
687 if(imass<0) return 0;
688 Float_t imassrec = TMath::Sqrt((e1+e2)*(e1+e2)-((px1+px2)*(px1+px2)+(py1+py2)*(py1+py2)+(pz1+pz2)*(pz1+pz2)));
689 return imassrec;
690}
691//________________________________________________________________________
692Double_t AliCFMuonResUpsilon::Rap(Float_t e, Float_t pz) const
693{
694// calculate rapidity
695 Float_t rap;
696 if(e!=pz){
697 rap = 0.5*TMath::Log((e+pz)/(e-pz));
698 return rap;
699 }
700 else{
701 rap = -200;
702 return rap;
703 }
704}
705//________________________________________________________________________
706Double_t AliCFMuonResUpsilon::CostCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
707Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2,
708Double_t Energy)
709{
710 TLorentzVector PMu1CM, PMu2CM, PProjCM, PTargCM, PDimuCM; // In the CM. frame
711 TLorentzVector PMu1Dimu, PMu2Dimu, PProjDimu, PTargDimu; // In the dimuon rest frame
712 TVector3 beta,zaxisCS;
713 Double_t mp=0.93827231;
714 //
715 // --- Fill the Lorentz vector for projectile and target in the CM frame
716 //
717 PProjCM.SetPxPyPzE(0.,0.,-Energy,TMath::Sqrt(Energy*Energy+mp*mp));
718 PTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(Energy*Energy+mp*mp));
719 //
720 // --- Get the muons parameters in the CM frame
721 //
722 PMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
723 PMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
724 //
725 // --- Obtain the dimuon parameters in the CM frame
726 //
727 PDimuCM=PMu1CM+PMu2CM;
728 //
729 // --- Translate the dimuon parameters in the dimuon rest frame
730 //
731 beta=(-1./PDimuCM.E())*PDimuCM.Vect();
732 PMu1Dimu=PMu1CM;
733 PMu2Dimu=PMu2CM;
734 PProjDimu=PProjCM;
735 PTargDimu=PTargCM;
736 PMu1Dimu.Boost(beta);
737 PMu2Dimu.Boost(beta);
738 PProjDimu.Boost(beta);
739 PTargDimu.Boost(beta);
740 //
741 // --- Determine the z axis for the CS angle
742 //
743 zaxisCS=(((PProjDimu.Vect()).Unit())-((PTargDimu.Vect()).Unit())).Unit();
744 //
745 // --- Determine the CS angle (angle between mu+ and the z axis defined above)
746 Double_t cost;
747
748 if(charge1>0) cost = zaxisCS.Dot((PMu1Dimu.Vect()).Unit());
749 else cost = zaxisCS.Dot((PMu2Dimu.Vect()).Unit());
750 return cost;
751}
752//________________________________________________________________________
753
754//________________________________________________________________________
755Double_t AliCFMuonResUpsilon::CostHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
756Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
757{
758 TLorentzVector PMu1CM, PMu2CM, PDimuCM; // In the CM frame
759 TLorentzVector PMu1Dimu, PMu2Dimu; // In the dimuon rest frame
760 TVector3 beta,zaxisCS;
761 //
762 // --- Get the muons parameters in the CM frame
763 //
764 PMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
765 PMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
766 //
767 // --- Obtain the dimuon parameters in the CM frame
768 //
769 PDimuCM=PMu1CM+PMu2CM;
770 //
771 // --- Translate the muon parameters in the dimuon rest frame
772 //
773 beta=(-1./PDimuCM.E())*PDimuCM.Vect();
774 PMu1Dimu=PMu1CM;
775 PMu2Dimu=PMu2CM;
776 PMu1Dimu.Boost(beta);
777 PMu2Dimu.Boost(beta);
778 //
779 // --- Determine the z axis for the calculation of the polarization angle (i.e. the direction of the dimuon in the CM system)
780 //
781 TVector3 zaxis;
782 zaxis=(PDimuCM.Vect()).Unit();
783
784 // --- Calculation of the polarization angle (angle between mu+ and the z axis defined above)
785 Double_t cost;
786 if(charge1>0) cost = zaxis.Dot((PMu1Dimu.Vect()).Unit());
787 else cost = zaxis.Dot((PMu2Dimu.Vect()).Unit());
788
789 return cost;
790}
791//________________________________________________________________________
792
793//________________________________________________________________________
794Double_t AliCFMuonResUpsilon::PhiCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
795Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2,
796Double_t Energy)
797{
798 TLorentzVector PMu1CM, PMu2CM, PProjCM, PTargCM, PDimuCM; // In the CM frame
799 TLorentzVector PMu1Dimu, PMu2Dimu, PProjDimu, PTargDimu; // In the dimuon rest frame
800 TVector3 beta,yaxisCS, xaxisCS, zaxisCS;
801 Double_t mp=0.93827231;
802 //
803 // --- Fill the Lorentz vector for projectile and target in the CM frame
804 //
805 PProjCM.SetPxPyPzE(0.,0.,-Energy,TMath::Sqrt(Energy*Energy+mp*mp));
806 PTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(Energy*Energy+mp*mp));
807 //
808 // --- Get the muons parameters in the CM frame
809 //
810 PMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
811 PMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
812 //
813 // --- Obtain the dimuon parameters in the CM frame
814 //
815 PDimuCM=PMu1CM+PMu2CM;
816 //
817 // --- Translate the dimuon parameters in the dimuon rest frame
818 //
819 beta=(-1./PDimuCM.E())*PDimuCM.Vect();
820 PMu1Dimu=PMu1CM;
821 PMu2Dimu=PMu2CM;
822 PProjDimu=PProjCM;
823 PTargDimu=PTargCM;
824 PMu1Dimu.Boost(beta);
825 PMu2Dimu.Boost(beta);
826 PProjDimu.Boost(beta);
827 PTargDimu.Boost(beta);
828 //
829 // --- Determine the z axis for the CS angle
830 //
831 zaxisCS=(((PProjDimu.Vect()).Unit())-((PTargDimu.Vect()).Unit())).Unit();
832 yaxisCS=(((PProjDimu.Vect()).Unit()).Cross((PTargDimu.Vect()).Unit())).Unit();
833 xaxisCS=(yaxisCS.Cross(zaxisCS)).Unit();
834
835 Double_t phi;
836 if(charge1>0) phi = TMath::ATan2((PMu1Dimu.Vect()).Dot(yaxisCS),((PMu1Dimu.Vect()).Dot(xaxisCS)));
837 else phi = TMath::ATan2((PMu2Dimu.Vect()).Dot(yaxisCS),((PMu2Dimu.Vect()).Dot(xaxisCS)));
838
839 return phi;
840}
841//________________________________________________________________________
842
843//________________________________________________________________________
844Double_t AliCFMuonResUpsilon::PhiHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
845Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2, Double_t Energy)
846{
847 TLorentzVector PMu1CM, PMu2CM, PProjCM, PTargCM, PDimuCM; // In the CM frame
848 TLorentzVector PMu1Dimu, PMu2Dimu, PProjDimu, PTargDimu; // In the dimuon rest frame
849 TVector3 beta,xaxis,yaxis,zaxis;
850 Double_t mp=0.93827231;
851
852 //
853 // --- Get the muons parameters in the CM frame
854 //
855 PMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
856 PMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
857 //
858 // --- Obtain the dimuon parameters in the CM frame
859 //
860 PDimuCM=PMu1CM+PMu2CM;
861 //
862 // --- Translate the muon parameters in the dimuon rest frame
863 //
864 zaxis=(PDimuCM.Vect()).Unit();
865
866 beta=(-1./PDimuCM.E())*PDimuCM.Vect();
867
868 PProjCM.SetPxPyPzE(0.,0.,-Energy,TMath::Sqrt(Energy*Energy+mp*mp));
869 PTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(Energy*Energy+mp*mp));
870
871 PProjDimu=PProjCM;
872 PTargDimu=PTargCM;
873
874 PProjDimu.Boost(beta);
875 PTargDimu.Boost(beta);
876
877 yaxis=((PProjDimu.Vect()).Cross(PTargDimu.Vect())).Unit();
878 xaxis=(yaxis.Cross(zaxis)).Unit();
879
880 PMu1Dimu=PMu1CM;
881 PMu2Dimu=PMu2CM;
882 PMu1Dimu.Boost(beta);
883 PMu2Dimu.Boost(beta);
884
885 Double_t phi;
886 if(charge1>0) phi = TMath::ATan2((PMu1Dimu.Vect()).Dot(yaxis),(PMu1Dimu.Vect()).Dot(xaxis));
887 else phi = TMath::ATan2((PMu2Dimu.Vect()).Dot(yaxis),(PMu2Dimu.Vect()).Dot(xaxis));
888
889 return phi;
890}
891
892//________________________________________________________________________
893void AliCFMuonResUpsilon::Terminate(Option_t *)
894{
895 // draw result of the Invariant mass MC and ESD
896
897 TH1D *h1 = dynamic_cast<TH1D*>(GetOutputData(1));
898 AliCFContainer *cont = dynamic_cast<AliCFContainer*> (GetOutputData(2));
899
900 TH1D *kmass = cont->ShowProjection(3,0);
901 TH1D *rmass = cont->ShowProjection(3,1);
902 TH1D *mmass = cont->ShowProjection(3,4);
903
904 TCanvas *c1 = new TCanvas("AliCFMuonResUpsilon","UPSILON Container",0,0,800,800);
905 c1->Divide(2,2);
906 c1->cd(1);
907 kmass->Draw("HIST");
908 c1->cd(2);
909 rmass->Draw("HIST");
910 c1->cd(3);
911 mmass->Draw("HIST");
912 c1->cd(4);
913 h1->Draw("HIST");
914}
915//________________________________________________________________________
916
917#endif