PWG1/global/AliAnalysisTaskCTau.cxx

   1 #include <TCanvas.h>
   2 #include <TTree.h>
   3 #include <TFile.h>
   4 #include <TH1F.h>
   5 #include <TPDGCode.h>
   6 #include <TDatabasePDG.h>
   7 #include <TParticlePDG.h>
   8 #include <TParticle.h>
   9 #include <TROOT.h>
  10
  11 #include "AliESDEvent.h"
  12 #include "AliESDv0.h"
  13 #include "AliESDcascade.h"
  14
  15 #include "AliMCEvent.h"
  16 #include "AliStack.h"
  17
  18 #include "AliAnalysisTaskCTau.h"
  19
  20 extern TROOT *gROOT;
  21
  22 //
  23 //  This is a little task for checking the c*tau of the strange particles
  24 //
  25
  26 AliAnalysisTaskCTau::AliAnalysisTaskCTau(const char *name) :
  27 AliAnalysisTaskSE(name),
  28 fOutput(0),
  29 fESD(0),
  30 fK0s(0),
  31 fK0sMC(0),
  32 fLambdas(0),
  33 fLambdasMC(0),
  34 fLambdaBars(0),
  35 fLambdaBarsMC(0),
  36 fXis(0),
  37 fXisMC(0),
  38 fMass(0),
  39 fMassMC(0)
  40 {
  41   // Constructor. Initialization of pointers
  42   DefineOutput(1, TList::Class());
  43 }
  44
  45 void AliAnalysisTaskCTau::UserCreateOutputObjects()
  46 {
  47   fOutput = new TList();
  48   fOutput->SetOwner();
  49
  50   fK0s     = new TH1F("K0s", "c\\tau for K^{0}_{s}", 100, 0., 25.);//10*2.684);
  51   fK0s->GetXaxis()->SetTitle("L\\frac{m}{p} [cm]");
  52   fOutput->Add(fK0s);
  53   fK0sMC   = new TH1F("K0sMC", "c\\tau for K^{0}_{s}, MC", 100, 0., 25.);
  54   fK0sMC->GetXaxis()->SetTitle("L\\frac{m}{p} [cm]");
  55   fOutput->Add(fK0sMC);
  56
  57   fLambdas = new TH1F("Lambdas","c\\tau for \\Lambda", 100, 0.,80.);//10*7.89);
  58   fLambdas->GetXaxis()->SetTitle("L\\frac{m}{p} [cm]");
  59    //fLambdas->Sumw2();
  60   fOutput->Add(fLambdas);
  61   fLambdasMC = new TH1F("LambdasMC","c\\tau for \\Lambda, MC", 100, 0.,80.);
  62   fLambdasMC->GetXaxis()->SetTitle("L\\frac{m}{p} [cm]");
  63    //fLambdasMC->Sumw2();
  64   fOutput->Add(fLambdasMC);
  65
  66   fLambdaBars = new TH1F("LambdaBars","c\\tau for \\bar{\\Lambda}",100,0.,80.);
  67   fLambdaBars->GetXaxis()->SetTitle("L\\frac{m}{p} [cm]");
  68    //fLambdaBars->Sumw2();
  69   fOutput->Add(fLambdaBars);
  70   fLambdaBarsMC=new TH1F("LambdaBarsMC","c\\tau for \\bar{\\Lambda}, MC",100,0.,80.);
  71   fLambdaBarsMC->GetXaxis()->SetTitle("L\\frac{m}{p} [cm]");
  72    //fLambdaBarsMC->Sumw2();
  73   fOutput->Add(fLambdaBarsMC);
  74
  75   fXis = new TH1F("Xis","c\\tau for \\Xi^{-} + \\Xi^{+}",50,0.,30.);//10*4.91)
  76   fXis->GetXaxis()->SetTitle("L\\frac{m}{p} [cm]");
  77   fOutput->Add(fXis);
  78   fXisMC = new TH1F("XisMC","c\\tau for \\Xi^{-} + \\Xi^{+}, MC",50,0.,30.);
  79   fXisMC->GetXaxis()->SetTitle("L\\frac{m}{p} [cm]");
  80   fOutput->Add(fXisMC);
  81
  82    //fMass = new TH1F("Mass", "Mass", 50, 0.4477, 0.5477);
  83   fMass = new TH1F("Mass", "Mass", 50, 1.07, 1.16);
  84    //fMass = new TH1F("Mass", "Mass", 50, 1.271, 1.371);
  85   fOutput->Add(fMass);
  86   fMassMC = new TH1F("MassMC", "Mass, MC", 50, 1.07, 1.16);
  87   fOutput->Add(fMassMC);
  88
  89 }
  90
  91 void AliAnalysisTaskCTau::UserExec(Option_t *)
  92 {
  93
  94   fESD=(AliESDEvent *)InputEvent();
  95
  96   if (!fESD) {
  97     Printf("ERROR: fESD not available");
  98     return;
  99   }
 100
 101   //Char_t *tname="CINT1B-ABCE-NOPF-ALL";
 102   //Int_t ok=fESD->IsTriggerClassFired(tname);
 103   //if (!ok) return;
 104
 105   const AliESDVertex *vtx=fESD->GetPrimaryVertexSPD();
 106   //if (!vtx->GetStatus()) {
 107   //   vtx=fESD->GetPrimaryVertexTracks();
 108   //   if (!vtx->GetStatus()) return;
 109   //}
 110   Double_t xv=vtx->GetXv(), yv=vtx->GetYv(), zv=vtx->GetZv();
 111
 112
 113   //+++++++ MC
 114   AliStack *stack = 0x0;
 115   AliMCEvent *mcEvent = MCEvent();
 116   if (mcEvent) {
 117      stack = mcEvent->Stack();
 118      if (!stack) {
 119         Printf("ERROR: stack not available");
 120         return;
 121      }
 122
 123      const AliVVertex *mcVtx=mcEvent->GetPrimaryVertex();
 124      Double_t mcXv=mcVtx->GetX(), mcYv=mcVtx->GetY(), mcZv=mcVtx->GetZ();
 125      Int_t ntrk=stack->GetNtrack();
 126      while (ntrk--) {
 127        TParticle *p0=stack->Particle(ntrk);
 128        Int_t code=p0->GetPdgCode();
 129        if (code != kK0Short)
 130        if (code != kLambda0)
 131        if (code != kLambda0Bar) continue;
 132
 133        Int_t plab=p0->GetFirstDaughter(), nlab=p0->GetLastDaughter();
 134        if (nlab==plab) continue;
 135        if (nlab<0) continue;
 136        if (plab<0) continue;
 137        TParticle *part = stack->Particle(plab);
 138        if (!part) continue;
 139        TParticlePDG *partPDG = part->GetPDG();
 140        if (!partPDG) continue;
 141        Double_t charge=partPDG->Charge();
 142        if (charge == 0.) continue;
 143        if (charge < 0.) {
 144           Int_t i=plab; plab=nlab; nlab=i;
 145        }
 146
 147        Double_t pz=p0->Pz(), pt=p0->Pt(), p=p0->P();
 148        if (pt<0.1) continue;
 149        if (TMath::Abs(pz/pt)>0.7) continue;
 150        if (p<0.9) continue;
 151
 152        Double_t x=p0->Vx(), y=p0->Vy(), z=p0->Vz();
 153        Double_t dx=mcXv-x, dy=mcYv-y, dz=mcZv-z;
 154        Double_t l=TMath::Sqrt(dx*dx + dy*dy + dz*dz);
 155
 156        if (l > 0.1) continue; // secondary V0
 157
 158        x=part->Vx(); y=part->Vy(); z=part->Vz();
 159        dx=mcXv-x; dy=mcYv-y; dz=mcZv-z;
 160        l=TMath::Sqrt(dx*dx + dy*dy + dz*dz);
 161
 162        if (code == kK0Short) {
 163           Double_t m=TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass();
 164           Double_t ct=l*m/p;
 165           fK0sMC->Fill(ct);
 166        }
 167        if (code == kLambda0) {
 168           Double_t m=TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass();
 169           Double_t ct=l*m/p;
 170           fLambdasMC->Fill(ct);
 171        }
 172        if (code == kLambda0Bar) {
 173          Double_t m=TDatabasePDG::Instance()->GetParticle(kLambda0Bar)->Mass();
 174          Double_t ct=l*m/p;
 175          fLambdaBarsMC->Fill(ct);
 176        }
 177      }
 178   }
 179
 180   Int_t nv0 = fESD->GetNumberOfV0s();
 181   while (nv0--) {
 182       AliESDv0 *v0=fESD->GetV0(nv0);
 183
 184       if (v0->GetOnFlyStatus()) continue;
 185
 186       Int_t nidx=TMath::Abs(v0->GetNindex());
 187       AliESDtrack *ntrack=fESD->GetTrack(nidx);
 188       if (!ntrack->IsOn(AliESDtrack::kTPCrefit)) continue;
 189       //if (!ntrack->HasPointOnITSLayer(5)) continue;
 190
 191       Int_t pidx=TMath::Abs(v0->GetPindex());
 192       AliESDtrack *ptrack=fESD->GetTrack(pidx);
 193       if (!ptrack->IsOn(AliESDtrack::kTPCrefit)) continue;
 194       //if (!ptrack->HasPointOnITSLayer(5)) continue;
 195
 196
 197       Double_t pz=v0->Pz(),pt=v0->Pt(),p=v0->P();
 198       if (pt<0.1) continue;
 199       if (TMath::Abs(pz/pt)>0.7) continue;
 200       if (p<0.9) continue;
 201
 202       Double_t cpa=v0->GetV0CosineOfPointingAngle();
 203       if (cpa<0.995) continue;
 204
 205
 206       Double_t x,y,z; v0->GetXYZ(x,y,z);
 207       Double_t dx=x-xv, dy=y-yv, dz=z-zv;
 208       //Double_t r=TMath::Sqrt(x*x + y*y);
 209       //if (r > 3.) continue;
 210       Double_t l=TMath::Sqrt(dx*dx + dy*dy + dz*dz);
 211
 212       v0->ChangeMassHypothesis(kK0Short);
 213       Double_t mass=v0->GetEffMass();
 214       Double_t m=TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass();
 215       if (TMath::Abs(m-mass)<0.016) {
 216          Double_t ct=l*m/p;
 217          fK0s->Fill(ct);
 218       }
 219
 220       //if (ntrack->GetP()<0.7) continue;
 221
 222       v0->ChangeMassHypothesis(kLambda0);
 223       mass=v0->GetEffMass();
 224       m=TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass();
 225       if (TMath::Abs(m-mass)<0.006) {
 226          Double_t ct=l*m/p;
 227          fLambdas->Fill(ct);
 228       }
 229
 230       v0->ChangeMassHypothesis(kLambda0Bar);
 231       mass=v0->GetEffMass();
 232       m=TDatabasePDG::Instance()->GetParticle(kLambda0Bar)->Mass();
 233       if (TMath::Abs(m-mass)<0.006) {
 234          Double_t ct=l*m/p;
 235          fLambdaBars->Fill(ct);
 236       }
 237       fMass->Fill(mass);
 238
 239
 240       //+++++++ MC
 241       Int_t nlab=TMath::Abs(ntrack->GetLabel());
 242       Int_t plab=TMath::Abs(ptrack->GetLabel());
 243
 244       if ((plab-nlab) != 1)
 245         if ((nlab-plab) != 1) continue;
 246
 247       if (!stack) continue;
 248       TParticle *np=stack->Particle(nlab);
 249       Int_t i0=np->GetFirstMother();
 250       TParticle *p0=stack->Particle(i0);
 251       if (!p0) continue;
 252       if (p0->GetPdgCode() == kLambda0Bar) fMassMC->Fill(mass);
 253
 254   }
 255
 256   Double_t kine0;
 257   Int_t ncs=fESD->GetNumberOfCascades();
 258   for (Int_t i=0; i<ncs; i++) {
 259       AliESDcascade *cs=fESD->GetCascade(i);
 260
 261       Int_t nidx=TMath::Abs(cs->GetNindex());
 262       AliESDtrack *ntrack=fESD->GetTrack(nidx);
 263       if (!ntrack->IsOn(AliESDtrack::kTPCrefit)) continue;
 264
 265       Int_t pidx=TMath::Abs(cs->GetPindex());
 266       AliESDtrack *ptrack=fESD->GetTrack(pidx);
 267       if (!ptrack->IsOn(AliESDtrack::kTPCrefit)) continue;
 268
 269       Int_t bidx=TMath::Abs(cs->GetBindex());
 270       AliESDtrack *btrack=fESD->GetTrack(bidx);
 271       if (!btrack->IsOn(AliESDtrack::kTPCrefit)) continue;
 272
 273       Int_t charge=cs->Charge();
 274       if (charge<0)
 275          cs->ChangeMassHypothesis(kine0,kXiMinus);
 276       else
 277          cs->ChangeMassHypothesis(kine0,kXiPlusBar);
 278
 279       Double_t mass=cs->GetEffMassXi();
 280
 281       if (mass>1.305)
 282       if (mass<1.335) {
 283         //TTree *tree=fChain->GetTree();
 284         //TDirectory *dir=tree->GetDirectory();
 285         //Int_t n=fESD->GetEventNumberInFile();
 286         //Printf("%d %d Xi found !!! %s %d\n", ncs, charge, dir->GetName(), n);
 287         Double_t p=cs->P();
 288         Double_t m=TDatabasePDG::Instance()->GetParticle(kXiMinus)->Mass();
 289         Double_t x,y,z; cs->GetXYZcascade(x,y,z);
 290         Double_t dx=x-xv, dy=y-yv, dz=z-zv;
 291         Double_t l=TMath::Sqrt(dx*dx + dy*dy + dz*dz);
 292
 293         if (p<0.1) continue;
 294
 295         Double_t ct=l*m/p;
 296         fXis->Fill(ct);
 297       }
 298
 299   }
 300
 301   PostData(1, fOutput);
 302 }
 303
 304 void AliAnalysisTaskCTau::Terminate(Option_t *)
 305 {
 306    // The Terminate() function is the last function to be called during
 307    // a query. It always runs on the client, it can be used to present
 308    // the results graphically or save the results to file.
 309
 310   fOutput=(TList*)GetOutputData(1);
 311   if (!fOutput) {
 312      Printf("ERROR: fOutput not available");
 313      return;
 314   }
 315
 316   fMass = dynamic_cast<TH1F*>(fOutput->FindObject("Mass")) ;
 317   if (!fMass) {
 318      Printf("ERROR: fMass not available");
 319      return;
 320   }
 321   fLambdas = dynamic_cast<TH1F*>(fOutput->FindObject("Lambdas")) ;
 322   if (!fLambdas) {
 323      Printf("ERROR: fLambdas not available");
 324      return;
 325   }
 326   fLambdaBars = dynamic_cast<TH1F*>(fOutput->FindObject("LambdaBars")) ;
 327   if (!fLambdaBars) {
 328      Printf("ERROR: fLambdaBars not available");
 329      return;
 330   }
 331   fXis = dynamic_cast<TH1F*>(fOutput->FindObject("Xis")) ;
 332   if (!fXis) {
 333      Printf("ERROR: fXis not available");
 334      return;
 335   }
 336
 337   if (!gROOT->IsBatch()) {
 338     TCanvas *c1 = new TCanvas("c1","Lambdas");
 339     c1->SetFillColor(10);
 340     c1->SetHighLightColor(10);
 341     fLambdas->DrawCopy("E") ;
 342     TCanvas *c2 = new TCanvas("c2","LambdaBar");
 343     c2->SetFillColor(10);
 344     c2->SetHighLightColor(10);
 345     fLambdaBars->DrawCopy("E") ;
 346   }
 347 }