changes for Vertex and Tracks classes
[u/mrichter/AliRoot.git] / PWGLF / 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 #include "AliESDpid.h"
50 #include "AliTOFT0maker.h"
51 #include "AliTOFcalib.h"
52 #include "AliCDBManager.h"
53 #include "AliESDtrackCuts.h"
54
55 ClassImp(AliAnalysisTaskSigma1385)
56
57 //________________________________________________________________________
58 AliAnalysisTaskSigma1385::AliAnalysisTaskSigma1385()
59    : AliAnalysisTaskSE(),
60      
61      //-------------------------------- For PID
62      
63      fisMC(0),
64      fIsMC(fisMC),
65      fCheckITS(kTRUE),
66      fCheckTPC(kTRUE),
67      fCheckTOF(kTRUE),
68      fUseGlobal(kTRUE),
69      fUseITSSA(kTRUE),
70      fMaxITSband(3.0),
71      fTPCpLimit(0.35),
72      fTPCpar(),
73      fMinTPCband(3.0),
74      fMaxTPCband(5.0),
75      fESDpid(0x0),
76      fTOFmaker(0x0),
77      fTOFcalib(0x0),
78      fTOFcalibrateESD(!fisMC),
79      fTOFcorrectTExp(kTRUE),
80      fTOFuseT0(kTRUE),
81      fTOFtuneMC(fisMC),
82      fTOFresolution(100.0),
83      fMinTOF(-2.5),
84      fMaxTOF(3.5),
85      fLastRun(-1),
86      fOkTrack(),
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
92
93      //--------------------------------
94
95 {
96    // Dummy Constructor
97    Int_t i;
98    for (i = 0; i < 5; i++) fOkTrack[i] = kFALSE;
99 }
100
101 //________________________________________________________________________
102 AliAnalysisTaskSigma1385::AliAnalysisTaskSigma1385(const char *name)
103    : AliAnalysisTaskSE(name),
104        
105
106      //-------------------------------- For PID
107
108      fisMC(0),
109      fIsMC(fisMC),
110      fCheckITS(kTRUE),
111      fCheckTPC(kTRUE),
112      fCheckTOF(kTRUE),
113      fUseGlobal(kTRUE),
114      fUseITSSA(kTRUE),
115      fMaxITSband(3.0),
116      fTPCpLimit(0.35),
117      fTPCpar(),
118      fMinTPCband(3.0),
119      fMaxTPCband(5.0),
120      fESDpid(0x0),
121      fTOFmaker(0x0),
122      fTOFcalib(0x0),
123      fTOFcalibrateESD(!fisMC),
124      fTOFcorrectTExp(kTRUE),
125      fTOFuseT0(kTRUE),
126      fTOFtuneMC(fisMC),
127      fTOFresolution(100.0),
128      fMinTOF(-2.5),
129      fMaxTOF(3.5),
130      fLastRun(-1),
131      fOkTrack(),
132      fAnalysisType("ESD"), fCollidingSystems(0), fDataType("REAL"), fListHistCascade(0), 
133      fHistEventMultiplicity(0), fHistEventMultiplicityRAVS(0), 
134      fNtuple1(0), fNtuple2(0), fNtuple3(0), fNtuple4(0)
135
136      //--------------------------------
137 {
138
139    // Output slot #0 writes into a TList container (Cascade)
140    DefineOutput(1, TList::Class());
141    
142    Int_t i;
143    for (i = 0; i < 5; i++) fOkTrack[i] = kFALSE;
144 }
145
146 //________________________________________________________________________
147 void AliAnalysisTaskSigma1385::UserCreateOutputObjects()
148 {
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
185 }// end UserCreateOutputObjects
186
187
188 //________________________________________________________________________
189 void AliAnalysisTaskSigma1385::UserExec(Option_t *)
190 {
191
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;
218    if (fDataType == "SIM") {stack = mcEvent->Stack(); fIsMC = 1; fisMC = 1;}
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());
228       if (lESDevent) 
229          ncascades = lESDevent->GetNumberOfCascades();
230       else {
231          Printf("ERROR: lESDevent not available \n");
232          return;
233       }
234    } else if (fAnalysisType == "AOD") {
235       lAODevent = dynamic_cast<AliAODEvent*>(InputEvent());
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");
244       return;
245    }
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
285    if (!lESDevent) {
286       Printf("ERROR: lESDevent not available \n");
287       return;
288    }
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
325
326    Double_t b = lESDevent->GetMagneticField();
327    Int_t trackNumber = lESDevent->GetNumberOfTracks();
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;
343       }
344
345
346    }
347
348    Double_t xPrimaryVertex = vtxT3D->GetX();
349    Double_t yPrimaryVertex = vtxT3D->GetY();
350    Double_t zPrimaryVertex = vtxT3D->GetZ();
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;
414       Float_t lambdaDCA = 0; // DCA between Lambda and Primary Vertex
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
427       if ( (TMath::Abs(pTrackXi->GetSign()) - TMath::Abs(nTrackXi->GetSign()) ) < 0.1) continue;
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
435       if ((TMath::Abs(v->GetEffMass() - massLambda)) < 0.2) {
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
444          strness = 1; lambdaMass = v->GetEffMass(); lambdaP = v->P(); lambdaPt = v->Pt(); lambdaDCA = v->GetD(xPrimaryVertex, yPrimaryVertex, zPrimaryVertex); v0cospointangle = v->GetV0CosineOfPointingAngle(); v0daughtersDCA = v->GetDcaV0Daughters();
445
446
447
448       }
449
450
451       v->ChangeMassHypothesis(kLambda0Bar); // the v0 must be Anti Lambda
452
453
454
455       if ((TMath::Abs(v->GetEffMass() - massLambda)) < 0.2) {
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
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();
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
487          //Bool_t *IsOkTrack = IsSelected(btrk); //for Alberto's PID
488
489          //Bool_t *okTrack = new Bool_t[5];
490
491          //for (Int_t k = 0; k < 5; k++) okTrack[k] = IsOkTrack[k];
492          
493          IsSelected(btrk);
494
495
496          Int_t bachCharge = btrk->Charge();
497          Float_t pionDCA = TMath::Abs(btrk->GetD(xPrimaryVertex, yPrimaryVertex, b)); // DCA between Bachelor and Primary Vertex
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
522          Short_t mcTrue=0;
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) ||
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;
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
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) ;
592
593          lambda.SetPxPyPzE(lLambdaMomX, lLambdaMomY, lLambdaMomZ, eLambda);
594          pion.SetPxPyPzE(lBachMomX, lBachMomY, lBachMomZ, ePion);
595
596          sigma = lambda + pion;
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
603          invmass = sigma.M();
604
605
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);
611
612          ncasc++;
613
614          //delete IsOkTrack;
615          //delete okTrack;
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
627 }
628
629 //________________________________________________________________________
630
631 //Bool_t *AliAnalysisTaskSigma1385::IsSelected(AliESDtrack *track)
632 void AliAnalysisTaskSigma1385::IsSelected(AliESDtrack *track)
633 {
634 //
635 //
636    //Bool_t okTrack[5];
637
638    //for (Int_t i = 0; i < 5; i++) okTrack[i] = kFALSE;
639    for (Int_t i = 0; i < 5; i++) fOkTrack[i] = kFALSE;
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);
660    isITSSA = ((!isTPC) && (track->IsOn(AliESDtrack::kITSrefit)) && (!track->IsOn(AliESDtrack::kITSpureSA))); //(!isTPC && (status & AliESDtrack::kITSrefit) != 0 && (status & AliESDtrack::kITSpureSA) == 0 && (status & AliESDtrack::kITSpid) != 0);
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");
667       //return okTrack;
668       return;
669
670    } else if (isTPC && !fUseGlobal) {
671       AliDebug(AliLog::kDebug + 2, "Global tracks not used. Rejected");
672       //return okTrack;
673       return;
674
675    } else if (isITSSA && !fUseITSSA) {
676       AliDebug(AliLog::kDebug + 2, "ITS standalone not used. Rejected");
677       //return okTrack;
678       return;
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;
728       }
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;
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          }
747       } else {
748          okITS = kTRUE;
749       }
750    } else {
751       // if we are here, the track is surely bad
752       okITS = kFALSE;
753       okTPC = kFALSE;
754       okTOF = kFALSE;
755    }
756
757
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;
769
770    //cout<<"##########################################"<<endl;
771    //cout<<"isITSSA "<<isITSSA<< " isTPC "<<isTPC<<endl;
772    //cout<<"##########################################"<<endl;
773
774
775    //return okTrack;
776 }
777
778
779
780
781
782
783
784
785
786
787
788 //________________________________________________________________________
789 void AliAnalysisTaskSigma1385::Terminate(Option_t *)
790 {
791    // Draw result to the screen
792    // Called once at the end of the query
793 }
794
795
796