PWGUD/DIFFRACTIVE/AliAnalysisTaskDDMeson.cxx

   1 /**************************************************************************
   2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15 //
  16 // Reconstruct mesons
  17 // in central diffractive event
  18 // in pp collisions
  19 //
  20 // Authors:
  21 //  Xianguo Lu <lu@physi.uni-heidelberg.de>
  22 //
  23
  24 #include <TH1I.h>
  25 #include <TH2I.h>
  26 #include <TH2D.h>
  27 #include <TList.h>
  28 #include <TLorentzVector.h>
  29 #include <TRandom3.h>
  30 #include <THnSparse.h>
  31
  32 #include "AliCDBManager.h"
  33 #include "AliCDBEntry.h"
  34 #include "AliESDInputHandler.h"
  35 #include "AliESDtrackCuts.h"
  36 #include "AliGeomManager.h"
  37 #include "AliITSAlignMille2Module.h"
  38 #include "AliITSsegmentationSPD.h"
  39 #include "AliKFParticle.h"
  40 #include "AliMultiplicity.h"
  41 #include "AliPID.h"
  42 #include "AliSPDUtils.h"
  43 #include "AliTriggerAnalysis.h"
  44
  45 #include "AliAnalysisTaskDDMeson.h"
  46
  47 void AliAnalysisTaskDDMeson::IniTask()
  48 {
  49   //
  50   //initialize values used in many places
  51   //
  52   fnmass = 512;
  53   fmass1 = 0.01 * fnmass;
  54
  55   fnptv = 10;
  56   fptv1 = 0.2 * fnptv;
  57
  58   fnpta = 51*2;
  59   fpta1 = 0.05 * fnpta;
  60
  61   fneta = 64;
  62   feta = 0.04 * fneta/2.;
  63
  64   fnnsel = 2;
  65   fnsel1 = 2+fnnsel;
  66
  67   fncts = 10;
  68   fnctlab = 20;
  69
  70   fCHECKVBA = 2;
  71   fCHECKVBC = 4;
  72
  73   //------------------
  74   fOpt.ToUpper();
  75 }
  76
  77 AliAnalysisTaskDDMeson::AliAnalysisTaskDDMeson(const TString opt):
  78   AliAnalysisTaskSE("DDMeson")
  79   , fOpt(opt)
  80   , fESD(0x0)
  81
  82   , fnmass(-999)
  83   , fmass1(-999)
  84   , fnptv(-999)
  85   , fptv1(-999)
  86   , fnpta(-999)
  87   , fpta1(-999)
  88   , fneta(-999)
  89   , feta(-999)
  90   , fnnsel(-999)
  91   , fnsel1(-999)
  92   , fncts(-999)
  93   , fnctlab(-999)
  94   , fCHECKVBA(-999)
  95   , fCHECKVBC(-999)
  96
  97   , fHBIT(0x0)
  98   , fBitcg(0)
  99
 100   , fRun(-999)
 101   , fat(0x0)
 102   , fct(0x0)
 103   , fbt(0x0)
 104   , fnt(0x0)
 105   , ftt(0x0)
 106
 107   , fv0ntrk(0x0)
 108   , frsntrk(0x0)
 109
 110   , fhps(0x0)
 111   , fhfo(0x0)
 112   , fhspd(0x0)
 113   , fhv0fmd(0x0)
 114   , fhpriv(0x0)
 115   , fhntrk(0x0)
 116
 117   , fHist(0x0)
 118   , fThnMass(0x0)
 119   , fThnDPt(0x0)
 120   , fThnDEta(0x0)
 121   , fThnKF(0x0)
 122 {
 123   //
 124   // Dummy constructor
 125   //
 126   //slot in TaskSE must start from 1
 127   DefineOutput(1, TList::Class());
 128   DefineOutput(2, THnSparse::Class());
 129   DefineOutput(3, THnSparse::Class());
 130   DefineOutput(4, THnSparse::Class());
 131   DefineOutput(5, THnSparse::Class());
 132
 133   IniTask();
 134 }
 135
 136 AliAnalysisTaskDDMeson::AliAnalysisTaskDDMeson(const AliAnalysisTaskDDMeson  &p):
 137   AliAnalysisTaskSE("DDMeson")
 138   , fOpt(p.fOpt)
 139   , fESD(p.fESD)
 140
 141   , fnmass(p.fnmass)
 142   , fmass1(p.fmass1)
 143   , fnptv(p.fnptv)
 144   , fptv1(p.fptv1)
 145   , fnpta(p.fnpta)
 146   , fpta1(p.fpta1)
 147   , fneta(p.fneta)
 148   , feta(p.feta)
 149   , fnnsel(p.fnnsel)
 150   , fnsel1(p.fnsel1)
 151   , fncts(p.fncts)
 152   , fnctlab(p.fnctlab)
 153   , fCHECKVBA(p.fCHECKVBA)
 154   , fCHECKVBC(p.fCHECKVBC)
 155
 156   , fHBIT(p.fHBIT)
 157   , fBitcg(p.fBitcg)
 158
 159   , fRun(p.fRun)
 160   , fat(p.fat)
 161   , fct(p.fct)
 162   , fbt(p.fbt)
 163   , fnt(p.fnt)
 164   , ftt(p.ftt)
 165
 166   , fv0ntrk(p.fv0ntrk)
 167   , frsntrk(p.frsntrk)
 168
 169   , fhps(p.fhps)
 170   , fhfo(p.fhfo)
 171   , fhspd(p.fhspd)
 172   , fhv0fmd(p.fhv0fmd)
 173   , fhpriv(p.fhpriv)
 174   , fhntrk(p.fhntrk)
 175
 176   , fHist(p.fHist)
 177   , fThnMass(p.fThnMass)
 178   , fThnDPt(p.fThnDPt)
 179   , fThnDEta(p.fThnDEta)
 180   , fThnKF(p.fThnKF)
 181 {
 182   //
 183   // Copy constructor
 184   //
 185
 186   IniTask();
 187 }
 188
 189 AliAnalysisTaskDDMeson & AliAnalysisTaskDDMeson::operator=(const AliAnalysisTaskDDMeson  &p)
 190 {
 191   //
 192   // overload =
 193   //
 194   if(&p == this) return *this;
 195   AliAnalysisTaskSE::operator=(p);
 196
 197   fESD = p.fESD;
 198
 199   fnmass = p.fnmass;
 200   fmass1 = p.fmass1;
 201   fnptv = p.fnptv;
 202   fptv1 = p.fptv1;
 203   fnpta = p.fnpta;
 204   fpta1 = p.fpta1;
 205   fneta = p.fneta;
 206   feta = p.feta;
 207   fnnsel = p.fnnsel;
 208   fnsel1 = p.fnsel1;
 209   fncts = p.fncts;
 210   fnctlab = p.fnctlab;
 211   fCHECKVBA = p.fCHECKVBA;
 212   fCHECKVBC = p.fCHECKVBC;
 213
 214   fHBIT = p.fHBIT;
 215   fBitcg = p.fBitcg;
 216
 217   fRun = p.fRun;
 218
 219   fat = p.fat;
 220   fct = p.fct;
 221   fbt = p.fbt;
 222   fnt = p.fnt;
 223   ftt = p.ftt;
 224
 225   fv0ntrk = p.fv0ntrk;
 226   frsntrk = p.frsntrk;
 227
 228   fhps = p.fhps;
 229   fhfo = p.fhfo;
 230   fhspd = p.fhspd;
 231   fhv0fmd = p.fhv0fmd;
 232   fhpriv = p.fhpriv;
 233   fhntrk = p.fhntrk;
 234
 235   fHist = p.fHist;
 236   fThnMass = p.fThnMass;
 237   fThnDPt = p.fThnDPt;
 238   fThnDEta = p.fThnDEta;
 239   fThnKF = p.fThnKF;
 240
 241   IniTask();
 242
 243   return *this;
 244 }
 245
 246 AliAnalysisTaskDDMeson::~AliAnalysisTaskDDMeson()
 247 {
 248   //
 249   //Destructor
 250   //
 251   if(fESD) delete fESD;
 252
 253   delete fHBIT;
 254   delete fat;
 255   delete fct;
 256   delete fbt;
 257   delete fnt;
 258   delete ftt;
 259
 260   delete fv0ntrk;
 261   delete frsntrk;
 262
 263   delete fhps;
 264   delete fhfo;
 265   delete fhspd;
 266   delete fhv0fmd;
 267   delete fhpriv;
 268   delete fhntrk;
 269
 270   delete fThnMass;
 271   delete fThnDPt;
 272   delete fThnDEta;
 273   delete fThnKF;
 274
 275   fHist->Clear();
 276   delete fHist;
 277
 278 }
 279
 280 void AliAnalysisTaskDDMeson::CheckRange(Double_t &ptv, Double_t &pta, Double_t &etaa
 281                                         , Double_t &mpi
 282                                         ) const
 283 {
 284   //
 285   //save over/under flow bins
 286   //
 287
 288   const Double_t eps = 1e-6;
 289   if(ptv>=fptv1) ptv = fptv1 - eps;
 290   if(pta>=fpta1) pta = fpta1 - eps;
 291
 292   if(etaa >= feta) etaa = feta - eps;
 293
 294   if(etaa <= -feta) etaa = -feta + eps;
 295
 296   if(mpi >= fmass1) mpi = fmass1 - eps;
 297
 298 }
 299
 300 void AliAnalysisTaskDDMeson::UserCreateOutputObjects()
 301 {
 302   //
 303   //CreateOutputObjects
 304   //
 305   //=======================================
 306   const Double_t kvz0=1, kvz1=5;
 307   const Int_t nkvz=(Int_t)(kvz1-kvz0);
 308
 309   const Double_t charge0=1, charge1=4;
 310   const Int_t ncharge=(Int_t)(charge1-charge0);
 311
 312   const Double_t mass0 = 0;
 313   //=======================================
 314   //total momentum and pt in lab frame:
 315   const Double_t p0=0;
 316   //=======================================
 317   const Double_t nsel0 = 2;
 318
 319   //=======================================
 320   const Double_t cts0=0, cts1=1;
 321
 322   //=======================================
 323   const Double_t ctlab0 = -1, ctlab1 = 1;
 324
 325   //=======================================
 326   //------- Mass
 327   const Int_t    nparPtMass = 7;
 328   const Int_t    nbinPtMass[]     ={fnnsel, nkvz, ncharge,  fnmass, fnptv, fncts, fnctlab};
 329   const Double_t xminPtMass[]     ={nsel0,  kvz0, charge0,   mass0,    p0,  cts0,  ctlab0};
 330   const Double_t xmaxPtMass[]     ={fnsel1, kvz1, charge1,  fmass1, fptv1,  cts1,  ctlab1};
 331   fThnMass = new THnSparseD("DDMeson_Mass","nsel, kvz, charge, Mpipi, ptv, cts, ctlab", nparPtMass, nbinPtMass, xminPtMass, xmaxPtMass);
 332
 333   //------- DPt
 334   const Int_t    nparDPt = 4;
 335   const Int_t    nbinDPt[] ={fnnsel, nkvz, ncharge,  fnpta};
 336   const Double_t xminDPt[] ={nsel0, kvz0, charge0,  p0};
 337   const Double_t xmaxDPt[] ={fnsel1, kvz1, charge1,  fpta1};
 338   fThnDPt = new THnSparseD("DDMeson_DPt","nsel, kvz, charge, pt1", nparDPt, nbinDPt, xminDPt, xmaxDPt);
 339   //------- DEta
 340   const Int_t    nparDEta = 4;
 341   const Int_t    nbinDEta[] ={fnnsel, nkvz, ncharge,  fneta};
 342   const Double_t xminDEta[] ={nsel0, kvz0, charge0,  -feta};
 343   const Double_t xmaxDEta[] ={fnsel1, kvz1, charge1,  feta};
 344   fThnDEta = new THnSparseD("DDMeson_DEta"," nsel, kvz, charge, eta1", nparDEta, nbinDEta, xminDEta, xmaxDEta);
 345
 346   //------- KF
 347   const Double_t mks = 0.494;
 348   const Double_t dks = 0.024;
 349   const Int_t nparKF = 4;
 350   const Int_t nbinKF[] = {nkvz, 200, 200, 200};
 351   const Double_t xminKF[] ={kvz0, mks-dks, mks-dks, 0};
 352   const Double_t xmaxKF[] ={kvz1, mks+dks, mks+dks, 10};
 353   fThnKF = new THnSparseD("DDMeson_KF","kvz, rawks, kfks, chi2", nparKF, nbinKF, xminKF, xmaxKF);
 354   //=======================================
 355   //=======================================
 356   fHBIT = new TH1I("HBIT","",32,1,33);
 357
 358   const Int_t runb0 = 114650; //runa1=114649, runb1=117630, runc0=117631, runc1=121526, rund0=121527, rund1=126460, rune0=126461, rune1 = 130930; //runf0=130931, all checked on elog
 359   const Int_t runf1 = 133900;//02OCT2010
 360   const Int_t run0 = runb0-1, run1 = runf1+1;
 361   const Int_t nrun = (Int_t)(run1-run0);
 362
 363   fat = new TH1I("at","",nrun,run0,run1);
 364   fct = new TH1I("ct","",nrun,run0,run1);
 365   fbt = new TH1I("bt","",nrun,run0,run1);
 366   fnt = new TH1I("nt","",nrun,run0,run1);
 367   ftt = new TH1I("tt","",nrun,run0,run1);
 368
 369   fv0ntrk = new TH2D("c_v0ntrk","",80,0,80, 4, 1, 5);//x: ntrk; y: GetV0
 370   frsntrk = new TH2D("c_rsntrk","",1000,0,1000, 200, 0,200);
 371
 372   fhps =    new TH1I("a000_ps",    "", 20, 0,20);
 373   fhfo =    new TH2I("a010_fo",    "", 50, 0, 50, 50, 0, 50);
 374   fhspd =   new TH1I("a011_spd",   "", 50, 0, 50);
 375   fhv0fmd = new TH2I("a020_v0fmd", "", 4, 0, 4, 4, 0, 4);
 376   fhpriv =  new TH1I("a100_priv",  "", 2, 0,2);
 377   fhntrk =  new TH1I("a110_ntrk",  "", (Int_t)fnsel1, 0, (Int_t)fnsel1);
 378
 379   //------
 380
 381   fHist = new TList;
 382
 383   //------>>> Very important!!!!!!!
 384   fHist->SetOwner();
 385
 386   fHist->Add(fat);
 387   fHist->Add(fct);
 388   fHist->Add(fbt);
 389   fHist->Add(fnt);
 390   fHist->Add(ftt);
 391
 392   fHist->Add(fv0ntrk);
 393   fHist->Add(frsntrk);
 394   fHist->Add(fhps);
 395   fHist->Add(fhfo);
 396   fHist->Add(fhspd);
 397   fHist->Add(fhv0fmd);
 398   fHist->Add(fhpriv);
 399   fHist->Add(fhntrk);
 400 }
 401
 402 void AliAnalysisTaskDDMeson::UserExec(Option_t *)
 403 {
 404   //
 405   //Execute
 406   //
 407   //----------------------------------- general protection
 408   if(!CheckESD())
 409     return;
 410
 411   //----------------------------------- trigger bits
 412   if(!CheckBit())
 413     return;
 414
 415   //----------------------------------- track cuts
 416   const Int_t ntrk0 = fESD->GetNumberOfTracks();
 417   const AliESDtrack * trks[ntrk0];
 418   for(Int_t ii=0; ii<ntrk0; ii++){
 419     trks[ii] = 0x0;
 420   }
 421   const Int_t nsel = CutESD(trks);
 422   if(nsel<2)
 423     return;
 424
 425   for(Int_t ii=0; ii<nsel; ii++){
 426     for(Int_t jj=ii+1; jj<nsel; jj++){
 427       const AliESDtrack * trkpair[2]={trks[ii], trks[jj]};
 428
 429       //----------------------------------- tracks initailization
 430       TRandom3 tmprd(0);
 431       if(tmprd.Rndm()>0.5){
 432         SwapTrack(trkpair);
 433       }
 434
 435       //------------------------------------------------------------------------
 436       const Int_t pipdg = 211;
 437       const AliKFParticle daughter0(*(trkpair[0]), pipdg);
 438       const AliKFParticle daughter1(*(trkpair[1]), pipdg);
 439
 440       const AliKFParticle parent(daughter0, daughter1);
 441       const Double_t kfchi2 = parent.GetChi2();
 442       const Double_t kfmass = parent.GetMass();
 443
 444       //------------------------------------------------------------------------
 445       const Int_t kvz = GetV0();
 446
 447       const AliESDtrack * trk1=trkpair[0];
 448       const AliESDtrack * trk2=trkpair[1];
 449
 450       const Double_t ch1 = trk1->GetSign();
 451       const Double_t ch2 = trk2->GetSign();
 452       const Int_t combch = GetCombCh(ch1, ch2);
 453
 454       Double_t pta  = trk1->Pt();
 455       Double_t etaa = trk1->Eta();
 456
 457       Double_t tmpp1[3], tmpp2[3];
 458       trk1->GetPxPyPz(tmpp1);
 459       trk2->GetPxPyPz(tmpp2);
 460
 461       //------------------------------------------------------------------------
 462       //------------------------------------------------------------------------
 463       const Double_t ctlab = GetCtlab(tmpp1, tmpp2);
 464
 465       Double_t cts = -999;
 466       const Double_t masspi = AliPID::ParticleMass(AliPID::kPion);
 467       const TLorentzVector sumv = GetKinematics(tmpp1, tmpp2, masspi, masspi, cts);
 468       Double_t pipi = sumv.Mag();
 469       Double_t ptv = sumv.Pt();
 470
 471       CheckRange(ptv, pta, etaa, pipi);
 472
 473       const Double_t varMass[]  ={nsel, kvz, combch
 474                                        , pipi, ptv, cts
 475                                        , ctlab
 476       };
 477       const Double_t varDPt[]        ={nsel, kvz, combch, pta};
 478       const Double_t varDEta[]       ={nsel, kvz, combch, etaa};
 479       fThnMass->Fill(varMass);
 480       fThnDPt->Fill(varDPt);
 481       fThnDEta->Fill(varDEta);
 482
 483       const Double_t varKF[] = {kvz, pipi, kfmass, kfchi2};
 484       if(nsel==2 && combch==3)
 485         fThnKF->Fill(varKF);
 486     }
 487   }
 488   //------------------------------------------------------------------------
 489   //------------------------------------------------------------------------
 490
 491   PostData(1, fHist);
 492   PostData(2, fThnMass);
 493   PostData(3, fThnDPt);
 494   PostData(4, fThnDEta);
 495   PostData(5, fThnKF);
 496 }
 497 //=======================================================================
 498 //=======================================================================
 499
 500 void AliAnalysisTaskDDMeson::Terminate(Option_t *)
 501 {
 502   //
 503   //fillbit before leaving
 504   //
 505   FillBit();
 506 }
 507
 508 void AliAnalysisTaskDDMeson::FillBit()
 509 {
 510   //
 511   //save the bit count in f*t
 512   //
 513   Double_t tot[5];
 514   CalcBit(fHBIT, tot);
 515
 516   const Int_t bin = fat->FindBin(fRun);
 517   fat->AddBinContent(bin, tot[0]);
 518   fct->AddBinContent(bin, tot[1]);
 519   fbt->AddBinContent(bin, tot[2]);
 520   fnt->AddBinContent(bin, tot[3]);
 521   ftt->AddBinContent(bin, tot[4]);
 522 }
 523
 524 void AliAnalysisTaskDDMeson::CalcBit(TH1I *hc, Double_t tot[]) const
 525 {
 526   //
 527   //calculate the bit count in f*t
 528   //
 529   if(hc->GetBinLowEdge(1)!=1){
 530     printf("xlulog hc defined incorrectly!! %f\n",hc->GetBinLowEdge(1));
 531     return;
 532   }
 533   if(hc->GetBinContent(0)){
 534     printf("\nxlulog hc !!!!!!!!!!!!!!!!!!! underflow!!!!!!!!!!!!!!!!! %f\n\n", hc->GetBinContent(0));
 535     return;
 536   }
 537   if(hc->GetBinContent(hc->GetNbinsX())){
 538     printf("\nxlulog hc !!!!!!!!!!!!!!!!! OVERflow!!!!!!!!!!!!!! %f\n\n", hc->GetBinContent(hc->GetNbinsX()));
 539     return;
 540   }
 541
 542
 543   Double_t atot=0, ctot=0, btot=0, ntot=0;
 544   for(Int_t ib=1; ib<=hc->GetNbinsX(); ib++){
 545     const Double_t cont = hc->GetBinContent(ib);
 546     if(       (ib & fCHECKVBA) && !(ib & fCHECKVBC) ){
 547       atot += cont;
 548     }
 549     else if( !(ib & fCHECKVBA) &&  (ib & fCHECKVBC) ){
 550       ctot += cont;
 551     }
 552     else if(  (ib & fCHECKVBA) &&  (ib & fCHECKVBC) ){
 553       btot += cont;
 554     }
 555     else{
 556       ntot += cont;
 557     }
 558   }
 559
 560   const Double_t ttot   = atot + ctot + btot + ntot;
 561
 562   tot[0]=atot;
 563   tot[1]=ctot;
 564   tot[2]=btot;
 565   tot[3]=ntot;
 566   tot[4]=ttot;
 567
 568   //---
 569   char fullname[5][20];
 570   strncpy(fullname[0],"V0A-ONLY",19);
 571   strncpy(fullname[1],"V0C-ONLY",19);
 572   strncpy(fullname[2],"V0A & V0C",19);
 573   strncpy(fullname[3],"!V0A & !V0C",19);
 574   strncpy(fullname[4],"TOTAL",19);
 575   for (Int_t ii=0;ii<5;ii++){
 576     fullname[ii][19] = '\0';
 577   }
 578
 579   for(Int_t ii=0;ii<5;ii++){
 580     printf("xlulog %s CTP-L0: %s\tTOT: %10.0f\n",fOpt.Data(), fullname[ii],tot[ii]);
 581   }
 582
 583   //---
 584
 585   const Int_t nb=hc->GetNbinsX();
 586   for(Int_t ib=0; ib<=nb+1; ib++){
 587     hc->SetBinContent(ib, 0);
 588   }
 589 }
 590
 591 //====================================================================================
 592
 593 Bool_t AliAnalysisTaskDDMeson::CheckESD()
 594 {
 595   //
 596   //general protection
 597   //
 598   const AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler *>(fInputHandler);
 599   if(esdH)
 600     fESD = esdH->GetEvent();
 601   if(!fESD){
 602     AliError("xlulog No ESD Event");
 603     return 0;
 604   }
 605
 606   if(fabs(fESD->GetMagneticField())<1){
 607     printf("xlulog strange Bfield! %f\n", fESD->GetMagneticField());
 608     return 0;
 609   }
 610
 611   const Int_t tmprun = fESD->GetRunNumber();
 612
 613   if(fRun!=tmprun){
 614     if(fRun>0){
 615       printf("xlulog run changed during exec!! %d %d\n", tmprun, fRun);
 616       FillBit();
 617     }
 618     fRun = tmprun;
 619     if(fOpt.Contains("SPD")){
 620       SPDLoadGeom();
 621     }
 622   }
 623
 624   return 1;
 625 }
 626
 627 void AliAnalysisTaskDDMeson::SPDLoadGeom() const
 628 {
 629   //method to get the gGeomanager
 630   // it is called at the CreatedOutputObject stage
 631   // to comply with the CAF environment
 632   AliCDBManager *man = AliCDBManager::Instance();
 633   TString cdbpath("local:///lustre/alice/alien/alice/data/2010/OCDB");
 634   man->SetDefaultStorage(cdbpath);
 635   man->SetRun(fRun);
 636   man->SetSpecificStorage("ITS/Align/Data",cdbpath);
 637   man->SetSpecificStorage("GRP/Geometry/Data",cdbpath);
 638
 639   AliCDBEntry* obj = man->Get(AliCDBPath("GRP", "Geometry", "Data"));
 640   if (obj != NULL) {
 641     AliGeomManager::SetGeometry((TGeoManager*)obj->GetObject());
 642     AliGeomManager::ApplyAlignObjsFromCDB("ITS");
 643     printf("xlulog DDMeson %s Geom Loaded for run %d !\n", fOpt.Data(), fRun);
 644   }
 645   else {
 646     printf("xlulog DDMeson %s was unable to load Geom for run %d !\n", fOpt.Data(), fRun);
 647   }
 648 }
 649
 650 Bool_t AliAnalysisTaskDDMeson::SPDLoc2Glo(const Int_t id, const Double_t *loc, Double_t *glo) const
 651 {
 652   //
 653   //SPDLoc2Glo
 654   //
 655
 656   static TGeoHMatrix mat;
 657   Int_t vid = AliITSAlignMille2Module::GetVolumeIDFromIndex(id);
 658   if (vid<0) {
 659     printf("xlulog Did not find module with such ID %d\n",id);
 660     return kFALSE;
 661   }
 662   AliITSAlignMille2Module::SensVolMatrix(vid,&mat);
 663   mat.LocalToMaster(loc,glo);
 664   return kTRUE;
 665 }
 666
 667 Bool_t AliAnalysisTaskDDMeson::CheckChipEta(const Int_t chipKey) const
 668 {
 669   //
 670   //CheckChipEta
 671   //
 672
 673   //no eta cut
 674   if(fOpt.Contains("SPD0"))
 675     return kTRUE;
 676
 677   const Double_t etacut = 1.0;
 678
 679   UInt_t module=999, offchip=999;
 680   AliSPDUtils::GetOfflineFromOfflineChipKey(chipKey,module,offchip);
 681   UInt_t hs = AliSPDUtils::GetOnlineHSFromOffline(module);
 682   if(hs<2) offchip = 4 - offchip; // inversion  in the inner layer...
 683
 684   const Int_t col[]={
 685     hs<2? 0 : 31,
 686     hs<2? 31 : 0,
 687     hs<2? 31 : 0,
 688     hs<2? 0 : 31};
 689   const Int_t aa[]={0, 0, 255, 255};
 690   const AliITSsegmentationSPD seg;
 691
 692   for(Int_t ic=0; ic<4; ic++){
 693     Float_t localchip[3]={0.,0.,0.};
 694     seg.DetToLocal(aa[ic],col[ic]+32*offchip,localchip[0],localchip[2]); // local coordinate of the chip center
 695     //printf("local coordinates %d %d: %f %f \n",chipKey, ic, localchip[0],localchip[2]);
 696     const Double_t local[3] = {localchip[0],localchip[1],localchip[2]};
 697     Double_t glochip[3]={0.,0.,0.};
 698     if(!SPDLoc2Glo(module,local,glochip)){
 699       return kFALSE;
 700     }
 701
 702     const TVector3 pos(glochip[0],glochip[1],glochip[2]);
 703     //pos.Print();
 704
 705     if(fabs(pos.Eta()) > etacut)
 706       return kFALSE;
 707   }
 708
 709   return kTRUE;
 710 }
 711
 712 void AliAnalysisTaskDDMeson::GetNFO(Int_t &ni, Int_t &no) const
 713 {
 714   //
 715   //GetNFO
 716   //
 717
 718   const AliMultiplicity *mult = fESD->GetMultiplicity();
 719   Int_t nFOinner=0;
 720   Int_t nFOouter=0;
 721   for(Int_t iChipKey=0; iChipKey < 1200; iChipKey++){
 722     if(mult->TestFastOrFiredChips(iChipKey) && CheckChipEta(iChipKey)){ // here you check if the FastOr bit is 1 or 0
 723       if(iChipKey<400)
 724         nFOinner++;  // here you count the FastOr bits in the inner layer
 725       else
 726         nFOouter++;  // here you count the FastOr bits in the outer layer
 727     }
 728   }
 729
 730   ni = nFOinner;
 731   no = nFOouter;
 732 }
 733
 734 Bool_t AliAnalysisTaskDDMeson::CheckBit()
 735 {
 736   //
 737   //get offline trigger
 738   //
 739   fhps->Fill(1);
 740
 741   AliTriggerAnalysis triggerAnalysis;
 742
 743   //hardware 1: hw 0: offline
 744   const Bool_t khw = kFALSE;
 745   const Bool_t v0A       = (triggerAnalysis.V0Trigger(fESD, AliTriggerAnalysis::kASide, khw) == AliTriggerAnalysis::kV0BB);
 746   const Bool_t v0C       = (triggerAnalysis.V0Trigger(fESD, AliTriggerAnalysis::kCSide, khw) == AliTriggerAnalysis::kV0BB);
 747
 748   /*
 749     const Bool_t v0ABG     = (triggerAnalysis.V0Trigger(fESD, AliTriggerAnalysis::kASide, khw) == AliTriggerAnalysis::kV0BG);
 750     const Bool_t v0CBG     = (triggerAnalysis.V0Trigger(fESD, AliTriggerAnalysis::kCSide, khw) == AliTriggerAnalysis::kV0BG);
 751     const Bool_t v0BG      = v0ABG || v0CBG;
 752     const Int_t fastOROffline = triggerAnalysis.SPDFiredChips(fESD, 0); // SPD number of chips from clusters (!)
 753
 754     //http://alisoft.cern.ch/viewvc/trunk/ANALYSIS/AliPhysicsSelection.cxx?view=markup&root=AliRoot
 755     //Bool_t mb1 = (fastOROffline > 0 || v0A || v0C) && (!v0BG);
 756    */
 757
 758   const Bool_t fmdA = triggerAnalysis.FMDTrigger(fESD, AliTriggerAnalysis::kASide);
 759   const Bool_t fmdC = triggerAnalysis.FMDTrigger(fESD, AliTriggerAnalysis::kCSide);
 760
 761   Bool_t bitA = kFALSE, bitC = kFALSE;
 762   if(fOpt.Contains("V0")){
 763     bitA = bitA || v0A;
 764     bitC = bitC || v0C;
 765   }
 766   if(fOpt.Contains("FMD")){
 767     bitA = bitA || fmdA;
 768     bitC = bitC || fmdC;
 769   }
 770
 771   //---------------------------------------------------------------------------------------------
 772   //---------------------------------------------------------------------------------------------
 773   //considering the efficiency per chip, nfo and fastORhw is required only to have non-0 count!!
 774
 775   //simulate SPD _hardware_ trigger in eta range
 776   if(fOpt.Contains("SPD")){
 777     Int_t ni=-999, no = -999;
 778     GetNFO(ni, no);
 779     if(!bitA && !bitC){
 780       fhfo->Fill(ni, no);
 781     }
 782     if(ni<2 || no<2)
 783       return 0;
 784   }
 785
 786   //simulate MB condition, since for double gap SPD is definitely fired, fastORhw>0 to reduce GA, GC, NG
 787   //http://alisoft.cern.ch/viewvc/trunk/ANALYSIS/AliTriggerAnalysis.cxx?view=markup&root=AliRoot
 788   //Int_t AliTriggerAnalysis::SPDFiredChips(const AliESDEvent* aEsd, Int_t origin, Bool_t fillHists, Int_t layer)
 789   // returns the number of fired chips in the SPD
 790   // origin = 0 --> aEsd->GetMultiplicity()->GetNumberOfFiredChips() (filled from clusters)
 791   // origin = 1 --> aEsd->GetMultiplicity()->TestFastOrFiredChips() (from hardware bits)
 792   const Int_t fastORhw = triggerAnalysis.SPDFiredChips(fESD, 1);
 793   fhspd->Fill(fastORhw);
 794   if(fastORhw<1)
 795     return 0;
 796
 797   //---------------------------------------------------------------------------------------------
 798   //---------------------------------------------------------------------------------------------
 799
 800   const Int_t iv0 = v0A + 2*v0C;
 801   const Int_t ifmd = fmdA + 2*fmdC;
 802   fhv0fmd->Fill(iv0, ifmd);
 803
 804   //---------------------------------------------------------------------------------------------
 805   //---------------------------------------------------------------------------------------------
 806
 807   //need a base number to distringuish from null entry; fHBit starts from 1; important!!
 808   fBitcg=1;
 809
 810   if(bitA)
 811     fBitcg += fCHECKVBA;
 812   if(bitC)
 813     fBitcg += fCHECKVBC;
 814
 815   fHBIT->Fill( fBitcg );
 816
 817   return 1;
 818 }
 819
 820 Int_t AliAnalysisTaskDDMeson::CutESD(const AliESDtrack *outtrk[])
 821 {
 822   //
 823   //track cuts
 824   //
 825
 826   //collision vertex cut
 827   //A cut in XY is implicitly done during the reconstruction by constraining the vertex to the beam diamond.
 828
 829   // Primary vertex
 830   Bool_t kpr0 = kTRUE;
 831   const AliESDVertex *vertex = fESD->GetPrimaryVertexTracks();
 832   if(vertex->GetNContributors()<1) {
 833     // SPD vertex
 834     vertex = fESD->GetPrimaryVertexSPD();
 835     if(vertex->GetNContributors()<1) {
 836       // NO GOOD VERTEX, SKIP EVENT
 837       kpr0 = kFALSE;
 838     }
 839   }
 840   const Bool_t kpriv = kpr0 && ( fabs(fESD->GetPrimaryVertex()->GetZ())<10 );
 841   fhpriv->Fill(kpriv);
 842   if(!kpriv)
 843     return 0;
 844
 845   AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts;//::GetStandardITSTPCTrackCuts2010(kTRUE);
 846
 847   //--->
 848   // TPC
 849   esdTrackCuts->SetMinNClustersTPC(70);
 850   esdTrackCuts->SetMaxChi2PerClusterTPC(4);
 851   esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
 852   esdTrackCuts->SetRequireTPCRefit(kTRUE);
 853   // ITS
 854   esdTrackCuts->SetRequireITSRefit(kTRUE);
 855   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 856                                          AliESDtrackCuts::kAny);
 857   // 7*(0.0026+0.0050/pt^1.01)
 858   esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
 859
 860   esdTrackCuts->SetMaxDCAToVertexZ(2);
 861   esdTrackCuts->SetDCAToVertex2D(kFALSE);
 862   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
 863   //---<
 864
 865   esdTrackCuts->SetEtaRange(-0.9, 0.9);
 866
 867   const TObjArray* seltrk = esdTrackCuts->GetAcceptedTracks(fESD);
 868   const Int_t n0=seltrk->GetEntries();
 869   const AliESDtrack * trks[n0];
 870   Int_t nsel=0;
 871   for(Int_t ii=0; ii<n0; ii++){
 872     trks[ii]=0x0;
 873     const AliESDtrack *esdtrack = (AliESDtrack *)(seltrk->At(ii));
 874
 875     /*
 876     if(!CutTrack(esdtrack))
 877       continue;
 878     */
 879
 880     trks[nsel++]=esdtrack;
 881   }
 882
 883   delete seltrk;
 884   delete esdTrackCuts;
 885
 886   fv0ntrk->Fill(nsel,GetV0());
 887
 888   fhntrk->Fill(nsel);
 889
 890   frsntrk->Fill(fESD->GetNumberOfTracks(), nsel);
 891
 892   if(nsel>=fnsel1)
 893     return 0;
 894
 895   for(Int_t ii=0; ii<nsel; ii++){
 896     outtrk[ii]=trks[ii];
 897   }
 898
 899   return nsel;
 900 }
 901
 902 /*
 903 Bool_t AliAnalysisTaskDDMeson::CutTrack(const AliESDtrack *  esdtrack) const
 904 {
 905   //
 906   //track cut
 907   //
 908
 909   return kTRUE;
 910 }
 911 */
 912
 913 //=================================================================================
 914
 915 void AliAnalysisTaskDDMeson::SwapTrack(const AliESDtrack * trks[]) const
 916 {
 917   //
 918   //swap two esdtracks
 919   //
 920   const AliESDtrack * tmpt = trks[0];
 921   trks[0]=trks[1];
 922   trks[1]=tmpt;
 923 }
 924
 925 Int_t AliAnalysisTaskDDMeson::GetV0() const
 926 {
 927   //
 928   //V0 bit combination
 929   //
 930
 931   const Bool_t ka = (fBitcg & fCHECKVBA);
 932   const Bool_t kc = (fBitcg & fCHECKVBC);
 933
 934
 935   //a1c0 a0c1 must be adjacent so that a cut for signal gas can be eaily applied
 936   const Int_t a0c0 = 1;
 937   const Int_t a1c0 = 2;
 938   const Int_t a0c1 = 3;
 939   const Int_t a1c1 = 4;
 940
 941   if(!ka && !kc)
 942     return a0c0;
 943   else{
 944     if(ka && kc)
 945       return a1c1;
 946     else if(ka)
 947       return a1c0;
 948     else
 949       return a0c1;
 950   }
 951 }
 952
 953 Int_t AliAnalysisTaskDDMeson::GetCombCh(const Double_t s1, const Double_t s2) const
 954 {
 955   //
 956   //get combination of charges
 957   //
 958   const Int_t combPP = 1;
 959   const Int_t combMM = 2;
 960   const Int_t combPM = 3;
 961
 962   if(s1*s2<0){
 963     return combPM;
 964   }
 965   else if(s1>0)
 966     return combPP;
 967   else
 968     return combMM;
 969 }
 970
 971 //==========================================================================
 972
 973 TLorentzVector AliAnalysisTaskDDMeson::GetKinematics(const Double_t *pa, const Double_t *pb, const Double_t ma, const Double_t mb, Double_t & cts) const
 974 {
 975   //
 976   //get kinematics, cts = cos(theta_{#star})
 977   //
 978   TLorentzVector va, vb, sumv;
 979   va.SetXYZM(pa[0], pa[1], pa[2], ma);
 980   vb.SetXYZM(pb[0], pb[1], pb[2], mb);
 981   sumv = va+vb;
 982
 983   const TVector3 bv = -sumv.BoostVector();
 984
 985   va.Boost(bv);
 986   vb.Boost(bv);
 987   //printf("test a %f %f %f -- %f %f\n", va.X(), va.Y(), va.Z(), va.Theta(), va.Phi());
 988   //printf("test b %f %f %f -- %f %f\n", vb.X(), vb.Y(), vb.Z(), vb.Theta(), vb.Phi());
 989
 990   cts = fabs(cos(va.Theta()));
 991
 992   return sumv;
 993 }
 994
 995 Double_t AliAnalysisTaskDDMeson::GetCtlab(const Double_t *pa, const Double_t *pb) const
 996 {
 997   //
 998   //cos(theat_lab)
 999   //
1000
1001   TVector3 va, vb;
1002   va.SetXYZ(pa[0], pa[1], pa[2]);
1003   vb.SetXYZ(pb[0], pb[1], pb[2]);
1004
1005   const TVector3 ua = va.Unit();
1006   const TVector3 ub = vb.Unit();
1007
1008   return ua.Dot(ub);
1009 }