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