PWGLF/STRANGENESS/Hypernuclei/AliAnalysisTaskHelium3PiMC.cxx

   1 /**************************************************************************
   2  * Contributors are not mentioned at all.                                 *
   3  *                                                                        *
   4  * Permission to use, copy, modify and distribute this software and its   *
   5  * documentation strictly for non-commercial purposes is hereby granted   *
   6  * without fee, provided that the above copyright notice appears in all   *
   7  * copies and that both the copyright notice and this permission notice   *
   8  * appear in the supporting documentation. The authors make no claims     *
   9  * about the suitability of this software for any purpose. It is          *
  10  * provided "as is" without express or implied warranty.                  *
  11  **************************************************************************/
  12 //
  13 //-----------------------------------------------------------------
  14 //                 AliAnalysisTaskHelium3Pi class
  15 //-----------------------------------------------------------------
  16
  17
  18 class TTree;
  19 class TParticle;
  20 class TVector3;
  21
  22 #include "AliAnalysisManager.h"
  23 #include <AliMCEventHandler.h>
  24 #include <AliMCEvent.h>
  25 #include <AliStack.h>
  26
  27 class AliESDVertex;
  28 class AliAODVertex;
  29 class AliESDv0;
  30 class AliAODv0;
  31 class AliCascadeVertexer;
  32
  33 #include <iostream>
  34 #include "AliAnalysisTaskSE.h"
  35 #include "TList.h"
  36 #include "TH1.h"
  37 #include "TH2.h"
  38 #include "TH3.h"
  39 #include "TNtuple.h"
  40 #include "TGraph.h"
  41 #include "TCutG.h"
  42 #include "TF1.h"
  43 #include "TCanvas.h"
  44 #include "TMath.h"
  45 #include "TChain.h"
  46 #include "Riostream.h"
  47 #include "AliLog.h"
  48 #include "AliCascadeVertexer.h"
  49 #include "AliESDEvent.h"
  50 #include "AliESDtrack.h"
  51 #include "AliExternalTrackParam.h"
  52 #include "AliAODEvent.h"
  53 #include "AliInputEventHandler.h"
  54 #include "AliESDcascade.h"
  55 #include "AliAODcascade.h"
  56 #include "AliAnalysisTaskHelium3PiMC.h"
  57 #include "AliESDtrackCuts.h"
  58 #include "AliCentrality.h"
  59 #include "TString.h"
  60 #include <TDatime.h>
  61 #include <TRandom3.h>
  62 using std::endl;
  63 using std::cout;
  64
  65 const Int_t AliAnalysisTaskHelium3PiMC::fgNrot = 15;
  66
  67
  68 ClassImp(AliAnalysisTaskHelium3PiMC)
  69
  70 //________________________________________________________________________
  71 AliAnalysisTaskHelium3PiMC::AliAnalysisTaskHelium3PiMC()
  72 : AliAnalysisTaskSE(),
  73   fAnalysisType("ESD"),
  74   fCollidingSystems(0),
  75   fDataType("SIM"),
  76   fListHistCascade(0),
  77   fHistEventMultiplicity(0),
  78   fHistTrackMultiplicity(0),
  79   fHistMCMultiplicityTracks(0),
  80   fHistMCEta(0),
  81   fHistMCPt(0),
  82   fHistMCTheta(0),
  83   fHistMCDecayPosition(0),
  84   fHistMCDecayRho(0),
  85   fhRigidityHevsMomPiMC(0),
  86   fhRigidityHevsMomPiRec(0),
  87   fhInvMassMC(0),
  88   fhInvMassMum(0),
  89   fhInvMassRec(0),
  90   fhInvMassRec1(0),
  91   fhInvMassRec2(0),
  92   fhInvMassRec3(0),
  93   fhInvMassRec4(0),
  94   fhInvMassRec5(0),
  95   fhInvMassRec6(0),
  96   fhInvMassRec7(0),
  97   fhHeMCRigidity(0),
  98   fhPioneMC(0),
  99   hBBTPCnoCuts(0),
 100   fhBBTPC(0),
 101   fhBBTPCNegativePions(0),
 102   fhBBTPCPositivePions(0),
 103   fhBBTPCHe3(0),
 104   fHistProvaDCA(0),
 105   fHistPercentileVsTrackNumber(0),
 106   fhHeDCAXY(0),
 107   fhHeDCAZ(0),
 108   fhPiDCAXY(0),
 109   fhPiDCAZ(0),
 110   hITSClusterMap(0),
 111   fNtuple1(0),
 112   fNtuple2(0)
 113
 114 {
 115   // Dummy Constructor(0);
 116 }
 117
 118 //________________________________________________________________________
 119 AliAnalysisTaskHelium3PiMC::AliAnalysisTaskHelium3PiMC(const char *name)
 120   : AliAnalysisTaskSE(name),
 121     fAnalysisType("ESD"),
 122     fCollidingSystems(0),
 123     fDataType("SIM"),
 124     fListHistCascade(0),
 125     fHistEventMultiplicity(0),
 126     fHistTrackMultiplicity(0),
 127     fHistMCMultiplicityTracks(0),
 128     fHistMCEta(0),
 129     fHistMCPt(0),
 130     fHistMCTheta(0),
 131     fHistMCDecayPosition(0),
 132     fHistMCDecayRho(0),
 133     fhRigidityHevsMomPiMC(0),
 134     fhRigidityHevsMomPiRec(0),
 135     fhInvMassMC(0),
 136     fhInvMassMum(0),
 137     fhInvMassRec(0),
 138     fhInvMassRec1(0),
 139     fhInvMassRec2(0),
 140     fhInvMassRec3(0),
 141     fhInvMassRec4(0),
 142     fhInvMassRec5(0),
 143     fhInvMassRec6(0),
 144     fhInvMassRec7(0),
 145     fhHeMCRigidity(0),
 146     fhPioneMC(0),
 147     hBBTPCnoCuts(0),
 148     fhBBTPC(0),
 149     fhBBTPCNegativePions(0),
 150     fhBBTPCPositivePions(0),
 151     fhBBTPCHe3(0),
 152     fHistProvaDCA(0),
 153     fHistPercentileVsTrackNumber(0),
 154     fhHeDCAXY(0),
 155     fhHeDCAZ(0),
 156     fhPiDCAXY(0),
 157     fhPiDCAZ(0),
 158     hITSClusterMap(0),
 159     fNtuple1(0),
 160     fNtuple2(0)
 161
 162
 163 {
 164   // Define input and output slots here
 165   // Input slot #0 works with a TChain
 166   //DefineInput(0, TChain::Class());
 167   // Output slot #0 writes into a TList container (Cascade)
 168   DefineOutput(1, TList::Class());
 169 }
 170 //_______________________________________________________
 171 AliAnalysisTaskHelium3PiMC::~AliAnalysisTaskHelium3PiMC()
 172 {
 173   // Destructor
 174   if (fListHistCascade) {
 175     delete fListHistCascade;
 176     fListHistCascade = 0;
 177   }
 178
 179 }
 180 //=================DEFINITION BETHE BLOCH==============================
 181
 182 Double_t AliAnalysisTaskHelium3PiMC::BetheBloch(Double_t betaGamma,Double_t charge,Bool_t isPbPb) {
 183
 184   Double_t kp1, kp2, kp3, kp4, kp5;
 185
 186   if(isPbPb){
 187
 188     //pass2 //to be checked
 189     kp1 = 5.2*charge*charge;
 190     kp2 = 8.98482806165147636e+00;
 191     kp3 = 1.54000000000000005e-05;
 192     kp4 = 2.30445734159456084e+00;
 193     kp5 = 2.25624744086878559e+00;
 194
 195
 196   }
 197
 198   else{
 199
 200     //pass2 // to be defined
 201     kp1 = 5.2*charge*charge;
 202     kp2 = 8.98482806165147636e+00;
 203     kp3 = 1.54000000000000005e-05;
 204     kp4 = 2.30445734159456084e+00;
 205     kp5 = 2.25624744086878559e+00;
 206
 207   }
 208
 209   Double_t beta = betaGamma / TMath::Sqrt(1.0 + betaGamma * betaGamma);
 210
 211   Double_t aa = TMath::Power(beta, kp4);
 212   Double_t bb = TMath::Power(1.0 / betaGamma, kp5);
 213
 214   bb = TMath::Log(kp3 + bb);
 215
 216   Double_t out = (kp2 - aa - bb) * kp1 / aa;
 217
 218   return out;
 219
 220 }
 221
 222 //==================DEFINITION OF OUTPUT OBJECTS==============================
 223
 224 void AliAnalysisTaskHelium3PiMC::UserCreateOutputObjects()
 225 {
 226   fListHistCascade = new TList();
 227   fListHistCascade->SetOwner();  // IMPORTANT!
 228
 229   if(! fHistEventMultiplicity ){
 230     fHistEventMultiplicity   = new TH1F( "fHistEventMultiplicity" , "Nb of Events" , 6 , -1, 5 );
 231     fHistEventMultiplicity->GetXaxis()->SetTitle("Event Type");
 232     fListHistCascade->Add(fHistEventMultiplicity);
 233   }
 234
 235   if(! fHistTrackMultiplicity ){
 236
 237     fHistTrackMultiplicity   = new TH1F( "fHistTrackMultiplicity" , "Nb of Tracks" , 25000,0, 25000 );
 238     fHistTrackMultiplicity->GetXaxis()->SetTitle("Number of tracks");
 239     fListHistCascade->Add(fHistTrackMultiplicity);
 240   }
 241
 242   if(! fHistMCMultiplicityTracks){
 243     fHistMCMultiplicityTracks =new TH1F("fHistMCMultiplicityTracks","MC Multiplicity Tracks",1000,0,1000);
 244     fHistMCMultiplicityTracks->GetXaxis()->SetTitle("MC Number of tracks");
 245     fListHistCascade->Add(fHistMCMultiplicityTracks);
 246   }
 247   if(!fHistMCEta ){
 248     fHistMCEta=new TH1F("fHistMCEta","MC eta",1000,-3,3);
 249     fHistMCEta->GetXaxis()->SetTitle("Injected Eta");
 250     fListHistCascade->Add(fHistMCEta);
 251   }
 252   if(! fHistMCPt){
 253     fHistMCPt =new TH1F("fHistMCPt","MC pt",1000,0,20);
 254     fHistMCPt->GetXaxis()->SetTitle("Injected Pt");
 255     fListHistCascade->Add(fHistMCPt);
 256   }
 257   if(!fHistMCTheta ){
 258     fHistMCTheta=new TH1F("fHistMCTheta","MC theta",1000,-6,6);
 259     fHistMCTheta->GetXaxis()->SetTitle("Injected Theta");
 260     fListHistCascade->Add(fHistMCTheta);
 261   }
 262   if(!fHistMCDecayPosition){
 263     fHistMCDecayPosition =new TH1F("fHistMCDecayPosition","MC Decay Position",10000,0,1000);
 264     fHistMCDecayPosition->GetXaxis()->SetTitle("Decay Position");
 265     fListHistCascade->Add(fHistMCDecayPosition);
 266   }
 267   if(!fHistMCDecayRho ){
 268     fHistMCDecayRho  =new TH1F("fHistMCDecayRho","MC decay position 3d",10000,0,1000);
 269     fHistMCDecayRho->GetXaxis()->SetTitle("Decay rho");
 270     fListHistCascade->Add(fHistMCDecayRho);
 271   }
 272
 273   if(!fhRigidityHevsMomPiMC ){
 274     fhRigidityHevsMomPiMC=new TH2F("fhRigidityHevsMomPiMC","Rigidity He vs Mom Pi MC",20,0,10,300,0,30);
 275     fhRigidityHevsMomPiMC->GetXaxis()->SetTitle("He3 Rigidity");
 276     fhRigidityHevsMomPiMC->GetYaxis()->SetTitle("Pi momentum");
 277     fListHistCascade->Add(fhRigidityHevsMomPiMC);
 278   }
 279
 280   if(! fhRigidityHevsMomPiRec){
 281     fhRigidityHevsMomPiRec=new TH2F("fhRigidityHevsMomPiRec","Rigidity He vs Mom Pi Rec",20,0,10,300,0,30);
 282     fhRigidityHevsMomPiRec->GetXaxis()->SetTitle("He3 Rigidity");
 283     fhRigidityHevsMomPiRec->GetYaxis()->SetTitle("Pi momentum");
 284     fListHistCascade->Add(fhRigidityHevsMomPiRec);
 285   }
 286
 287   if(!fhInvMassMC){
 288     fhInvMassMC=new TH1F("fhInvMassMC","fhInvMassMC",800,2.,6.);
 289     fhInvMassMC->GetXaxis()->SetTitle("(He3,#pi) InvMass");
 290     fListHistCascade->Add(fhInvMassMC);
 291   }
 292
 293   if(!fhInvMassMum){
 294     fhInvMassMum=new TH1F("fhInvMassMum","fhInvMassMum",800,2.,6.);
 295     fhInvMassMum->GetXaxis()->SetTitle("(He3,#pi) InvMass");
 296     fListHistCascade->Add(fhInvMassMum);
 297   }
 298
 299   if(!fhInvMassRec){
 300     fhInvMassRec=new TH1F("fhInvMassRec","fhInvMassRec",800,2.,6.);
 301     fhInvMassRec->GetXaxis()->SetTitle("(He3,#pi) InvMass");
 302     fListHistCascade->Add(fhInvMassRec);
 303   }
 304
 305   if(!fhInvMassRec1){
 306     fhInvMassRec1=new TH1F("fhInvMassRec1","No Altri tagli",800,2.,6.);
 307     fhInvMassRec1->GetXaxis()->SetTitle("(He3,#pi) InvMass");
 308     fListHistCascade->Add(fhInvMassRec1);
 309   }
 310   if(!fhInvMassRec2){
 311     fhInvMassRec2=new TH1F("fhInvMassRec2","DCA pi > 0.1",800,2.,6.);
 312     fhInvMassRec2->GetXaxis()->SetTitle("(He3,#pi) InvMass");
 313     fListHistCascade->Add(fhInvMassRec2);
 314   }
 315   if(!fhInvMassRec3){
 316     fhInvMassRec3=new TH1F("fhInvMassRec3","DCA He > 0.05",800,2.,6.);
 317     fhInvMassRec3->GetXaxis()->SetTitle("(He3,#pi) InvMass");
 318     fListHistCascade->Add(fhInvMassRec3);
 319   }
 320   if(!fhInvMassRec4){
 321     fhInvMassRec4=new TH1F("fhInvMassRec4","DCA tracks < 1 cm",800,2.,6.);
 322     fhInvMassRec4->GetXaxis()->SetTitle("(He3,#pi) InvMass");
 323     fListHistCascade->Add(fhInvMassRec4);
 324   }
 325   if(!fhInvMassRec5){
 326     fhInvMassRec5=new TH1F("fhInvMassRec5","Condizione xn+xp",800,2.,6.);
 327     fhInvMassRec5->GetXaxis()->SetTitle("(He3,#pi) InvMass");
 328     fListHistCascade->Add(fhInvMassRec5);
 329   }
 330   if(!fhInvMassRec6){
 331     fhInvMassRec6=new TH1F("fhInvMassRec6","Ho fatto V0 ",800,2.,6.);
 332     fhInvMassRec6->GetXaxis()->SetTitle("(He3,#pi) InvMass");
 333     fListHistCascade->Add(fhInvMassRec6);
 334   }
 335   if(!fhInvMassRec7){
 336     fhInvMassRec7=new TH1F("fhInvMassRec7","V0+Taglio CPA",800,2.,6.);
 337     fhInvMassRec7->GetXaxis()->SetTitle("(He3,#pi) InvMass");
 338     fListHistCascade->Add(fhInvMassRec7);
 339   }
 340
 341   if(!fhHeMCRigidity ){
 342     fhHeMCRigidity=new TH1F("fhHeMCRigidity","He3 rigidity distribution",200,0,20);
 343     fhHeMCRigidity->GetXaxis()->SetTitle("He3 rigidity");
 344     fListHistCascade->Add(fhHeMCRigidity);
 345   }
 346   if(!fhPioneMC ){
 347     fhPioneMC=new TH1F("hPioneMC","Pion mom distribution",200,0,50);
 348     fhPioneMC->GetXaxis()->SetTitle("Pion momentum");
 349     fListHistCascade->Add(fhPioneMC);
 350   }
 351
 352   if(!hBBTPCnoCuts ){
 353     hBBTPCnoCuts=new TH2F("hBBTPCnoCuts","scatterPlot TPC no cuts",2000,-10,10,1000,0,3000);
 354     hBBTPCnoCuts->GetXaxis()->SetTitle("p/Z (GeV/#it{c})");
 355     hBBTPCnoCuts->GetYaxis()->SetTitle("TPC Signal (a.u)");
 356     fListHistCascade->Add(hBBTPCnoCuts);
 357   }
 358   if(!fhBBTPC ){
 359     fhBBTPC=new TH2F("fhBBTPC","scatterPlot TPC",2000,-10,10,1000,0,3000);
 360     fhBBTPC->GetXaxis()->SetTitle("p/Z (GeV/#it{c})");
 361     fhBBTPC->GetYaxis()->SetTitle("TPC Signal (a.u)");
 362     fListHistCascade->Add(fhBBTPC);
 363   }
 364   if(!fhBBTPCNegativePions ){
 365     fhBBTPCNegativePions=new TH2F("fhBBTPCNegativePions","scatterPlot Neg Pions",2000,-10,10,1000,0,3000);
 366     fhBBTPCNegativePions->GetXaxis()->SetTitle("p/Z (GeV/#it{c})");
 367     fhBBTPCNegativePions->GetYaxis()->SetTitle("TPC Signal (a.u)");
 368     fListHistCascade->Add(fhBBTPCNegativePions);
 369   }
 370   if(!fhBBTPCPositivePions ){
 371     fhBBTPCPositivePions=new TH2F("fhBBTPCPositivePions","scatterPlot Pos Pions",2000,-10,10,1000,0,3000);
 372     fhBBTPCPositivePions->GetXaxis()->SetTitle("p/Z (GeV/#it{c})");
 373     fhBBTPCPositivePions->GetYaxis()->SetTitle("TPC Signal (a.u)");
 374     fListHistCascade->Add(fhBBTPCPositivePions);
 375   }
 376   if(!fhBBTPCHe3 ){
 377     fhBBTPCHe3=new TH2F("fhBBTPCHe3","scatterPlot TPC -  He3",2000,-10,10,1000,0,3000);
 378     fhBBTPCHe3->GetXaxis()->SetTitle("p/Z (GeV/#it{c})");
 379     fhBBTPCHe3->GetYaxis()->SetTitle("TPC Signal (a.u)");
 380     fListHistCascade->Add(fhBBTPCHe3);
 381   }
 382   if(!fHistProvaDCA ){
 383     fHistProvaDCA=new TH2F("fHistProvaDCA","fHistProvaDCA",1000,-50,50,1000,0,100);
 384     fHistProvaDCA->GetXaxis()->SetTitle("xn+xp");
 385     fHistProvaDCA->GetYaxis()->SetTitle("dca tracks");
 386     fListHistCascade->Add(fHistProvaDCA);
 387   }
 388
 389   if(!hITSClusterMap ){
 390     hITSClusterMap=new TH1F("hITSClusterMap","hITSClusterMap",65,-1,64);
 391     fListHistCascade->Add(hITSClusterMap);
 392   }
 393
 394   if(!fHistPercentileVsTrackNumber){
 395     fHistPercentileVsTrackNumber=new TH2F("fHistPercentileVsTrackNumber","fHistPercentileVsTrackNumber",120,-3,117,2500,0,25000);
 396     fHistPercentileVsTrackNumber->GetXaxis()->SetTitle("Percentile");
 397     fHistPercentileVsTrackNumber->GetYaxis()->SetTitle("Tracks Number");
 398     fListHistCascade->Add(fHistPercentileVsTrackNumber);
 399   }
 400
 401   if(!fhHeDCAXY){
 402     fhHeDCAXY=new TH1F("fhHeDCAXY","fhHeDCAXY",800,-4,4);
 403     fListHistCascade->Add(fhHeDCAXY);
 404   }
 405   if(!fhHeDCAZ){
 406     fhHeDCAZ=new TH1F("fhHeDCAZ","fhHeDCAZ",800,-30,30);
 407     fListHistCascade->Add(fhHeDCAZ);
 408   }
 409   if(!fhPiDCAXY){
 410     fhPiDCAXY=new TH1F("fhPiDCAXY","fhPiDCAXY",800,-4,4);
 411     fListHistCascade->Add(fhPiDCAXY);
 412   }
 413   if(!fhPiDCAZ){
 414     fhPiDCAZ=new TH1F("fhPiDCAZ","fhPiDCAZ",800,-30,30);
 415     fListHistCascade->Add(fhPiDCAZ);
 416   }
 417
 418   if(! fNtuple1 ) {
 419     fNtuple1 = new TNtuple("fNtuple1","Ntuple1","runNumber:evNumber:TrackNumber:percentile:xPrimaryVertex:yPrimaryVertex:zPrimaryVertex:xSecondaryVertex:ySecondaryVertex:zSecondaryVertex:dcaTracks:CosPointingAngle:DCAV0toPrimaryVertex:HeSign:HepInTPC:HeTPCsignal:DcaHeToPrimVertex:HeEta:momHex:momHey:momHez:momHeAtSVx:momHeAtSVy:momHeAtSVz:HeTPCNcls:HeimpactXY:HeimpactZ:isTOFHe:HeBeta:HeITSClusterMap:IsHeITSRefit:PionSign:PionpInTPC:PionTPCsignal:DcaPionToPrimVertex:PionEta:momPionx:momPiony:momPionz:momNegPionAtSVx:momNegPionAtSVy:momNegPionAtSVz:PionTPCNcls:PionimpactXY:PionimpactZ:isTOFPion:PionBeta:PionITSClusterMap:IsPiITSRefit:PDGCodeNeg:PDCCodePos:motherPDGNeg:motherPDGPos:labelPi:labelHe:mumidNeg:mumidPos");
 420
 421     fListHistCascade->Add(fNtuple1);
 422   }
 423
 424   if(! fNtuple2 ) {
 425
 426     fNtuple2 = new TNtuple("fNtuple2","Ntuple2","run:event:iMC:Centrality:PVx:PVy:PVz:PDGcodeMum:MotherIndex:SVxD0:SVyD0:SVzD0:SVxD1:SVyD1:SVzD1:SV3d:EtaMum:YMum:ThetaMum:PhiMum:PxMum:PyMum:PzMum:PdgDaughter0:PdgDaughter1:PxD0:PyD0:PzD0:PxD1:PyD1:PzD1");
 427
 428     fListHistCascade->Add(fNtuple2);
 429   }
 430
 431   PostData(1,fListHistCascade);
 432
 433 }// end UserCreateOutputObjects
 434
 435
 436
 437 //====================== USER EXEC ========================
 438
 439 void AliAnalysisTaskHelium3PiMC::UserExec(Option_t *)
 440 {
 441   //_______________________________________________________________________
 442
 443   //!*********************!//
 444   //!  Define variables   !//
 445   //!*********************!//
 446   Float_t vett1[60];
 447   for(Int_t i=0;i<60;i++) vett1[i]=0;
 448
 449   Float_t vett2[40];
 450   for(Int_t i=0;i<40;i++) vett2[i]=0;
 451
 452
 453   Double_t  pinTPC=0.,TPCSignal=0.;
 454   Double_t xPrimaryVertex=0.,yPrimaryVertex=0.,zPrimaryVertex=0.;
 455
 456   ULong_t  status=0;
 457   ULong_t  statusT=0;
 458   ULong_t  statusPi=0;
 459
 460   Bool_t   isTPC=kFALSE,isTOFHe3=kFALSE,isTOFPi=kFALSE;
 461
 462   Double_t fPos[3]={0.,0.,0.};
 463   Double_t runNumber=0.;
 464   Double_t evNumber=0.;
 465
 466   Int_t id0           = 0, id1          = 0;
 467   Double_t mcDecayPosXD0 = 0, mcDecayPosYD0 = 0, mcDecayPosR = 0, mcDecayPosZD0 =0, mcDecayPosRho=0.;
 468   Double_t mcDecayPosXD1 = 0, mcDecayPosYD1 = 0, mcDecayPosZD1 =0;
 469
 470   Double_t lEtaCurrentPart =0., lPtCurrentPart = 0.,lThetaCurrentPart = 0., lPhiCurrentPart = 0.;
 471   Int_t iCurrentMother = 0;
 472   Double_t mcPosX = 0., mcPosY = 0.,mcPosZ = 0.;
 473
 474   Double_t lPdgCurrentDaughter0 = 0, lPdgCurrentDaughter1= 0., /*lPdgCurrentMother=0.,*/lPdgCurrentDaughter =0;
 475
 476   Double_t PxD0 = 0, PyD0 = 0,PzD0 = 0;
 477   Double_t PxD1 = 0, PyD1 = 0,PzD1 = 0;
 478
 479   Int_t lNbMCPart           = 0;
 480   Int_t lPdgcodeCurrentPart = 0;
 481   //!----------------------------------------------------------------
 482
 483   //! A set of very loose parameters for cuts
 484
 485   Double_t fgChi2max=33.;     //! max chi2
 486   Double_t fgDNmin=0.05;      //! min imp parameter for the 1st daughter = 500um
 487   Double_t fgDCAmax=1.;       //! max DCA between the daughter tracks in cm
 488   Double_t fgCPAmin=0.9;      //! min cosine of V0's pointing angle
 489   Double_t fgRmin=0.1;        //! min radius of the fiducial volume = 1 mm
 490   Double_t fgRmax=200.;       //! max radius of the fiducial volume = 2 m
 491
 492   //------------------------------------------
 493   // create pointer to event
 494
 495   AliVEvent *event = InputEvent();
 496   if (!event) { Printf("ERROR: Could not retrieve event"); return; }
 497
 498
 499
 500 //   AliVEvent *lESDevent = InputEvent();
 501 //   if (!lESDevent) {
 502 //     Printf("ERROR: Could not retrieve event");
 503 //     return;
 504 //  }
 505
 506   Info("AliAnalysisTaskHelium3PiMC","Starting UserExec");
 507
 508   SetDataType("SIM");
 509   AliStack *stack = 0;
 510   if(fDataType == "SIM") {
 511
 512     // Main loop
 513     // Called for EACH event
 514     AliMCEvent *mcEvent = MCEvent();
 515     if (!mcEvent) {
 516       Printf("ERROR: Could not retrieve MC event");
 517       return;
 518     }
 519
 520     Printf("MC particles: %d", mcEvent->GetNumberOfTracks());
 521
 522     // set up a stack for use in check for primary/stable particles
 523     stack = mcEvent->Stack();
 524     if( !stack ) { Printf( "Stack not available"); return; }
 525   }
 526
 527   else{
 528     Printf( "This Task Works Only on Simulation");
 529     return;
 530   }
 531   AliESDEvent *lESDevent = 0;
 532
 533   //********************************** Connect to the InputEvent ******//
 534
 535   //Int_t TrackNumber = 0;
 536   if(fAnalysisType == "ESD"){
 537     lESDevent = dynamic_cast<AliESDEvent*>(event);
 538     if (!lESDevent) {
 539       Printf("ERROR: lESDevent not available \n");
 540       return;
 541     }
 542   }
 543
 544   else{
 545     Printf("This Analysis Works Only for ESD\n");
 546     return;
 547   }
 548
 549   //*****************//
 550   //*   Centrality  *//
 551   //*****************//
 552
 553   AliCentrality *centrality = (AliCentrality*)lESDevent->GetCentrality();
 554
 555   Float_t percentile=centrality->GetCentralityPercentile("V0M");
 556
 557   //------------------------------
 558
 559   runNumber = lESDevent->GetRunNumber();
 560   evNumber =lESDevent->GetEventNumberInFile();
 561
 562   //---------------------
 563
 564   //  Int_t primary = stack->GetNprimary();
 565   Int_t label =-1;
 566
 567   lNbMCPart    = stack->GetNtrack();
 568
 569   fHistMCMultiplicityTracks->Fill(lNbMCPart);     //histo
 570
 571   TArrayD MomPionsMC(lNbMCPart);        //Neg pions
 572   Int_t nPionsMC=0;
 573   TArrayD MomHeMC(lNbMCPart);        //helium3
 574   Int_t nHeMC=0;
 575
 576   //------ Trimomento pion
 577   TArrayD PxPionsMC(lNbMCPart);
 578   Int_t nPxPionsMC=0;
 579   TArrayD PyPionsMC(lNbMCPart);
 580   Int_t nPyPionsMC=0;
 581   TArrayD PzPionsMC(lNbMCPart);
 582   Int_t nPzPionsMC=0;
 583   //------ Trimomento He
 584   TArrayD PxHeMC(lNbMCPart);
 585   Int_t nPxHeMC=0;
 586   TArrayD PyHeMC(lNbMCPart);
 587   Int_t nPyHeMC=0;
 588   TArrayD PzHeMC(lNbMCPart);
 589   Int_t nPzHeMC=0;
 590
 591   //Mass Definition
 592
 593   Double_t        Helium3Mass = 2.80839;
 594   Double_t        PionMass = 0.13957;
 595
 596   TLorentzVector  vPionMC,vHeliumMC,vSumMC;
 597   TLorentzVector  vPionMum,vHeliumMum,vSumMum;
 598   TLorentzVector  vPionRec,vHeliumRec,vSumRec;
 599   Bool_t isTwoBody=kFALSE;
 600
 601   for (Int_t iMC=0; iMC<stack->GetNtrack(); iMC++)
 602     {
 603       TParticle *p0 = stack->Particle(iMC);
 604
 605       if (!p0) {
 606         //Printf("ERROR: particle with label %d not found in stack (mc loop)", iMc);
 607         continue;
 608       }
 609
 610       lPdgcodeCurrentPart  = p0->GetPdgCode();
 611
 612       if(lPdgcodeCurrentPart == 1000020030 || lPdgcodeCurrentPart == -1000020030 ){
 613
 614         MomHeMC[nHeMC++]=p0->P();
 615
 616         PxHeMC[nPxHeMC++]=p0->Px();
 617         PyHeMC[nPyHeMC++]=p0->Py();
 618         PzHeMC[nPzHeMC++]=p0->Pz();
 619
 620         fhHeMCRigidity->Fill(p0->P()/2);
 621       }
 622
 623       if(lPdgcodeCurrentPart == 211 || lPdgcodeCurrentPart == -211 ){
 624
 625         MomPionsMC[nPionsMC++]=p0->P();
 626
 627         PxPionsMC[nPxPionsMC++]=p0->Px();
 628         PyPionsMC[nPyPionsMC++]=p0->Py();
 629         PzPionsMC[nPzPionsMC++]=p0->Pz();
 630
 631         fhPioneMC->Fill(p0->P());
 632       }
 633
 634       if ( lPdgcodeCurrentPart == 1010010030 || lPdgcodeCurrentPart == -1010010030 ){
 635
 636         lEtaCurrentPart   = p0->Eta();
 637         lPtCurrentPart    = p0->Pt();
 638         lThetaCurrentPart = p0->Theta();
 639         lPhiCurrentPart   = p0->Phi();
 640         iCurrentMother    = p0->GetFirstMother();
 641
 642         fHistMCEta->Fill(lEtaCurrentPart);
 643         fHistMCPt->Fill(lPtCurrentPart);
 644         fHistMCTheta->Fill(lThetaCurrentPart);
 645
 646         //      if (iCurrentMother == -1){lPdgCurrentMother=0; } else {lPdgCurrentMother = stack->Particle(iCurrentMother)->GetPdgCode();}
 647
 648         mcPosX = p0->Vx();
 649         mcPosY = p0->Vy();
 650         mcPosZ = p0->Vz();
 651
 652         isTwoBody=kFALSE;
 653
 654         for(Int_t i=p0->GetFirstDaughter(); i<= p0->GetLastDaughter(); i++){
 655           TParticle *pDaughter = stack->Particle(i);
 656           lPdgCurrentDaughter= pDaughter->GetPdgCode();
 657           cout<<lPdgCurrentDaughter<<endl;
 658           if(lPdgCurrentDaughter == 1000020030 || lPdgCurrentDaughter ==-1000020030 ){
 659             isTwoBody=kTRUE;
 660
 661           }
 662         }
 663
 664         if(isTwoBody){
 665
 666           for(Int_t i=p0->GetFirstDaughter(); i<= p0->GetLastDaughter(); i++){
 667
 668             TParticle *pDaughter = stack->Particle(i);
 669
 670             lPdgCurrentDaughter= pDaughter->GetPdgCode();
 671
 672             if(lPdgCurrentDaughter == 211 || lPdgCurrentDaughter == -211 ){
 673               id0 = i;
 674             }
 675
 676             if(lPdgCurrentDaughter == 1000020030 || lPdgCurrentDaughter == -1000020030 ){
 677               id1 = i;
 678             }
 679           }
 680
 681           TParticle *pDaughter0 = stack->Particle(id0);
 682           TParticle *pDaughter1 = stack->Particle(id1);
 683           lPdgCurrentDaughter0 = pDaughter0->GetPdgCode();
 684           lPdgCurrentDaughter1 = pDaughter1->GetPdgCode();
 685
 686           // Decay Radius and Production Radius
 687
 688           if ( id0 <= lNbMCPart && id0 > 0 && id1 <= lNbMCPart && id1 > 0) {
 689
 690             lPdgCurrentDaughter0 = pDaughter0->GetPdgCode();
 691             lPdgCurrentDaughter1 = pDaughter1->GetPdgCode();
 692
 693             PxD0 = pDaughter0->Px();
 694             PyD0 = pDaughter0->Py();
 695             PzD0 = pDaughter0->Pz();
 696
 697             PxD1 = pDaughter1->Px();
 698             PyD1 = pDaughter1->Py();
 699             PzD1 = pDaughter1->Pz();
 700
 701             mcDecayPosXD0 = pDaughter0->Vx();
 702             mcDecayPosYD0 = pDaughter0->Vy();
 703             mcDecayPosZD0 = pDaughter0->Vz();
 704
 705             mcDecayPosXD1 = pDaughter0->Vx();
 706             mcDecayPosYD1 = pDaughter0->Vy();
 707             mcDecayPosZD1 = pDaughter0->Vz();
 708
 709             mcDecayPosR = TMath::Sqrt(mcDecayPosXD0*mcDecayPosXD0+mcDecayPosYD0*mcDecayPosYD0);
 710             fHistMCDecayPosition->Fill(mcDecayPosR);
 711
 712             mcDecayPosRho = TMath::Sqrt(mcDecayPosXD0*mcDecayPosXD0+mcDecayPosYD0*mcDecayPosYD0+mcDecayPosZD0*mcDecayPosZD0);
 713             fHistMCDecayRho->Fill(mcDecayPosRho);
 714
 715             //---- Initial mass Test
 716
 717             vHeliumMum.SetXYZM(PxD1,PyD1,PzD1,Helium3Mass);
 718             vPionMum.SetXYZM(PxD0,PyD0,PzD0,PionMass);
 719             vSumMum=vHeliumMum+vPionMum;
 720
 721             fhInvMassMum->Fill(vSumMum.M());
 722
 723             //Ntupla hyper triton
 724
 725             vett2[0]=(Float_t)lESDevent->GetRunNumber();
 726             vett2[1]=(Float_t)lESDevent->GetEventNumberInFile();
 727             vett2[2]=(Float_t)iMC;
 728             vett2[3]=(Float_t)percentile;
 729             vett2[4]=(Float_t)mcPosX;
 730             vett2[5]=(Float_t)mcPosY;
 731             vett2[6]=(Float_t)mcPosZ;
 732             vett2[7]=(Float_t)lPdgcodeCurrentPart;
 733             vett2[8]=(Float_t)iCurrentMother;
 734             vett2[9]=(Float_t)mcDecayPosXD0;
 735             vett2[10]=(Float_t)mcDecayPosYD0;
 736             vett2[11]=(Float_t)mcDecayPosZD0;
 737             vett2[12]=(Float_t)mcDecayPosXD1;
 738             vett2[13]=(Float_t)mcDecayPosYD1;
 739             vett2[14]=(Float_t)mcDecayPosZD1;
 740             vett2[15]=(Float_t)mcDecayPosRho;
 741             vett2[16]=(Float_t)lEtaCurrentPart;
 742             vett2[17]=(Float_t)p0->Y();
 743             vett2[18]=(Float_t)lThetaCurrentPart;
 744             vett2[19]=(Float_t)lPhiCurrentPart;
 745             vett2[20]=(Float_t)p0->Px();
 746             vett2[21]=(Float_t)p0->Py();
 747             vett2[22]=(Float_t)p0->Pz();
 748             vett2[23]=(Float_t)lPdgCurrentDaughter0;
 749             vett2[24]=(Float_t)lPdgCurrentDaughter1;
 750             vett2[25]=(Float_t)PxD0; //pion
 751             vett2[26]=(Float_t)PyD0;
 752             vett2[27]=(Float_t)PzD0;
 753             vett2[28]=(Float_t)PxD1; //He3
 754             vett2[29]=(Float_t)PyD1;
 755             vett2[30]=(Float_t)PzD1;
 756
 757             fNtuple2->Fill(vett2);
 758
 759           }//if check daughters index
 760         }//is twobody
 761       } // if mother
 762     } // Kinetic Track loop ends here
 763
 764   // Loop phase - space
 765
 766   Double_t HeMomMC =0;
 767   Double_t PionMomMC=0;
 768   Double_t PxHeMc=0, PyHeMc=0,PzHeMc=0;
 769   Double_t PxPionMc=0, PyPionMc=0,PzPionMc=0;
 770
 771   for(Int_t l=0; l < nHeMC; l++){
 772
 773     HeMomMC=MomHeMC[l];
 774
 775     PxHeMc=PxHeMC[l];
 776     PyHeMc=PyHeMC[l];
 777     PzHeMc=PzHeMC[l];
 778
 779     for(Int_t k=0; k < nPionsMC; k++){
 780
 781       PionMomMC=MomPionsMC[k];
 782
 783       PxPionMc=PxPionsMC[k];
 784       PyPionMc=PyPionsMC[k];
 785       PzPionMc=PzPionsMC[k];
 786
 787       fhRigidityHevsMomPiMC->Fill(HeMomMC/2,PionMomMC);
 788
 789       vHeliumMC.SetXYZM(PxHeMc,PyHeMc,PzHeMc,Helium3Mass);
 790       vPionMC.SetXYZM(PxPionMc,PyPionMc,PzPionMc,PionMass);
 791       vSumMC=vHeliumMC+vPionMC;
 792
 793       fhInvMassMC->Fill(vSumMC.M());
 794
 795     }
 796
 797   } // end loop phase space
 798
 799   //-------------- RECONSTRUCTION -------------------
 800
 801   fHistEventMultiplicity->Fill(0);
 802
 803   Double_t lMagneticField=lESDevent->GetMagneticField();
 804
 805   Int_t TrackNumber = -1;
 806
 807   // ANALISYS reconstructed tracks
 808
 809   // Primary vertex cut
 810
 811   const AliESDVertex *vtx = lESDevent->GetPrimaryVertexTracks();
 812
 813   if(vtx->GetNContributors()<1) {
 814
 815     // SPD vertex cut
 816     vtx = lESDevent->GetPrimaryVertexSPD();
 817
 818     if(vtx->GetNContributors()<1) {
 819       Info("AliAnalysisTaskHelium3PiMC","No good vertex, skip event");
 820       return; // NO GOOD VERTEX, SKIP EVENT
 821     }
 822   }
 823
 824   fHistEventMultiplicity->Fill(1); // analyzed events with PV
 825
 826   xPrimaryVertex=vtx->GetXv();
 827   yPrimaryVertex=vtx->GetYv();
 828   zPrimaryVertex=vtx->GetZv();
 829
 830   TrackNumber = lESDevent->GetNumberOfTracks();
 831   fHistTrackMultiplicity->Fill(TrackNumber); //tracce per evento
 832
 833   fHistPercentileVsTrackNumber->Fill(percentile,TrackNumber);
 834
 835   if (TrackNumber<2) return;
 836
 837   fHistEventMultiplicity->Fill(2);
 838
 839   //Find Pair candidates
 840
 841   TArrayI PionsTPC(TrackNumber);        //Neg pions
 842   Int_t nPionsTPC=0;
 843
 844   TArrayI HeTPC(TrackNumber);        //helium3
 845   Int_t nHeTPC=0;
 846
 847   // find pairs candidates phase daughter tracks rec
 848
 849   TArrayD MomPionsRec(TrackNumber);        //Neg pions
 850   Int_t nPionsRec=0;
 851
 852   TArrayD MomHeRec(TrackNumber);        //helium3
 853   Int_t nHeRec=0;
 854
 855   //------ Trimomento pion
 856   TArrayD PxPionsRec(TrackNumber);
 857   Int_t nPxPionsRec=0;
 858   TArrayD PyPionsRec(TrackNumber);
 859   Int_t nPyPionsRec=0;
 860   TArrayD PzPionsRec(TrackNumber);
 861   Int_t nPzPionsRec=0;
 862
 863   //------ Trimomento He
 864   TArrayD PxHeRec(TrackNumber);
 865   Int_t nPxHeRec=0;
 866   TArrayD PyHeRec(TrackNumber);
 867   Int_t nPyHeRec=0;
 868   TArrayD PzHeRec(TrackNumber);
 869   Int_t nPzHeRec=0;
 870
 871   Float_t impactXY=-999, impactZ=-999;
 872   Float_t impactXYpi=-999, impactZpi=-999;
 873
 874   Int_t PDGCodePos=0;
 875   Int_t PDGCodeNeg=0;
 876   Int_t motherPDGNeg=0;
 877   Int_t motherPDGPos=0;
 878
 879   //!   SELECTIONS:
 880   //! - No ITSpureSA
 881   //! - ITSrefit
 882   //! - TPCrefit
 883   //! - ITSmap !=0
 884
 885   // ******************* Track Cuts Definitions ********************//
 886
 887   //! ITS
 888
 889   AliESDtrackCuts* esdtrackCutsITS = new AliESDtrackCuts("esdtrackCutsITS");
 890   esdtrackCutsITS->SetRequireITSStandAlone(kFALSE);
 891   esdtrackCutsITS->SetRequireITSPureStandAlone(kFALSE);
 892
 893   //! TPC
 894
 895   Int_t    minclsTPC=60;
 896   Double_t maxchi2perTPCcl=5.;
 897
 898   AliESDtrackCuts* esdtrackCutsTPC = new AliESDtrackCuts("esdtrackCutsTPC");
 899   esdtrackCutsTPC->SetRequireTPCRefit(kTRUE);
 900   esdtrackCutsTPC->SetAcceptKinkDaughters(kFALSE);
 901   esdtrackCutsTPC->SetMinNClustersTPC(minclsTPC);
 902   esdtrackCutsTPC->SetMaxChi2PerClusterTPC(maxchi2perTPCcl);
 903
 904   //*************************************************************
 905
 906   for (Int_t j=0; j<TrackNumber; j++) { //loop on tracks
 907
 908     AliESDtrack *esdtrack=lESDevent->GetTrack(j);
 909
 910     if(!esdtrack) {
 911       AliError(Form("ERROR: Could not retrieve esdtrack %d",j));
 912       continue;
 913     }
 914
 915     hBBTPCnoCuts->Fill(esdtrack->GetSign()*esdtrack->P(),esdtrack->GetTPCsignal());
 916
 917     // ************** Track cuts ****************
 918
 919     status  = (ULong_t)esdtrack->GetStatus();
 920
 921     isTPC   = (((status) & (AliESDtrack::kTPCin))  != 0);
 922
 923     Bool_t IsTrackAcceptedTPC =  esdtrackCutsTPC->AcceptTrack(esdtrack);
 924     Bool_t IsTrackAcceptedITS =  esdtrackCutsITS->AcceptTrack(esdtrack);
 925
 926     if (!(IsTrackAcceptedTPC && IsTrackAcceptedITS)) continue;
 927
 928     //----------------------------------------------
 929
 930     //****** Cuts from  AliV0Vertex.cxx *************
 931
 932     Double_t d=esdtrack->GetD(xPrimaryVertex,yPrimaryVertex,lMagneticField);
 933     //    if (TMath::Abs(d)<fgDPmin) continue;
 934     if (TMath::Abs(d)>fgRmax) continue;
 935
 936     //---- (Usefull) Stuff
 937
 938     TPCSignal=esdtrack->GetTPCsignal();
 939
 940     if (TPCSignal<10)continue;
 941
 942     if(!isTPC)continue;
 943
 944     if(!esdtrack->GetTPCInnerParam())continue;
 945
 946     AliExternalTrackParam trackIn(*esdtrack->GetInnerParam());
 947     pinTPC = trackIn.GetP();
 948
 949     fhBBTPC->Fill(pinTPC*esdtrack->GetSign(),TPCSignal);
 950
 951     d=esdtrack->GetD(xPrimaryVertex,yPrimaryVertex,lMagneticField);
 952     //    if (TMath::Abs(d)<fgDPmin) continue;
 953     if (TMath::Abs(d)>fgRmax) continue;
 954
 955     label = TMath::Abs(esdtrack->GetLabel());
 956
 957     if (label>=10000000) {
 958       // Underlying event. 10000000 is the
 959       // value of fkMASKSTEP in AliRunDigitizer
 960       cout <<"Strange, there should be no underlying event"<<endl;
 961     }
 962
 963     else {
 964       if (label>=lNbMCPart) {
 965         cout <<"Strange, label outside the range"<< endl;
 966         continue;
 967       }
 968     }
 969
 970     TParticle * part = stack->Particle(label);
 971
 972     Int_t PDGCode=part->GetPdgCode();
 973
 974     if(PDGCode==-211)
 975       fhBBTPCNegativePions->Fill(esdtrack->GetSign()*esdtrack->P(),esdtrack->GetTPCsignal());
 976
 977     if(PDGCode==+211)
 978       fhBBTPCPositivePions->Fill(esdtrack->GetSign()*esdtrack->P(),esdtrack->GetTPCsignal());
 979
 980
 981     //    if(PDGCode == 211){
 982
 983     if(PDGCode==-211 || PDGCode==+211){
 984
 985       PionsTPC[nPionsTPC++]=j;
 986
 987       esdtrack->GetImpactParameters(impactXY, impactZ);
 988       fhPiDCAXY->Fill(impactXY);
 989       fhPiDCAZ->Fill(impactZ);
 990
 991       MomPionsRec[nPionsRec++]=esdtrack->P();
 992
 993       PxPionsRec[nPxPionsRec++]=esdtrack->Px();
 994       PyPionsRec[nPyPionsRec++]=esdtrack->Py();
 995       PzPionsRec[nPzPionsRec++]=esdtrack->Pz();
 996
 997     }
 998
 999     if(PDGCode==1000020030 ||PDGCode==-1000020030 ){
1000
1001
1002       HeTPC[nHeTPC++]=j;
1003
1004       fhBBTPCHe3->Fill(esdtrack->GetSign()*esdtrack->P(),esdtrack->GetTPCsignal());
1005
1006       esdtrack->GetImpactParameters(impactXY, impactZ);
1007       fhHeDCAXY->Fill(impactXY);
1008       fhHeDCAZ->Fill(impactZ);
1009
1010       MomHeRec[nHeRec++]=esdtrack->P();
1011
1012       PxHeRec[nPxHeRec++]=esdtrack->Px();
1013       PyHeRec[nPyHeRec++]=esdtrack->Py();
1014       PzHeRec[nPzHeRec++]=esdtrack->Pz();
1015
1016     }
1017
1018   }  //     ! track
1019
1020   //-------------- LOOP pairs 1 -------------
1021   // Fill phase space and inva mass before cuts
1022
1023   Double_t HeMomRec =0;
1024   Double_t PionMomRec=0;
1025   Double_t PxHeReco=0, PyHeReco=0,PzHeReco=0;
1026   Double_t PxPionReco=0, PyPionReco=0,PzPionReco=0;
1027
1028   for(Int_t l=0; l < nHeRec; l++){
1029
1030     HeMomRec=MomHeRec[l];
1031
1032     PxHeReco=PxHeRec[l];
1033     PyHeReco=PyHeRec[l];
1034     PzHeReco=PzHeRec[l];
1035
1036     for(Int_t k=0; k < nPionsRec; k++){
1037
1038       PionMomRec=MomPionsRec[k];
1039
1040       PxPionReco=PxPionsRec[k];
1041       PyPionReco=PyPionsRec[k];
1042       PzPionReco=PzPionsRec[k];
1043
1044       fhRigidityHevsMomPiRec->Fill(HeMomRec,PionMomRec);
1045
1046       vHeliumRec.SetXYZM(2*PxHeReco,2*PyHeReco,2*PzHeReco,Helium3Mass);
1047       vPionRec.SetXYZM(PxPionReco,PyPionReco,PzPionReco,PionMass);
1048       vSumRec=vHeliumRec+vPionRec;
1049
1050       fhInvMassRec->Fill(vSumRec.M());
1051     }
1052
1053   } // fine loop phase space
1054
1055   //--------------- LOOP PAIRS ----------------------//
1056
1057   Double_t        DcaHeToPrimVertex=0;
1058   Double_t        DcaPionToPrimVertex=0;
1059
1060   impactXY=-999, impactZ=-999;
1061   impactXYpi=-999, impactZpi=-999;
1062
1063
1064   // Vettori per il PxPyPz
1065
1066   Double_t momPionVett[3];
1067   for(Int_t i=0;i<3;i++)momPionVett[i]=0;
1068
1069   Double_t momHeVett[3];
1070   for(Int_t i=0;i<3;i++)momHeVett[i]=0;
1071
1072   //At SV
1073
1074   Double_t momPionVettAt[3];
1075   for(Int_t i=0;i<3;i++)momPionVettAt[i]=0;
1076
1077   Double_t momHeVettAt[3];
1078   for(Int_t i=0;i<3;i++)momHeVettAt[i]=0;
1079
1080   Bool_t   IsHeITSRefit,IsPiITSRefit ;
1081
1082   //----------- My 2nd Vertex Finder
1083
1084   for (Int_t k=0; k < nPionsTPC; k++) {                           //! Pions Loop
1085
1086     DcaPionToPrimVertex=0.;
1087     DcaHeToPrimVertex=0;
1088
1089     Int_t PionIdx=PionsTPC[k];
1090
1091     AliESDtrack  *PionTrack=lESDevent->GetTrack(PionIdx);
1092
1093     statusPi = (ULong_t)PionTrack->GetStatus();
1094     IsPiITSRefit = ((statusPi) & (AliESDtrack::kITSrefit));
1095
1096     Int_t labelPi = TMath::Abs(PionTrack->GetLabel());
1097     TParticle * partNeg = stack->Particle(labelPi);
1098     PDGCodeNeg=partNeg->GetPdgCode();
1099
1100     Int_t mumidNeg = partNeg->GetFirstMother();
1101     if(mumidNeg>-1){
1102       TParticle *motherNeg=(TParticle*)stack->Particle(mumidNeg);
1103       motherPDGNeg = motherNeg->GetPdgCode();
1104     }
1105
1106     if (PionTrack)
1107       DcaPionToPrimVertex = TMath::Abs(PionTrack->GetD(xPrimaryVertex, yPrimaryVertex,lMagneticField)); //OK
1108
1109     if(DcaPionToPrimVertex<0.1)continue;
1110
1111     AliExternalTrackParam trackInPion(*PionTrack);
1112
1113     for (Int_t i=0; i<nHeTPC; i++){                               //! Helium Loop
1114
1115       Int_t HeIdx=HeTPC[i];
1116
1117       AliESDtrack  *HeTrack=lESDevent->GetTrack(HeIdx);
1118
1119       statusT= (ULong_t)HeTrack->GetStatus();
1120       IsHeITSRefit = ((statusT) & (AliESDtrack::kITSrefit));
1121
1122       Int_t labelHe = TMath::Abs(HeTrack->GetLabel());
1123       TParticle * partPos = stack->Particle(labelHe);
1124       PDGCodePos=partPos->GetPdgCode();
1125
1126       Int_t mumidPos = partPos->GetFirstMother();
1127       if(mumidPos>-1){
1128         TParticle *motherPos=(TParticle*)stack->Particle(mumidPos);
1129         motherPDGPos = motherPos->GetPdgCode();
1130       }
1131
1132       if (HeTrack)
1133         DcaHeToPrimVertex = TMath::Abs(HeTrack->GetD(xPrimaryVertex, yPrimaryVertex,lMagneticField)); //OK
1134
1135       AliExternalTrackParam trackInHe(*HeTrack);
1136
1137       HeTrack->PxPyPz(momHeVett);
1138       PionTrack->PxPyPz(momPionVett);
1139
1140       vHeliumRec.SetXYZM(2*momHeVett[0],2*momHeVett[1],2*momHeVett[2],Helium3Mass);
1141       vPionRec.SetXYZM(momPionVett[0],momPionVett[1],momPionVett[2],PionMass);
1142       vSumRec=vHeliumRec+vPionRec;
1143
1144       fhInvMassRec1->Fill(vSumRec.M());
1145
1146       fhInvMassRec2->Fill(vSumRec.M());
1147
1148       if ( DcaPionToPrimVertex < fgDNmin)                //OK
1149         if ( DcaHeToPrimVertex < fgDNmin) continue;    //OK
1150
1151       fhInvMassRec3->Fill(vSumRec.M());
1152
1153       Double_t xn, xp;
1154       Double_t dca=0.;
1155
1156       dca= PionTrack->GetDCA(HeTrack,lMagneticField,xn,xp); //!distance between two tracks (Neg to Pos)
1157       fHistProvaDCA->Fill(xn-xp,dca);
1158       if (dca > fgDCAmax) continue;
1159
1160       fhInvMassRec4->Fill(vSumRec.M());
1161
1162       if ((xn+xp) > 2*fgRmax) continue;
1163       if ((xn+xp) < 2*fgRmin) continue;
1164       fhInvMassRec5->Fill(vSumRec.M());
1165
1166       //CORREZIONE da AliV0Vertex
1167
1168       Bool_t corrected=kFALSE;
1169       if ((trackInPion.GetX() > 3.) && (xn < 3.)) {
1170         //correct for the beam pipe material
1171         corrected=kTRUE;
1172       }
1173       if ((trackInHe.GetX() > 3.) && (xp < 3.)) {
1174         //correct for the beam pipe material
1175         corrected=kTRUE;
1176       }
1177       if (corrected) {
1178         dca=trackInPion.GetDCA(&trackInHe,lMagneticField,xn,xp);
1179         if (dca > fgDCAmax) continue;
1180         if ((xn+xp) > 2*fgRmax) continue;
1181         if ((xn+xp) < 2*fgRmin) continue;
1182       }
1183
1184       //=============================================//
1185       // Make  a  "V0" with Tracks                   //
1186       //=============================================//
1187
1188       trackInPion.PropagateTo(xn,lMagneticField);
1189       trackInHe.PropagateTo(xp,lMagneticField);
1190
1191       AliESDv0 vertex(trackInPion,PionIdx,trackInHe,HeIdx);
1192       if (vertex.GetChi2V0() > fgChi2max) continue;
1193       fhInvMassRec6->Fill(vSumRec.M());
1194
1195       Float_t CosPointingAngle=vertex.GetV0CosineOfPointingAngle(xPrimaryVertex,yPrimaryVertex,zPrimaryVertex); //PointingAngle
1196       if (CosPointingAngle < fgCPAmin) continue;
1197
1198       fhInvMassRec7->Fill(vSumRec.M());
1199
1200       vertex.SetDcaV0Daughters(dca);
1201       vertex.SetV0CosineOfPointingAngle(CosPointingAngle);
1202
1203       fPos[0]=vertex.Xv();
1204       fPos[1]=vertex.Yv();
1205       fPos[2]=vertex.Zv();
1206
1207
1208
1209       Double_t raggio=TMath::Sqrt(fPos[0]*fPos[0]+fPos[1]*fPos[1]+fPos[2]*fPos[2]);
1210       HeTrack->GetPxPyPzAt(raggio,lMagneticField,momHeVettAt);
1211       PionTrack->GetPxPyPzAt(raggio,lMagneticField,momPionVettAt);
1212
1213       //------------------------------------------------------------------------//
1214
1215       HeTrack->GetImpactParameters(impactXY, impactZ);
1216
1217       PionTrack->GetImpactParameters(impactXYpi, impactZpi);
1218
1219       Float_t timeTOFHe= HeTrack->GetTOFsignal();                 // ps
1220       Float_t trackLenghtTOFHe= HeTrack->GetIntegratedLength();  // cm
1221
1222       Float_t timeTOFPi= PionTrack->GetTOFsignal();                 // ps
1223       Float_t trackLenghtTOFPi= PionTrack->GetIntegratedLength();  // cm
1224
1225       //----------------------------------------------------------------------//
1226
1227       vett1[0]=(Float_t)runNumber;
1228       vett1[1]=(Float_t)evNumber;
1229       vett1[2]=(Float_t)lNbMCPart;
1230       vett1[3]=(Float_t)percentile;
1231       vett1[4]=(Float_t)xPrimaryVertex; //PRIMARY
1232       vett1[5]=(Float_t)yPrimaryVertex;
1233       vett1[6]=(Float_t)zPrimaryVertex;
1234       vett1[7]=(Float_t)fPos[0]; //SECONDARY
1235       vett1[8]=(Float_t)fPos[1];
1236       vett1[9]=(Float_t)fPos[2];
1237       vett1[10]=(Float_t)dca;           //between 2 tracks
1238       vett1[11]=(Float_t)CosPointingAngle;          //cosPointingAngle da V0
1239       vett1[12]=(Float_t)vertex.GetD(xPrimaryVertex,yPrimaryVertex,zPrimaryVertex);
1240       vett1[13]=(Float_t)HeTrack->GetSign(); //He
1241       vett1[14]=(Float_t)trackInHe.GetP();
1242       vett1[15]=(Float_t)HeTrack->GetTPCsignal();
1243       vett1[16]=(Float_t)DcaHeToPrimVertex;
1244       vett1[17]=(Float_t)HeTrack->Eta();
1245       vett1[18]=(Float_t)momHeVett[0];
1246       vett1[19]=(Float_t)momHeVett[1];
1247       vett1[20]=(Float_t)momHeVett[2];
1248       vett1[21]=(Float_t)momHeVettAt[0];
1249       vett1[22]=(Float_t)momHeVettAt[1];
1250       vett1[23]=(Float_t)momHeVettAt[2];
1251       vett1[24]=(Float_t)HeTrack->GetTPCNcls();
1252       vett1[25]=(Float_t)impactXY;
1253       vett1[26]=(Float_t)impactZ;
1254       vett1[27]=(Float_t)isTOFHe3;
1255       vett1[28]=(Float_t)(trackLenghtTOFHe/timeTOFHe)/2.99792458e-2;
1256       vett1[29]=(Float_t)HeTrack->GetITSClusterMap();
1257       vett1[30]=(Float_t)IsHeITSRefit;
1258       vett1[31]=(Float_t)PionTrack->GetSign(); //Pion
1259       vett1[32]=(Float_t)trackInPion.GetP();
1260       vett1[33]=(Float_t)PionTrack->GetTPCsignal();
1261       vett1[34]=(Float_t)DcaPionToPrimVertex;
1262       vett1[35]=(Float_t)PionTrack->Eta();
1263       vett1[36]=(Float_t)momPionVett[0];
1264       vett1[37]=(Float_t)momPionVett[1];
1265       vett1[38]=(Float_t)momPionVett[2];
1266       vett1[39]=(Float_t)momPionVettAt[0];
1267       vett1[40]=(Float_t)momPionVettAt[1];
1268       vett1[41]=(Float_t)momPionVettAt[2];
1269       vett1[42]=(Float_t)PionTrack->GetTPCNcls();
1270       vett1[43]=(Float_t)impactXYpi;
1271       vett1[44]=(Float_t)impactZpi;
1272       vett1[45]=(Float_t)isTOFPi;
1273       vett1[46]=(Float_t)(trackLenghtTOFPi/timeTOFPi)/2.99792458e-2;
1274       vett1[47]=(Float_t)PionTrack->GetITSClusterMap();
1275       vett1[48]=(Float_t)IsPiITSRefit;
1276       vett1[49]=(Float_t)PDGCodeNeg;
1277       vett1[50]=(Float_t)PDGCodePos;
1278       vett1[51]=(Float_t)motherPDGNeg;
1279       vett1[52]=(Float_t)motherPDGPos;
1280       vett1[53]=(Float_t)labelPi;
1281       vett1[54]=(Float_t)labelHe;
1282       vett1[55]=(Float_t)mumidNeg;
1283       vett1[56]=(Float_t)mumidPos;
1284
1285       fNtuple1->Fill(vett1);
1286
1287     }// positive TPC
1288
1289   } //negative tpc
1290
1291   PostData(1,fListHistCascade);
1292
1293 }// end userexec
1294
1295
1296 //________________________________________________________________________
1297 //
1298 void AliAnalysisTaskHelium3PiMC::Terminate(Option_t *)
1299 {
1300   //  Draw result to the screen
1301   //   Called once at the end of the query
1302 }
1303