PWGUD/DIFFRACTIVE/xsAndTwoProng/AliCDMesonUtils.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 // AliCDMesonUtils
  17 // for
  18 // AliAnalysisTaskCDMeson
  19 //
  20 //  Author:
  21 //  Xianguo Lu <lu@physi.uni-heidelberg.de>
  22 //  continued by
  23 //  Felix Reidt <Felix.Reidt@cern.ch>
  24
  25 #include <TH1.h>
  26 #include <TH2.h>
  27 #include <TGeoMatrix.h>
  28 #include <THnSparse.h>
  29 #include <TLorentzVector.h>
  30 #include <TRandom3.h>
  31 #include <TParticle.h>
  32 #include <TObjArray.h>
  33
  34 #include "AliCDBManager.h"
  35 #include "AliCDBEntry.h"
  36 #include "AliCDBStorage.h"
  37 #include "AliESDEvent.h"
  38 #include "AliPIDResponse.h"
  39 #include "AliVTrack.h"
  40 #include "AliVParticle.h"
  41 #include "AliESDtrack.h"
  42 #include "AliESDtrackCuts.h"
  43 #include "AliESDFMD.h"
  44 #include "AliESDVZERO.h"
  45 #include "AliGeomManager.h"
  46 #include "AliITSAlignMille2Module.h"
  47 #include "AliITSsegmentationSPD.h"
  48 #include "AliMultiplicity.h"
  49 #include "AliPIDResponse.h"
  50 #include "AliSPDUtils.h"
  51 #include "AliTriggerAnalysis.h"
  52 #include "AliVZEROTriggerData.h"
  53 #include "AliVZEROCalibData.h"
  54
  55 #include "AliAODTracklets.h"
  56 #include "AliAODEvent.h"
  57
  58 #include "AliCDMesonBase.h"
  59 #include "AliCDMesonTracks.h"
  60
  61 #include "AliCDMesonUtils.h"
  62
  63
  64 //==============================================================================
  65 //------------------------------------------------------------------------------
  66 void AliCDMesonUtils::GetMassPtCtsOA(Int_t pid,
  67                                      const TVector3* momenta[], Double_t & mass,
  68                                      Double_t &pt, Double_t &cts, Double_t &oa)
  69 {
  70         //
  71         // Get Mass Pt Cts OA
  72         //
  73
  74         Double_t tmpp1[3], tmpp2[3];
  75         momenta[0]->GetXYZ(tmpp1);
  76         momenta[1]->GetXYZ(tmpp2);
  77
  78         Double_t masstrk = -999;
  79         if(pid==AliCDMesonBase::kBinPionE || pid==AliCDMesonBase::kBinPion ||
  80            pid==AliCDMesonBase::kBinSinglePion)
  81                 masstrk = AliPID::ParticleMass(AliPID::kPion);
  82         else if(pid==AliCDMesonBase::kBinKaonE || pid==AliCDMesonBase::kBinKaon ||
  83                 pid==AliCDMesonBase::kBinSingleKaon)
  84                 masstrk = AliPID::ParticleMass(AliPID::kKaon);
  85         else if(pid==AliCDMesonBase::kBinProtonE || pid==AliCDMesonBase::kBinProton ||
  86                 pid==AliCDMesonBase::kBinSingleProton)
  87                 masstrk = AliPID::ParticleMass(AliPID::kProton);
  88         else if(pid==AliCDMesonBase::kBinElectronE ||
  89                 pid==AliCDMesonBase::kBinElectron ||
  90                 pid==AliCDMesonBase::kBinSingleElectron)
  91                 masstrk = AliPID::ParticleMass(AliPID::kElectron);
  92         else
  93                 masstrk = AliPID::ParticleMass(AliPID::kPion);
  94
  95         const TLorentzVector sumv = GetKinematics(tmpp1, tmpp2, masstrk, masstrk,
  96                                                   cts);
  97         mass = sumv.M();
  98         pt = sumv.Pt();
  99
 100         oa = GetOA(tmpp1, tmpp2);
 101 }
 102
 103
 104 //------------------------------------------------------------------------------
 105 void AliCDMesonUtils::GetPWAinfo(Int_t pid, const AliVTrack *trks[],
 106                                  Float_t& theta, Float_t& phi, Float_t& mass,
 107                                  Float_t momentum[])
 108 {
 109         // transforms the coordiante system to the helicity frame and determines
 110         // invariant mass and 3-vector of the two track system
 111         //
 112
 113         Double_t tmpp1[3], tmpp2[3]; // 3-vectors of the daughter tracks
 114         trks[0]->GetPxPyPz(tmpp1);
 115         trks[1]->GetPxPyPz(tmpp2);
 116
 117         // determine mass of the daugther tracks
 118         Double_t masstrk = -999;
 119         if(pid==AliCDMesonBase::kBinPion || pid==AliCDMesonBase::kBinSinglePion)
 120                 masstrk = AliPID::ParticleMass(AliPID::kPion);
 121         else if(pid==AliCDMesonBase::kBinKaon || pid==AliCDMesonBase::kBinSingleKaon)
 122                 masstrk = AliPID::ParticleMass(AliPID::kKaon);
 123         else if(pid==AliCDMesonBase::kBinProton ||
 124                 pid==AliCDMesonBase::kBinSingleProton)
 125                 masstrk = AliPID::ParticleMass(AliPID::kProton);
 126         else if(pid==AliCDMesonBase::kBinElectron ||
 127                 pid==AliCDMesonBase::kBinSingleElectron)
 128                 masstrk = AliPID::ParticleMass(AliPID::kElectron);
 129         else
 130                 masstrk = AliPID::ParticleMass(AliPID::kPion);
 131
 132         TLorentzVector va, vb; // 4-vectors of the two tracks
 133         va.SetXYZM(tmpp1[0], tmpp1[1], tmpp1[2], masstrk);
 134         vb.SetXYZM(tmpp2[0], tmpp2[1], tmpp2[2], masstrk);
 135
 136         TLorentzVector sumv = va+vb; // 4-vector of the pair
 137         TVector3 sumXZ(sumv.X(), 0, sumv.Z()); // its projection to the xz-plane
 138
 139         // mass and momentum
 140         mass = sumv.M();
 141         momentum[0] = sumv.X();
 142         momentum[1] = sumv.Y();
 143         momentum[2] = sumv.Z();
 144
 145         // coordinate axis vectors
 146         const TVector3 x(1.,0,0);
 147         const TVector3 y(0,1.,0);
 148         const TVector3 z(0,0,1.);
 149
 150         const Double_t alpha = -z.Angle(sumXZ);
 151
 152         sumXZ.Rotate(alpha, y);
 153         sumv.Rotate(alpha, y);
 154         va.Rotate(alpha, y);
 155         vb.Rotate(alpha, y);
 156
 157         const Double_t beta = z.Angle(sumv.Vect());
 158         sumv.Rotate(beta, x);
 159         va.Rotate(beta, x);
 160         vb.Rotate(beta, x);
 161
 162
 163         const TVector3 bv = -sumv.BoostVector();
 164         va.Boost(bv);
 165         vb.Boost(bv);
 166
 167         // angles
 168         theta = va.Theta();
 169         phi = va.Phi();
 170 }
 171
 172 //------------------------------------------------------------------------------
 173 Int_t AliCDMesonUtils::GetEventType(const AliESDEvent *ESDEvent)
 174 {
 175         // checks of which type a event is:
 176         // beam-beam interaction (I), beam-gas (A/C), empty (E)
 177         //
 178
 179         TString firedTriggerClasses = ESDEvent->GetFiredTriggerClasses();
 180
 181         if (firedTriggerClasses.Contains("CINT1A-ABCE-NOPF-ALL")) { // A
 182                 return AliCDMesonBase::kBinEventA;
 183         }
 184         if (firedTriggerClasses.Contains("CINT1C-ABCE-NOPF-ALL")) { // C
 185                 return AliCDMesonBase::kBinEventC;
 186         }
 187         if (firedTriggerClasses.Contains("CINT1B-ABCE-NOPF-ALL")) { // I
 188                 return AliCDMesonBase::kBinEventI;
 189         }
 190         if (firedTriggerClasses.Contains("CINT1-E-NOPF-ALL")) { // E
 191                 return AliCDMesonBase::kBinEventE;
 192         }
 193         if (firedTriggerClasses.Contains("CDG5-E")) { // E
 194                 return AliCDMesonBase::kBinEventE;
 195         }
 196         if (firedTriggerClasses.Contains("CDG5-I")) { // I
 197                 return AliCDMesonBase::kBinEventI;
 198         }
 199         if (firedTriggerClasses.Contains("CDG5-AC")) { // AC
 200                 return AliCDMesonBase::kBinEventAC;
 201         }
 202         return AliCDMesonBase::kBinEventUnknown;
 203 }
 204
 205
 206 //------------------------------------------------------------------------------
 207 void AliCDMesonUtils::SwapTrack(const AliVTrack* trks[])
 208 {
 209         //
 210         // swap two esdtracks, needed since only the properties of one a stored and
 211         // AliRoot sorts them
 212         //
 213
 214         TRandom3 tmprd(0);
 215         if(tmprd.Rndm()>0.5)
 216                 return;
 217
 218         const AliVTrack* tmpt = trks[0];
 219         trks[0]=trks[1];
 220         trks[1]=tmpt;
 221 }
 222
 223
 224 //------------------------------------------------------------------------------
 225 Int_t AliCDMesonUtils::GetCombCh(const AliVTrack * trks[])
 226 {
 227         //
 228         // get combination of charges
 229         //
 230
 231         const Double_t ch1 = trks[0]->Charge();
 232         const Double_t ch2 = trks[1]->Charge();
 233
 234         if(ch1*ch2<0){
 235                 return AliCDMesonBase::kBinPM;
 236         }
 237         else
 238                 return AliCDMesonBase::kBinPPMM;
 239 }
 240
 241
 242 //------------------------------------------------------------------------------
 243 Int_t AliCDMesonUtils::GetCombPID(AliPIDResponse* pid, const AliVTrack* trks[],
 244                                   Int_t mode, TH2* comb2trkPID /*= 0x0 */)
 245 {
 246         //
 247         // Get PID for a two track system (data)
 248         //
 249
 250         if (!pid){
 251                 return AliCDMesonBase::kBinPIDUnknown;
 252         }
 253
 254         // get PID for the single tracks
 255         Int_t kpid[2];
 256         for(Int_t ii=0; ii<2; ii++){
 257                 kpid[ii] = GetPID(pid, trks[ii], mode);
 258         }
 259
 260         // PID QA histogram
 261         if (comb2trkPID) {
 262                 comb2trkPID->Fill(kpid[0], kpid[1]);
 263         }
 264
 265         return CombinePID(kpid);
 266 }
 267
 268
 269 //------------------------------------------------------------------------------
 270 Int_t AliCDMesonUtils::GetCombPID(const TParticle* particles[],
 271                                   TH2 *comb2trkPID /* = 0x0 */)
 272 {
 273         //
 274         // Get PID for a two track system (MC)
 275         //
 276
 277         // get PID for the single tracks
 278         Int_t kpid[2];
 279         for(Int_t ii=0; ii<2; ii++){
 280                 kpid[ii] = GetPID(particles[ii]->GetPdgCode());
 281         }
 282
 283         // PID QA histogram
 284         if (comb2trkPID) {
 285                 comb2trkPID->Fill(kpid[0], kpid[1]);
 286         }
 287
 288         return CombinePID(kpid);
 289 }
 290
 291
 292 //------------------------------------------------------------------------------
 293 void AliCDMesonUtils::GetSPDTrackletMult(const AliESDEvent *ESDEvent,
 294                                          Int_t& sum, Int_t& forwardA,
 295                                          Int_t& forwardC, Int_t& central)
 296 {
 297         // obtain the multiplicity for eta < -.9 (forwardC), eta > 0.9 (fowardA), in
 298         // the central barrel (central) and its sum from SPD tracklets with the
 299         // AliMultiplicity class
 300         // code simular to PWGPP/ITS/AliAnalysisTaskSPD.cxx
 301
 302         sum = forwardA = forwardC = central = 0; // initialize values
 303
 304         const AliMultiplicity *mult = ESDEvent->GetMultiplicity();
 305         for (Int_t iTracklet = 0; iTracklet < mult->GetNumberOfTracklets();
 306              iTracklet++) {
 307                 float_t eta = mult->GetEta(iTracklet);
 308                 if (eta < -0.9) {
 309                         forwardC++;
 310                 }
 311                 else if (eta < 0.9) {
 312                         central++;
 313                 }
 314                 else {
 315                         forwardA++;
 316                 }
 317                 sum++;
 318         }
 319 }
 320
 321
 322 //------------------------------------------------------------------------------
 323 Bool_t AliCDMesonUtils::CutEvent(const AliESDEvent *ESDEvent, TH1 *hspd,
 324                                  TH1 *hpriv, TH1 *hpriVtxPos, TH1 *hpriVtxDist,
 325                                  TH2 *hfo, TH1* hfochans, Int_t &kfo,
 326                                  Int_t &nip, Int_t &nop, TH1 *hpriVtxX,
 327                                  TH1 *hpriVtxY, TH1 *hpriVtxZ)
 328 {
 329         //
 330         // CutEvent
 331         //
 332
 333         AliTriggerAnalysis triggerAnalysis;
 334         /*
 335         //printf("AliCDMesonUtils active triggers: %s\n",
 336         //       ESDEvent->GetHeader()->GetActiveTriggerInputs().Data());
 337         //printf("AliCDMesonUtils fired triggers: %s\n",
 338         //       ESDEvent->GetHeader()->GetFiredTriggerInputs().Data());
 339         //*/
 340         //http://alisoft.cern.ch/viewvc/trunk/ANALYSIS/AliTriggerAnalysis.cxx?view=markup&root=AliRoot
 341         //Int_t AliTriggerAnalysis::SPDFiredChips(const AliESDEvent* aEsd, Int_t origin, Bool_t fillHists, Int_t layer)
 342         // returns the number of fired chips in the SPD
 343         // origin = 0 --> aEsd->GetMultiplicity()->GetNumberOfFiredChips() (filled from clusters)
 344         // origin = 1 --> aEsd->GetMultiplicity()->TestFastOrFiredChips() (from hardware bits)
 345
 346
 347         // SPD FastOR cut
 348         const Int_t fastORhw = triggerAnalysis.SPDFiredChips(ESDEvent, 1);
 349         if (hspd) hspd->Fill(fastORhw);
 350         /*
 351         if(fastORhw<1)
 352                 return kFALSE;
 353         */
 354
 355         if (hfochans) {
 356                 const AliMultiplicity *mult = ESDEvent->GetMultiplicity();
 357
 358                 for(Int_t iChipKey=0; iChipKey < 1200; iChipKey++){
 359                         if (mult->TestFastOrFiredChips(iChipKey)) {
 360                                 hfochans->Fill((Double_t)iChipKey);
 361                         }
 362                 }
 363         }
 364
 365         if(kfo){ // spd trigger study
 366                 Int_t nfoctr[10];
 367                 GetNFO(ESDEvent, "[1.0]", nfoctr, 0x0, 0x0);
 368                 nip = nfoctr[kInnerPixel];
 369                 nop = nfoctr[kOuterPixel];
 370
 371                 if(AliCDMesonBase::GetGapBin("V0", GetV0(ESDEvent)) ==
 372                     AliCDMesonBase::kBinDG) {
 373                         hfo->Fill(nip, nop);
 374                 }
 375
 376                 if(nop<2)
 377                         kfo = 1;
 378                 else
 379                         kfo = nop;
 380
 381                 if(kfo>=10)
 382                         kfo=9;
 383         }
 384
 385         // collision vertex cut
 386         // A cut in XY is implicitly done during the reconstruction by constraining
 387         // the vertex to the beam diamond.
 388
 389         // Primary vertex
 390         Bool_t kpr0 = kTRUE;
 391         const AliESDVertex *vertex = ESDEvent->GetPrimaryVertexTracks();
 392         if(vertex->GetNContributors()<1) {
 393                 // SPD vertex
 394                 vertex = ESDEvent->GetPrimaryVertexSPD();
 395                 if(vertex->GetNContributors()<1) {
 396                         // NO VERTEX, SKIP EVENT
 397                         kpr0 = kFALSE;
 398                 }
 399         }
 400         const Bool_t kpriv = kpr0 && (fabs(vertex->GetZ()) < 10.);
 401         // 10 is the common value, unit: cm
 402         if (hpriv) hpriv->Fill(kpriv);
 403         if (hpriVtxDist) hpriVtxDist->Fill(vertex->GetZ());
 404         if (hpriVtxPos && kpr0) hpriVtxPos->Fill(!kpriv);
 405         if(!kpriv)
 406                 return kFALSE;
 407
 408         if(hpriVtxX) hpriVtxX->Fill(vertex->GetX());
 409         if(hpriVtxY) hpriVtxY->Fill(vertex->GetY());
 410         if(hpriVtxZ) hpriVtxZ->Fill(vertex->GetZ());
 411         return kTRUE;
 412 }
 413
 414 //------------------------------------------------------------------------------
 415 Bool_t AliCDMesonUtils::CutEvent(const AliAODEvent *AODEvent, TH1 *hpriv,
 416                                  TH1 *hpriVtxX, TH1 *hpriVtxY, TH1 *hpriVtxZ,
 417                                  TH1 *hpriVtxPos, TH1 *hpriVtxDist)
 418 {
 419         //
 420         // Cut Event for AOD Events, to be combined with the ESD Track Cut
 421         //
 422
 423         // TODO: no idea about fast or yet, to be thought of
 424
 425         // Primary vertex
 426         Bool_t kpr0 = kTRUE;
 427         const AliAODVertex *vertex = AODEvent->GetPrimaryVertex();
 428         if(vertex->GetNContributors()<1) {
 429                 // SPD vertex
 430                 vertex = AODEvent->GetPrimaryVertexSPD();
 431                 if(vertex->GetNContributors()<1) {
 432                         // NO GOOD VERTEX, SKIP EVENT
 433                         kpr0 = kFALSE;
 434                 }
 435         }
 436         const Bool_t kpriv = kpr0 && (fabs(vertex->GetZ())<10.);
 437         // 10 is the common value, unit: cm
 438         hpriv->Fill(kpriv);
 439         if (hpriVtxDist) hpriVtxDist->Fill(vertex->GetZ());
 440         if (hpriVtxPos && kpr0) hpriVtxPos->Fill(!kpriv);
 441         if(!kpriv)
 442                 return kFALSE;
 443
 444         if(hpriVtxX) hpriVtxX->Fill(vertex->GetX());
 445         if(hpriVtxY) hpriVtxY->Fill(vertex->GetY());
 446         if(hpriVtxZ) hpriVtxZ->Fill(vertex->GetZ());
 447         return kTRUE;
 448 }
 449
 450
 451 //------------------------------------------------------------------------------
 452 void AliCDMesonUtils::DoVZEROStudy(const AliESDEvent *ESDEvent,
 453                                    TObjArray* hists, Int_t run)
 454 {
 455         //
 456         //
 457         // IMPORTANT: the order of the histograms here and in
 458         // AliCDMesonBase::GetHistVZEROStudies(..) has to match
 459
 460         const AliESDVZERO* esdV0 = ESDEvent->GetVZEROData();
 461         if (!esdV0) {
 462                 Printf("ERROR: esd V0  not available");
 463                 return;
 464         }
 465
 466         // determine trigger decision
 467         AliTriggerAnalysis triggerAnalysis;
 468         //const Bool_t khw = kFALSE;
 469
 470         //Float_t v0a = (triggerAnalysis.V0Trigger(ESDEvent, AliTriggerAnalysis::kASide, khw) ==
 471         //             AliTriggerAnalysis::kV0BB) ? 1. : 0.;
 472         //Float_t v0c = (triggerAnalysis.V0Trigger(ESDEvent, AliTriggerAnalysis::kCSide, khw) ==
 473         //             AliTriggerAnalysis::kV0BB) ? 1. : 0.;
 474
 475         // obtain OCDB objects
 476         AliCDBManager *man = AliCDBManager::Instance();
 477         TString cdbpath;
 478         if (man->IsDefaultStorageSet()) {
 479                 const AliCDBStorage *dsto = man->GetDefaultStorage();
 480                 cdbpath = TString(dsto->GetBaseFolder());
 481         }
 482         else { //should not be used!
 483                 man->SetDefaultStorage(gSystem->Getenv("TRAIN_CDB_PATH"));
 484                 cdbpath = TString(gSystem->Getenv("TRAIN_CDB_PATH"));
 485         }
 486         man->SetSpecificStorage("VZERO/Trigger/Data",cdbpath);
 487         man->SetSpecificStorage("VZERO/Calib/Data",cdbpath);
 488         man->SetRun(run);
 489
 490         AliCDBEntry *ent1 = man->Get("VZERO/Trigger/Data");
 491         if (!ent1) {
 492                 printf("AliCDMesonUtils failed loading VZERO trigger entry from OCDB\n");
 493                 return;
 494         }
 495         AliVZEROTriggerData *trigData = (AliVZEROTriggerData*)ent1->GetObject();
 496         if (!trigData) {
 497                 printf("AliCDMesonUtils failed loading VZERO trigger data from OCDB\n");
 498                 return;
 499         }
 500
 501         AliCDBEntry *ent2 = man->Get("VZERO/Calib/Data");
 502         if (!ent2) {
 503                 printf("AliCDMesonUtils failed loading VZERO calib entry from OCDB\n");
 504                 return;
 505         }
 506         AliVZEROCalibData *calData = (AliVZEROCalibData*)ent2->GetObject();
 507         if (!calData) {
 508                 printf("AliCDMesonUtils failed loading VZERO calibration data from OCDB\n");
 509                 return;
 510         }
 511
 512         // fill histograms
 513         Int_t pmtHist = 0;
 514         Int_t iPMT = 0;
 515         for (iPMT = 0; iPMT < 64; ++iPMT) {
 516                 Int_t board   = AliVZEROCalibData::GetBoardNumber(iPMT);
 517                 Int_t channel = AliVZEROCalibData::GetFEEChannelNumber(iPMT);
 518                 ((TH2*)hists->At(iPMT+pmtHist))->Fill(esdV0->GetAdc(iPMT),
 519                                                       trigData->GetPedestalCut(0, board, channel));
 520         //                                      calData->GetDiscriThr(iPMT));
 521         }
 522         pmtHist = iPMT;
 523         for (iPMT = 0; iPMT < 64; ++iPMT) {
 524                 ((TH2*)hists->At(iPMT+pmtHist))->Fill(esdV0->GetAdc(iPMT),
 525                                                       esdV0->GetMultiplicity(iPMT));
 526         }
 527         /*
 528         // not used in 2010 for pp
 529         ((TH2*)hists->At(pmtHist++))->Fill((Float_t)esdV0->GetTriggerChargeA(), v0a);
 530         ((TH2*)hists->At(pmtHist++))->Fill((Float_t)esdV0->GetTriggerChargeC(), v0c);
 531         */
 532 }
 533
 534
 535 //------------------------------------------------------------------------------
 536 Int_t AliCDMesonUtils::GetGapConfig(const AliESDEvent *ESDEvent,
 537                                     TH2 *hitMapSPDinner, TH2 *hitMapSPDouter,
 538                                     TH2 *hitMapSPDtrklt, TH2 *hitMapFMDa,
 539                                     TH2 *hitMapFMDc, TH1 **fmdSums,
 540                                     TH2 *TPCGapDCAaSide, TH2 *TPCGapDCAcSide)
 541 {
 542         //
 543         // GetGapConfigAndTracks
 544         //
 545         // retrieves the gap configuration of a track and returns it as
 546         // an bit vector
 547         // kBaseLine ensures, that this event is valid
 548         // + is equivalent to | in this case
 549         Int_t gapConfig = AliCDMesonBase::kBitBaseLine + GetV0(ESDEvent)
 550                 + GetFMD(ESDEvent, hitMapFMDa, hitMapFMDc, fmdSums)
 551                 + GetSPD(ESDEvent, hitMapSPDinner, hitMapSPDouter, hitMapSPDtrklt);
 552         if (gapConfig == AliCDMesonBase::kBitBaseLine) {
 553                 gapConfig += GetTPC(ESDEvent, TPCGapDCAaSide, TPCGapDCAcSide);
 554         }
 555         else {
 556                 gapConfig += GetTPC(ESDEvent, 0x0, 0x0);
 557         }
 558         if (GetFastORmultiplicity(ESDEvent) > 0) {
 559                 gapConfig += AliCDMesonBase::kBitCentAct;
 560         }
 561         return gapConfig; // + GetZDC(ESDEvent);
 562 }
 563
 564
 565 //------------------------------------------------------------------------------
 566 void AliCDMesonUtils::FillEtaPhiMap(const AliVEvent *event,
 567                                     const AliCDMesonTracks* tracks, TH2 *map,
 568                                     TH2 *map_c)
 569 {
 570         //
 571         // Fills the eta phi information about all tracks and about the tracks surving
 572         // the track cuts (CutTrack) into two separate histograms
 573         //
 574
 575         // all tracks
 576         for (Int_t itrk = 0; itrk < event->GetNumberOfTracks(); ++itrk) {
 577                 if (AliVParticle *trk = event->GetTrack(itrk)) {
 578                         if (map) {
 579                                 map->Fill(trk->Eta(), trk->Phi());
 580                         }
 581                 }
 582         }
 583
 584         // tracks that survived the cuts
 585         for (Int_t itrk = 0; itrk < tracks->GetCombinedTracks(); ++itrk) {
 586                 if (map_c) {
 587                         map_c->Fill(tracks->GetTrack(itrk)->Eta(), tracks->GetTrack(itrk)->Phi());
 588                 }
 589         }
 590 }
 591
 592
 593 //------------------------------------------------------------------------------
 594 void AliCDMesonUtils::GetMultFMD(const AliESDEvent *ESDEvent, Int_t& fmdA,
 595                                  Int_t& fmdC, Float_t  *fmdSums)
 596 {
 597         //
 598         // Multiplicity seen by FMD
 599         //
 600         // WARNING: this function is only working with a modified AliRoot so far
 601
 602 #ifdef STD_ALIROOT
 603         fmdA = FMDHitCombinations(ESDEvent, 0);
 604         fmdC = FMDHitCombinations(ESDEvent, 1);
 605 #else
 606         AliTriggerAnalysis triggerAnalysis;
 607         triggerAnalysis.SetFMDThreshold(0.3, 0.5); // parameters got from FMD
 608         triggerAnalysis.FMDTrigger(ESDEvent, AliTriggerAnalysis::kASide, &fmdA);
 609         triggerAnalysis.FMDTrigger(ESDEvent, AliTriggerAnalysis::kCSide, &fmdC);
 610 #endif
 611
 612         if (fmdSums) {
 613                 const AliESDFMD* fmdData =
 614                         (const_cast<AliESDEvent*>(ESDEvent))->GetFMDData();
 615                 for (UShort_t det=1; det<=3; det++) {
 616                         Int_t nRings = (det == 1 ? 1 : 2);
 617                         for (UShort_t ir = 0; ir < nRings; ir++) {
 618                                 Char_t   ring = (ir == 0 ? 'I' : 'O');
 619                                 UShort_t nsec = (ir == 0 ? 20  : 40);
 620                                 UShort_t nstr = (ir == 0 ? 512 : 256);
 621                                 for (UShort_t sec =0; sec < nsec;  sec++) {
 622                                         for (UShort_t strip = 0; strip < nstr; strip++) {
 623                                                 Float_t mult = fmdData->Multiplicity(det,ring,sec,strip);
 624                                                 if (mult == AliESDFMD::kInvalidMult) continue;
 625
 626                                                 if (det == 1 && ring == 'I') fmdSums[0] += mult;
 627                                                 else if (det == 2 && ring == 'I') fmdSums[1] += mult;
 628                                                 else if (det == 2 && ring == 'O') fmdSums[2] += mult;
 629                                                 else if (det == 3 && ring == 'I') fmdSums[3] += mult;
 630                                                 else if (det == 3 && ring == 'O') fmdSums[4] += mult;
 631                                         }
 632                                 }
 633                         }
 634                 }
 635         }
 636 }
 637
 638
 639 //------------------------------------------------------------------------------
 640 void AliCDMesonUtils::GetMultSPD(const AliESDEvent * ESDEvent, Int_t& spdIA,
 641                                  Int_t& spdIC, Int_t& spdOA, Int_t& spdOC)
 642 {
 643         //
 644         // Retrieves the multiplicity seen by the SPD FastOR
 645         //
 646
 647         Int_t nfoctr[10];
 648         GetNFO(ESDEvent, "]0.9[", nfoctr, 0x0, 0x0);
 649
 650         spdIA = nfoctr[kIPA];
 651         spdIC = nfoctr[kIPC];
 652         spdOA = nfoctr[kOPA];
 653         spdOC = nfoctr[kOPC];
 654 }
 655
 656
 657 //==============================================================================
 658 //------------------------------------------------------------------------------
 659 Int_t AliCDMesonUtils::GetV0(const AliESDEvent * ESDEvent)
 660 {
 661         //
 662         //GetV0
 663         //
 664
 665         AliTriggerAnalysis triggerAnalysis;
 666         const Bool_t khw = kFALSE;
 667         const Bool_t v0A =
 668                 (triggerAnalysis.V0Trigger(ESDEvent, AliTriggerAnalysis::kASide, khw) ==
 669                  AliTriggerAnalysis::kV0BB);
 670         const Bool_t v0C =
 671                 (triggerAnalysis.V0Trigger(ESDEvent, AliTriggerAnalysis::kCSide, khw) ==
 672                  AliTriggerAnalysis::kV0BB);
 673
 674         return v0A * AliCDMesonBase::kBitV0A + v0C * AliCDMesonBase::kBitV0C;
 675 }
 676
 677
 678 //------------------------------------------------------------------------------
 679 Int_t AliCDMesonUtils::GetFMD(const AliESDEvent *ESDEvent, TH2 *hitMapFMDa,
 680                               TH2 *hitMapFMDc, TH1 **fmdSums)
 681 {
 682         //
 683         // GetFMD
 684         //
 685
 686         AliTriggerAnalysis triggerAnalysis;
 687         triggerAnalysis.SetFMDThreshold(0.3, 0.5); // parameters got from FMD
 688         const Bool_t fmdA =
 689                 triggerAnalysis.FMDTrigger(ESDEvent, AliTriggerAnalysis::kASide);
 690         const Bool_t fmdC =
 691                 triggerAnalysis.FMDTrigger(ESDEvent, AliTriggerAnalysis::kCSide);
 692
 693         //printf("FR - GetFMD\n");
 694
 695         // prepartions for a charge summation algorithm
 696         Bool_t hitMaps = (Bool_t)(hitMapFMDa && hitMapFMDc);
 697         Bool_t calcSum = fmdSums ?
 698                 (fmdSums[0] && fmdSums[1] && fmdSums[2] && fmdSums[3] && fmdSums[4]) :
 699                 kFALSE; // are the histograms defined?
 700         Float_t sum[] = { 0., 0., 0., 0., 0. };
 701         // summed multiplicity in the FMD detectors
 702
 703         // hit map generation
 704         if (hitMaps || calcSum) {
 705                 const AliESDFMD* fmdData =
 706                         (const_cast<AliESDEvent*>(ESDEvent))->GetFMDData();
 707
 708                 for (UShort_t det=1; det<=3;det++) {
 709                         Int_t nRings = (det == 1 ? 1 : 2);
 710                         for (UShort_t ir = 0; ir < nRings; ir++) {
 711                                 Char_t   ring = (ir == 0 ? 'I' : 'O');
 712                                 UShort_t nsec = (ir == 0 ? 20  : 40);
 713                                 UShort_t nstr = (ir == 0 ? 512 : 256);
 714                                 for (UShort_t sec =0; sec < nsec;  sec++) {
 715                                         for (UShort_t strip = 0; strip < nstr; strip++) {
 716                                                 Float_t mult = fmdData->Multiplicity(det,ring,sec,strip);
 717                                                 if (mult == AliESDFMD::kInvalidMult) continue;
 718
 719                                                 if (calcSum) {
 720                                                         if (det == 1 && ring == 'I') sum[0] += mult;
 721                                                         else if (det == 2 && ring == 'I') sum[1] += mult;
 722                                                         else if (det == 2 && ring == 'O') sum[2] += mult;
 723                                                         else if (det == 3 && ring == 'I') sum[3] += mult;
 724                                                         else if (det == 3 && ring == 'O') sum[4] += mult;
 725                                                 }
 726
 727                                                 if (hitMaps) { // care about hit map specific information
 728                                                         const Float_t eta = fmdData->Eta(det,ring,sec,strip);
 729                                                         const Float_t phi =
 730                                                                 fmdData->Phi(det,ring,sec,strip) / 180. * TMath::Pi();
 731                                                         //printf("FR - GetFMD: %f %f %f\n", eta, phi, mult);
 732                                                         if (eta != AliESDFMD::kInvalidEta) {
 733                                                                 if ((-3.5 < eta) && (eta < -1.5)) {
 734                                                                         hitMapFMDc->Fill(eta, phi, mult); // use mult as weight
 735                                                                 }
 736                                                                 else if ((1.5 < eta) && (eta < 5.5)) {
 737                                                                         hitMapFMDa->Fill(eta, phi, mult); // use mult as weight
 738                                                                 }
 739                                                         }
 740                                                 }
 741                                         }
 742                                 }
 743                         }
 744                 }
 745         }
 746
 747         if (calcSum) {
 748                 //printf("DEBUG -- SUM(%f,%f,%f,%f,%f)\n", sum[0], sum[1], sum[2], sum[3],
 749                 //       sum[4]);
 750                 for (UInt_t i = 0; i < 5; i++) { //
 751                         fmdSums[i]->Fill(sum[i]);
 752                 }
 753         }
 754
 755         return fmdA * AliCDMesonBase::kBitFMDA + fmdC * AliCDMesonBase::kBitFMDC;
 756 }
 757
 758
 759 //------------------------------------------------------------------------------
 760 Int_t AliCDMesonUtils::GetSPD(const AliESDEvent *ESDEvent, TH2 *hitMapSPDinner,
 761                               TH2 *hitMapSPDouter, TH2 *hitMapSPDtrklt)
 762 {
 763         //
 764         // GetSPD
 765         //
 766
 767         Int_t nfoctr[10];
 768         GetNFO(ESDEvent, "]0.9[", nfoctr, hitMapSPDinner, hitMapSPDouter);
 769         // get multiplicity from fastOR and fill corresponding hit maps
 770
 771         if (hitMapSPDtrklt) FillSPDtrkltMap(ESDEvent, hitMapSPDtrklt);
 772         // fill tracklet hit map
 773
 774         const Int_t ipA = nfoctr[kIPA]; // inner layer A side
 775         const Int_t ipC = nfoctr[kIPC]; // inner layer C side
 776         const Int_t opA = nfoctr[kOPA]; // outer layer A side
 777         const Int_t opC = nfoctr[kOPC]; // outer layer C side
 778
 779         const Bool_t spdA = ipA + opA; // A side hit?
 780         const Bool_t spdC = ipC + opC; // C side hit?
 781
 782         return spdA * AliCDMesonBase::kBitSPDA + spdC * AliCDMesonBase::kBitSPDC;
 783 }
 784
 785
 786 //------------------------------------------------------------------------------
 787 Int_t AliCDMesonUtils::GetTPC(const AliESDEvent * ESDEvent, TH2 *TPCGapDCAaSide,
 788                               TH2 *TPCGapDCAcSide)
 789 {
 790         //
 791         //GetTPC
 792         //
 793
 794         const Double_t etacut = 0.9;
 795         Int_t nA = 0;
 796         Int_t nC = 0;
 797
 798         AliESDtrackCuts cuts;
 799         cuts.SetMaxDCAToVertexXY(0.1);
 800         cuts.SetMaxDCAToVertexZ(2.);
 801
 802         for(Int_t itrack = 0; itrack < ESDEvent->GetNumberOfTracks(); itrack++){
 803                 const AliESDtrack* esdtrack = ESDEvent->GetTrack(itrack);
 804                 Float_t b[2];
 805                 Float_t bCov[3];
 806                 esdtrack->GetImpactParameters(b,bCov);
 807
 808                 if (bCov[0]<=0 || bCov[2]<=0) {
 809                         printf("AliCDMesonUtils - Estimated b resolution lower or equal zero!\n");
 810                         bCov[0]=0;
 811                         bCov[2]=0;
 812                 }
 813
 814                 Float_t dcaToVertexXY = b[0];
 815                 Float_t dcaToVertexZ = b[1];
 816                 if (esdtrack->Eta() > etacut) {
 817                         if (cuts.AcceptTrack(esdtrack)) nA++;
 818                         if (TPCGapDCAaSide) {
 819                                 TPCGapDCAaSide->Fill(dcaToVertexXY, dcaToVertexZ);
 820                         }
 821                 }
 822                 else if (esdtrack->Eta() < -etacut) {
 823                         if (cuts.AcceptTrack(esdtrack)) nC++;
 824                         if (TPCGapDCAcSide) {
 825                                 TPCGapDCAcSide->Fill(dcaToVertexXY, dcaToVertexZ);
 826                         }
 827                 }
 828         }
 829
 830         const Bool_t tpcA = nA;
 831         const Bool_t tpcC = nC;
 832
 833         return tpcA * AliCDMesonBase::kBitTPCA + tpcC * AliCDMesonBase::kBitTPCC;
 834 }
 835
 836
 837 //------------------------------------------------------------------------------
 838 Int_t AliCDMesonUtils::GetZDC(const AliESDEvent * ESDEvent)
 839 {
 840         //
 841         //GetZDC
 842         //
 843
 844         const Int_t qa = ESDEvent->GetESDZDC()->GetESDQuality();
 845         Bool_t zdcA = kFALSE, zdcC = kFALSE;
 846         for(Int_t ii=0; ii<6; ii++){
 847                 if(qa & (1<<ii)){
 848                         if(ii<4) zdcA = kTRUE;
 849                         else zdcC = kTRUE;
 850                 }
 851         }
 852
 853         return zdcA * AliCDMesonBase::kBitZDCA + zdcC * AliCDMesonBase::kBitZDCC;
 854 }
 855
 856
 857 //==============================================================================
 858 //------------------------------------------------------------------------------
 859 #ifdef STD_ALIROOT
 860 Int_t AliCDMesonUtils::FMDHitCombinations(const AliESDEvent* ESDEvent,
 861                                           Int_t side)
 862 {
 863         //
 864         // copy of the FMDHitCombinations function originating from AliTriggerAnalysis
 865         //
 866         // side == 0 -> A side, side > 0 -> C side
 867
 868         // workaround for AliESDEvent::GetFMDData is not const!
 869         const AliESDFMD* fmdData =
 870                 (const_cast<AliESDEvent*>(ESDEvent))->GetFMDData();
 871         if (!fmdData)
 872                 {
 873                         puts("AliESDFMD not available");
 874                         return -1;
 875                 }
 876
 877         Int_t detFrom = (side == 0) ? 1 : 3;
 878         Int_t detTo   = (side == 0) ? 2 : 3;
 879
 880         Int_t triggers = 0;
 881         Float_t totalMult = 0;
 882         for (UShort_t det=detFrom;det<=detTo;det++) {
 883                 Int_t nRings = (det == 1 ? 1 : 2);
 884                 for (UShort_t ir = 0; ir < nRings; ir++) {
 885                         Char_t   ring = (ir == 0 ? 'I' : 'O');
 886                         UShort_t nsec = (ir == 0 ? 20  : 40);
 887                         UShort_t nstr = (ir == 0 ? 512 : 256);
 888                         for (UShort_t sec =0; sec < nsec;  sec++) {
 889                                 for (UShort_t strip = 0; strip < nstr; strip++) {
 890                                         Float_t mult = fmdData->Multiplicity(det,ring,sec,strip);
 891                                         if (mult == AliESDFMD::kInvalidMult) continue;
 892
 893                                         if (mult > 0.3) // fFMDLowCut
 894                                                 totalMult = totalMult + mult;
 895                                         else {
 896                                                 if (totalMult > 0.5) // fFMDHitCut
 897                                                         triggers++;
 898                                         }
 899                                 }
 900                         }
 901                 }
 902         }
 903         return triggers;
 904 }
 905 #endif
 906
 907
 908 //==============================================================================
 909 //------------------------------------------------------------------------------
 910 void AliCDMesonUtils::SPDLoadGeom(Int_t run)
 911 {
 912         // method to get the gGeomanager
 913         // it is called at the CreatedOutputObject stage
 914         // to comply with the CAF environment
 915
 916         AliCDBManager *man = AliCDBManager::Instance();
 917
 918         TString cdbpath;
 919         if (man->IsDefaultStorageSet()) {
 920                 const AliCDBStorage *dsto = man->GetDefaultStorage();
 921                 cdbpath = TString(dsto->GetBaseFolder());
 922                 //printf("default was set to: %s\n", cdbpath.Data());
 923         }
 924         else { //should not be used!
 925                 // man->SetDefaultStorage("alien://folder=/alice/data/2010/OCDB");
 926                 // would be needed on grid
 927                 man->SetDefaultStorage(gSystem->Getenv("TRAIN_CDB_PATH"));
 928                 cdbpath = TString(gSystem->Getenv("TRAIN_CDB_PATH"));
 929         }
 930
 931         man->SetSpecificStorage("ITS/Align/Data",cdbpath);
 932         man->SetSpecificStorage("GRP/Geometry/Data",cdbpath);
 933         man->SetRun(run);
 934
 935         AliCDBEntry* obj = man->Get(AliCDBPath("GRP", "Geometry", "Data"));
 936         if (!obj) {
 937                 printf("AliCDMesonUtils failed loading geometry object\n");
 938                 return;
 939         }
 940         AliGeomManager::SetGeometry((TGeoManager*)obj->GetObject());
 941         AliGeomManager::ApplyAlignObjsFromCDB("ITS");
 942 }
 943
 944 //------------------------------------------------------------------------------
 945 Bool_t AliCDMesonUtils::SPDLoc2Glo(Int_t id, const Double_t *loc,
 946                                    Double_t *glo)
 947 {
 948         //
 949         //SPDLoc2Glo, do not touch
 950         //
 951
 952         static TGeoHMatrix mat;
 953         Int_t vid = AliITSAlignMille2Module::GetVolumeIDFromIndex(id);
 954         if (vid<0) {
 955                 printf("AliCDMesonUtils Did not find module with such ID %d\n",id);
 956                 return kFALSE;
 957         }
 958         AliITSAlignMille2Module::SensVolMatrix(vid,&mat);
 959         mat.LocalToMaster(loc,glo);
 960         return kTRUE;
 961 }
 962
 963
 964 //------------------------------------------------------------------------------
 965 Int_t AliCDMesonUtils::CheckChipEta(Int_t chipKey, TString scut,
 966                                     const Double_t vtxPos[],
 967                                     TH2 *hitMapSPDinner, TH2 *hitMapSPDouter)
 968 {
 969         //
 970         //CheckChipEta
 971         //
 972
 973         // retrieves the position in eta for a given chip and applies the cut
 974         // results:
 975         // 0 <= out of range
 976         // -1 <= negative pseudo-rapidity position, in range (C-Side)
 977         // 1 <= positive pseudo-rapidity position, in range (A-Side)
 978         //
 979         // scut: "[0.9" or "]0.9", only 3 digits for the value!!
 980
 981
 982         const Bool_t kincl = (scut[0] == '[');
 983         const TString cutval = scut(1,3);
 984         const Double_t etacut = fabs(cutval.Atof());
 985
 986         //no eta cut, save time
 987         if(kincl && etacut>=2)
 988                 return kTRUE;
 989
 990         Int_t etaside = 1;
 991         //------------------------------- NOT TO TOUCH ------------------------>>
 992         UInt_t module=999, offchip=999;
 993         AliSPDUtils::GetOfflineFromOfflineChipKey(chipKey,module,offchip);
 994         UInt_t hs = AliSPDUtils::GetOnlineHSFromOffline(module);
 995         if(hs<2) offchip = 4 - offchip; // inversion  in the inner layer...
 996
 997         const Int_t col[]={
 998                 hs<2? 0 : 31,
 999                 hs<2? 31 : 0,
1000                 hs<2? 31 : 0,
1001                 hs<2? 0 : 31};
1002         const Int_t aa[]={0, 0, 255, 255};
1003         const AliITSsegmentationSPD seg;
1004
1005         for(Int_t ic=0; ic<4; ic++){
1006                 Float_t localchip[3]={0.,0.,0.};
1007                 seg.DetToLocal(aa[ic],col[ic]+32*offchip,localchip[0],localchip[2]);
1008                 // local coordinate of the chip center
1009                 //printf("local coordinates %d %d: %f %f \n",chipKey, ic, localchip[0],localchip[2]);
1010                 const Double_t local[3] = {localchip[0],localchip[1],localchip[2]};
1011                 Double_t glochip[3]={0.,0.,0.};
1012                 if(!SPDLoc2Glo(module,local,glochip)){
1013                         return kFALSE;
1014                 }
1015
1016                 //-------------------------------------------------------------------<<
1017
1018                 const TVector3 pos(glochip[0]-vtxPos[0], glochip[1]-vtxPos[1],
1019                                    glochip[2]-vtxPos[2]);
1020                 //pos.Print();
1021
1022                 if (chipKey < 400) { // inner SPD layer
1023                         if (hitMapSPDinner) {
1024                                 Double_t phi = pos.Phi(); // output in the range -Pi +Pi
1025                                 if (phi < 0.) phi += TMath::TwoPi(); // remap to the interval [0, TwoPi)
1026                                 const Double_t eta = pos.Eta();
1027                                 hitMapSPDinner->Fill(eta, phi);
1028                         }
1029                 }
1030                 else {
1031                         if (hitMapSPDouter) { // outer SPD layer
1032                                 Double_t phi = pos.Phi(); // output in the range -Pi +Pi
1033                                 if (phi < 0.) phi += TMath::TwoPi(); // remap to the interval [0, TwoPi)
1034                                 const Double_t eta = pos.Eta();
1035                                 hitMapSPDouter->Fill(eta, phi);
1036                         }
1037                 }
1038
1039                 if( kincl && fabs(pos.Eta()) > etacut)
1040                         return kFALSE;
1041
1042                 if(!kincl){
1043                         if(fabs(pos.Eta()) < etacut)
1044                                 return kFALSE;
1045                         else if(pos.Eta()<0)
1046                                 etaside = -1;
1047                         else
1048                                 etaside = 1;
1049                 }
1050         }
1051
1052         return etaside;
1053 }
1054
1055
1056 //------------------------------------------------------------------------------
1057 void AliCDMesonUtils::GetNFO(const AliESDEvent *ESDEvent, TString etacut,
1058                              Int_t ctr[], TH2 *hitMapSPDinner,
1059                              TH2 *hitMapSPDouter)
1060 {
1061         //
1062         // GetNFO
1063         //
1064         // analyzes the SPD fastOR for a given eta range and returns
1065         // an array with the number of hits in:
1066
1067         Int_t ninner=0; // inner layer
1068         Int_t nouter=0; // outer layer
1069         Int_t ipA = 0; // inner layer A side
1070         Int_t ipC = 0; // inner layer C side
1071         Int_t opA = 0; // outer layer A side
1072         Int_t opC = 0; // outer layer C side
1073
1074         const AliMultiplicity *mult = ESDEvent->GetMultiplicity();
1075
1076         // position of the primary vertex
1077         Double_t tmp[3] = { 0., 0., 0. };
1078         ESDEvent->GetPrimaryVertex()->GetXYZ(tmp);
1079         Double_t vtxPos[3] = { tmp[0], tmp[1], tmp[2] };
1080
1081
1082         for(Int_t iChipKey=0; iChipKey < 1200; iChipKey++){
1083                 if(mult->TestFastOrFiredChips(iChipKey)){
1084                         // here you check if the FastOr bit is 1 or 0
1085                         const Int_t iseta = CheckChipEta(iChipKey, etacut, vtxPos, hitMapSPDinner,
1086                                                          hitMapSPDouter);
1087                         if(iseta==0)
1088                                 continue;
1089
1090                         if(iChipKey<400) {
1091                                 ninner++;  // here you count the FastOr bits in the inner layer
1092                                 if(iseta>0)
1093                                         ipA ++;
1094                                 else
1095                                         ipC ++;
1096                         }
1097                         else {
1098                                 nouter++;  // here you count the FastOr bits in the outer layer
1099                                 if(iseta>0)
1100                                         opA ++;
1101                                 else
1102                                         opC ++;
1103                         }
1104                 }
1105         }
1106
1107         ctr[kInnerPixel]= ninner;
1108         ctr[kOuterPixel]= nouter;
1109         ctr[kIPA]= ipA;
1110         ctr[kIPC]= ipC;
1111         ctr[kOPA]= opA;
1112         ctr[kOPC]= opC;
1113
1114         return;
1115 }
1116
1117
1118 //--------------------------------------------------------------------------
1119 Int_t AliCDMesonUtils::GetFastORmultiplicity(const AliESDEvent* ESDEvent)
1120 {
1121         // determine the number of fired fastOR chips in both layers within
1122         // -0.9 < eta < 0.9
1123         //
1124
1125         const AliMultiplicity *mult = ESDEvent->GetMultiplicity();
1126
1127         // position of the primary vertex
1128         Double_t tmp[3] = { 0., 0., 0. };
1129         ESDEvent->GetPrimaryVertex()->GetXYZ(tmp);
1130         Double_t vtxPos[3] = { tmp[0], tmp[1], tmp[2] };
1131
1132         Int_t multiplicity = 0;
1133
1134         for (Int_t iChipKey=0; iChipKey < 1200; iChipKey++) {
1135                 if(mult->TestFastOrFiredChips(iChipKey)){
1136                         // here you check if the FastOr bit is 1 or 0
1137                         const Int_t iseta = CheckChipEta(iChipKey, "[0.9]", vtxPos, 0x0, 0x0);
1138                         if(iseta==0)
1139                                 continue;
1140                         else
1141                                 ++multiplicity;
1142                 }
1143         }
1144
1145         return multiplicity;
1146 }
1147
1148
1149 //--------------------------------------------------------------------------
1150 void AliCDMesonUtils::FillSPDtrkltMap(const AliVEvent* event,
1151                                       TH2 *hitMapSPDtrklt)
1152 {
1153         //
1154         // fill eta phi map of SPD tracklets
1155         //
1156
1157         if (hitMapSPDtrklt) {
1158                 if (TString(event->ClassName()) == "AliESDEvent") {
1159                         const AliMultiplicity *mult = ((AliESDEvent*)event)->GetMultiplicity();
1160                         for (Int_t iTracklet = 0; iTracklet < mult->GetNumberOfTracklets();
1161                              iTracklet++) {
1162                                 Double_t eta = mult->GetEta(iTracklet);
1163                                 Double_t phi = mult->GetPhi(iTracklet);
1164                                 hitMapSPDtrklt->Fill(eta, phi);
1165                         }
1166                 }
1167                 else if (TString(event->ClassName()) == "AliAODEvent") {
1168                         const AliAODTracklets* mult = ((AliAODEvent*)event)->GetTracklets();
1169                         for (Int_t iTracklet = 0; iTracklet < mult->GetNumberOfTracklets();
1170                              iTracklet++) {
1171                                 Double_t eta = -TMath::Log(TMath::Tan(mult->GetTheta(iTracklet)/2.));
1172                                 Double_t phi = mult->GetPhi(iTracklet);
1173                                 hitMapSPDtrklt->Fill(eta, phi);
1174                         }
1175                 }
1176         }
1177 }
1178
1179
1180 //==========================================================================
1181 Int_t AliCDMesonUtils::GetPID(AliPIDResponse *pid, const AliVTrack *trk,
1182                               Int_t mode /* = 0 */)
1183 {
1184         // determines PID for ESDs and AODs
1185         //
1186         //
1187
1188         if (!pid) return AliCDMesonBase::kBinPIDUnknown; // no pid available
1189
1190         Double_t tpcpion = -999., tpckaon = -999., tpcproton = -999.,
1191                 tpcelectron = -999.;
1192         Double_t tofpion = -999., tofkaon = -999., tofproton = -999.,
1193                 tofelectron = -999.;
1194         Double_t itspion = -999., itskaon = -999., itsproton = -999.,
1195                 itselectron = -999.;
1196
1197
1198         // check whether the track was pure ITS standalone track (has not left TPC)
1199         const Bool_t kits = !(trk->GetStatus() & AliESDtrack::kTPCout);
1200
1201         if (kits) { // do ITS pid
1202                 const Double_t nin = 3.;
1203
1204                 itspion = pid->NumberOfSigmasITS( trk, AliPID::kPion );
1205                 itskaon = pid->NumberOfSigmasITS( trk, AliPID::kKaon );
1206                 itsproton = pid->NumberOfSigmasITS( trk, AliPID::kProton );
1207                 itselectron = pid->NumberOfSigmasITS( trk, AliPID::kElectron );
1208
1209                 if(fabs(itspion) < nin) // inclusion only
1210                         return AliCDMesonBase::kBinPion;
1211                 else if(fabs(itskaon) < nin)
1212                         return AliCDMesonBase::kBinKaon;
1213                 else if(fabs(itsproton) < nin)
1214                         return AliCDMesonBase::kBinProton;
1215                 else if(fabs(itselectron) < nin)
1216                         return AliCDMesonBase::kBinElectron;
1217                 else{
1218                         return AliCDMesonBase::kBinPIDUnknown;
1219                 }
1220         }
1221
1222         tpcpion     = pid->NumberOfSigmasTPC( trk, AliPID::kPion );
1223         tpckaon     = pid->NumberOfSigmasTPC( trk, AliPID::kKaon );
1224         tpcproton   = pid->NumberOfSigmasTPC( trk, AliPID::kProton );
1225         tpcelectron = pid->NumberOfSigmasTPC( trk, AliPID::kElectron );
1226
1227         // derive information, whether tof pid is available
1228         const Bool_t ka = !(trk->GetStatus() & AliESDtrack::kTOFmismatch);
1229         const Bool_t kb =  (trk->GetStatus() & AliESDtrack::kTOFpid);
1230         const Bool_t ktof = ka && kb;
1231
1232         // retrieve TOF information if available
1233         if(ktof){
1234                 tofpion = pid->NumberOfSigmasTOF(trk, AliPID::kPion);
1235                 tofkaon = pid->NumberOfSigmasTOF(trk, AliPID::kKaon);
1236                 tofproton = pid->NumberOfSigmasTOF(trk, AliPID::kProton);
1237                 tofelectron = pid->NumberOfSigmasTOF(trk, AliPID::kElectron);
1238         }
1239
1240         if (mode == 0) {
1241                 // 2010 cuts for PID
1242                 const Double_t nin = 3.;
1243
1244                 // inclusion + exclusion (either TPC or TOF)
1245                 if((fabs(tpcpion) < nin && fabs(tpckaon) > nin && fabs(tpcproton) > nin
1246                     && fabs(tpcelectron) > nin) ||
1247                    (fabs(tofpion) < nin && fabs(tofkaon) > nin && fabs(tofproton) > nin
1248                     && fabs(tofelectron) > nin))
1249                         return AliCDMesonBase::kBinPionE;
1250                 else if((fabs(tpcpion) > nin && fabs(tpckaon) < nin && fabs(tpcproton) > nin
1251                          && fabs(tpcelectron) > nin) ||
1252                         (fabs(tofpion) > nin && fabs(tofkaon) < nin && fabs(tofproton) > nin
1253                          && fabs(tofelectron) > nin))
1254                         return AliCDMesonBase::kBinKaonE;
1255                 else if((fabs(tpcpion) > nin && fabs(tpckaon) > nin && fabs(tpcproton) < nin
1256                          && fabs(tpcelectron) > nin) ||
1257                         (fabs(tofpion) > nin && fabs(tofkaon) > nin && fabs(tofproton) < nin
1258                          && fabs(tofelectron) > nin))
1259                         return AliCDMesonBase::kBinProtonE;
1260                 else if((fabs(tpcpion) > nin && fabs(tpckaon) > nin && fabs(tpcproton) > nin
1261                          && fabs(tpcelectron) < nin) ||
1262                         (fabs(tofpion) > nin && fabs(tofkaon) > nin && fabs(tofproton) > nin
1263                          && fabs(tofelectron) < nin))
1264                         return AliCDMesonBase::kBinElectronE;
1265                 else if(fabs(tpcpion) < nin && fabs(tofpion) < nin) // inclusion (TPC + TOF)
1266                         return AliCDMesonBase::kBinPion;
1267                 else if(fabs(tpckaon) < nin && fabs(tofkaon) < nin)
1268                         return AliCDMesonBase::kBinKaon;
1269                 else if(fabs(tpcproton) < nin && fabs(tofproton) < nin)
1270                         return AliCDMesonBase::kBinProton;
1271                 else if(fabs(tpcelectron) < nin && fabs(tofelectron) < nin)
1272                         return AliCDMesonBase::kBinElectron;
1273                 else{
1274                         return AliCDMesonBase::kBinPIDUnknown;
1275                 }
1276         }
1277         else if (mode == 1) {
1278                 // 2011 cuts for PID in LHC11f
1279                 // TPC: [-3,5] sigma (pion)
1280                 // TOF: 3 sigma for all,
1281                 // only Pion is tuned!
1282                 // ONLY INCLUSION CUTS NO EXCLUSION
1283                 const Double_t nin = 3.;
1284
1285                 if(tpcpion < 4. && tpcpion > -2. && fabs(tofpion) < -3. )
1286                         return AliCDMesonBase::kBinPion;
1287                 else if(fabs(tpckaon) < nin && fabs(tofkaon) < nin)
1288                         return AliCDMesonBase::kBinKaon;
1289                 else if(fabs(tpcproton) < nin &&  fabs(tofproton) < nin)
1290                         return AliCDMesonBase::kBinProton;
1291                 else if(fabs(tpcelectron) < nin &&  fabs(tofelectron) < nin)
1292                         return AliCDMesonBase::kBinElectron;
1293                 else{
1294                         return AliCDMesonBase::kBinPIDUnknown;
1295                 }
1296         }
1297         return AliCDMesonBase::kBinPIDUnknown;
1298 }
1299
1300
1301 //==========================================================================
1302 Int_t AliCDMesonUtils::GetPID(Int_t pdgCode)
1303 {
1304         //
1305         // determine particle type based on PDG code
1306         //
1307
1308         if (TMath::Abs(pdgCode) == 211) return AliCDMesonBase::kBinPionE;
1309         else if (TMath::Abs(pdgCode) == 321) return AliCDMesonBase::kBinKaonE;
1310         else if (TMath::Abs(pdgCode) == 2212) return AliCDMesonBase::kBinProtonE;
1311         else if (TMath::Abs(pdgCode) == 11) return AliCDMesonBase::kBinElectronE;
1312         else return AliCDMesonBase::kBinPIDUnknown;
1313 }
1314
1315
1316 //------------------------------------------------------------------------------
1317 Int_t AliCDMesonUtils::CombinePID(const Int_t pid[])
1318 {
1319         //
1320         // combine the PID result
1321         //
1322
1323         // determine return value
1324         if (pid[0] == pid[1]) { // same result for both tracks
1325                 return pid[0];
1326         }
1327         // one track identified with exclusion the other only without
1328         else if ((pid[0] == AliCDMesonBase::kBinPionE &&
1329                   pid[1] == AliCDMesonBase::kBinPion) ||
1330                  (pid[1] == AliCDMesonBase::kBinPionE &&
1331                   pid[0] == AliCDMesonBase::kBinPion)) {
1332                 return AliCDMesonBase::kBinPion;
1333         }
1334         else if ((pid[0] == AliCDMesonBase::kBinKaonE &&
1335                   pid[1] == AliCDMesonBase::kBinKaon) ||
1336                  (pid[1] == AliCDMesonBase::kBinKaonE &&
1337                   pid[0] == AliCDMesonBase::kBinKaon)) {
1338                 return AliCDMesonBase::kBinKaon;
1339         }
1340         else if ((pid[0] == AliCDMesonBase::kBinProtonE &&
1341                   pid[1] == AliCDMesonBase::kBinProton) ||
1342                  (pid[1] == AliCDMesonBase::kBinProtonE &&
1343                   pid[0] == AliCDMesonBase::kBinProton)) {
1344                 return AliCDMesonBase::kBinProton;
1345         }
1346         else if ((pid[0] == AliCDMesonBase::kBinElectronE &&
1347                   pid[1] == AliCDMesonBase::kBinElectron) ||
1348                  (pid[1] == AliCDMesonBase::kBinElectronE &&
1349                   pid[0] == AliCDMesonBase::kBinElectron)) {
1350                 return AliCDMesonBase::kBinElectron;
1351         }
1352         // one track identified and one not
1353         else if (((pid[0] == AliCDMesonBase::kBinPionE ||
1354                    pid[0] == AliCDMesonBase::kBinPion) &&
1355                   pid[1] == AliCDMesonBase::kBinPIDUnknown) ||
1356                  (pid[1] == AliCDMesonBase::kBinPion &&
1357                   pid[0] == AliCDMesonBase::kBinPIDUnknown)) {
1358                 return AliCDMesonBase::kBinSinglePion;
1359         }
1360         else if (((pid[0] == AliCDMesonBase::kBinKaonE ||
1361                    pid[0] == AliCDMesonBase::kBinKaon)&&
1362                   pid[1] == AliCDMesonBase::kBinPIDUnknown) ||
1363                  (pid[1] == AliCDMesonBase::kBinKaon &&
1364                   pid[0] == AliCDMesonBase::kBinPIDUnknown)) {
1365                 return AliCDMesonBase::kBinSingleKaon;
1366         }
1367         else if (((pid[0] == AliCDMesonBase::kBinProtonE ||
1368                    pid[0] == AliCDMesonBase::kBinProton) &&
1369                   pid[1] == AliCDMesonBase::kBinPIDUnknown) ||
1370                  (pid[1] == AliCDMesonBase::kBinProton &&
1371                   pid[0] == AliCDMesonBase::kBinPIDUnknown)) {
1372                 return AliCDMesonBase::kBinSingleProton;
1373         }
1374         else if (((pid[0] == AliCDMesonBase::kBinElectronE ||
1375                    pid[0] == AliCDMesonBase::kBinElectron) &&
1376                   pid[1] == AliCDMesonBase::kBinPIDUnknown) ||
1377                  (pid[1] == AliCDMesonBase::kBinElectron &&
1378                   pid[0] == AliCDMesonBase::kBinPIDUnknown)) {
1379                 return AliCDMesonBase::kBinSingleElectron;
1380         }
1381         else
1382                 return AliCDMesonBase::kBinPIDUnknown;
1383 }
1384
1385
1386 //------------------------------------------------------------------------------
1387 TLorentzVector AliCDMesonUtils::GetKinematics(const Double_t *pa,
1388                                               const Double_t *pb,
1389                                               Double_t ma,
1390                                               Double_t mb, Double_t& cts)
1391 {
1392         //
1393         //get kinematics, cts = cos(theta_{#star})
1394         //
1395
1396         TLorentzVector va, vb;
1397         va.SetXYZM(pa[0], pa[1], pa[2], ma);
1398         vb.SetXYZM(pb[0], pb[1], pb[2], mb);
1399         const TLorentzVector sumv = va+vb;
1400
1401         const TVector3 bv = -sumv.BoostVector();
1402
1403         va.Boost(bv);
1404         vb.Boost(bv);
1405
1406         // 3-vectors in the restframe of the mother particle
1407         const TVector3 pra = va.Vect();
1408         const TVector3 prb = vb.Vect();
1409         const TVector3 diff = pra - prb; // their difference
1410
1411         cts = (diff.Mag2()-pra.Mag2()-prb.Mag2()) / (-2.*pra.Mag()*prb.Mag());
1412         // cosine theta star, calculated according to the law-of-cosines
1413
1414         return sumv;
1415 }
1416
1417
1418 //------------------------------------------------------------------------------
1419 Double_t AliCDMesonUtils::GetOA(const Double_t *pa, const Double_t *pb)
1420 {
1421         //
1422         //cosOpeningAngle
1423         //
1424
1425         TVector3 va, vb;
1426         va.SetXYZ(pa[0], pa[1], pa[2]);
1427         vb.SetXYZ(pb[0], pb[1], pb[2]);
1428
1429         const TVector3 ua = va.Unit();
1430         const TVector3 ub = vb.Unit();
1431
1432         return ua.Dot(ub);
1433 }