changes for Vertex and Tracks classes
[u/mrichter/AliRoot.git] / PWGLF / RESONANCES / extra / AliAnalysisTaskSigma1385.cxx
CommitLineData
35e49ca5 1/**************************************************************************
2 * Authors : Massimo Venaruzzo (massimo.venaruzzo@ts.infn.it) *
7356f978 3 * Enrico Fragiacomo (enrico.fragiacomo@ts.infn.it) *
35e49ca5 4 * Contributors are mentioned in the code where appropriate. *
5 * *
6 * Permission to use, copy, modify and distribute this software and its *
7 * documentation strictly for non-commercial purposes is hereby granted *
8 * without fee, provided that the above copyright notice appears in all *
9 * copies and that both the copyright notice and this permission notice *
10 * appear in the supporting documentation. The authors make no claims *
11 * about the suitability of this software for any purpose. It is *
12 * provided "as is" without express or implied warranty. *
13 **************************************************************************/
14
15//-----------------------------------------------------------------
16// AliAnalysisTaskSigma1385 class
17//-----------------------------------------------------------------
18
19class TTree;
20class TParticle;
21class TVector3;
22
23#include "AliAnalysisManager.h"
24#include <AliMCEventHandler.h>
25#include <AliMCEvent.h>
26#include <AliStack.h>
27
28class AliESDVertex;
29class AliESDv0;
30class AliAODv0;
31
32#include <iostream>
33
34#include "TList.h"
35#include "TH1.h"
36#include "TNtuple.h"
37#include "TGraph.h"
38#include "TCanvas.h"
39#include "TMath.h"
40#include "TChain.h"
41#include "AliLog.h"
42#include "AliESDEvent.h"
43#include "AliAODEvent.h"
44#include "AliCascadeVertexer.h"
45#include "AliESDcascade.h"
46#include "AliAODcascade.h"
47#include "AliAnalysisTaskSigma1385.h"
48#include "AliESDtrackCuts.h"
6e92ab0d 49#include "AliESDpid.h"
50#include "AliTOFT0maker.h"
51#include "AliTOFcalib.h"
52#include "AliCDBManager.h"
53#include "AliESDtrackCuts.h"
35e49ca5 54
55ClassImp(AliAnalysisTaskSigma1385)
7356f978 56
35e49ca5 57//________________________________________________________________________
7356f978 58AliAnalysisTaskSigma1385::AliAnalysisTaskSigma1385()
6e92ab0d 59 : AliAnalysisTaskSE(),
60
7356f978 61 //-------------------------------- For PID
6e92ab0d 62
63 fisMC(0),
64 fIsMC(fisMC),
7356f978 65 fCheckITS(kTRUE),
66 fCheckTPC(kTRUE),
67 fCheckTOF(kTRUE),
68 fUseGlobal(kTRUE),
69 fUseITSSA(kTRUE),
70 fMaxITSband(3.0),
71 fTPCpLimit(0.35),
a61c64bc 72 fTPCpar(),
7356f978 73 fMinTPCband(3.0),
74 fMaxTPCband(5.0),
75 fESDpid(0x0),
76 fTOFmaker(0x0),
77 fTOFcalib(0x0),
6e92ab0d 78 fTOFcalibrateESD(!fisMC),
7356f978 79 fTOFcorrectTExp(kTRUE),
80 fTOFuseT0(kTRUE),
6e92ab0d 81 fTOFtuneMC(fisMC),
7356f978 82 fTOFresolution(100.0),
83 fMinTOF(-2.5),
84 fMaxTOF(3.5),
6e92ab0d 85 fLastRun(-1),
a61c64bc 86 fOkTrack(),
6e92ab0d 87 fAnalysisType("ESD"), fCollidingSystems(0), fDataType("REAL"), fListHistCascade(0),
88 fHistEventMultiplicity(0), fHistEventMultiplicityRAVS(0),
89 fNtuple1(0), fNtuple2(0), fNtuple3(0), fNtuple4(0)
90
91
7356f978 92
93 //--------------------------------
94
35e49ca5 95{
7356f978 96 // Dummy Constructor
7d460484 97 Int_t i;
98 for (i = 0; i < 5; i++) fOkTrack[i] = kFALSE;
35e49ca5 99}
100
101//________________________________________________________________________
7356f978 102AliAnalysisTaskSigma1385::AliAnalysisTaskSigma1385(const char *name)
6e92ab0d 103 : AliAnalysisTaskSE(name),
104
7356f978 105
106 //-------------------------------- For PID
107
6e92ab0d 108 fisMC(0),
109 fIsMC(fisMC),
7356f978 110 fCheckITS(kTRUE),
111 fCheckTPC(kTRUE),
112 fCheckTOF(kTRUE),
113 fUseGlobal(kTRUE),
114 fUseITSSA(kTRUE),
115 fMaxITSband(3.0),
116 fTPCpLimit(0.35),
a61c64bc 117 fTPCpar(),
7356f978 118 fMinTPCband(3.0),
119 fMaxTPCband(5.0),
120 fESDpid(0x0),
121 fTOFmaker(0x0),
122 fTOFcalib(0x0),
6e92ab0d 123 fTOFcalibrateESD(!fisMC),
7356f978 124 fTOFcorrectTExp(kTRUE),
125 fTOFuseT0(kTRUE),
6e92ab0d 126 fTOFtuneMC(fisMC),
7356f978 127 fTOFresolution(100.0),
128 fMinTOF(-2.5),
129 fMaxTOF(3.5),
6e92ab0d 130 fLastRun(-1),
a61c64bc 131 fOkTrack(),
6e92ab0d 132 fAnalysisType("ESD"), fCollidingSystems(0), fDataType("REAL"), fListHistCascade(0),
133 fHistEventMultiplicity(0), fHistEventMultiplicityRAVS(0),
134 fNtuple1(0), fNtuple2(0), fNtuple3(0), fNtuple4(0)
7356f978 135
136 //--------------------------------
35e49ca5 137{
138
7356f978 139 // Output slot #0 writes into a TList container (Cascade)
140 DefineOutput(1, TList::Class());
7d460484 141
142 Int_t i;
143 for (i = 0; i < 5; i++) fOkTrack[i] = kFALSE;
35e49ca5 144}
145
146//________________________________________________________________________
147void AliAnalysisTaskSigma1385::UserCreateOutputObjects()
148{
7356f978 149 fListHistCascade = new TList();
150
151 if (! fHistEventMultiplicity) {
152 fHistEventMultiplicity = new TH1F("fHistEventMultiplicity" , "Nb of Events" , 4, -1.0, 3.0);
153 fListHistCascade->Add(fHistEventMultiplicity);
154 }
155
156 if (! fHistEventMultiplicityRAVS) {
157 fHistEventMultiplicityRAVS = new TH1F("fHistEventMultiplicityRAVS" , "Nb of Events Rejected After Vertex selection" , 4, -1.0, 3.0);
158 fListHistCascade->Add(fHistEventMultiplicityRAVS);
159 }
160
161 if (! fNtuple1) {
162 fNtuple1 = new TNtuple("fNtuple1", "Ntuple1", "TrkNmb");
163 fNtuple1->SetDirectory(0);
164 fListHistCascade->Add(fNtuple1);
165 }
166
167 if (! fNtuple2) {
168 fNtuple2 = new TNtuple("fNtuple2", "Ntuple2", "s:dcal:lCosPoinAn:lDaugDCA:lambdap:lambdapt:lambdamass");
169 fNtuple2->SetDirectory(0);
170 fListHistCascade->Add(fNtuple2);
171 }
172
173 if (! fNtuple3) {
174 fNtuple3 = new TNtuple("fNtuple3", "Ntuple3", "c:dcapi:ppi:ptpi:bachphi:bachtheta:okPiTPC:okPiTOF");
175 fNtuple3->SetDirectory(0);
176 fListHistCascade->Add(fNtuple3);
177 }
178
179 if (! fNtuple4) {
180 fNtuple4 = new TNtuple("fNtuple4", "Ntuple4", "dca:mc:phi:theta:eta:y:pt:p:opang:invmass");
181 fListHistCascade->Add(fNtuple4);
182 }
183
184
35e49ca5 185}// end UserCreateOutputObjects
186
187
188//________________________________________________________________________
7356f978 189void AliAnalysisTaskSigma1385::UserExec(Option_t *)
35e49ca5 190{
191
7356f978 192 // Main loop
193 // Called for each event
194
195
196 Info("AliAnalysisTaskSigma1385", "Starting UserExec");
197
198 AliMCEventHandler* eventHandler;
199 AliMCEvent* mcEvent = 0;
200
201 if (fDataType == "SIM") {
202
203 eventHandler = dynamic_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
204 if (!eventHandler) {
205 Printf("ERROR: Could not retrieve MC event handler");
206 return;
207 }
208
209 mcEvent = eventHandler->MCEvent();
210 if (!mcEvent) {
211 Printf("ERROR: Could not retrieve MC event");
212 return;
213 }
214
215 }
216
217 AliStack* stack = 0;
6e92ab0d 218 if (fDataType == "SIM") {stack = mcEvent->Stack(); fIsMC = 1; fisMC = 1;}
7356f978 219
220 AliESDEvent *lESDevent = 0x0;
221 AliAODEvent *lAODevent = 0x0;
222
223
224 // Connect to the InputEvent
225 Int_t ncascades = -1;
226 if (fAnalysisType == "ESD") {
227 lESDevent = dynamic_cast<AliESDEvent*>(InputEvent());
a50065d4 228 if (lESDevent)
229 ncascades = lESDevent->GetNumberOfCascades();
230 else {
231 Printf("ERROR: lESDevent not available \n");
232 return;
233 }
4d921955 234 } else if (fAnalysisType == "AOD") {
7356f978 235 lAODevent = dynamic_cast<AliAODEvent*>(InputEvent());
a50065d4 236 if (lAODevent)
237 ncascades = lAODevent->GetNumberOfCascades();
238 else {
239 Printf("ERROR: lAODevent not available \n");
240 return;
241 }
242 } else {
243 Printf("Input type undefined \n");
4d921955 244 return;
245 }
7356f978 246
247
248 //------------------------------for PID
249
250 SetCheckITS(kTRUE);
251 SetCheckTPC(kTRUE);
252 SetCheckTOF(kTRUE);
253
254// ----> set TPC range for PID and calibration
255 SetTPCrange(5.0, 3.0);
256 SetTPCpLimit(0.35);
257
258 // ----> set ITS range for PID
259 SetITSband(4.0);
260
261 // ----> set TPC calibration
262 if (fDataType == "SIM") SetTPCpar(2.15898 / 50.0, 1.75295E1, 3.40030E-9, 1.96178, 3.91720);
263 else SetTPCpar(1.41543 / 50.0, 2.63394E1, 5.0411E-11, 2.12543, 4.88663);
264
265 // ----> set the TOF calibration depending on type of input (sim/data)
266 SetTOFcorrectTExp(kTRUE);
267 SetTOFuseT0(kTRUE);
268 SetTOFresolution(100.0);
269 if (fDataType == "SIM") {
270 SetTOFcalibrateESD(kFALSE);
271 SetTOFtuneMC(kTRUE);
272 } else {
273
274 SetTOFcalibrateESD(kTRUE);
275 SetTOFtuneMC(kFALSE);
276 }
277
278
279 if (!fESDpid) {
280 fESDpid = new AliESDpid;
281 fESDpid->GetTPCResponse().SetBetheBlochParameters(fTPCpar[0], fTPCpar[1], fTPCpar[2], fTPCpar[3], fTPCpar[4]);
282 }
283
284 // initialize DB to current run
4f926d1d 285 if (!lESDevent) {
286 Printf("ERROR: lESDevent not available \n");
287 return;
288 }
7356f978 289 Int_t run = lESDevent->GetRunNumber();
290 if (run != fLastRun) {
291 //cout << "Run = " << run << " -- LAST = " << fLastRun << endl;
292 fLastRun = run;
293
294 // setup TOF maker & calibration
295 if (!fTOFcalib) fTOFcalib = new AliTOFcalib;
296 fTOFcalib->SetCorrectTExp(fTOFcorrectTExp);
297 if (!fTOFmaker) fTOFmaker = new AliTOFT0maker(fESDpid, fTOFcalib);
298 fTOFmaker->SetTimeResolution(fTOFresolution);
299
300 AliCDBManager *cdb = AliCDBManager::Instance();
301 cdb->ClearCache(); // suggestion by Annalisa
302 cdb->Clear(); // suggestion by Annalisa
303 cdb->SetDefaultStorage("raw://");
304 cdb->SetRun(run);
305 fTOFcalib->SetCorrectTExp(fTOFcorrectTExp);
306 fTOFcalib->Init();
307 }
308
309 // if required, calibrate the TOF t0 maker with current event
310 if (fTOFcalibrateESD) fTOFcalib->CalibrateESD(lESDevent);
311 if (fTOFtuneMC) fTOFmaker->TuneForMC(lESDevent);
312 if (fTOFuseT0) {
313 fTOFmaker->ComputeT0TOF(lESDevent);
314 fTOFmaker->ApplyT0TOF(lESDevent);
315 fESDpid->MakePID(lESDevent, kFALSE, 0.);
316 }
317
318
319 //--------------------------------------------------------
320
321
322 fHistEventMultiplicity->Fill(1);
323
324 //Some Quantities to characterize the event
35e49ca5 325
7356f978 326 Double_t b = lESDevent->GetMagneticField();
6e92ab0d 327 Int_t trackNumber = lESDevent->GetNumberOfTracks();
7356f978 328
329
330 //---------------------------
331 // new part from AliCascadeVertexer (use vertexer for reconstructing Sigma(1385)
332 //
333 //const AliESDVertex *vtxT3D=lESDevent->GetPrimaryVertex();
334 const AliESDVertex *vtxT3D = lESDevent->GetPrimaryVertexTracks();
335 if (vtxT3D->GetNContributors() < 1) {
336 // SPD vertex
337 vtxT3D = lESDevent->GetPrimaryVertexSPD();
338 if (vtxT3D->GetNContributors() < 1) {
339
340
341 fHistEventMultiplicityRAVS->Fill(1);
342 return;
35e49ca5 343 }
7356f978 344
345
346 }
347
e690d4d0 348 Double_t xPrimaryVertex = vtxT3D->GetX();
349 Double_t yPrimaryVertex = vtxT3D->GetY();
350 Double_t zPrimaryVertex = vtxT3D->GetZ();
7356f978 351
352 if (zPrimaryVertex > 10 || zPrimaryVertex < -10) return;
353
354 //-------------------------------------------------------
355 // New Part about tracks global selection criteria
356 //-------------------------------------------------------
357
358 AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts");
359
360 esdTrackCuts->SetAcceptKinkDaughters(0); // 0 = kFalse
361 esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
362 esdTrackCuts->SetMaxChi2PerClusterTPC(4);
363 esdTrackCuts->SetMinNClustersTPC(70);
364
365
366 //-------------------------------------------------------
367 // loops over V0s
368 //stores relevant V0s in an array
369 Int_t nV0 = (Int_t)lESDevent->GetNumberOfV0s();
370
371 TObjArray vtcs(nV0);
372 for (Int_t i = 0; i < nV0; i++) {
373 AliESDv0 *v = lESDevent->GetV0(i);
374 if (v->GetOnFlyStatus()) continue; // if kTRUE, then this V0 is recontructed
375
376 vtcs.AddLast(v);
377 }
378 nV0 = vtcs.GetEntriesFast();
379
380 //-------------------------------------------------------
381 // loops over bachelor tracks
382 // stores relevant tracks in another array
383
384 Int_t nentr = (Int_t)lESDevent->GetNumberOfTracks();
385 TArrayI trk(nentr); Int_t ntr = 0;
386
387 for (Int_t i = 0; i < nentr; i++) {
388 AliESDtrack *esdtr = lESDevent->GetTrack(i);
389
390 if (!(esdtr->GetStatus() & AliESDtrack::kITSrefit)) continue;
391 if (!(esdtr->GetStatus() & AliESDtrack::kTPCrefit)) continue;
392
393
394 if (!(esdTrackCuts->AcceptTrack(esdtr))) continue;
395
396 trk[ntr++] = i;
397 }
398
399
400
401 //-----------------------------------------------------
402 // nested loops over V0s and bachelors
403 Float_t massLambda = 1.11568;
404 Int_t ncasc = 0;
405 AliCascadeVertexer *cascvert = new AliCascadeVertexer();
406 for (Int_t i = 0; i < nV0; i++) { //loop on V0s
407
408 // Lambda
409 AliESDv0 *v = (AliESDv0*)vtcs.UncheckedAt(i);
410 Int_t strness = 0;
411 Double_t lambdaMass = 0;
412 Float_t lambdaP = 0;
413 Float_t lambdaPt = 0;
6e92ab0d 414 Float_t lambdaDCA = 0; // DCA between Lambda and Primary Vertex
7356f978 415 Double_t v0cospointangle = 0;
416 Float_t v0daughtersDCA = 0;
417
418
419 // Lambda quality cuts
420 UInt_t lIdxPosXi = (UInt_t) TMath::Abs(v->GetPindex());
421 UInt_t lIdxNegXi = (UInt_t) TMath::Abs(v->GetNindex());
422
423 AliESDtrack *pTrackXi = lESDevent->GetTrack(lIdxPosXi);
424 AliESDtrack *nTrackXi = lESDevent->GetTrack(lIdxNegXi);
425
426 // Filter like-sign V0
6e92ab0d 427 if ( (TMath::Abs(pTrackXi->GetSign()) - TMath::Abs(nTrackXi->GetSign()) ) < 0.1) continue;
7356f978 428
429 // WARNING: the following selections cannot be done for AOD yet...
430
431
432 v->ChangeMassHypothesis(kLambda0); // the v0 must be Lambda
433
434
6e92ab0d 435 if ((TMath::Abs(v->GetEffMass() - massLambda)) < 0.2) {
7356f978 436
437 if (!(pTrackXi->GetStatus() & AliESDtrack::kTPCrefit)) continue;
438 if (!(nTrackXi->GetStatus() & AliESDtrack::kTPCrefit)) continue;
439
440
441 if ((pTrackXi->GetTPCNcls()) < 70) continue;
442 if ((nTrackXi->GetTPCNcls()) < 70) continue;
443
6e92ab0d 444 strness = 1; lambdaMass = v->GetEffMass(); lambdaP = v->P(); lambdaPt = v->Pt(); lambdaDCA = v->GetD(xPrimaryVertex, yPrimaryVertex, zPrimaryVertex); v0cospointangle = v->GetV0CosineOfPointingAngle(); v0daughtersDCA = v->GetDcaV0Daughters();
7356f978 445
446
447
448 }
449
450
451 v->ChangeMassHypothesis(kLambda0Bar); // the v0 must be Anti Lambda
452
453
454
6e92ab0d 455 if ((TMath::Abs(v->GetEffMass() - massLambda)) < 0.2) {
7356f978 456
457 if (!(pTrackXi->GetStatus() & AliESDtrack::kTPCrefit)) continue;
458 if (!(nTrackXi->GetStatus() & AliESDtrack::kTPCrefit)) continue;
459
460
461 if ((pTrackXi->GetTPCNcls()) < 70) continue;
462 if ((nTrackXi->GetTPCNcls()) < 70) continue;
463
6e92ab0d 464 Int_t temp = strness + 1; strness = -1 * temp; lambdaMass = v->GetEffMass(); lambdaP = v->P(); lambdaPt = v->Pt(); lambdaDCA = v->GetD(xPrimaryVertex, yPrimaryVertex, zPrimaryVertex); v0cospointangle = v->GetV0CosineOfPointingAngle(); v0daughtersDCA = v->GetDcaV0Daughters();
7356f978 465
466
467 }
468
469 if (strness == 0) continue;
470
471
472 for (Int_t j = 0; j < ntr; j++) { //loop on tracks
473
474 // Pion bachelor
475 Int_t bidx = trk[j];
476 Int_t bachTPCcls = 0;
477
478
479 if ((bidx == v->GetIndex(0)) || (bidx == v->GetIndex(1))) continue; // bach and V0's daughter's must be different!
480
481 AliESDtrack *btrk = lESDevent->GetTrack(bidx);
482
483 if (!(btrk->GetStatus() & AliESDtrack::kTPCrefit)) continue;
484
485 if ((bachTPCcls = btrk->GetTPCNcls()) < 70) continue;
486
7d460484 487 //Bool_t *IsOkTrack = IsSelected(btrk); //for Alberto's PID
7356f978 488
7d460484 489 //Bool_t *okTrack = new Bool_t[5];
7356f978 490
7d460484 491 //for (Int_t k = 0; k < 5; k++) okTrack[k] = IsOkTrack[k];
492
493 IsSelected(btrk);
7356f978 494
495
6e92ab0d 496 Int_t bachCharge = btrk->Charge();
497 Float_t pionDCA = TMath::Abs(btrk->GetD(xPrimaryVertex, yPrimaryVertex, b)); // DCA between Bachelor and Primary Vertex
7356f978 498
499 Double_t bachphi = btrk->Phi();
500 Double_t bachtheta = btrk->Theta();
501 Double_t bachmass = 0.13957;
502 Double_t bachp = btrk->GetP();
503 Double_t bachpt = btrk->Pt();
504
505
506 // Distance between Lambda and Pion bachelor
507 AliESDv0 v0(*v), *pv0 = &v0;
508 AliExternalTrackParam bt(*btrk), *pbt = &bt;
509 Double_t dca = cascvert->PropagateToDCA(pv0, pbt, b); // distance bwn V0 and bach (in cm)
510 if (dca > 10.0) continue; // Note: was 0.1! Enlarged->further filter in second pass analysis
511
512
513 AliESDcascade cascade(*pv0, *pbt, bidx); //constucts a sigma1385 candidate
514 AliESDcascade *xi = &cascade;
515
516
517 UInt_t lBachIdx = (UInt_t) TMath::Abs(xi->GetBindex());
518 AliESDtrack *bachTrackXi = lESDevent->GetTrack(lBachIdx);
519
520
521
6e92ab0d 522 Short_t mcTrue=0;
7356f978 523 if (fDataType == "SIM") {
524
525 Int_t pLabel = TMath::Abs(pTrackXi->GetLabel());
526 Int_t nLabel = TMath::Abs(nTrackXi->GetLabel());
527 Int_t bachLabel = TMath::Abs(bachTrackXi->GetLabel());
528
529 Int_t motherpLabel = stack->Particle(pLabel)->GetFirstMother();
530 Int_t mothernLabel = stack->Particle(nLabel)->GetFirstMother();
531 Int_t motherbachLabel = stack->Particle(bachLabel)->GetFirstMother();
532
533
534 //cout<< "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
535 //cout<< "plabel " << pLabel << " nlabel " << nLabel << " mother p " << motherpLabel << " mother n " << mothernLabel << " bachlabel " << bachLabel << " mother bach " << motherbachLabel << endl;
536 //cout<< "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
537
538 if (motherpLabel > -1 && mothernLabel > -1) {
539 TParticle *plambda = stack->Particle(motherpLabel);
540 Int_t grandmother = plambda->GetFirstMother();
541
542 motherpLabel = TMath::Abs(motherpLabel);
543 mothernLabel = TMath::Abs(mothernLabel);
544 //motherbachLabel = TMath::Abs(motherbachLabel);
545
546 //cout<< "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
547 //cout<< "plabel " << pLabel << " nlabel " << nLabel << " mother p " << motherpLabel << " mother n " << mothernLabel << " mother lambda " << grandmother << " bachlabel " << bachLabel << " mother bach " << motherbachLabel << endl;
548 //cout<< "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
549
550 if (motherbachLabel > -1) {
551
552
553 //cout<< "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
554 //cout<< "mother lambda " << grandmother << " mother bach " << motherbachLabel << " PDG Code "<< stack->Particle(grandmother)->GetPdgCode() << endl;
555 //cout<< "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
556
557 if ((motherpLabel == mothernLabel) && (grandmother == motherbachLabel) &&
558 ((TMath::Abs(stack->Particle(grandmother)->GetPdgCode()) == 3114) ||
6e92ab0d 559 (TMath::Abs(stack->Particle(grandmother)->GetPdgCode()) == 3224))) mcTrue = 1;
560
561 if( (motherpLabel == mothernLabel) && (grandmother == motherbachLabel) &&
562 (TMath::Abs(stack->Particle(grandmother)->GetPdgCode())==3312) ) mcTrue = 2;
7356f978 563
564 }
565
566 }
567
568 }
569
570 Double_t lBachMomX = 0., lBachMomY = 0., lBachMomZ = 0.;
571 Double_t lPMom[3] = {0., 0., 0.,};
572 Double_t lNMom[3] = {0., 0., 0.,};
573 Float_t lLambdaMomX = 0., lLambdaMomY = 0., lLambdaMomZ = 0.;
574 xi->GetBPxPyPz(lBachMomX, lBachMomY, lBachMomZ);
575 xi->GetPPxPyPz(lPMom[0], lPMom[1], lPMom[2]);
576 xi->GetNPxPyPz(lNMom[0], lNMom[1], lNMom[2]);
577 lLambdaMomX = lPMom[0] + lNMom[0];
578 lLambdaMomY = lPMom[1] + lNMom[1];
579 lLambdaMomZ = lPMom[2] + lNMom[2];
580
581 Float_t lRapXi = xi->RapXi();
582 Float_t lEta = xi->Eta();
583 Float_t lTheta = xi->Theta();
584 Float_t lPhi = xi->Phi();
585 Float_t lPt = xi->Pt();
586 Float_t lP = xi->P();
587
588 // Support variables for invariant mass calculation
6e92ab0d 589 TLorentzVector lambda, pion, sigma;
590 Double_t eLambda = TMath::Sqrt(lLambdaMomX * lLambdaMomX + lLambdaMomY * lLambdaMomY + lLambdaMomZ * lLambdaMomZ + lambdaMass * lambdaMass);
591 Double_t ePion = TMath::Sqrt(lBachMomX * lBachMomX + lBachMomY * lBachMomY + lBachMomZ * lBachMomZ + bachmass * bachmass) ;
7356f978 592
6e92ab0d 593 lambda.SetPxPyPzE(lLambdaMomX, lLambdaMomY, lLambdaMomZ, eLambda);
594 pion.SetPxPyPzE(lBachMomX, lBachMomY, lBachMomZ, ePion);
7356f978 595
6e92ab0d 596 sigma = lambda + pion;
7356f978 597
598 Double_t openingangle, invmass = 0;
599
600 openingangle = TMath::ACos((lBachMomX * lLambdaMomX + lBachMomY * lLambdaMomY + lBachMomZ * lLambdaMomZ) / ((TMath::Sqrt(lBachMomX * lBachMomX + lBachMomY * lBachMomY + lBachMomZ * lBachMomZ)) * (TMath::Sqrt(lLambdaMomX * lLambdaMomX + lLambdaMomY * lLambdaMomY +
601 lLambdaMomZ * lLambdaMomZ))));
602
6e92ab0d 603 invmass = sigma.M();
7356f978 604
605
6e92ab0d 606 fNtuple1->Fill(trackNumber);
607 fNtuple2->Fill((1.*strness), lambdaDCA, v0cospointangle, v0daughtersDCA , lambdaP, lambdaPt, lambdaMass);
608 //fNtuple3->Fill((1.*bachCharge), pionDCA, bachp, bachpt, bachphi, bachtheta, okTrack[1], okTrack[2]);
609 fNtuple3->Fill((1.*bachCharge), pionDCA, bachp, bachpt, bachphi, bachtheta, fOkTrack[1], fOkTrack[2]);
610 fNtuple4->Fill(dca, (1.*mcTrue), lPhi, lTheta, lEta, lRapXi, lPt, lP, openingangle, invmass);
7356f978 611
612 ncasc++;
613
7d460484 614 //delete IsOkTrack;
615 //delete okTrack;
7356f978 616
617
618 } // end loop tracks
619 } // end loop V0s
620
621 Info("AliAnalysisTaskSigma1385", "Number of reconstructed Sigma(1385): %d", ncasc);
622
623
624 // Post output data.
625 PostData(1, fListHistCascade);
626
35e49ca5 627}
628
629//________________________________________________________________________
630
7d460484 631//Bool_t *AliAnalysisTaskSigma1385::IsSelected(AliESDtrack *track)
632void AliAnalysisTaskSigma1385::IsSelected(AliESDtrack *track)
35e49ca5 633{
634//
635//
7d460484 636 //Bool_t okTrack[5];
7356f978 637
7d460484 638 //for (Int_t i = 0; i < 5; i++) okTrack[i] = kFALSE;
639 for (Int_t i = 0; i < 5; i++) fOkTrack[i] = kFALSE;
7356f978 640
641
642 AliITSPIDResponse itsrsp(fIsMC);
643
644
645 ULong_t status;
646 Int_t k, nITS;
647 Double_t times[10], tpcNSigma, tpcMaxNSigma, itsSignal, itsNSigma, mom, tofTime, tofSigma, tofRef, tofRel;
648 Bool_t okTOF = kFALSE;
649 Bool_t okTPC = kFALSE;
650 Bool_t okITS = kFALSE;
651 Bool_t isTPC = kFALSE;
652 Bool_t isITSSA = kFALSE;
653 Bool_t isTOF = kFALSE;
654 UChar_t itsCluMap;
655
656 // get commonly used variables
657 status = (ULong_t)track->GetStatus();
658 mom = track->P();
659 isTPC = ((status & AliESDtrack::kTPCin) != 0);
a50065d4 660 isITSSA = ((!isTPC) && (track->IsOn(AliESDtrack::kITSrefit)) && (!track->IsOn(AliESDtrack::kITSpureSA))); //(!isTPC && (status & AliESDtrack::kITSrefit) != 0 && (status & AliESDtrack::kITSpureSA) == 0 && (status & AliESDtrack::kITSpid) != 0);
7356f978 661 isTOF = (((status & AliESDtrack::kTOFout) != 0) && ((status & AliESDtrack::kTIME) != 0) /* && mom > TMath::Max(b1, b2)*/);
662
663
664 // check if the track type matches what is required
665 if (!isTPC && !isITSSA) {
666 AliDebug(AliLog::kDebug + 2, "Track is not either a TPC track or a ITS standalone. Rejected");
7d460484 667 //return okTrack;
668 return;
7356f978 669
670 } else if (isTPC && !fUseGlobal) {
671 AliDebug(AliLog::kDebug + 2, "Global tracks not used. Rejected");
7d460484 672 //return okTrack;
673 return;
7356f978 674
675 } else if (isITSSA && !fUseITSSA) {
676 AliDebug(AliLog::kDebug + 2, "ITS standalone not used. Rejected");
7d460484 677 //return okTrack;
678 return;
7356f978 679 }
680
681 // does a preliminary check on TOF values, if necessary
682 // then, if the reference time or TOF signal are meaningless
683 // even if the 'isTOF' flag is true, switch it to false
684 if (isTOF) {
685 track->GetIntegratedTimes(times);
686 tofTime = (Double_t)track->GetTOFsignal();
687 tofSigma = fTOFmaker->GetExpectedSigma(mom, times[AliPID::kPion], AliPID::ParticleMass(AliPID::kPion));
688 tofRef = times[AliPID::kPion];
689 if (tofRef <= 0.0 && tofSigma <= 0.0) isTOF = kFALSE;
690 }
691
692
693
694 // check TPC dE/dx
695 if (isTPC) { // this branch is entered by all global tracks
696 // check TPC dE/dx:
697 if (fCheckTPC) {
698 tpcNSigma = TMath::Abs(fESDpid->NumberOfSigmasTPC(track, AliPID::kPion));
699 if (track->GetInnerParam()->P() > fTPCpLimit) tpcMaxNSigma = fMinTPCband; else tpcMaxNSigma = fMaxTPCband;
700 okTPC = (tpcNSigma <= tpcMaxNSigma);
701 AliDebug(AliLog::kDebug + 2, Form("TPC nsigma = %f -- max = %f -- cut %s", tpcNSigma, tpcMaxNSigma, (okTPC ? "passed" : "failed")));
702 } else {
703 // if TPC is not checked, it is as if all tracks do pass the cut
704 okTPC = kTRUE;
705 AliDebug(AliLog::kDebug + 2, "TPC not checked, track accepted");
706 }
707
708 // check TOF (only if flags are OK)
709 if (fCheckTOF) {
710 if (isTOF) {
711 // TOF can be checked only when track is matched there
712 track->GetIntegratedTimes(times);
713 tofTime = (Double_t)track->GetTOFsignal();
714 tofSigma = fTOFmaker->GetExpectedSigma(mom, times[AliPID::kPion], AliPID::ParticleMass(AliPID::kPion));
715 tofRef = times[AliPID::kPion];
716
717 tofRel = (tofTime - tofRef) / tofSigma;
718 okTOF = ((tofRel >= fMinTOF) && (tofRel <= fMaxTOF));
719 AliDebug(AliLog::kDebug + 2, Form("TOF nsigma = %f -- range = %f %f -- cut %s", tofRel, fMinTOF, fMaxTOF, (okTOF ? "passed" : "failed")));
720 } else {
721 // if TOF is not matched, the answer depends on TPC:
722 // - if TPC is required, track is checked only there and TOF simply ignored
723 // - if TPC is not required, track is rejected when TOF does not match it, if TOF check is required
724 if (fCheckTPC) okTOF = kTRUE; else okTOF = kFALSE;
725 }
726 } else {
727 okTOF = kTRUE;
35e49ca5 728 }
7356f978 729
730 } else if (isITSSA) { // this branch is entered by all ITS standalone tracks
731 // check dE/dx only if this is required, otherwise ITS standalone are just added but not checked for PID
732 if (fCheckITS) {
733 itsSignal = track->GetITSsignal();
734 itsCluMap = track->GetITSClusterMap();
735 nITS = 0;
736 for (k = 2; k < 6; k++) if (itsCluMap & (1 << k)) ++nITS;
7d460484 737 //if (nITS < 3) return kFALSE;
738 //itsNSigma = itsrsp.GetNumberOfSigmas(mom, itsSignal, AliPID::kPion, nITS, kTRUE);
739 //okITS = ((TMath::Abs(itsNSigma)) <= fMaxITSband);
740 if (nITS < 3)
741 okITS = kFALSE;
742 else {
743 itsNSigma = itsrsp.GetNumberOfSigmas(mom, itsSignal, AliPID::kPion, nITS, kTRUE);
744 okITS = ((TMath::Abs(itsNSigma)) <= fMaxITSband);
745 AliDebug(AliLog::kDebug + 2, Form("ITS nsigma = %f -- max = %f -- cut %s", itsNSigma, fMaxITSband, (okITS ? "passed" : "failed")));
746 }
7356f978 747 } else {
748 okITS = kTRUE;
35e49ca5 749 }
7356f978 750 } else {
751 // if we are here, the track is surely bad
752 okITS = kFALSE;
753 okTPC = kFALSE;
754 okTOF = kFALSE;
35e49ca5 755 }
7356f978 756
757
7d460484 758 //okTrack[0] = okITS;
759 //okTrack[1] = okTPC;
760 //okTrack[2] = okTOF;
761 //okTrack[3] = isTPC;
762 //okTrack[4] = isITSSA;
763
764 fOkTrack[0] = okITS;
765 fOkTrack[1] = okTPC;
766 fOkTrack[2] = okTOF;
767 fOkTrack[3] = isTPC;
768 fOkTrack[4] = isITSSA;
7356f978 769
770 //cout<<"##########################################"<<endl;
771 //cout<<"isITSSA "<<isITSSA<< " isTPC "<<isTPC<<endl;
772 //cout<<"##########################################"<<endl;
773
774
7d460484 775 //return okTrack;
35e49ca5 776}
777
778
779
780
781
782
783
784
785
786
787
788//________________________________________________________________________
7356f978 789void AliAnalysisTaskSigma1385::Terminate(Option_t *)
35e49ca5 790{
7356f978 791 // Draw result to the screen
792 // Called once at the end of the query
35e49ca5 793}
794
795
796