1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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
22 //-----------------------------------------------------------------------
23 // Author : R. Vernet, Consorzio Cometa - Catania (it)
24 //-----------------------------------------------------------------------
25 // Modification done by X. Lopez - LPC Clermont (fr)
26 //-----------------------------------------------------------------------
29 #ifndef ALICFMUONRESTASK1_CXX
30 #define ALICFMUONRESTASK1_CXX
32 #include "AliCFMuonResTask1.h"
33 #include "AliHeader.h"
34 #include "AliESDHeader.h"
36 #include "TParticle.h"
37 #include "TLorentzVector.h"
39 #include "AliMCEvent.h"
40 #include "AliAnalysisManager.h"
41 #include "AliESDEvent.h"
42 #include "AliAODEvent.h"
43 #include "AliCFManager.h"
44 #include "AliCFCutBase.h"
45 #include "AliCFContainer.h"
47 #include "AliESDtrack.h"
49 #include "AliESDMuonTrack.h"
50 #include "AliESDtrack.h"
51 #include "AliESDInputHandler.h"
54 ClassImp(AliCFMuonResTask1)
56 //__________________________________________________________________________
57 AliCFMuonResTask1::AliCFMuonResTask1() :
61 fHistEventsProcessed(0x0),
68 //___________________________________________________________________________
69 AliCFMuonResTask1::AliCFMuonResTask1(const Char_t* name) :
70 AliAnalysisTaskSE(name),
74 fHistEventsProcessed(0x0),
78 // Constructor. Initialization of Inputs and Outputs
80 Info("AliCFMuonResTask1","Calling Constructor");
82 fHistEventsProcessed = new TH1I("fHistEventsProcessed","",1,0,1) ;
84 DefineOutput(1,TH1I::Class());
85 DefineOutput(2,AliCFContainer::Class());
89 //___________________________________________________________________________
90 AliCFMuonResTask1& AliCFMuonResTask1::operator=(const AliCFMuonResTask1& c)
93 // Assignment operator
96 AliAnalysisTaskSE::operator=(c) ;
97 fReadAODData = c.fReadAODData ;
98 fCFManager = c.fCFManager;
99 fQAHistList = c.fQAHistList ;
100 fHistEventsProcessed = c.fHistEventsProcessed;
106 //___________________________________________________________________________
107 AliCFMuonResTask1::AliCFMuonResTask1(const AliCFMuonResTask1& c) :
108 AliAnalysisTaskSE(c),
109 fReadAODData(c.fReadAODData),
110 fCFManager(c.fCFManager),
111 fQAHistList(c.fQAHistList),
112 fHistEventsProcessed(c.fHistEventsProcessed),
120 //___________________________________________________________________________
121 AliCFMuonResTask1::~AliCFMuonResTask1() {
125 Info("~AliCFMuonResTask1","Calling Destructor");
126 if (fCFManager) delete fCFManager ;
127 if (fHistEventsProcessed) delete fHistEventsProcessed ;
128 if (fQAHistList) {fQAHistList->Clear(); delete fQAHistList;}
131 //_________________________________________________
132 void AliCFMuonResTask1::UserExec(Option_t *)
135 // Main loop function
138 Info("UserExec","") ;
140 Error("UserExec","NO MC EVENT FOUND!");
145 Double_t containerInput[13] ;
146 Double_t beamEnergy=7000;
152 fCFManager->SetMCEventInfo(fMCEvent);
153 AliStack* stack = fMCEvent->Stack();
155 // loop on the MC event
156 for (Int_t ipart=0; ipart<fMCEvent->GetNumberOfTracks(); ipart++) {
157 AliMCParticle *mcPart = (AliMCParticle*) fMCEvent->GetTrack(ipart);
159 TParticle *part = mcPart->Particle();
160 TParticle *part0 = mcPart->Particle();
161 TParticle *part1 = mcPart->Particle();
163 // Selection of the resonance
164 if (!fCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,mcPart)) continue;
167 Float_t e = part->Energy();
168 Float_t pz = part->Pz();
169 //Float_t py = part->Py();
170 //Float_t px = part->Px();
171 Float_t rapmc = Rap(e,pz);
172 //Float_t mass = part->GetCalcMass();
176 Int_t p0 = part->GetDaughter(0);
177 part0 = stack->Particle(p0);
178 // selection of the rapidity for first muon
179 AliMCParticle *mcpart0 = new AliMCParticle(part0);
180 if (!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mcpart0)) continue;
181 Int_t pdg0 = part0->GetPdgCode();
183 Int_t p1 = part->GetDaughter(1);
184 part1 = stack->Particle(p1);
185 // selection of the rapidity for second muon
186 AliMCParticle *mcpart1 = new AliMCParticle(part1);
187 if (!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mcpart1)) continue;
188 Int_t pdg1 = part1->GetPdgCode();
191 Float_t e0=-999, pz0=-999, py0=-999, px0=-999, phi0=-999, rapmc0=-999;
192 Float_t e1=-999, pz1=-999, py1=-999, px1=-999, phi1=-999, rapmc1=-999;
193 Double_t charge0=-999, charge1=-999;
195 // ordering sign: first = mu-
197 e0 = part0->Energy();
204 charge0 = part0->GetPDG()->Charge()/3;
206 e1 = part1->Energy();
213 charge1 = part1->GetPDG()->Charge()/3;
217 e1 = part0->Energy();
224 charge1 = part0->GetPDG()->Charge()/3;
226 e0 = part1->Energy();
233 charge0 = part1->GetPDG()->Charge()/3;
237 if(pdg0==13 || pdg1==13) {
239 Float_t pmc = TMath::Sqrt((px0+px1)*(px0+px1)+(py0+py1)*(py0+py1)+
240 (pz0+pz1)*(pz0+pz1));
241 Float_t ptmc = TMath::Sqrt((px0+px1)*(px0+px1)+(py0+py1)*(py0+py1));
242 Float_t imassmc = Imass(e0,px0,py0,pz0,e1,px1,py1,pz1);
243 //Float_t rapmc_check=Rap((e0+e1),(pz0+pz1));
245 Double_t costCS = CostCS((Double_t) px0,(Double_t) py0,(Double_t) pz0,(Double_t) e0, (Double_t)charge0,(Double_t) px1,(Double_t) py1,(Double_t) pz1,(Double_t) e1, beamEnergy);
246 Double_t costHE = CostHE((Double_t) px0,(Double_t) py0,(Double_t) pz0,(Double_t) e0, (Double_t)charge0,(Double_t) px1,(Double_t) py1,(Double_t) pz1,(Double_t) e1);
247 Double_t phiCS = PhiCS((Double_t) px0,(Double_t) py0,(Double_t) pz0,(Double_t) e0, (Double_t)charge0,(Double_t) px1,(Double_t) py1,(Double_t) pz1,(Double_t) e1, beamEnergy);
248 Double_t phiHE = PhiHE((Double_t) px0,(Double_t) py0,(Double_t) pz0,(Double_t) e0, (Double_t)charge0,(Double_t) px1,(Double_t) py1,(Double_t) pz1,(Double_t) e1, beamEnergy);
250 containerInput[0] = fNevt+0.5 ;
251 containerInput[1] = rapmc0 ;
252 containerInput[2] = phi0 ;
253 containerInput[3] = rapmc1 ;
254 containerInput[4] = phi1 ;
255 containerInput[5] = imassmc ;
256 containerInput[6] = rapmc ;
257 containerInput[7] = ptmc;
258 containerInput[8] = pmc;
259 containerInput[9] = costCS;
260 containerInput[10] = phiCS;
261 containerInput[11] = costHE;
262 containerInput[12] = phiHE;
264 // fill the container at the first step
265 fCFManager->GetParticleContainer()->Fill(containerInput,0);
274 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>
275 (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
276 fESD = esdH->GetEvent();
277 Int_t mult1 = fESD->GetNumberOfMuonTracks() ;
280 for (Int_t j = 0; j<mult1; j++) {
281 AliESDMuonTrack* mu1 = new AliESDMuonTrack(*(fESD->GetMuonTrack(j)));
282 Float_t zr1 = mu1->Charge();
285 Float_t pxr1 = mu1->Px();
286 Float_t pyr1 = mu1->Py();
287 Float_t pzr1 = mu1->Pz();
288 Float_t phir1 = mu1->Phi();
290 Float_t er1 = mu1->E();
291 Float_t rap1=Rap(er1,pzr1);
292 // selection of the rapidity for first muon
293 if (!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mu1)) continue;
295 for (Int_t jj = 0; jj<mult1; jj++) {
296 AliESDMuonTrack* mu2 = new AliESDMuonTrack(*(fESD->GetMuonTrack(jj)));
297 Float_t zr2 = mu2->Charge();
300 Float_t pxr2 = mu2->Px();
301 Float_t pyr2 = mu2->Py();
302 Float_t pzr2 = mu2->Pz();
303 Float_t phir2 = mu2->Phi();
304 phir2 = Phideg(phir2);
305 Float_t er2 = mu2->E();
306 Float_t rap2=Rap(er2,pzr2);
307 // selection of the rapidity for second muon
308 if (!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mu2)) continue;
310 Float_t prec = TMath::Sqrt((pxr1+pxr2)*(pxr1+pxr2)+(pyr1+pyr2)*(pyr1+pyr2)+
311 (pzr1+pzr2)*(pzr1+pzr2));
312 Float_t ptrec = TMath::Sqrt((pxr1+pxr2)*(pxr1+pxr2)+(pyr1+pyr2)*(pyr1+pyr2));
313 Float_t raprec=Rap((er1+er2),(pzr1+pzr2));
314 Float_t imassrec = Imass(er1,pxr1,pyr1,pzr1,er2,pxr2,pyr2,pzr2);
316 Double_t costCSrec = CostCS((Double_t) pxr1,(Double_t) pyr1,(Double_t)pzr1,(Double_t) er1, (Double_t)zr1,(Double_t) pxr2,(Double_t) pyr2,(Double_t)pzr2,(Double_t) er2, beamEnergy);
317 Double_t costHErec = CostHE((Double_t) pxr1,(Double_t) pyr1,(Double_t)pzr1,(Double_t) er1, (Double_t)zr1,(Double_t) pxr2,(Double_t) pyr2,(Double_t)pzr2,(Double_t) er2);
318 Double_t phiCSrec = PhiCS((Double_t) pxr1,(Double_t) pyr1,(Double_t)pzr1,(Double_t) er1, (Double_t)zr1,(Double_t) pxr2,(Double_t) pyr2,(Double_t)pzr2,(Double_t) er2, beamEnergy);
319 Double_t phiHErec = PhiHE((Double_t) pxr1,(Double_t) pyr1,(Double_t)pzr1,(Double_t) er1, (Double_t)zr1,(Double_t) pxr2,(Double_t) pyr2,(Double_t)pzr2,(Double_t) er2, beamEnergy);
321 containerInput[0] = fNevt+0.5 ;
322 containerInput[1] = rap1 ;
323 containerInput[2] = phir1 ;
324 containerInput[3] = rap2 ;
325 containerInput[4] = phir2 ;
326 containerInput[5] = imassrec ;
327 containerInput[6] = raprec ;
328 containerInput[7] = ptrec;
329 containerInput[8] = prec;
330 containerInput[9] = costCSrec;
331 containerInput[10] = phiCSrec;
332 containerInput[11] = costHErec;
333 containerInput[12] = phiHErec;
335 // fill the container at the second step
336 fCFManager->GetParticleContainer()->Fill(containerInput,1);
344 fHistEventsProcessed->Fill(0);
345 PostData(1,fHistEventsProcessed) ;
346 PostData(2,fCFManager->GetParticleContainer()) ;
348 //________________________________________________________________________
349 Float_t AliCFMuonResTask1::Imass(Float_t e1, Float_t px1, Float_t py1, Float_t pz1,
350 Float_t e2, Float_t px2, Float_t py2, Float_t pz2) const
352 // invariant mass calculation
353 Float_t imassrec = TMath::Sqrt((e1+e2)*(e1+e2)-((px1+px2)*(px1+px2)+
354 (py1+py2)*(py1+py2)+(pz1+pz2)*(pz1+pz2)));
357 //________________________________________________________________________
358 Float_t AliCFMuonResTask1::Rap(Float_t e, Float_t pz) const
360 // calculate rapidity
362 if(TMath::Abs(e-pz)>1e-7){
363 rap = 0.5*TMath::Log((e+pz)/(e-pz));
371 //________________________________________________________________________
372 Float_t AliCFMuonResTask1::Phideg(Float_t phi) const
374 // calculate Phi in range [-180,180]
377 phideg = phi-TMath::Pi();
378 phideg = phideg*57.296;
381 //________________________________________________________________________
382 Double_t AliCFMuonResTask1::CostCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
383 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2,
387 TLorentzVector pMu1CM, pMu2CM, pProjCM, pTargCM, pDimuCM; // In the CM. frame
388 TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame
389 TVector3 beta,zaxisCS;
390 Double_t mp=0.93827231;
392 // --- Fill the Lorentz vector for projectile and target in the CM frame
394 pProjCM.SetPxPyPzE(0.,0.,-Energy,TMath::Sqrt(Energy*Energy+mp*mp));
395 pTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(Energy*Energy+mp*mp));
397 // --- Get the muons parameters in the CM frame
399 pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
400 pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
402 // --- Obtain the dimuon parameters in the CM frame
404 pDimuCM=pMu1CM+pMu2CM;
406 // --- Translate the dimuon parameters in the dimuon rest frame
408 beta=(-1./pDimuCM.E())*pDimuCM.Vect();
413 pMu1Dimu.Boost(beta);
414 pMu2Dimu.Boost(beta);
415 pProjDimu.Boost(beta);
416 pTargDimu.Boost(beta);
418 // --- Determine the z axis for the CS angle
420 zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit();
422 // --- Determine the CS angle (angle between mu+ and the z axis defined above)
425 if(charge1>0) cost = zaxisCS.Dot((pMu1Dimu.Vect()).Unit());
426 else cost = zaxisCS.Dot((pMu2Dimu.Vect()).Unit());
429 //________________________________________________________________________
431 //________________________________________________________________________
432 Double_t AliCFMuonResTask1::CostHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
433 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
436 TLorentzVector pMu1CM, pMu2CM, pDimuCM; // In the CM frame
437 TLorentzVector pMu1Dimu, pMu2Dimu; // In the dimuon rest frame
438 TVector3 beta,zaxisCS;
440 // --- Get the muons parameters in the CM frame
442 pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
443 pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
445 // --- Obtain the dimuon parameters in the CM frame
447 pDimuCM=pMu1CM+pMu2CM;
449 // --- Translate the muon parameters in the dimuon rest frame
451 beta=(-1./pDimuCM.E())*pDimuCM.Vect();
454 pMu1Dimu.Boost(beta);
455 pMu2Dimu.Boost(beta);
457 // --- Determine the z axis for the calculation of the polarization angle (i.e. the direction of the dimuon in the CM system)
460 zaxis=(pDimuCM.Vect()).Unit();
462 // --- Calculation of the polarization angle (angle between mu+ and the z axis defined above)
464 if(charge1>0) cost = zaxis.Dot((pMu1Dimu.Vect()).Unit());
465 else cost = zaxis.Dot((pMu2Dimu.Vect()).Unit());
469 //________________________________________________________________________
471 //________________________________________________________________________
472 Double_t AliCFMuonResTask1::PhiCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
473 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2,
477 TLorentzVector pMu1CM, pMu2CM, pProjCM, pTargCM, pDimuCM; // In the CM frame
478 TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame
479 TVector3 beta,yaxisCS, xaxisCS, zaxisCS;
480 Double_t mp=0.93827231;
482 // --- Fill the Lorentz vector for projectile and target in the CM frame
484 pProjCM.SetPxPyPzE(0.,0.,-Energy,TMath::Sqrt(Energy*Energy+mp*mp));
485 pTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(Energy*Energy+mp*mp));
487 // --- Get the muons parameters in the CM frame
489 pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
490 pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
492 // --- Obtain the dimuon parameters in the CM frame
494 pDimuCM=pMu1CM+pMu2CM;
496 // --- Translate the dimuon parameters in the dimuon rest frame
498 beta=(-1./pDimuCM.E())*pDimuCM.Vect();
503 pMu1Dimu.Boost(beta);
504 pMu2Dimu.Boost(beta);
505 pProjDimu.Boost(beta);
506 pTargDimu.Boost(beta);
508 // --- Determine the z axis for the CS angle
510 zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit();
511 yaxisCS=(((pProjDimu.Vect()).Unit()).Cross((pTargDimu.Vect()).Unit())).Unit();
512 xaxisCS=(yaxisCS.Cross(zaxisCS)).Unit();
515 if(charge1>0) phi = TMath::ATan2((pMu1Dimu.Vect()).Dot(yaxisCS),((pMu1Dimu.Vect()).Dot(xaxisCS)));
516 else phi = TMath::ATan2((pMu2Dimu.Vect()).Dot(yaxisCS),((pMu2Dimu.Vect()).Dot(xaxisCS)));
520 //________________________________________________________________________
522 //________________________________________________________________________
523 Double_t AliCFMuonResTask1::PhiHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
524 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2, Double_t Energy)
527 TLorentzVector pMu1CM, pMu2CM, pProjCM, pTargCM, pDimuCM; // In the CM frame
528 TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame
529 TVector3 beta,xaxis,yaxis,zaxis;
530 Double_t mp=0.93827231;
533 // --- Get the muons parameters in the CM frame
535 pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
536 pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
538 // --- Obtain the dimuon parameters in the CM frame
540 pDimuCM=pMu1CM+pMu2CM;
542 // --- Translate the muon parameters in the dimuon rest frame
544 zaxis=(pDimuCM.Vect()).Unit();
546 beta=(-1./pDimuCM.E())*pDimuCM.Vect();
548 pProjCM.SetPxPyPzE(0.,0.,-Energy,TMath::Sqrt(Energy*Energy+mp*mp));
549 pTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(Energy*Energy+mp*mp));
554 pProjDimu.Boost(beta);
555 pTargDimu.Boost(beta);
557 yaxis=((pProjDimu.Vect()).Cross(pTargDimu.Vect())).Unit();
558 xaxis=(yaxis.Cross(zaxis)).Unit();
562 pMu1Dimu.Boost(beta);
563 pMu2Dimu.Boost(beta);
566 if(charge1>0) phi = TMath::ATan2((pMu1Dimu.Vect()).Dot(yaxis),(pMu1Dimu.Vect()).Dot(xaxis));
567 else phi = TMath::ATan2((pMu2Dimu.Vect()).Dot(yaxis),(pMu2Dimu.Vect()).Dot(xaxis));
572 //________________________________________________________________________
573 void AliCFMuonResTask1::Terminate(Option_t *)
575 // draw result of the Invariant mass MC and ESD
577 AliCFContainer *cont = dynamic_cast<AliCFContainer*> (GetOutputData(2));
579 TH1D *kmass = cont->ShowProjection(5,0);
580 TH1D *rmass = cont->ShowProjection(5,1);
581 TH1D *kcostCS = cont->ShowProjection(9,0);
582 TH1D *rcostCS = cont->ShowProjection(9,1);
584 TCanvas *c1 = new TCanvas("AliCFMuonResTask1","JPSI MC & ESD",10,10,510,510);
591 kcostCS->Draw("HIST");
593 rcostCS->Draw("HIST");
595 //________________________________________________________________________