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