PWGLF/STRANGENESS/Hypernuclei/AliAnalysisTaskHelium3Pi.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 "AliAnalysisTaskHelium3Pi.h"
  57 #include "AliESDtrackCuts.h"
  58 #include "AliCentrality.h"
  59 #include "TString.h"
  60 #include <TDatime.h>
  61 #include <TRandom3.h>
  62 #include <TLorentzVector.h>
  63 //#include <AliVTrack.h>
  64
  65 const Int_t AliAnalysisTaskHelium3Pi::fgNrot = 15;
  66
  67 ClassImp(AliAnalysisTaskHelium3Pi)
  68
  69 //________________________________________________________________________
  70 AliAnalysisTaskHelium3Pi::AliAnalysisTaskHelium3Pi()
  71 : AliAnalysisTaskSE(),
  72   fAnalysisType("ESD"),
  73   fCollidingSystems(0),
  74   fESDtrackCuts(0),
  75   fDataType("REAL"),
  76   fListHist(0),
  77   fHistEventMultiplicity(0),
  78   fHistTrackMultiplicity(0),
  79   fHistTrackMultiplicityCent(0),
  80   fHistTrackMultiplicitySemiCent(0),
  81   fHistTrackMultiplicityMB(0),
  82   fHistTrackMultiplicityPVCent(0),
  83   fHistTrackMultiplicityPVSemiCent(0),
  84   fHistTrackMultiplicityPVMB(0),
  85   fhBB(0),
  86   fhTOF(0),
  87   fhMassTOF(0),
  88   fhBBPions(0),
  89   fhBBHe(0),
  90   fhNaPos(0),
  91   fhNaNeg(0),
  92   fBetavsTPCsignalPos(0),
  93   fBetavsTPCsignalNeg(0),
  94   fNtuple1(0),
  95   trunNumber(0),
  96   tbunchcross(0),
  97   torbit(0),
  98   tperiod(0),
  99   teventtype(0),
 100   tTrackNumber(0),
 101   tpercentile(0),
 102   txPrimaryVertex(0),
 103   tyPrimaryVertex(0),
 104   tzPrimaryVertex(0),
 105   txSecondaryVertex(0),
 106   tySecondaryVertex(0),
 107   tzSecondaryVertex(0),
 108   tdcaTracks(0),
 109   tCosPointingAngle(0),
 110   tDCAV0toPrimaryVertex(0),
 111   tHeSign(0),
 112   tHepInTPC(0),
 113   tHeTPCsignal(0),
 114   tDcaHeToPrimVertex(0),
 115   tHeEta(0),
 116   tmomHex(0),
 117   tmomHey(0),
 118   tmomHez(0),
 119   tmomHeAtSVx(0),
 120   tmomHeAtSVy(0),
 121   tmomHeAtSVz(0),
 122   tHeTPCNcls(0),
 123   tHeimpactXY(0),
 124   tHeimpactZ(0),
 125   tHeITSClusterMap(0),
 126   tIsHeITSRefit(0),
 127   tPionSign(0),
 128   tPionpInTPC(0),
 129   tPionTPCsignal(0),
 130   tDcaPionToPrimVertex(0),
 131   tPionEta(0),
 132   tmomPionx(0),
 133   tmomPiony(0),
 134   tmomPionz(0),
 135   tmomNegPionAtSVx(0),
 136   tmomNegPionAtSVy(0),
 137   tmomNegPionAtSVz(0),
 138   tPionTPCNcls(0),
 139   tPionimpactXY(0),
 140   tPionimpactZ(0),
 141   tPionITSClusterMap(0),
 142   tIsPiITSRefit(0),
 143   txn(0),
 144   txp(0),
 145   tchi2He(0),
 146   tchi2Pi(0),
 147   fNtuple4(0),
 148   tHelrunNumber(0),
 149   tHelBCNumber(0),
 150   tHelOrbitNumber(0),
 151   tHelPeriodNumber(0),
 152   tHeleventtype(0),
 153   tHelisHeITSrefit(0),
 154   tHelpercentile(0),
 155   tHelSign(0),
 156   tHelpinTPC(0),
 157   tHelGetTPCsignal(0),
 158   tHelPx(0),
 159   tHelPy(0),
 160   tHelPz(0),
 161   tHelEta(0),
 162   tHelisTOF(0),
 163   tHelpoutTPC(0),
 164   tHeltimeTOF(0),
 165   tHeltrackLenghtTOF(0),
 166   tHelimpactXY(0),
 167   tHelimpactZ(0),
 168   tHelmapITS(0),
 169   tHelTPCNcls(0),
 170   tHelTRDsignal(0),
 171   tHelxPrimaryVertex(0),
 172   tHelyPrimaryVertex(0),
 173   tHelzPrimaryVertex(0),
 174   tHelchi2PerClusterTPC(0),
 175   fPIDResponse(0)
 176
 177 {
 178   // Dummy Constructor
 179
 180   //  printf("Dummy Constructor");
 181
 182   fESDtrackCuts = new AliESDtrackCuts("fESDtrackCuts");
 183   fESDtrackCuts->SetRequireITSStandAlone(kFALSE);
 184   fESDtrackCuts->SetRequireITSPureStandAlone(kFALSE);
 185
 186   fESDtrackCuts->SetRequireTPCRefit(kTRUE);
 187   fESDtrackCuts->SetAcceptKinkDaughters(kFALSE);
 188   fESDtrackCuts->SetMinNClustersTPC(60);
 189   fESDtrackCuts->SetMaxChi2PerClusterTPC(5);
 190   fESDtrackCuts->SetEtaRange(-0.9,0.9);
 191
 192 }
 193
 194 //________________________________________________________________________
 195 AliAnalysisTaskHelium3Pi::AliAnalysisTaskHelium3Pi(const char *name)
 196 : AliAnalysisTaskSE(name),
 197   fAnalysisType("ESD"),
 198   fCollidingSystems(0),
 199   fESDtrackCuts(0),
 200   fDataType("REAL"),
 201   fListHist(0),
 202   fHistEventMultiplicity(0),
 203   fHistTrackMultiplicity(0),
 204   fHistTrackMultiplicityCent(0),
 205   fHistTrackMultiplicitySemiCent(0),
 206   fHistTrackMultiplicityMB(0),
 207   fHistTrackMultiplicityPVCent(0),
 208   fHistTrackMultiplicityPVSemiCent(0),
 209   fHistTrackMultiplicityPVMB(0),
 210   fhBB(0),
 211   fhTOF(0),
 212   fhMassTOF(0),
 213   fhBBPions(0),
 214   fhBBHe(0),
 215   fhNaPos(0),
 216   fhNaNeg(0),
 217   fBetavsTPCsignalPos(0),
 218   fBetavsTPCsignalNeg(0),
 219   fNtuple1(0),
 220   trunNumber(0),
 221   tbunchcross(0),
 222   torbit(0),
 223   tperiod(0),
 224   teventtype(0),
 225   tTrackNumber(0),
 226   tpercentile(0),
 227   txPrimaryVertex(0),
 228   tyPrimaryVertex(0),
 229   tzPrimaryVertex(0),
 230   txSecondaryVertex(0),
 231   tySecondaryVertex(0),
 232   tzSecondaryVertex(0),
 233   tdcaTracks(0),
 234   tCosPointingAngle(0),
 235   tDCAV0toPrimaryVertex(0),
 236   tHeSign(0),
 237   tHepInTPC(0),
 238   tHeTPCsignal(0),
 239   tDcaHeToPrimVertex(0),
 240   tHeEta(0),
 241   tmomHex(0),
 242   tmomHey(0),
 243   tmomHez(0),
 244   tmomHeAtSVx(0),
 245   tmomHeAtSVy(0),
 246   tmomHeAtSVz(0),
 247   tHeTPCNcls(0),
 248   tHeimpactXY(0),
 249   tHeimpactZ(0),
 250   tHeITSClusterMap(0),
 251   tIsHeITSRefit(0),
 252   tPionSign(0),
 253   tPionpInTPC(0),
 254   tPionTPCsignal(0),
 255   tDcaPionToPrimVertex(0),
 256   tPionEta(0),
 257   tmomPionx(0),
 258   tmomPiony(0),
 259   tmomPionz(0),
 260   tmomNegPionAtSVx(0),
 261   tmomNegPionAtSVy(0),
 262   tmomNegPionAtSVz(0),
 263   tPionTPCNcls(0),
 264   tPionimpactXY(0),
 265   tPionimpactZ(0),
 266   tPionITSClusterMap(0),
 267   tIsPiITSRefit(0),
 268   txn(0),
 269   txp(0),
 270   tchi2He(0),
 271   tchi2Pi(0),
 272   fNtuple4(0),
 273   tHelrunNumber(0),
 274   tHelBCNumber(0),
 275   tHelOrbitNumber(0),
 276   tHelPeriodNumber(0),
 277   tHeleventtype(0),
 278   tHelisHeITSrefit(0),
 279   tHelpercentile(0),
 280   tHelSign(0),
 281   tHelpinTPC(0),
 282   tHelGetTPCsignal(0),
 283   tHelPx(0),
 284   tHelPy(0),
 285   tHelPz(0),
 286   tHelEta(0),
 287   tHelisTOF(0),
 288   tHelpoutTPC(0),
 289   tHeltimeTOF(0),
 290   tHeltrackLenghtTOF(0),
 291   tHelimpactXY(0),
 292   tHelimpactZ(0),
 293   tHelmapITS(0),
 294   tHelTPCNcls(0),
 295   tHelTRDsignal(0),
 296   tHelxPrimaryVertex(0),
 297   tHelyPrimaryVertex(0),
 298   tHelzPrimaryVertex(0),
 299   tHelchi2PerClusterTPC(0),
 300   fPIDResponse(0)
 301 {
 302
 303   // Define input and output slots here
 304   // Input slot #0 works with a TChain
 305   //DefineInput(0, TChain::Class());
 306   // Output slot #0 writes into a TList container ()
 307
 308   DefineInput(0, TChain::Class());
 309
 310   DefineOutput(1, TList::Class());
 311   DefineOutput(2, TTree::Class());
 312   DefineOutput(3, TTree::Class());
 313
 314   fESDtrackCuts = new AliESDtrackCuts("fESDtrackCuts");
 315   fESDtrackCuts->SetRequireITSStandAlone(kFALSE);
 316   fESDtrackCuts->SetRequireITSPureStandAlone(kFALSE);
 317
 318   fESDtrackCuts->SetRequireTPCRefit(kTRUE);
 319   fESDtrackCuts->SetAcceptKinkDaughters(kFALSE);
 320   fESDtrackCuts->SetMinNClustersTPC(60);
 321   fESDtrackCuts->SetMaxChi2PerClusterTPC(5);
 322   fESDtrackCuts->SetEtaRange(-0.9,0.9);
 323
 324 }
 325 //_______________________________________________________
 326 AliAnalysisTaskHelium3Pi::~AliAnalysisTaskHelium3Pi()
 327 {
 328   // Destructor
 329   if (fListHist) {
 330     delete fListHist;
 331     fListHist = 0;
 332   }
 333
 334   if (fESDtrackCuts) delete fESDtrackCuts;
 335   if(fNtuple1) delete fNtuple1;
 336   if(fNtuple4) delete fNtuple4;
 337 }
 338 //=================DEFINITION BETHE BLOCH==============================
 339
 340 Double_t AliAnalysisTaskHelium3Pi::BetheBloch(Double_t betaGamma,Double_t charge,Bool_t isPbPb) {
 341
 342   Double_t kp1, kp2, kp3, kp4, kp5;
 343
 344   if(isPbPb){
 345
 346     //    pass1 2011
 347     kp1 = 4.7*charge*charge;
 348     kp2 = 8.98482806165147636e+00;
 349     kp3 = 1.54000000000000005e-05;
 350     kp4 = 2.30445734159456084e+00;
 351     kp5 = 2.25624744086878559e+00;
 352
 353   }
 354
 355   else{
 356
 357     // to be defined ...
 358     //pass1 2011
 359     kp1 = 4.7*charge*charge;
 360     kp2 = 8.98482806165147636e+00;
 361     kp3 = 1.54000000000000005e-05;
 362     kp4 = 2.30445734159456084e+00;
 363     kp5 = 2.25624744086878559e+00;
 364
 365   }
 366
 367   Double_t beta = betaGamma / TMath::Sqrt(1.0 + betaGamma * betaGamma);
 368
 369   Double_t aa = TMath::Power(beta, kp4);
 370   Double_t bb = TMath::Power(1.0 / betaGamma, kp5);
 371
 372   bb = TMath::Log(kp3 + bb);
 373
 374   Double_t out = (kp2 - aa - bb) * kp1 / aa;
 375
 376   return out;
 377
 378 }
 379
 380 //==================DEFINITION OF OUTPUT OBJECTS==============================
 381
 382 void AliAnalysisTaskHelium3Pi::UserCreateOutputObjects()
 383 {
 384
 385   fListHist = new TList();
 386   fListHist->SetOwner();  // IMPORTANT!
 387
 388   if(! fHistEventMultiplicity ){
 389     fHistEventMultiplicity   = new TH1F( "fHistEventMultiplicity" , "Nb of Events" , 10 , 0, 10);
 390     fHistEventMultiplicity->GetXaxis()->SetBinLabel(1,"All Events");
 391     fHistEventMultiplicity->GetXaxis()->SetBinLabel(2,"Events w/PV");
 392     fHistEventMultiplicity->GetXaxis()->SetBinLabel(3,"Events w/|Vz|<10cm");
 393     fHistEventMultiplicity->GetXaxis()->SetBinLabel(4,"Central Events");
 394     fHistEventMultiplicity->GetXaxis()->SetBinLabel(5,"SemiCentral Events");
 395     fHistEventMultiplicity->GetXaxis()->SetBinLabel(6,"MB Events");
 396     fHistEventMultiplicity->GetXaxis()->SetBinLabel(7,"Central Events  w/|Vz|<10cm");
 397     fHistEventMultiplicity->GetXaxis()->SetBinLabel(8,"SemiCentral Events  w/|Vz|<10cm");
 398     fHistEventMultiplicity->GetXaxis()->SetBinLabel(9,"MB Events   w/|Vz|<10cm");
 399
 400     fListHist->Add(fHistEventMultiplicity);
 401   }
 402
 403   if(! fHistTrackMultiplicity ){
 404     fHistTrackMultiplicity   = new TH2F( "fHistTrackMultiplicity" , "Nb of Tracks", 2500,0, 25000,210,-1,104);
 405     fHistTrackMultiplicity->GetXaxis()->SetTitle("Number of tracks");
 406     fHistTrackMultiplicity->GetYaxis()->SetTitle("Percentile");
 407     fListHist->Add(fHistTrackMultiplicity);
 408   }
 409
 410   if(! fHistTrackMultiplicityCent ){
 411     fHistTrackMultiplicityCent   = new TH2F( "fHistTrackMultiplicityCent", "Nb of Tracks Central Events", 2500,0, 25000,210,-1,104 );
 412     fHistTrackMultiplicityCent->GetXaxis()->SetTitle("Number of tracks");
 413     fHistTrackMultiplicityCent->GetYaxis()->SetTitle("Percentile");
 414     fListHist->Add(fHistTrackMultiplicityCent);
 415   }
 416
 417   if(! fHistTrackMultiplicitySemiCent ){
 418     fHistTrackMultiplicitySemiCent   = new TH2F( "fHistTrackMultiplicitySemiCent" , "Nb of Tracks SemiCentral Events", 2500,0, 25000 ,210,-1,104);
 419     fHistTrackMultiplicitySemiCent->GetXaxis()->SetTitle("Number of tracks");
 420     fHistTrackMultiplicitySemiCent->GetYaxis()->SetTitle("Percentile");
 421     fListHist->Add(fHistTrackMultiplicitySemiCent);
 422   }
 423
 424   if(! fHistTrackMultiplicityMB ){
 425     fHistTrackMultiplicityMB   = new TH2F( "fHistTrackMultiplicityMB" , "Nb of Tracks MBral Events", 2500,0, 25000,210,-1,104 );
 426     fHistTrackMultiplicityMB->GetXaxis()->SetTitle("Number of tracks");
 427     fHistTrackMultiplicityMB->GetYaxis()->SetTitle("Percentile");
 428     fListHist->Add(fHistTrackMultiplicityMB);
 429   }
 430
 431   if(! fHistTrackMultiplicityPVCent ){
 432     fHistTrackMultiplicityPVCent   = new TH2F( "fHistTrackMultiplicityPVCent" , "Nb of Tracks Central Events", 2500,0, 25000,210,-1,104 );
 433     fHistTrackMultiplicityPVCent->GetXaxis()->SetTitle("Number of tracks");
 434     fHistTrackMultiplicityPVCent->GetYaxis()->SetTitle("Percentile");
 435     fListHist->Add(fHistTrackMultiplicityPVCent);
 436   }
 437
 438   if(! fHistTrackMultiplicityPVSemiCent ){
 439     fHistTrackMultiplicityPVSemiCent   = new TH2F( "fHistTrackMultiplicityPVSemiCent" , "Nb of Tracks SemiCentral Events", 2500,0, 25000 ,210,-1,104);
 440     fHistTrackMultiplicityPVSemiCent->GetXaxis()->SetTitle("Number of tracks");
 441     fHistTrackMultiplicityPVSemiCent->GetYaxis()->SetTitle("Percentile");
 442     fListHist->Add(fHistTrackMultiplicityPVSemiCent);
 443   }
 444
 445   if(! fHistTrackMultiplicityPVMB ){
 446     fHistTrackMultiplicityPVMB   = new TH2F( "fHistTrackMultiplicityPVMB" , "Nb of Tracks MBral Events", 2500,0, 25000,210,-1,104 );
 447     fHistTrackMultiplicityPVMB->GetXaxis()->SetTitle("Number of tracks");
 448     fHistTrackMultiplicityPVMB->GetYaxis()->SetTitle("Percentile");
 449     fListHist->Add(fHistTrackMultiplicityPVMB);
 450   }
 451
 452   if(! fhBB ){
 453     fhBB = new TH2F( "fhBB" , "BetheBlochTPC" , 120,-6,6,150,0,1500);
 454     fhBB->GetXaxis()->SetTitle("p/z (GeV/#it{c})");
 455     fhBB->GetYaxis()->SetTitle("TPC Signal");
 456     fListHist->Add(fhBB);
 457   }
 458
 459   if(! fhTOF ){
 460     fhTOF = new TH2F( "fhTOF" , "Scatter Plot TOF" , 120,-6,6,100,0,1.2);
 461     fhTOF->GetXaxis()->SetTitle("p/z (GeV/#it{c})");
 462     fhTOF->GetYaxis()->SetTitle("#beta");
 463     fListHist->Add(fhTOF);
 464   }
 465
 466   if(! fhMassTOF){
 467     fhMassTOF=new TH1F ("fhMassTOF","Particle Mass - TOF", 300,0 ,5);
 468     fhMassTOF->GetXaxis()->SetTitle("Mass (GeV/#it{c}^{2})");
 469     fListHist->Add(fhMassTOF);
 470   }
 471
 472   if(! fhBBPions ){
 473     fhBBPions = new TH2F( "fhBBPions" , "Bethe-Bloch TPC Pions" , 120,-6,6,150,0,1500);
 474     fhBBPions->GetXaxis()->SetTitle("p/z (GeV/#it{c})");
 475     fhBBPions->GetYaxis()->SetTitle("TPC Signal");
 476     fListHist->Add(fhBBPions);
 477   }
 478
 479   if(! fhBBHe ){
 480     fhBBHe = new TH2F( "fhBBHe" , "Bethe-Bloch TPC He" , 120,-6,6,150,0,1500);
 481     fhBBHe->GetXaxis()->SetTitle("p/z (GeV/#it{c})");
 482     fhBBHe->GetYaxis()->SetTitle("TPC Signal");
 483     fListHist->Add(fhBBHe);
 484   }
 485
 486   if(! fhNaPos ){
 487     fhNaPos = new TH2F( "fhNaPos" , "Distribution Pos" , 500,0,5,40,-10,10);
 488     fhNaPos->GetXaxis()->SetTitle("p/z (GeV/#it{c})");
 489     fhNaPos->GetYaxis()->SetTitle("(TPCSignal-bbtheo)/bbtheo (He)");
 490     fListHist->Add(fhNaPos);
 491   }
 492
 493   if(! fhNaNeg ){
 494     fhNaNeg = new TH2F( "fhNaNeg" , "Distribution Neg" , 500,0,5,40,-10,10);
 495     fhNaNeg->GetXaxis()->SetTitle("p/z (GeV/#it{c})");
 496     fhNaNeg->GetYaxis()->SetTitle("(TPCSignal-bbtheo)/bbtheo (He)");
 497     fListHist->Add(fhNaNeg);
 498   }
 499
 500   if(! fBetavsTPCsignalPos ){
 501     fBetavsTPCsignalPos = new TH2F("fBetavsTPCsignalPos","fBetavsTPCsignalPos",100,0,1.2,150,0,1500);
 502     fBetavsTPCsignalPos->GetXaxis()->SetTitle("#beta");
 503     fBetavsTPCsignalPos->GetYaxis()->SetTitle("TPC Signal");
 504     fListHist->Add(fBetavsTPCsignalPos);
 505   }
 506
 507   if(! fBetavsTPCsignalNeg ){
 508     fBetavsTPCsignalNeg = new TH2F("fBetavsTPCsignalNeg","fBetavsTPCsignalNeg",100,0,1.2,150,0,1500);
 509     fBetavsTPCsignalNeg->GetXaxis()->SetTitle("#beta");
 510     fBetavsTPCsignalNeg->GetYaxis()->SetTitle("TPC Signal");
 511     fListHist->Add(fBetavsTPCsignalNeg);
 512   }
 513
 514
 515
 516   if(! fNtuple1 ) {
 517
 518     fNtuple1 = new TTree("fNtuple1","fNtuple1");
 519
 520     fNtuple1->Branch("trunNumber"           ,&trunNumber           ,"trunNumber/F");
 521     fNtuple1->Branch("tbunchcross"          ,&tbunchcross          ,"tbunchcross/F");
 522     fNtuple1->Branch("torbit"               ,&torbit               ,"torbit/F");
 523     fNtuple1->Branch("tperiod"              ,&tperiod              ,"tperiod/F");
 524     fNtuple1->Branch("teventtype"           ,&teventtype           ,"teventtype/F");
 525     fNtuple1->Branch("tTrackNumber"         ,&tTrackNumber         ,"tTrackNumber/F");
 526     fNtuple1->Branch("tpercentile"          ,&tpercentile          ,"tpercentile/F") ;
 527     fNtuple1->Branch("txPrimaryVertex"      ,&txPrimaryVertex      ,"txPrimaryVertex/F");
 528     fNtuple1->Branch("tyPrimaryVertex"      ,&tyPrimaryVertex      ,"tyPrimaryVertex/F");
 529     fNtuple1->Branch("tzPrimaryVertex"      ,&tzPrimaryVertex      ,"tzPrimaryVertex/F");
 530     fNtuple1->Branch("txSecondaryVertex"    ,&txSecondaryVertex    ,"txSecondaryVertex/F");
 531     fNtuple1->Branch("tySecondaryVertex"    ,&tySecondaryVertex    ,"tySecondaryVertex/F");
 532     fNtuple1->Branch("tzSecondaryVertex"    ,&tzSecondaryVertex    ,"tzSecondaryVertex/F");
 533     fNtuple1->Branch("tdcaTracks"           ,&tdcaTracks           ,"tdcaTracks/F");
 534     fNtuple1->Branch("tCosPointingAngle"    ,&tCosPointingAngle    ,"tCosPointingAngle/F");
 535     fNtuple1->Branch("tDCAV0toPrimaryVertex",&tDCAV0toPrimaryVertex,"tDCAV0toPrimaryVertex/F");
 536     fNtuple1->Branch("tHeSign"              ,&tHeSign              ,"tHeSign/F");
 537     fNtuple1->Branch("tHepInTPC"            ,&tHepInTPC            ,"tHepInTPC/F");
 538     fNtuple1->Branch("tHeTPCsignal"         ,&tHeTPCsignal         ,"tHeTPCsignal/F");
 539     fNtuple1->Branch("tDcaHeToPrimVertex"   ,&tDcaHeToPrimVertex   ,"tDcaHeToPrimVertex/F");
 540     fNtuple1->Branch("tHeEta"               ,&tHeEta               ,"tHeEta/F");
 541     fNtuple1->Branch("tmomHex"              ,&tmomHex              ,"tmomHex/F");
 542     fNtuple1->Branch("tmomHey"              ,&tmomHey              ,"tmomHey/F");
 543     fNtuple1->Branch("tmomHez"              ,&tmomHez              ,"tmomHez/F");
 544     fNtuple1->Branch("tmomHeAtSVx"          ,&tmomHeAtSVx          ,"tmomHeAtSVx/F");
 545     fNtuple1->Branch("tmomHeAtSVy"          ,&tmomHeAtSVy          ,"tmomHeAtSVy/F");
 546     fNtuple1->Branch("tmomHeAtSVz"          ,&tmomHeAtSVz          ,"tmomHeAtSVz/F");
 547     fNtuple1->Branch("tHeTPCNcls"           ,&tHeTPCNcls           ,"tHeTPCNcls/F");
 548     fNtuple1->Branch("tHeimpactXY"          ,&tHeimpactXY          ,"tHeimpactXY/F");
 549     fNtuple1->Branch("tHeimpactZ"           ,&tHeimpactZ           ,"tHeimpactZ/F");
 550     fNtuple1->Branch("tHeITSClusterMap"     ,&tHeITSClusterMap     ,"tHeITSClusterMap/F");
 551     fNtuple1->Branch("tIsHeITSRefit"        ,&tIsHeITSRefit        ,"tIsHeITSRefit/F");
 552     fNtuple1->Branch("tPionSign"            ,&tPionSign            ,"tPionSign/F");
 553     fNtuple1->Branch("tPionpInTPC"          ,&tPionpInTPC          ,"tPionpInTPC/F");
 554     fNtuple1->Branch("tPionTPCsignal"       ,&tPionTPCsignal       ,"tPionTPCsignal/F");
 555     fNtuple1->Branch("tDcaPionToPrimVertex" ,&tDcaPionToPrimVertex ,"tDcaPionToPrimVertex/F");
 556     fNtuple1->Branch("tPionEta"             ,&tPionEta             ,"tPionEta/F");
 557     fNtuple1->Branch("tmomPionx"            ,&tmomPionx            ,"tmomPionx/F");
 558     fNtuple1->Branch("tmomPiony"            ,&tmomPiony            ,"tmomPiony/F");
 559     fNtuple1->Branch("tmomPionz"            ,&tmomPionz            ,"tmomPionz/F");
 560     fNtuple1->Branch("tmomNegPionAtSVx"     ,&tmomNegPionAtSVx     ,"tmomNegPionAtSVx/F");
 561     fNtuple1->Branch("tmomNegPionAtSVy"     ,&tmomNegPionAtSVy     ,"tmomNegPionAtSVy/F");
 562     fNtuple1->Branch("tmomNegPionAtSVz"     ,&tmomNegPionAtSVz     ,"tmomNegPionAtSVz/F");
 563     fNtuple1->Branch("tPionTPCNcls"         ,&tPionTPCNcls         ,"tPionTPCNcls/F");
 564     fNtuple1->Branch("tPionimpactXY"        ,&tPionimpactXY        ,"tPionimpactXY/F");
 565     fNtuple1->Branch("tPionimpactZ"         ,&tPionimpactZ         ,"tPionimpactZ/F");
 566     fNtuple1->Branch("tPionITSClusterMap"   ,&tPionITSClusterMap   ,"tPionITSClusterMap/F");
 567     fNtuple1->Branch("tIsPiITSRefit"        ,&tIsPiITSRefit        ,"tIsPiITSRefit/F");
 568     fNtuple1->Branch("txn"                  ,&txn                  ,"txn/F");
 569     fNtuple1->Branch("txp"                  ,&txp                  ,"txp/F");
 570     fNtuple1->Branch("tchi2He"              ,&tchi2He              ,"tchi2He/F");
 571     fNtuple1->Branch("tchi2Pi"              ,&tchi2Pi              ,"tchi2Pi/F");
 572
 573   }
 574
 575   if(! fNtuple4 ) {
 576
 577     fNtuple4 = new TTree("fNtuple4","fNtuple4");
 578
 579     fNtuple4->Branch("tHelrunNumber"        ,&tHelrunNumber        ,"tHelrunNumber/F");
 580     fNtuple4->Branch("tHelBCNumber"         ,&tHelBCNumber         ,"tHelBCNumber/F");
 581     fNtuple4->Branch("tHelOrbitNumber"      ,&tHelOrbitNumber      ,"tHelOrbitNumber/F");
 582     fNtuple4->Branch("tHelPeriodNumber"     ,&tHelPeriodNumber     ,"tHelPeriodNumber/F");
 583     fNtuple4->Branch("tHeleventtype"        ,&tHeleventtype        ,"tHeleventtype/F");
 584     fNtuple4->Branch("tHelisHeITSrefit"     ,&tHelisHeITSrefit     ,"tHelisHeITSrefit/F");
 585     fNtuple4->Branch("tHelpercentile"       ,&tHelpercentile       ,"tHelpercentile/F");
 586     fNtuple4->Branch("tHelSign"             ,&tHelSign             ,"tHelSign/F");
 587     fNtuple4->Branch("tHelpinTPC"           ,&tHelpinTPC           ,"tHelpinTPC/F");
 588     fNtuple4->Branch("tHelGetTPCsignal"     ,&tHelGetTPCsignal     ,"tHelGetTPCsignal/F");
 589     fNtuple4->Branch("tHelPx"               ,&tHelPx               ,"tHelPx/F");
 590     fNtuple4->Branch("tHelPy"               ,&tHelPy               ,"tHelPy/F");
 591     fNtuple4->Branch("tHelPz"               ,&tHelPz               ,"tHelPz/F");
 592     fNtuple4->Branch("tHelEta"              ,&tHelEta              ,"tHelEta/F");
 593     fNtuple4->Branch("tHelisTOF"            ,&tHelisTOF            ,"tHelisTOF/F");
 594     fNtuple4->Branch("tHelpoutTPC"          ,&tHelpoutTPC          ,"tHelpoutTPC/F");
 595     fNtuple4->Branch("tHeltimeTOF"          ,&tHeltimeTOF          ,"tHeltimeTOF/F");
 596     fNtuple4->Branch("tHeltrackLenghtTOF"   ,&tHeltrackLenghtTOF   ,"tHeltrackLenghtTOF/F");
 597     fNtuple4->Branch("tHelimpactXY"         ,&tHelimpactXY         ,"tHelimpactXY/F");
 598     fNtuple4->Branch("tHelimpactZ"          ,&tHelimpactZ          ,"tHelimpactZ/F");
 599     fNtuple4->Branch("tHelmapITS"           ,&tHelmapITS           ,"tHelmapITS/F");
 600     fNtuple4->Branch("tHelTPCNcls"          ,&tHelTPCNcls          ,"tHelTPCNcls/F");
 601     fNtuple4->Branch("tHelTRDsignal"        ,&tHelTRDsignal        ,"tHelTRDsignal/F");
 602     fNtuple4->Branch("tHelxPrimaryVertex"   ,&tHelxPrimaryVertex   ,"tHelxPrimaryVertex/F");
 603     fNtuple4->Branch("tHelyPrimaryVertex"   ,&tHelyPrimaryVertex   ,"tHelyPrimaryVertex/F");
 604     fNtuple4->Branch("tHelzPrimaryVertex"   ,&tHelzPrimaryVertex   ,"tHelzPrimaryVertex/F");
 605     fNtuple4->Branch("tHelchi2PerClusterTPC",&tHelchi2PerClusterTPC,"tHelchi2PerClusterTPC/F");
 606
 607
 608   }
 609
 610   PostData(1,  fListHist);
 611   PostData(2,  fNtuple1);
 612   PostData(3,  fNtuple4);
 613 }// end UserCreateOutputObjects
 614
 615
 616 //====================== USER EXEC ========================
 617
 618 void AliAnalysisTaskHelium3Pi::UserExec(Option_t *)
 619 {
 620   //_______________________________________________________________________
 621
 622   //!*********************!//
 623   //!  Define variables   !//
 624   //!*********************!//
 625
 626   Double_t pinTPC=0.,poutTPC=0.,TPCSignal=0.;
 627   Double_t xPrimaryVertex=0.,yPrimaryVertex=0.,zPrimaryVertex=0.;
 628   Double_t massTOF=0.,timeTOF=0.,trackLenghtTOF=0.,betaTOF=0.;
 629
 630   ULong_t  status=0;
 631   //  ULong_t  statusT=0;
 632   ULong_t  statusPi=0;
 633
 634   Bool_t   isTPC=kFALSE,isTOF=kFALSE,IsHeITSRefit=kFALSE,IsPiITSRefit=kFALSE ;
 635
 636   Float_t nSigmaNegPion=0.;
 637
 638   Double_t cutNSigma = 3;
 639   Double_t bbtheoM=0.,bbtheo=0.;
 640   Double_t zNathashaNeg=0;
 641   Double_t zNathashaPos=0;
 642   Double_t fPos[3]={0.,0.,0.};
 643   Double_t runNumber=0.;
 644   //  Double_t evNumber=0.;
 645
 646   Double_t BCNumber=0.;
 647   Double_t OrbitNumber=0.;
 648   Double_t PeriodNumber=0.;
 649
 650   Double_t        Helium3Mass = 2.80839;
 651   Double_t        PionMass    = 0.13957;
 652   // TLORENTZ vectors
 653
 654   TLorentzVector  vPion,vHelium,vSum;
 655
 656   //!----------------------------------------------------------------
 657
 658   //! A set of very loose parameters for cuts
 659
 660   Double_t fgChi2max=33.;     //! max chi2
 661   Double_t fgDNmin=0.05;      //! min imp parameter for the 1st daughter = 500um
 662   Double_t fgDCAmax=1.0;      //! max DCA between the daughter tracks in cm
 663   Double_t fgCPAmin=0.99;     //! min cosine of V0's pointing angle
 664   //  Double_t fgRmin=0.2;    //! min radius of the fiducial volume //original
 665   Double_t fgRmin=0.1;        //! min radius of the fiducial volume = 1 mm
 666   Double_t fgRmax=200.;       //! max radius of the fiducial volume = 2 m
 667
 668   //------------------------------------------
 669
 670   // Main loop
 671   // Called for EACH event
 672
 673   AliVEvent *event = InputEvent();
 674   if (!event) { Printf("ERROR: Could not retrieve event"); return; }
 675
 676   Info("AliAnalysisTaskHelium3Pi","Starting UserExec");
 677
 678   SetDataType("REAL");
 679
 680   // create pointer to event
 681   AliESDEvent* lESDevent = dynamic_cast<AliESDEvent*>(event);
 682   if (!lESDevent) {
 683     AliError("Cannot get the ESD event");
 684     return;
 685   }
 686
 687   fHistEventMultiplicity->Fill(0);
 688
 689   Double_t lMagneticField=lESDevent->GetMagneticField();
 690   Int_t TrackNumber = -1;
 691
 692
 693   //*****************//
 694   //*   Centrality  *//
 695   //*****************//
 696
 697   AliCentrality *centrality = lESDevent->GetCentrality();
 698   Float_t percentile=centrality->GetCentralityPercentile("V0M");
 699
 700   TrackNumber = lESDevent->GetNumberOfTracks();
 701   if (TrackNumber<2) return;
 702
 703   fHistTrackMultiplicity->Fill(TrackNumber,percentile); //tracce per evento
 704
 705   //****************************************
 706
 707   // PID
 708
 709   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 710   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 711   fPIDResponse=inputHandler->GetPIDResponse(); // data member di tipo "const AliPIDResponse *fPIDResponse;"
 712   //  cout<<"fPIDResponse "<<fPIDResponse<<endl;
 713   //===========================================
 714
 715   Int_t eventtype=-99;
 716
 717   Bool_t isSelectedCentral     = (inputHandler->IsEventSelected() & AliVEvent::kCentral);
 718   Bool_t isSelectedSemiCentral = (inputHandler->IsEventSelected() & AliVEvent::kSemiCentral);
 719   Bool_t isSelectedMB          = (inputHandler->IsEventSelected() & AliVEvent::kMB);
 720
 721   if(isSelectedCentral){
 722     fHistEventMultiplicity->Fill(3);
 723     fHistTrackMultiplicityCent->Fill(TrackNumber,percentile);
 724     eventtype=1;
 725   }
 726
 727   if(isSelectedSemiCentral){
 728     fHistEventMultiplicity->Fill(4);
 729     fHistTrackMultiplicitySemiCent->Fill(TrackNumber,percentile);
 730     eventtype=2;
 731   }
 732
 733   if(isSelectedMB){
 734     fHistEventMultiplicity->Fill(5);
 735     fHistTrackMultiplicityMB->Fill(TrackNumber,percentile);
 736     eventtype=3;
 737   }
 738
 739   if(isSelectedCentral || isSelectedSemiCentral || isSelectedMB){
 740
 741     // ANALISYS
 742
 743     // Primary vertex cut
 744
 745     const AliESDVertex *vtx = lESDevent->GetPrimaryVertexTracks();
 746
 747     if(vtx->GetNContributors()<1) {
 748
 749       // SPD vertex cut
 750       vtx = lESDevent->GetPrimaryVertexSPD();
 751
 752       if(vtx->GetNContributors()<1) {
 753         Info("AliAnalysisTaskHelium3Pi","No good vertex, skip event");
 754         return; // NO GOOD VERTEX, SKIP EVENT
 755       }
 756     }
 757
 758     fHistEventMultiplicity->Fill(1); // analyzed events with PV
 759
 760     xPrimaryVertex=vtx->GetXv();
 761     yPrimaryVertex=vtx->GetYv();
 762     zPrimaryVertex=vtx->GetZv();
 763
 764     if(TMath::Abs(zPrimaryVertex)>10) return;
 765
 766     if(eventtype==1){
 767       fHistTrackMultiplicityPVCent->Fill(TrackNumber,percentile);
 768       fHistEventMultiplicity->Fill(6);
 769     }
 770
 771     if(eventtype==2){
 772       fHistTrackMultiplicityPVSemiCent->Fill(TrackNumber,percentile);
 773       fHistEventMultiplicity->Fill(7);
 774     }
 775
 776     if(eventtype==3){
 777       fHistTrackMultiplicityPVMB->Fill(TrackNumber,percentile);
 778       fHistEventMultiplicity->Fill(8);
 779     }
 780
 781
 782     fHistEventMultiplicity->Fill(2);
 783
 784     //Find Pair candidates
 785
 786     TArrayI PionsTPC(TrackNumber);        //Neg pions
 787     Int_t nPionsTPC=0;
 788
 789     TArrayI HeTPC(TrackNumber);        //helium3
 790     Int_t nHeTPC=0;
 791
 792     const Double_t speedOfLight =  TMath::C()*1E2*1E-12; // cm/ps
 793
 794     Float_t impactXY=-999, impactZ=-999;
 795     Float_t impactXYpi=-999, impactZpi=-999;
 796
 797
 798     //*************************************************************
 799
 800     runNumber = lESDevent->GetRunNumber();
 801     BCNumber    = lESDevent->GetBunchCrossNumber();
 802     OrbitNumber = lESDevent->GetOrbitNumber();
 803     PeriodNumber= lESDevent->GetPeriodNumber();
 804
 805     //*************************************************************
 806
 807     for (Int_t j=0; j<TrackNumber; j++) { //loop on tracks
 808
 809       AliESDtrack *esdtrack=lESDevent->GetTrack(j);
 810       //      AliVTrack*  esdtrack= (AliVTrack *) fEvent->GetTrack(iT);
 811
 812
 813       if(!esdtrack) {
 814         AliError(Form("ERROR: Could not retrieve esdtrack %d",j));
 815         continue;
 816       }
 817
 818       // ************** Track cuts ****************
 819
 820       if (!fESDtrackCuts->AcceptTrack(esdtrack)) continue;
 821
 822
 823       status  = (ULong_t)esdtrack->GetStatus();
 824       isTPC   = (((status) & AliESDtrack::kTPCin)  != 0);
 825       isTOF   = ((((status) & AliESDtrack::kTOFout) != 0) && (((status) & AliESDtrack::kTIME) != 0));
 826
 827
 828       UInt_t mapITS=esdtrack->GetITSClusterMap();
 829
 830       //----------------------------------------------
 831
 832       //****** Cuts from  AliV0Vertex.cxx *************
 833
 834       Double_t d=esdtrack->GetD(xPrimaryVertex,yPrimaryVertex,lMagneticField);
 835       //    if (TMath::Abs(d)<fgDPmin) continue;
 836       if (TMath::Abs(d)>fgRmax) continue;
 837
 838       //---- (Usefull) Stuff
 839
 840       TPCSignal=esdtrack->GetTPCsignal();
 841
 842       if (TPCSignal<10)continue;
 843       if (TPCSignal>1000)continue;
 844
 845       if(!isTPC)continue;
 846       if(!esdtrack->GetTPCInnerParam())continue;
 847
 848       AliExternalTrackParam trackIn(*esdtrack->GetInnerParam());
 849       pinTPC= trackIn.GetP();
 850
 851       //pinTPC= esdtrack->GetTPCMomentum();
 852
 853       poutTPC=pinTPC;
 854
 855
 856       if((status) & (AliESDtrack::kITSrefit!=0)){
 857         fhBB->Fill(pinTPC*esdtrack->GetSign(),TPCSignal);
 858       }
 859
 860       timeTOF=esdtrack->GetTOFsignal();                 // ps
 861       trackLenghtTOF= esdtrack->GetIntegratedLength();  // cm
 862
 863       if(isTOF){
 864
 865         if(!esdtrack->GetOuterParam())continue;
 866
 867         AliExternalTrackParam trackOut(*esdtrack->GetOuterParam());
 868
 869         poutTPC = trackOut.GetP();
 870
 871         betaTOF= (trackLenghtTOF/timeTOF)/2.99792458e-2;
 872
 873         fhTOF->Fill(poutTPC*esdtrack->GetSign(),betaTOF);
 874
 875         Double_t mass2=(poutTPC*poutTPC)*((((speedOfLight*speedOfLight)*(timeTOF*timeTOF))-(trackLenghtTOF*trackLenghtTOF))/(trackLenghtTOF*trackLenghtTOF));
 876         if(mass2>0) massTOF=TMath::Sqrt(mass2);
 877         fhMassTOF->Fill(massTOF);
 878
 879         if(esdtrack->GetSign() < 0.)fBetavsTPCsignalNeg->Fill(betaTOF,TPCSignal);
 880         if(esdtrack->GetSign() > 0.)fBetavsTPCsignalPos->Fill(betaTOF,TPCSignal);
 881
 882       }
 883
 884       //pass2
 885
 886       // bbtheo =BetheBloch((2*pinTPC)/3.,2,kTRUE);    //! OK
 887       // bbtheoM=(1 - 0.08*5)*bbtheo;                  //! OK
 888       // bbtheoP=(1 + 0.08*5)*bbtheo;                  //! OK
 889
 890
 891       bbtheo = fPIDResponse->NumberOfSigmas((AliPIDResponse::EDetector)0,esdtrack,(AliPID::EParticleType) 7);
 892
 893       if(esdtrack->GetSign()<0){
 894         zNathashaNeg=bbtheo;//(TPCSignal-bbtheo)/bbtheo;
 895         //      cout<<"BBtheo 1 :"<<zNathashaNeg<<endl;
 896         fhNaNeg->Fill(pinTPC,zNathashaNeg);
 897       }
 898
 899       if(esdtrack->GetSign() > 0.){
 900         zNathashaPos=bbtheo;//(TPCSignal-bbtheo)/bbtheo;
 901         fhNaPos->Fill(pinTPC,zNathashaPos);
 902       }
 903
 904       nSigmaNegPion=TMath::Abs(fPIDResponse->NumberOfSigmasTPC(esdtrack,(AliPID::EParticleType) 2));
 905       //2 is pion
 906
 907       if ( (nSigmaNegPion < cutNSigma)){
 908
 909         //      cout<<"Nsigma pi: "<<nSigmaNegPion<<endl;
 910
 911         fhBBPions->Fill(pinTPC*esdtrack->GetSign(),TPCSignal);
 912
 913         if(pinTPC<3.){
 914           PionsTPC[nPionsTPC++]=j;
 915         }
 916       }
 917
 918       //      nSigmaNegPion=(fPIDResponse->NumberOfSigmasTPC(esdtrack,(AliPID::EParticleType) 2));
 919
 920       bbtheoM = TMath::Abs((fPIDResponse->NumberOfSigmasTPC(esdtrack,(AliPID::EParticleType) 7)));
 921
 922       //      if( TPCSignal > bbtheoM ) {
 923       //      if( bbtheoM > -3.) {
 924       if( bbtheoM < 3.) {
 925
 926         if(pinTPC>0.6){
 927
 928           fhBBHe->Fill(pinTPC*esdtrack->GetSign(),TPCSignal);
 929           HeTPC[nHeTPC++]=j;
 930
 931           Bool_t isHeITSrefit=((status) & (AliESDtrack::kITSrefit));
 932
 933           esdtrack->GetImpactParameters(impactXY, impactZ);
 934
 935           Int_t  fIdxInt[200]; //dummy array
 936           Int_t nClustersTPC = esdtrack->GetTPCclusters(fIdxInt);
 937
 938           Float_t chi2PerClusterTPC = esdtrack->GetTPCchi2()/(Float_t)(nClustersTPC);
 939
 940           tHelrunNumber         =(Float_t)runNumber;
 941           tHelBCNumber          =(Float_t)BCNumber;
 942           tHelOrbitNumber       =(Float_t)OrbitNumber;
 943           tHelPeriodNumber      =(Float_t)PeriodNumber;
 944           tHeleventtype         =(Float_t)eventtype;
 945           tHelisHeITSrefit      =(Float_t)isHeITSrefit;
 946           tHelpercentile        =(Float_t)percentile;
 947           tHelSign              =(Float_t)esdtrack->GetSign();
 948           tHelpinTPC            =(Float_t)pinTPC;
 949           tHelGetTPCsignal      =(Float_t)esdtrack->GetTPCsignal();
 950           tHelPx                =(Float_t)esdtrack->Px();
 951           tHelPy                =(Float_t)esdtrack->Py();
 952           tHelPz                =(Float_t)esdtrack->Pz();
 953           tHelEta               =(Float_t)esdtrack->Eta();
 954           tHelisTOF             =(Float_t)isTOF;
 955           tHelpoutTPC           =(Float_t)poutTPC;
 956           tHeltimeTOF           =(Float_t)timeTOF;
 957           tHeltrackLenghtTOF    =(Float_t)trackLenghtTOF;
 958           tHelimpactXY          =(Float_t)impactXY;
 959           tHelimpactZ           =(Float_t)impactZ;
 960           tHelmapITS            =(Float_t)mapITS;
 961           tHelTPCNcls           =(Float_t)esdtrack->GetTPCNcls();
 962           tHelTRDsignal         =(Float_t)esdtrack->GetTRDsignal();
 963           tHelxPrimaryVertex    =(Float_t)xPrimaryVertex;
 964           tHelyPrimaryVertex    =(Float_t)yPrimaryVertex;
 965           tHelzPrimaryVertex    =(Float_t)zPrimaryVertex;
 966           tHelchi2PerClusterTPC =(Float_t)chi2PerClusterTPC;
 967
 968           fNtuple4->Fill();
 969         }
 970       }
 971     }  //! track
 972
 973     PionsTPC.Set(nPionsTPC);
 974     HeTPC.Set(nHeTPC);
 975
 976     Double_t        DcaHeToPrimVertex=0;
 977     Double_t        DcaPionToPrimVertex=0;
 978
 979     impactXY=-999, impactZ=-999;
 980     impactXYpi=-999, impactZpi=-999;
 981
 982     // Track
 983
 984     // Vettors for il PxPyPz
 985
 986     Double_t momPionVett[3];
 987     for(Int_t i=0;i<3;i++)momPionVett[i]=0;
 988
 989     Double_t momHeVett[3];
 990     for(Int_t i=0;i<3;i++)momHeVett[i]=0;
 991
 992     //At SV
 993
 994     Double_t momPionVettAt[3];
 995     for(Int_t i=0;i<3;i++)momPionVettAt[i]=0;
 996
 997     Double_t momHeVettAt[3];
 998     for(Int_t i=0;i<3;i++)momHeVettAt[i]=0;
 999
1000     //---------------   LOOP PAIRS   ----------------
1001
1002     for (Int_t k=0; k < nPionsTPC; k++) {                           //! Pions Loop
1003
1004       DcaPionToPrimVertex=0.;
1005       DcaHeToPrimVertex=0;
1006
1007       Int_t PionIdx=PionsTPC[k];
1008
1009       AliESDtrack  *PionTrack=lESDevent->GetTrack(PionIdx);
1010
1011       statusPi = (ULong_t)PionTrack->GetStatus();
1012       //      isTOFPi  = ((((statusPi) & (AliESDtrack::kTOFout)) != 0) && (((statusPi) & (AliESDtrack::kTIME)) != 0));
1013       IsPiITSRefit = ((statusPi) & (AliESDtrack::kITSrefit));
1014
1015       if (PionTrack)
1016         DcaPionToPrimVertex = TMath::Abs(PionTrack->GetD(xPrimaryVertex, yPrimaryVertex,lMagneticField)); //OK
1017
1018       if(DcaPionToPrimVertex<0.2)continue;
1019
1020       AliExternalTrackParam trackInPion(*PionTrack);
1021
1022       for (Int_t i=0; i<nHeTPC; i++){                               //! Helium Loop
1023
1024         Int_t HeIdx=HeTPC[i];
1025
1026          AliESDtrack  *HeTrack=lESDevent->GetTrack(HeIdx);
1027
1028         //      statusT= (ULong_t)HeTrack->GetStatus();
1029         //      isTOFHe   = (((statusT & AliESDtrack::kTOFout) != 0) && ((statusT & AliESDtrack::kTIME) != 0));
1030         IsHeITSRefit = (status & AliESDtrack::kITSrefit);
1031
1032         if (HeTrack)
1033           DcaHeToPrimVertex = TMath::Abs(HeTrack->GetD(xPrimaryVertex, yPrimaryVertex,lMagneticField)); //OK
1034
1035         AliExternalTrackParam trackInHe(*HeTrack);
1036
1037         if ( DcaPionToPrimVertex < fgDNmin)                //OK
1038           if ( DcaHeToPrimVertex < fgDNmin) continue;    //OK
1039
1040         Double_t xn, xp;
1041         Double_t dca=0.;
1042
1043         dca= PionTrack->GetDCA(HeTrack,lMagneticField,xn,xp); //!dca (Neg to Pos)
1044
1045         if (dca > fgDCAmax) continue;
1046         if ((xn+xp) > 2*fgRmax) continue;
1047         if ((xn+xp) < 2*fgRmin) continue;
1048
1049         //CORRECTION from AliV0Vertex
1050
1051         Bool_t corrected=kFALSE;
1052         if ((trackInPion.GetX() > 3.) && (xn < 3.)) {
1053           //correct for the beam pipe material
1054           corrected=kTRUE;
1055         }
1056         if ((trackInHe.GetX() > 3.) && (xp < 3.)) {
1057           //correct for the beam pipe material
1058           corrected=kTRUE;
1059         }
1060         if (corrected) {
1061           dca=trackInPion.GetDCA(&trackInHe,lMagneticField,xn,xp);
1062           if (dca > fgDCAmax) continue;
1063           if ((xn+xp) > 2*fgRmax) continue;
1064           if ((xn+xp) < 2*fgRmin) continue;
1065         }
1066
1067         //=============================================//
1068         // Make "V0" with found tracks                 //
1069         //=============================================//
1070
1071         trackInPion.PropagateTo(xn,lMagneticField);
1072         trackInHe.PropagateTo(xp,lMagneticField);
1073
1074         AliESDv0 vertex(trackInPion,PionIdx,trackInHe,HeIdx);
1075         if (vertex.GetChi2V0() > fgChi2max) continue;
1076
1077         Float_t CosPointingAngle=vertex.GetV0CosineOfPointingAngle(xPrimaryVertex,yPrimaryVertex,zPrimaryVertex); //PointingAngle
1078         if (CosPointingAngle < fgCPAmin) continue;
1079
1080         vertex.SetDcaV0Daughters(dca);
1081         vertex.SetV0CosineOfPointingAngle(CosPointingAngle);
1082
1083         fPos[0]=vertex.Xv();
1084         fPos[1]=vertex.Yv();
1085         fPos[2]=vertex.Zv();
1086
1087         HeTrack->PxPyPz(momHeVett);
1088         PionTrack->PxPyPz(momPionVett);
1089
1090         Double_t raggio=TMath::Sqrt(fPos[0]*fPos[0]+fPos[1]*fPos[1]+fPos[2]*fPos[2]);
1091         HeTrack->GetPxPyPzAt(raggio,lMagneticField,momHeVettAt);
1092         PionTrack->GetPxPyPzAt(raggio,lMagneticField,momPionVettAt);
1093
1094         //------------------------------------------------------------------------//
1095
1096         HeTrack->GetImpactParameters(impactXY, impactZ);
1097
1098         PionTrack->GetImpactParameters(impactXYpi, impactZpi);
1099
1100         if(vertex.GetD(xPrimaryVertex,yPrimaryVertex,zPrimaryVertex)>3) continue;
1101
1102         //salvo solo fino a 3.1 GeV/c2
1103
1104         vHelium.SetXYZM(2*momHeVettAt[0],2*momHeVettAt[1],2*momHeVettAt[2],Helium3Mass);
1105         vPion.SetXYZM(momPionVettAt[0],momPionVettAt[1],momPionVettAt[2],PionMass);
1106         vSum=vHelium+vPion;
1107
1108         if(vSum.M()>3.2)
1109           continue;
1110
1111         Int_t  fIdxInt[200]; //dummy array
1112         Int_t nClustersTPCHe = HeTrack->GetTPCclusters(fIdxInt);
1113         Int_t nClustersTPCPi = PionTrack->GetTPCclusters(fIdxInt);
1114
1115         //----------------------------------------------------------------------//
1116
1117         trunNumber              =(Float_t)runNumber;
1118         tbunchcross             =(Float_t)BCNumber;
1119         torbit                  =(Float_t)OrbitNumber;
1120         tperiod                 =(Float_t)PeriodNumber;
1121         teventtype              =(Float_t)eventtype;
1122         tTrackNumber            =(Float_t)TrackNumber;
1123         tpercentile             =(Float_t)percentile;
1124         txPrimaryVertex         =(Float_t)xPrimaryVertex; //PRIMARY
1125         tyPrimaryVertex         =(Float_t)yPrimaryVertex;
1126         tzPrimaryVertex         =(Float_t)zPrimaryVertex;
1127         txSecondaryVertex       =(Float_t)fPos[0]; //SECONDARY
1128         tySecondaryVertex       =(Float_t)fPos[1];
1129         tzSecondaryVertex       =(Float_t)fPos[2];
1130         tdcaTracks              =(Float_t)dca;           //between 2 tracks
1131         tCosPointingAngle       =(Float_t)CosPointingAngle;          //cosPointingAngle da V0
1132         tDCAV0toPrimaryVertex   =(Float_t)vertex.GetD(xPrimaryVertex,yPrimaryVertex,zPrimaryVertex);
1133         tHeSign                 =(Float_t)HeTrack->GetSign(); //He
1134         tHepInTPC               =(Float_t)trackInHe.GetP();
1135         tHeTPCsignal            =(Float_t)HeTrack->GetTPCsignal();
1136         tDcaHeToPrimVertex      =(Float_t)DcaHeToPrimVertex;
1137         tHeEta                  =(Float_t)HeTrack->Eta();
1138         tmomHex                 =(Float_t)momHeVett[0];
1139         tmomHey                 =(Float_t)momHeVett[1];
1140         tmomHez                 =(Float_t)momHeVett[2];
1141         tmomHeAtSVx             =(Float_t)momHeVettAt[0];
1142         tmomHeAtSVy             =(Float_t)momHeVettAt[1];
1143         tmomHeAtSVz             =(Float_t)momHeVettAt[2];
1144         tHeTPCNcls              =(Float_t)HeTrack->GetTPCNcls();
1145         tHeimpactXY             =(Float_t)impactXY;
1146         tHeimpactZ              =(Float_t)impactZ;
1147         tHeITSClusterMap        =(Float_t)HeTrack->GetITSClusterMap();
1148         tIsHeITSRefit           =(Float_t)IsHeITSRefit;
1149         tPionSign               =(Float_t)PionTrack->GetSign(); //Pion
1150         tPionpInTPC             =(Float_t)trackInPion.GetP();
1151         tPionTPCsignal          =(Float_t)PionTrack->GetTPCsignal();
1152         tDcaPionToPrimVertex    =(Float_t)DcaPionToPrimVertex;
1153         tPionEta                =(Float_t)PionTrack->Eta();
1154         tmomPionx               =(Float_t)momPionVett[0];
1155         tmomPiony               =(Float_t)momPionVett[1];
1156         tmomPionz               =(Float_t)momPionVett[2];
1157         tmomNegPionAtSVx        =(Float_t)momPionVettAt[0];
1158         tmomNegPionAtSVy        =(Float_t)momPionVettAt[1];
1159         tmomNegPionAtSVz        =(Float_t)momPionVettAt[2];
1160         tPionTPCNcls            =(Float_t)PionTrack->GetTPCNcls();
1161         tPionimpactXY           =(Float_t)impactXYpi;
1162         tPionimpactZ            =(Float_t)impactZpi;
1163         tPionITSClusterMap      =(Float_t)PionTrack->GetITSClusterMap();
1164         tIsPiITSRefit           =(Float_t)IsPiITSRefit;
1165         txn                     =(Float_t)xn;
1166         txp                     =(Float_t)xp;
1167         tchi2He                 =(Float_t)HeTrack->GetTPCchi2()/(Float_t)(nClustersTPCHe);
1168         tchi2Pi                 =(Float_t)PionTrack->GetTPCchi2()/(Float_t)(nClustersTPCPi);
1169
1170
1171         fNtuple1->Fill();
1172         vertex.Delete();
1173       }// positive TPC
1174
1175     } //negative tpc
1176
1177   }
1178
1179   PostData(1,fListHist);
1180   PostData(2,fNtuple1);
1181   PostData(3,fNtuple4);
1182
1183 } //end userexec
1184
1185
1186 //________________________________________________________________________
1187
1188 void AliAnalysisTaskHelium3Pi::Terminate(Option_t *)
1189 {
1190   // Draw result to the screen
1191   // Called once at the end of the query
1192 }
1193
1194