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