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