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   Double_t ITSsample[4];
 453   for(Int_t i=0;i<4;i++)ITSsample[i]=0;
 454
 455   Double_t ITSsamplePos[4];
 456   for(Int_t i=0;i<4;i++)ITSsamplePos[i]=0;
 457
 458   Double_t ITSsampleNeg[4];
 459   for(Int_t i=0;i<4;i++)ITSsampleNeg[i]=0;
 460
 461   Double_t  pinTPC=0.,poutTPC=0.,TPCSignal=0.;
 462   Double_t xPrimaryVertex=0.,yPrimaryVertex=0.,zPrimaryVertex=0.;
 463
 464   ULong_t  status=0;
 465   ULong_t  statusT=0;
 466   ULong_t  statusPi=0;
 467
 468   Bool_t   isTPC=kFALSE,isTOF=kFALSE,isTOFHe3=kFALSE,isTOFPi=kFALSE;
 469
 470   Double_t fPos[3]={0.,0.,0.};
 471   Double_t runNumber=0.;
 472   Double_t evNumber=0.;
 473
 474   Int_t id0           = 0, id1          = 0;
 475   Double_t mcDecayPosXD0 = 0, mcDecayPosYD0 = 0, mcDecayPosR = 0, mcDecayPosZD0 =0, mcDecayPosRho=0.;
 476   Double_t mcDecayPosXD1 = 0, mcDecayPosYD1 = 0, mcDecayPosZD1 =0;
 477
 478   Double_t lEtaCurrentPart =0., lPtCurrentPart = 0.,lThetaCurrentPart = 0., lPhiCurrentPart = 0.;
 479   Int_t iCurrentMother = 0;
 480   Double_t mcPosX = 0., mcPosY = 0.,mcPosZ = 0., mcPosR = 0.;
 481
 482   Double_t lPdgCurrentDaughter0 = 0, lPdgCurrentDaughter1= 0., lPdgCurrentMother=0.,lPdgCurrentDaughter =0;
 483
 484   Double_t PxD0 = 0, PyD0 = 0,PzD0 = 0;
 485   Double_t PxD1 = 0, PyD1 = 0,PzD1 = 0;
 486
 487   Int_t lNbMCPrimary        = 0;
 488   Int_t lNbMCPart           = 0;
 489   Int_t lPdgcodeCurrentPart = 0;
 490   //!----------------------------------------------------------------
 491
 492   //! A set of very loose parameters for cuts
 493
 494   Double_t fgChi2max=33.;     //! max chi2
 495   Double_t fgDNmin=0.05;      //! min imp parameter for the 1st daughter = 500um
 496   Double_t fgDCAmax=1.;       //! max DCA between the daughter tracks in cm
 497   Double_t fgCPAmin=0.9;      //! min cosine of V0's pointing angle
 498   Double_t fgRmin=0.1;        //! min radius of the fiducial volume = 1 mm
 499   Double_t fgRmax=200.;       //! max radius of the fiducial volume = 2 m
 500
 501   //------------------------------------------
 502   // create pointer to event
 503
 504   AliVEvent *event = InputEvent();
 505   if (!event) { Printf("ERROR: Could not retrieve event"); return; }
 506
 507
 508
 509 //   AliVEvent *lESDevent = InputEvent();
 510 //   if (!lESDevent) {
 511 //     Printf("ERROR: Could not retrieve event");
 512 //     return;
 513 //  }
 514
 515   Info("AliAnalysisTaskHelium3PiMC","Starting UserExec");
 516
 517   SetDataType("SIM");
 518   AliStack *stack=0;
 519
 520   if(fDataType == "SIM") {
 521
 522     // Main loop
 523     // Called for EACH event
 524     AliMCEvent *mcEvent = MCEvent();
 525     if (!mcEvent) {
 526       Printf("ERROR: Could not retrieve MC event");
 527       return;
 528     }
 529
 530     Printf("MC particles: %d", mcEvent->GetNumberOfTracks());
 531
 532     // set up a stack for use in check for primary/stable particles
 533     stack = mcEvent->Stack();
 534     if( !stack ) { Printf( "Stack not available"); return; }
 535   }
 536
 537
 538   AliESDEvent *lESDevent = 0x0;
 539
 540   //********************************** Connect to the InputEvent ******//
 541
 542   //Int_t TrackNumber = 0;
 543   if(fAnalysisType == "ESD"){
 544     lESDevent = dynamic_cast<AliESDEvent*>(event);
 545     if (!lESDevent) {
 546       Printf("ERROR: lESDevent not available \n");
 547       return;
 548     }
 549   }
 550
 551   //*****************//
 552   //*   Centrality  *//
 553   //*****************//
 554
 555   AliCentrality *centrality = lESDevent->GetCentrality();
 556
 557   Float_t percentile=centrality->GetCentralityPercentile("V0M");
 558
 559   //------------------------------
 560
 561   runNumber = lESDevent->GetRunNumber();
 562   evNumber =lESDevent->GetEventNumberInFile();
 563
 564   //---------------------
 565
 566   //  Int_t primary = stack->GetNprimary();
 567   Int_t label =-1;
 568
 569   lNbMCPrimary = stack->GetNprimary();
 570   lNbMCPart    = stack->GetNtrack();
 571
 572   fHistMCMultiplicityTracks->Fill(lNbMCPart);     //histo
 573
 574   TArrayD MomPionsMC(lNbMCPart);        //Neg pions
 575   Int_t nPionsMC=0;
 576   TArrayD MomHeMC(lNbMCPart);        //helium3
 577   Int_t nHeMC=0;
 578
 579   //------ Trimomento pion
 580   TArrayD PxPionsMC(lNbMCPart);
 581   Int_t nPxPionsMC=0;
 582   TArrayD PyPionsMC(lNbMCPart);
 583   Int_t nPyPionsMC=0;
 584   TArrayD PzPionsMC(lNbMCPart);
 585   Int_t nPzPionsMC=0;
 586   //------ Trimomento He
 587   TArrayD PxHeMC(lNbMCPart);
 588   Int_t nPxHeMC=0;
 589   TArrayD PyHeMC(lNbMCPart);
 590   Int_t nPyHeMC=0;
 591   TArrayD PzHeMC(lNbMCPart);
 592   Int_t nPzHeMC=0;
 593
 594   //Mass Definition
 595
 596   Double_t        Helium3Mass = 2.80839;
 597   Double_t        PionMass = 0.13957;
 598
 599   TLorentzVector  vPionMC,vHeliumMC,vSumMC;
 600   TLorentzVector  vPionMum,vHeliumMum,vSumMum;
 601   TLorentzVector  vPionRec,vHeliumRec,vSumRec;
 602   Bool_t isTwoBody=kFALSE;
 603
 604   for (Int_t iMC=0; iMC<stack->GetNtrack(); iMC++)
 605     {
 606       TParticle *p0 = stack->Particle(iMC);
 607
 608       if (!p0) {
 609         //Printf("ERROR: particle with label %d not found in stack (mc loop)", iMc);
 610         continue;
 611       }
 612
 613       lPdgcodeCurrentPart  = p0->GetPdgCode();
 614
 615       if(lPdgcodeCurrentPart == 1000020030 || lPdgcodeCurrentPart == -1000020030 ){
 616
 617         MomHeMC[nHeMC++]=p0->P();
 618
 619         PxHeMC[nPxHeMC++]=p0->Px();
 620         PyHeMC[nPyHeMC++]=p0->Py();
 621         PzHeMC[nPzHeMC++]=p0->Pz();
 622
 623         fhHeMCRigidity->Fill(p0->P()/2);
 624       }
 625
 626       if(lPdgcodeCurrentPart == 211 || lPdgcodeCurrentPart == -211 ){
 627
 628         MomPionsMC[nPionsMC++]=p0->P();
 629
 630         PxPionsMC[nPxPionsMC++]=p0->Px();
 631         PyPionsMC[nPyPionsMC++]=p0->Py();
 632         PzPionsMC[nPzPionsMC++]=p0->Pz();
 633
 634         fhPioneMC->Fill(p0->P());
 635       }
 636
 637       if ( lPdgcodeCurrentPart == 1010010030 || lPdgcodeCurrentPart == -1010010030 ){
 638
 639         lEtaCurrentPart   = p0->Eta();
 640         lPtCurrentPart    = p0->Pt();
 641         lThetaCurrentPart = p0->Theta();
 642         lPhiCurrentPart   = p0->Phi();
 643         iCurrentMother    = p0->GetFirstMother();
 644
 645         fHistMCEta->Fill(lEtaCurrentPart);
 646         fHistMCPt->Fill(lPtCurrentPart);
 647         fHistMCTheta->Fill(lThetaCurrentPart);
 648
 649         if (iCurrentMother == -1){lPdgCurrentMother=0; } else {lPdgCurrentMother = stack->Particle(iCurrentMother)->GetPdgCode();}
 650
 651         mcPosX = p0->Vx();
 652         mcPosY = p0->Vy();
 653         mcPosZ = p0->Vz();
 654         mcPosR = TMath::Sqrt(mcPosX*mcPosX+mcPosY*mcPosY);
 655
 656         isTwoBody=kFALSE;
 657
 658         for(Int_t i=p0->GetFirstDaughter(); i<= p0->GetLastDaughter(); i++){
 659           TParticle *pDaughter = stack->Particle(i);
 660           lPdgCurrentDaughter= pDaughter->GetPdgCode();
 661           cout<<lPdgCurrentDaughter<<endl;
 662           if(lPdgCurrentDaughter == 1000020030 || lPdgCurrentDaughter ==-1000020030 ){
 663             isTwoBody=kTRUE;
 664
 665           }
 666         }
 667
 668         if(isTwoBody){
 669
 670           for(Int_t i=p0->GetFirstDaughter(); i<= p0->GetLastDaughter(); i++){
 671
 672             TParticle *pDaughter = stack->Particle(i);
 673
 674             lPdgCurrentDaughter= pDaughter->GetPdgCode();
 675
 676             if(lPdgCurrentDaughter == 211 || lPdgCurrentDaughter == -211 ){
 677               id0 = i;
 678             }
 679
 680             if(lPdgCurrentDaughter == 1000020030 || lPdgCurrentDaughter == -1000020030 ){
 681               id1 = i;
 682             }
 683           }
 684
 685           TParticle *pDaughter0 = stack->Particle(id0);
 686           TParticle *pDaughter1 = stack->Particle(id1);
 687           lPdgCurrentDaughter0 = pDaughter0->GetPdgCode();
 688           lPdgCurrentDaughter1 = pDaughter1->GetPdgCode();
 689
 690           // Decay Radius and Production Radius
 691
 692           if ( id0 <= lNbMCPart && id0 > 0 && id1 <= lNbMCPart && id1 > 0) {
 693
 694             lPdgCurrentDaughter0 = pDaughter0->GetPdgCode();
 695             lPdgCurrentDaughter1 = pDaughter1->GetPdgCode();
 696
 697             PxD0 = pDaughter0->Px();
 698             PyD0 = pDaughter0->Py();
 699             PzD0 = pDaughter0->Pz();
 700
 701             PxD1 = pDaughter1->Px();
 702             PyD1 = pDaughter1->Py();
 703             PzD1 = pDaughter1->Pz();
 704
 705             mcDecayPosXD0 = pDaughter0->Vx();
 706             mcDecayPosYD0 = pDaughter0->Vy();
 707             mcDecayPosZD0 = pDaughter0->Vz();
 708
 709             mcDecayPosXD1 = pDaughter0->Vx();
 710             mcDecayPosYD1 = pDaughter0->Vy();
 711             mcDecayPosZD1 = pDaughter0->Vz();
 712
 713             mcDecayPosR = TMath::Sqrt(mcDecayPosXD0*mcDecayPosXD0+mcDecayPosYD0*mcDecayPosYD0);
 714             fHistMCDecayPosition->Fill(mcDecayPosR);
 715
 716             mcDecayPosRho = TMath::Sqrt(mcDecayPosXD0*mcDecayPosXD0+mcDecayPosYD0*mcDecayPosYD0+mcDecayPosZD0*mcDecayPosZD0);
 717             fHistMCDecayRho->Fill(mcDecayPosRho);
 718
 719             //---- Initial mass Test
 720
 721             vHeliumMum.SetXYZM(PxD1,PyD1,PzD1,Helium3Mass);
 722             vPionMum.SetXYZM(PxD0,PyD0,PzD0,PionMass);
 723             vSumMum=vHeliumMum+vPionMum;
 724
 725             fhInvMassMum->Fill(vSumMum.M());
 726
 727             //Ntupla hyper triton
 728
 729             vett2[0]=(Float_t)lESDevent->GetRunNumber();
 730             vett2[1]=(Float_t)lESDevent->GetEventNumberInFile();
 731             vett2[2]=(Float_t)iMC;
 732             vett2[3]=(Float_t)percentile;
 733             vett2[4]=(Float_t)mcPosX;
 734             vett2[5]=(Float_t)mcPosY;
 735             vett2[6]=(Float_t)mcPosZ;
 736             vett2[7]=(Float_t)lPdgcodeCurrentPart;
 737             vett2[8]=(Float_t)iCurrentMother;
 738             vett2[9]=(Float_t)mcDecayPosXD0;
 739             vett2[10]=(Float_t)mcDecayPosYD0;
 740             vett2[11]=(Float_t)mcDecayPosZD0;
 741             vett2[12]=(Float_t)mcDecayPosXD1;
 742             vett2[13]=(Float_t)mcDecayPosYD1;
 743             vett2[14]=(Float_t)mcDecayPosZD1;
 744             vett2[15]=(Float_t)mcDecayPosRho;
 745             vett2[16]=(Float_t)lEtaCurrentPart;
 746             vett2[17]=(Float_t)p0->Y();
 747             vett2[18]=(Float_t)lThetaCurrentPart;
 748             vett2[19]=(Float_t)lPhiCurrentPart;
 749             vett2[20]=(Float_t)p0->Px();
 750             vett2[21]=(Float_t)p0->Py();
 751             vett2[22]=(Float_t)p0->Pz();
 752             vett2[23]=(Float_t)lPdgCurrentDaughter0;
 753             vett2[24]=(Float_t)lPdgCurrentDaughter1;
 754             vett2[25]=(Float_t)PxD0; //pion
 755             vett2[26]=(Float_t)PyD0;
 756             vett2[27]=(Float_t)PzD0;
 757             vett2[28]=(Float_t)PxD1; //He3
 758             vett2[29]=(Float_t)PyD1;
 759             vett2[30]=(Float_t)PzD1;
 760
 761             fNtuple2->Fill(vett2);
 762
 763           }//if check daughters index
 764         }//is twobody
 765       } // if mother
 766     } // Kinetic Track loop ends here
 767
 768   // Loop phase - space
 769
 770   Double_t HeMomMC =0;
 771   Double_t PionMomMC=0;
 772   Double_t PxHeMc=0, PyHeMc=0,PzHeMc=0;
 773   Double_t PxPionMc=0, PyPionMc=0,PzPionMc=0;
 774
 775   for(Int_t l=0; l < nHeMC; l++){
 776
 777     HeMomMC=MomHeMC[l];
 778
 779     PxHeMc=PxHeMC[l];
 780     PyHeMc=PyHeMC[l];
 781     PzHeMc=PzHeMC[l];
 782
 783     for(Int_t k=0; k < nPionsMC; k++){
 784
 785       PionMomMC=MomPionsMC[k];
 786
 787       PxPionMc=PxPionsMC[k];
 788       PyPionMc=PyPionsMC[k];
 789       PzPionMc=PzPionsMC[k];
 790
 791       fhRigidityHevsMomPiMC->Fill(HeMomMC/2,PionMomMC);
 792
 793       vHeliumMC.SetXYZM(PxHeMc,PyHeMc,PzHeMc,Helium3Mass);
 794       vPionMC.SetXYZM(PxPionMc,PyPionMc,PzPionMc,PionMass);
 795       vSumMC=vHeliumMC+vPionMC;
 796
 797       fhInvMassMC->Fill(vSumMC.M());
 798
 799     }
 800
 801   } // end loop phase space
 802
 803   //-------------- RECONSTRUCTION -------------------
 804
 805   fHistEventMultiplicity->Fill(0);
 806
 807   Double_t lMagneticField=lESDevent->GetMagneticField();
 808
 809   Int_t TrackNumber = -1;
 810
 811   // ANALISYS reconstructed tracks
 812
 813   // Primary vertex cut
 814
 815   const AliESDVertex *vtx = lESDevent->GetPrimaryVertexTracks();
 816
 817   if(vtx->GetNContributors()<1) {
 818
 819     // SPD vertex cut
 820     vtx = lESDevent->GetPrimaryVertexSPD();
 821
 822     if(vtx->GetNContributors()<1) {
 823       Info("AliAnalysisTaskHelium3PiMC","No good vertex, skip event");
 824       return; // NO GOOD VERTEX, SKIP EVENT
 825     }
 826   }
 827
 828   fHistEventMultiplicity->Fill(1); // analyzed events with PV
 829
 830   xPrimaryVertex=vtx->GetXv();
 831   yPrimaryVertex=vtx->GetYv();
 832   zPrimaryVertex=vtx->GetZv();
 833
 834   TrackNumber = lESDevent->GetNumberOfTracks();
 835   fHistTrackMultiplicity->Fill(TrackNumber); //tracce per evento
 836
 837   fHistPercentileVsTrackNumber->Fill(percentile,TrackNumber);
 838
 839   if (TrackNumber<2) return;
 840
 841   fHistEventMultiplicity->Fill(2);
 842
 843   //Find Pair candidates
 844
 845   TArrayI PionsTPC(TrackNumber);        //Neg pions
 846   Int_t nPionsTPC=0;
 847
 848   TArrayI HeTPC(TrackNumber);        //helium3
 849   Int_t nHeTPC=0;
 850
 851   // find pairs candidates phase daughter tracks rec
 852
 853   TArrayD MomPionsRec(TrackNumber);        //Neg pions
 854   Int_t nPionsRec=0;
 855
 856   TArrayD MomHeRec(TrackNumber);        //helium3
 857   Int_t nHeRec=0;
 858
 859   //------ Trimomento pion
 860   TArrayD PxPionsRec(TrackNumber);
 861   Int_t nPxPionsRec=0;
 862   TArrayD PyPionsRec(TrackNumber);
 863   Int_t nPyPionsRec=0;
 864   TArrayD PzPionsRec(TrackNumber);
 865   Int_t nPzPionsRec=0;
 866
 867   //------ Trimomento He
 868   TArrayD PxHeRec(TrackNumber);
 869   Int_t nPxHeRec=0;
 870   TArrayD PyHeRec(TrackNumber);
 871   Int_t nPyHeRec=0;
 872   TArrayD PzHeRec(TrackNumber);
 873   Int_t nPzHeRec=0;
 874
 875   Float_t impactXY=-999, impactZ=-999;
 876   Float_t impactXYpi=-999, impactZpi=-999;
 877
 878   Int_t PDGCodePos=0;
 879   Int_t PDGCodeNeg=0;
 880   Int_t motherPDGNeg=0;
 881   Int_t motherPDGPos=0;
 882
 883   Int_t mumpdg=-100;
 884
 885   //!   SELECTIONS:
 886   //! - No ITSpureSA
 887   //! - ITSrefit
 888   //! - TPCrefit
 889   //! - ITSmap !=0
 890
 891   // ******************* Track Cuts Definitions ********************//
 892
 893   //! ITS
 894
 895   AliESDtrackCuts* esdtrackCutsITS = new AliESDtrackCuts("esdtrackCutsITS");
 896   esdtrackCutsITS->SetRequireITSStandAlone(kFALSE);
 897   esdtrackCutsITS->SetRequireITSPureStandAlone(kFALSE);
 898
 899   //! TPC
 900
 901   Int_t    minclsTPC=60;
 902   Double_t maxchi2perTPCcl=5.;
 903
 904   AliESDtrackCuts* esdtrackCutsTPC = new AliESDtrackCuts("esdtrackCutsTPC");
 905   esdtrackCutsTPC->SetRequireTPCRefit(kTRUE);
 906   esdtrackCutsTPC->SetAcceptKinkDaughters(kFALSE);
 907   esdtrackCutsTPC->SetMinNClustersTPC(minclsTPC);
 908   esdtrackCutsTPC->SetMaxChi2PerClusterTPC(maxchi2perTPCcl);
 909
 910   //*************************************************************
 911
 912   for (Int_t j=0; j<TrackNumber; j++) { //loop on tracks
 913
 914     AliESDtrack *esdtrack=lESDevent->GetTrack(j);
 915
 916     if(!esdtrack) {
 917       AliError(Form("ERROR: Could not retrieve esdtrack %d",j));
 918       continue;
 919     }
 920
 921     hBBTPCnoCuts->Fill(esdtrack->GetSign()*esdtrack->P(),esdtrack->GetTPCsignal());
 922
 923     // ************** Track cuts ****************
 924
 925     status  = (ULong_t)esdtrack->GetStatus();
 926
 927     isTPC   = (((status) & (AliESDtrack::kTPCin))  != 0);
 928     isTOF   = ((((status) & (AliESDtrack::kTOFout)) != 0) && (((status) & (AliESDtrack::kTIME)) != 0));
 929
 930     Bool_t IsTrackAcceptedTPC =  esdtrackCutsTPC->AcceptTrack(esdtrack);
 931     Bool_t IsTrackAcceptedITS =  esdtrackCutsITS->AcceptTrack(esdtrack);
 932
 933     if (!(IsTrackAcceptedTPC && IsTrackAcceptedITS)) continue;
 934
 935     //----------------------------------------------
 936
 937     //****** Cuts from  AliV0Vertex.cxx *************
 938
 939     Double_t d=esdtrack->GetD(xPrimaryVertex,yPrimaryVertex,lMagneticField);
 940     //    if (TMath::Abs(d)<fgDPmin) continue;
 941     if (TMath::Abs(d)>fgRmax) continue;
 942
 943     //---- (Usefull) Stuff
 944
 945     TPCSignal=esdtrack->GetTPCsignal();
 946
 947     if (TPCSignal<10)continue;
 948
 949     if(!isTPC)continue;
 950
 951     if(!esdtrack->GetTPCInnerParam())continue;
 952
 953     AliExternalTrackParam trackIn(*esdtrack->GetInnerParam());
 954     pinTPC = trackIn.GetP();
 955
 956     poutTPC=pinTPC;
 957
 958     fhBBTPC->Fill(pinTPC*esdtrack->GetSign(),TPCSignal);
 959
 960     d=esdtrack->GetD(xPrimaryVertex,yPrimaryVertex,lMagneticField);
 961     //    if (TMath::Abs(d)<fgDPmin) continue;
 962     if (TMath::Abs(d)>fgRmax) continue;
 963
 964     label = TMath::Abs(esdtrack->GetLabel());
 965
 966     if (label>=10000000) {
 967       // Underlying event. 10000000 is the
 968       // value of fkMASKSTEP in AliRunDigitizer
 969       cout <<"Strange, there should be no underlying event"<<endl;
 970     }
 971
 972     else {
 973       if (label>=lNbMCPart) {
 974         cout <<"Strange, label outside the range"<< endl;
 975         continue;
 976       }
 977     }
 978
 979     TParticle * part = stack->Particle(label);
 980
 981     Int_t PDGCode=part->GetPdgCode();
 982     Int_t mumid = part->GetFirstMother();
 983
 984     if(mumid>-1){
 985       TParticle *mother=(TParticle*)stack->Particle(mumid);
 986       mumpdg = mother->GetPdgCode();
 987     }
 988
 989     if(PDGCode==-211)
 990       fhBBTPCNegativePions->Fill(esdtrack->GetSign()*esdtrack->P(),esdtrack->GetTPCsignal());
 991
 992     if(PDGCode==+211)
 993       fhBBTPCPositivePions->Fill(esdtrack->GetSign()*esdtrack->P(),esdtrack->GetTPCsignal());
 994
 995
 996     //    if(PDGCode == 211){
 997
 998     if(PDGCode==-211 || PDGCode==+211){
 999
1000       PionsTPC[nPionsTPC++]=j;
1001
1002       esdtrack->GetImpactParameters(impactXY, impactZ);
1003       fhPiDCAXY->Fill(impactXY);
1004       fhPiDCAZ->Fill(impactZ);
1005
1006       MomPionsRec[nPionsRec++]=esdtrack->P();
1007
1008       PxPionsRec[nPxPionsRec++]=esdtrack->Px();
1009       PyPionsRec[nPyPionsRec++]=esdtrack->Py();
1010       PzPionsRec[nPzPionsRec++]=esdtrack->Pz();
1011
1012     }
1013
1014     if(PDGCode==1000020030 ||PDGCode==-1000020030 ){
1015
1016
1017       HeTPC[nHeTPC++]=j;
1018
1019       fhBBTPCHe3->Fill(esdtrack->GetSign()*esdtrack->P(),esdtrack->GetTPCsignal());
1020
1021       esdtrack->GetImpactParameters(impactXY, impactZ);
1022       fhHeDCAXY->Fill(impactXY);
1023       fhHeDCAZ->Fill(impactZ);
1024
1025       MomHeRec[nHeRec++]=esdtrack->P();
1026
1027       PxHeRec[nPxHeRec++]=esdtrack->Px();
1028       PyHeRec[nPyHeRec++]=esdtrack->Py();
1029       PzHeRec[nPzHeRec++]=esdtrack->Pz();
1030
1031     }
1032
1033   }  //     ! track
1034
1035   //-------------- LOOP pairs 1 -------------
1036   // Fill phase space and inva mass before cuts
1037
1038   Double_t HeMomRec =0;
1039   Double_t PionMomRec=0;
1040   Double_t PxHeReco=0, PyHeReco=0,PzHeReco=0;
1041   Double_t PxPionReco=0, PyPionReco=0,PzPionReco=0;
1042
1043   for(Int_t l=0; l < nHeRec; l++){
1044
1045     HeMomRec=MomHeRec[l];
1046
1047     PxHeReco=PxHeRec[l];
1048     PyHeReco=PyHeRec[l];
1049     PzHeReco=PzHeRec[l];
1050
1051     for(Int_t k=0; k < nPionsRec; k++){
1052
1053       PionMomRec=MomPionsRec[k];
1054
1055       PxPionReco=PxPionsRec[k];
1056       PyPionReco=PyPionsRec[k];
1057       PzPionReco=PzPionsRec[k];
1058
1059       fhRigidityHevsMomPiRec->Fill(HeMomRec,PionMomRec);
1060
1061       vHeliumRec.SetXYZM(2*PxHeReco,2*PyHeReco,2*PzHeReco,Helium3Mass);
1062       vPionRec.SetXYZM(PxPionReco,PyPionReco,PzPionReco,PionMass);
1063       vSumRec=vHeliumRec+vPionRec;
1064
1065       fhInvMassRec->Fill(vSumRec.M());
1066     }
1067
1068   } // fine loop phase space
1069
1070   //--------------- LOOP PAIRS ----------------------//
1071
1072   Double_t        DcaHeToPrimVertex=0;
1073   Double_t        DcaPionToPrimVertex=0;
1074
1075   impactXY=-999, impactZ=-999;
1076   impactXYpi=-999, impactZpi=-999;
1077
1078   // Track
1079
1080   AliESDtrack  *PionTrack = 0x0;
1081   AliESDtrack  *HeTrack = 0x0;
1082
1083   // Vettori per il PxPyPz
1084
1085   Double_t momPionVett[3];
1086   for(Int_t i=0;i<3;i++)momPionVett[i]=0;
1087
1088   Double_t momHeVett[3];
1089   for(Int_t i=0;i<3;i++)momHeVett[i]=0;
1090
1091   //At SV
1092
1093   Double_t momPionVettAt[3];
1094   for(Int_t i=0;i<3;i++)momPionVettAt[i]=0;
1095
1096   Double_t momHeVettAt[3];
1097   for(Int_t i=0;i<3;i++)momHeVettAt[i]=0;
1098
1099   Bool_t   IsHeITSRefit,IsPiITSRefit ;
1100
1101   //----------- My 2nd Vertex Finder
1102
1103   for (Int_t k=0; k < nPionsTPC; k++) {                           //! Pions Loop
1104
1105     DcaPionToPrimVertex=0.;
1106     DcaHeToPrimVertex=0;
1107
1108     Int_t PionIdx=PionsTPC[k];
1109
1110     PionTrack=lESDevent->GetTrack(PionIdx);
1111
1112     statusPi = (ULong_t)PionTrack->GetStatus();
1113     IsPiITSRefit = ((statusPi) & (AliESDtrack::kITSrefit));
1114
1115     Int_t labelPi = TMath::Abs(PionTrack->GetLabel());
1116     TParticle * partNeg = stack->Particle(labelPi);
1117     PDGCodeNeg=partNeg->GetPdgCode();
1118
1119     Int_t mumidNeg = partNeg->GetFirstMother();
1120     if(mumidNeg>-1){
1121       TParticle *motherNeg=(TParticle*)stack->Particle(mumidNeg);
1122       motherPDGNeg = motherNeg->GetPdgCode();
1123     }
1124
1125     if (PionTrack)
1126       DcaPionToPrimVertex = TMath::Abs(PionTrack->GetD(xPrimaryVertex, yPrimaryVertex,lMagneticField)); //OK
1127
1128     if(DcaPionToPrimVertex<0.1)continue;
1129
1130     AliExternalTrackParam trackInPion(*PionTrack);
1131
1132     for (Int_t i=0; i<nHeTPC; i++){                               //! Helium Loop
1133
1134       Int_t HeIdx=HeTPC[i];
1135
1136       HeTrack=lESDevent->GetTrack(HeIdx);
1137
1138       statusT= (ULong_t)HeTrack->GetStatus();
1139       IsHeITSRefit = ((statusT) & (AliESDtrack::kITSrefit));
1140
1141       Int_t labelHe = TMath::Abs(HeTrack->GetLabel());
1142       TParticle * partPos = stack->Particle(labelHe);
1143       PDGCodePos=partPos->GetPdgCode();
1144
1145       Int_t mumidPos = partPos->GetFirstMother();
1146       if(mumidPos>-1){
1147         TParticle *motherPos=(TParticle*)stack->Particle(mumidPos);
1148         motherPDGPos = motherPos->GetPdgCode();
1149       }
1150
1151       if (HeTrack)
1152         DcaHeToPrimVertex = TMath::Abs(HeTrack->GetD(xPrimaryVertex, yPrimaryVertex,lMagneticField)); //OK
1153
1154       AliExternalTrackParam trackInHe(*HeTrack);
1155
1156       HeTrack->PxPyPz(momHeVett);
1157       PionTrack->PxPyPz(momPionVett);
1158
1159       vHeliumRec.SetXYZM(2*momHeVett[0],2*momHeVett[1],2*momHeVett[2],Helium3Mass);
1160       vPionRec.SetXYZM(momPionVett[0],momPionVett[1],momPionVett[2],PionMass);
1161       vSumRec=vHeliumRec+vPionRec;
1162
1163       fhInvMassRec1->Fill(vSumRec.M());
1164
1165       fhInvMassRec2->Fill(vSumRec.M());
1166
1167       if ( DcaPionToPrimVertex < fgDNmin)                //OK
1168         if ( DcaHeToPrimVertex < fgDNmin) continue;    //OK
1169
1170       fhInvMassRec3->Fill(vSumRec.M());
1171
1172       Double_t xn, xp;
1173       Double_t dca=0.;
1174
1175       dca= PionTrack->GetDCA(HeTrack,lMagneticField,xn,xp); //!distance between two tracks (Neg to Pos)
1176       fHistProvaDCA->Fill(xn-xp,dca);
1177       if (dca > fgDCAmax) continue;
1178
1179       fhInvMassRec4->Fill(vSumRec.M());
1180
1181       if ((xn+xp) > 2*fgRmax) continue;
1182       if ((xn+xp) < 2*fgRmin) continue;
1183       fhInvMassRec5->Fill(vSumRec.M());
1184
1185       //CORREZIONE da AliV0Vertex
1186
1187       Bool_t corrected=kFALSE;
1188       if ((trackInPion.GetX() > 3.) && (xn < 3.)) {
1189         //correct for the beam pipe material
1190         corrected=kTRUE;
1191       }
1192       if ((trackInHe.GetX() > 3.) && (xp < 3.)) {
1193         //correct for the beam pipe material
1194         corrected=kTRUE;
1195       }
1196       if (corrected) {
1197         dca=trackInPion.GetDCA(&trackInHe,lMagneticField,xn,xp);
1198         if (dca > fgDCAmax) continue;
1199         if ((xn+xp) > 2*fgRmax) continue;
1200         if ((xn+xp) < 2*fgRmin) continue;
1201       }
1202
1203       //=============================================//
1204       // Make  a  "V0" with Tracks                   //
1205       //=============================================//
1206
1207       trackInPion.PropagateTo(xn,lMagneticField);
1208       trackInHe.PropagateTo(xp,lMagneticField);
1209
1210       AliESDv0 vertex(trackInPion,PionIdx,trackInHe,HeIdx);
1211       if (vertex.GetChi2V0() > fgChi2max) continue;
1212       fhInvMassRec6->Fill(vSumRec.M());
1213
1214       Float_t CosPointingAngle=vertex.GetV0CosineOfPointingAngle(xPrimaryVertex,yPrimaryVertex,zPrimaryVertex); //PointingAngle
1215       if (CosPointingAngle < fgCPAmin) continue;
1216
1217       fhInvMassRec7->Fill(vSumRec.M());
1218
1219       vertex.SetDcaV0Daughters(dca);
1220       vertex.SetV0CosineOfPointingAngle(CosPointingAngle);
1221
1222       fPos[0]=vertex.Xv();
1223       fPos[1]=vertex.Yv();
1224       fPos[2]=vertex.Zv();
1225
1226
1227
1228       Double_t raggio=TMath::Sqrt(fPos[0]*fPos[0]+fPos[1]*fPos[1]+fPos[2]*fPos[2]);
1229       HeTrack->GetPxPyPzAt(raggio,lMagneticField,momHeVettAt);
1230       PionTrack->GetPxPyPzAt(raggio,lMagneticField,momPionVettAt);
1231
1232       //------------------------------------------------------------------------//
1233
1234       HeTrack->GetImpactParameters(impactXY, impactZ);
1235
1236       PionTrack->GetImpactParameters(impactXYpi, impactZpi);
1237
1238       Float_t timeTOFHe= HeTrack->GetTOFsignal();                 // ps
1239       Float_t trackLenghtTOFHe= HeTrack->GetIntegratedLength();  // cm
1240
1241       Float_t timeTOFPi= PionTrack->GetTOFsignal();                 // ps
1242       Float_t trackLenghtTOFPi= PionTrack->GetIntegratedLength();  // cm
1243
1244       //----------------------------------------------------------------------//
1245
1246       vett1[0]=(Float_t)runNumber;
1247       vett1[1]=(Float_t)evNumber;
1248       vett1[2]=(Float_t)lNbMCPart;
1249       vett1[3]=(Float_t)percentile;
1250       vett1[4]=(Float_t)xPrimaryVertex; //PRIMARY
1251       vett1[5]=(Float_t)yPrimaryVertex;
1252       vett1[6]=(Float_t)zPrimaryVertex;
1253       vett1[7]=(Float_t)fPos[0]; //SECONDARY
1254       vett1[8]=(Float_t)fPos[1];
1255       vett1[9]=(Float_t)fPos[2];
1256       vett1[10]=(Float_t)dca;           //between 2 tracks
1257       vett1[11]=(Float_t)CosPointingAngle;          //cosPointingAngle da V0
1258       vett1[12]=(Float_t)vertex.GetD(xPrimaryVertex,yPrimaryVertex,zPrimaryVertex);
1259       vett1[13]=(Float_t)HeTrack->GetSign(); //He
1260       vett1[14]=(Float_t)trackInHe.GetP();
1261       vett1[15]=(Float_t)HeTrack->GetTPCsignal();
1262       vett1[16]=(Float_t)DcaHeToPrimVertex;
1263       vett1[17]=(Float_t)HeTrack->Eta();
1264       vett1[18]=(Float_t)momHeVett[0];
1265       vett1[19]=(Float_t)momHeVett[1];
1266       vett1[20]=(Float_t)momHeVett[2];
1267       vett1[21]=(Float_t)momHeVettAt[0];
1268       vett1[22]=(Float_t)momHeVettAt[1];
1269       vett1[23]=(Float_t)momHeVettAt[2];
1270       vett1[24]=(Float_t)HeTrack->GetTPCNcls();
1271       vett1[25]=(Float_t)impactXY;
1272       vett1[26]=(Float_t)impactZ;
1273       vett1[27]=(Float_t)isTOFHe3;
1274       vett1[28]=(Float_t)(trackLenghtTOFHe/timeTOFHe)/2.99792458e-2;
1275       vett1[29]=(Float_t)HeTrack->GetITSClusterMap();
1276       vett1[30]=(Float_t)IsHeITSRefit;
1277       vett1[31]=(Float_t)PionTrack->GetSign(); //Pion
1278       vett1[32]=(Float_t)trackInPion.GetP();
1279       vett1[33]=(Float_t)PionTrack->GetTPCsignal();
1280       vett1[34]=(Float_t)DcaPionToPrimVertex;
1281       vett1[35]=(Float_t)PionTrack->Eta();
1282       vett1[36]=(Float_t)momPionVett[0];
1283       vett1[37]=(Float_t)momPionVett[1];
1284       vett1[38]=(Float_t)momPionVett[2];
1285       vett1[39]=(Float_t)momPionVettAt[0];
1286       vett1[40]=(Float_t)momPionVettAt[1];
1287       vett1[41]=(Float_t)momPionVettAt[2];
1288       vett1[42]=(Float_t)PionTrack->GetTPCNcls();
1289       vett1[43]=(Float_t)impactXYpi;
1290       vett1[44]=(Float_t)impactZpi;
1291       vett1[45]=(Float_t)isTOFPi;
1292       vett1[46]=(Float_t)(trackLenghtTOFPi/timeTOFPi)/2.99792458e-2;
1293       vett1[47]=(Float_t)PionTrack->GetITSClusterMap();
1294       vett1[48]=(Float_t)IsPiITSRefit;
1295       vett1[49]=(Float_t)PDGCodeNeg;
1296       vett1[50]=(Float_t)PDGCodePos;
1297       vett1[51]=(Float_t)motherPDGNeg;
1298       vett1[52]=(Float_t)motherPDGPos;
1299       vett1[53]=(Float_t)labelPi;
1300       vett1[54]=(Float_t)labelHe;
1301       vett1[55]=(Float_t)mumidNeg;
1302       vett1[56]=(Float_t)mumidPos;
1303
1304       fNtuple1->Fill(vett1);
1305
1306     }// positive TPC
1307
1308   } //negative tpc
1309
1310   PostData(1,fListHistCascade);
1311
1312 }// end userexec
1313
1314
1315 //________________________________________________________________________
1316 //
1317 void AliAnalysisTaskHelium3PiMC::Terminate(Option_t *)
1318 {
1319   //  Draw result to the screen
1320   //   Called once at the end of the query
1321 }
1322