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