PWGLF/SPECTRA/ChargedHadrons/dNdPt/AlidNdPtCutAnalysis.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 // AlidNdPtCutAnalysis class.
  17 //
  18 // a. functionality:
  19 // - fills generic cut histograms
  20 // - generates cuts (selection criteria)
  21 //
  22 // b. data members:
  23 // - generic cut histograms
  24 // - control histograms
  25 //
  26 // Author: J.Otwinowski 04/11/2008
  27 //------------------------------------------------------------------------------
  28 #include "TH1.h"
  29 #include "TH2.h"
  30
  31 #include "AliHeader.h"
  32 #include "AliGenEventHeader.h"
  33 #include "AliInputEventHandler.h"
  34 #include "AliStack.h"
  35 #include "AliESDEvent.h"
  36 #include "AliMCEvent.h"
  37 #include "AliESDtrackCuts.h"
  38 #include "AliLog.h"
  39 #include "AliTracker.h"
  40
  41 #include "AlidNdPtEventCuts.h"
  42 #include "AlidNdPtAcceptanceCuts.h"
  43 #include "AlidNdPtBackgroundCuts.h"
  44 #include "AlidNdPtAnalysis.h"
  45 #include "AliPhysicsSelection.h"
  46
  47 #include "AliPWG0Helper.h"
  48 #include "AlidNdPtHelper.h"
  49 #include "AlidNdPtCutAnalysis.h"
  50
  51 using namespace std;
  52
  53 ClassImp(AlidNdPtCutAnalysis)
  54
  55 //_____________________________________________________________________________
  56   AlidNdPtCutAnalysis::AlidNdPtCutAnalysis(): AlidNdPt(),
  57   fAnalysisFolder(0),
  58   fEventCount(0),
  59   fRecEventHist(0),
  60   fMCEventHist(0),
  61   fRecMCEventHist(0),
  62   fRecMCTrackHist(0)
  63 {
  64   // default constructor
  65   Init();
  66 }
  67
  68 //_____________________________________________________________________________
  69 AlidNdPtCutAnalysis::AlidNdPtCutAnalysis(Char_t* name, Char_t* title): AlidNdPt(name,title),
  70   fAnalysisFolder(0),
  71   fEventCount(0),
  72   fRecEventHist(0),
  73   fMCEventHist(0),
  74   fRecMCEventHist(0),
  75   fRecMCTrackHist(0)
  76 {
  77   // constructor
  78   Init();
  79 }
  80
  81 //_____________________________________________________________________________
  82 AlidNdPtCutAnalysis::~AlidNdPtCutAnalysis() {
  83   //
  84   if(fEventCount) delete fEventCount; fEventCount=0;
  85   if(fRecEventHist) delete fRecEventHist; fRecEventHist=0;
  86   if(fMCEventHist) delete fMCEventHist; fMCEventHist=0;
  87   if(fRecMCEventHist) delete fRecMCEventHist; fRecMCEventHist=0;
  88   if(fRecMCTrackHist) delete fRecMCTrackHist; fRecMCTrackHist=0;
  89
  90   if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
  91 }
  92
  93 //_____________________________________________________________________________
  94 void AlidNdPtCutAnalysis::Init(){
  95   //
  96   // Init histograms
  97   //
  98   /*
  99   const Int_t ptNbins = 56;
 100   const Double_t ptMin = 0.;
 101   const Double_t ptMax = 16.;
 102   Double_t binsPt[ptNbins+1] = {0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.2,2.4,2.6,2.8,3.0,3.2,3.4,3.6,3.8,4.0,4.5,5.0,5.5,6.0,6.5,7.0,7.5,8.0,9.0,10.0,11.0,12.0,13.0,14.0,15.0,16.0};
 103   */
 104   // set pt bins
 105   const Int_t ptNbins = 50;
 106   const Double_t ptMin = 1.e-2, ptMax = 50.;
 107   Double_t *binsPt = CreateLogAxis(ptNbins,ptMin,ptMax);
 108
 109   //
 110   Int_t binsEventCount[2]={2,2};
 111   Double_t minEventCount[2]={0,0};
 112   Double_t maxEventCount[2]={2,2};
 113   fEventCount = new THnSparseF("fEventCount","trig vs trig+vertex",2,binsEventCount,minEventCount,maxEventCount);
 114   fEventCount->GetAxis(0)->SetTitle("trig");
 115   fEventCount->GetAxis(1)->SetTitle("trig+vert");
 116   fEventCount->Sumw2();
 117
 118   //Xv:Yv:Zv:ResZv:Mult
 119   Double_t kFact = 1.0;
 120
 121   Int_t binsRecEventHist[5]={80,80,100,80,150};
 122   Double_t minRecEventHist[5]={-3.*kFact,-3.*kFact,-35.,0.,0.};
 123   Double_t maxRecEventHist[5]={3.*kFact,3.*kFact,35.,10.,150.};
 124   fRecEventHist = new THnSparseF("fRecEventHist","Xv:Yv:Zv:ResZv:Mult",5,binsRecEventHist,minRecEventHist,maxRecEventHist);
 125   fRecEventHist->GetAxis(0)->SetTitle("Xv (cm)");
 126   fRecEventHist->GetAxis(1)->SetTitle("Yv (cm)");
 127   fRecEventHist->GetAxis(2)->SetTitle("Zv (cm)");
 128   fRecEventHist->GetAxis(3)->SetTitle("ResZv (cm)");
 129   fRecEventHist->GetAxis(4)->SetTitle("Mult");
 130   fRecEventHist->Sumw2();
 131
 132   //Xv:Yv:Zv
 133   Int_t binsMCEventHist[3]={80,80,100};
 134   Double_t minMCEventHist[3]={-0.1,-0.1,-35.};
 135   Double_t maxMCEventHist[3]={0.1,0.1,35.};
 136   fMCEventHist = new THnSparseF("fMCEventHist","mcXv:mcYv:mcZv",3,binsMCEventHist,minMCEventHist,maxMCEventHist);
 137   fMCEventHist->GetAxis(0)->SetTitle("mcXv (cm)");
 138   fMCEventHist->GetAxis(1)->SetTitle("mcYv (cm)");
 139   fMCEventHist->GetAxis(2)->SetTitle("mcZv (cm)");
 140   fMCEventHist->Sumw2();
 141
 142   //Xv-mcXv:Yv-mcYv:Zv-mcZv:Mult
 143   Int_t binsRecMCEventHist[4]={100,100,100,150};
 144   Double_t minRecMCEventHist[4]={-1.0*kFact,-1.0*kFact,-1.0*kFact,0.};
 145   Double_t maxRecMCEventHist[4]={1.0*kFact,1.0*kFact,1.0*kFact,150.};
 146   fRecMCEventHist = new THnSparseF("fRecMCEventHist","Xv-mcXv:Yv-mcYv:Zv-mcZv:Mult",4,binsRecMCEventHist,minRecMCEventHist,maxRecMCEventHist);
 147   fRecMCEventHist->GetAxis(0)->SetTitle("Xv-mcXv (cm)");
 148   fRecMCEventHist->GetAxis(1)->SetTitle("Yv-mcYv (cm)");
 149   fRecMCEventHist->GetAxis(2)->SetTitle("Zv-mcZv (cm)");
 150   fRecMCEventHist->GetAxis(3)->SetTitle("Mult");
 151   fRecMCEventHist->Sumw2();
 152
 153   //
 154   // THnSparse track histograms
 155   //
 156
 157  //nCrossRows:chi2PerClust:nCrossRows/nFindableClust:fracSharedClust:DCAy:DCAz:eta:phi:pt:hasStrangeMother:isFromMaterial:isPrim:charge
 158   Int_t binsRecMCTrackHist[13]=  { 160,  10,  20,  20, 50,  50,   20,  90,             ptNbins, 2,  2,  2,  3  };
 159   Double_t minRecMCTrackHist[13]={ 0.,   0.,  0.,  0., -0.5,-0.5,-1.0, 0.,             ptMin,   0., 0., 0.,-1. };
 160   Double_t maxRecMCTrackHist[13]={ 160., 10., 1.,  1., 0.5, 0.5,  1.0, 2.*TMath::Pi(), ptMax,   2., 2., 2., 2. };
 161
 162   fRecMCTrackHist = new THnSparseF("fRecMCTrackHist","nCrossRows:chi2PerClust:nCrossRows/nFindableClust:fracSharedClust:DCAy:DCAz:eta:phi:pt:hasStrangeMother:isFromMaterial:isPrim:charge",13,binsRecMCTrackHist,minRecMCTrackHist,maxRecMCTrackHist);
 163   fRecMCTrackHist->SetBinEdges(8,binsPt);
 164   fRecMCTrackHist->GetAxis(0)->SetTitle("nCrossRows");
 165   fRecMCTrackHist->GetAxis(1)->SetTitle("chi2PerClust");
 166   fRecMCTrackHist->GetAxis(2)->SetTitle("nCrossRows/nFindableClust");
 167   fRecMCTrackHist->GetAxis(3)->SetTitle("fracSharedClust");
 168   fRecMCTrackHist->GetAxis(4)->SetTitle("DCAy (cm)");
 169   fRecMCTrackHist->GetAxis(5)->SetTitle("DCAz (cm)");
 170   fRecMCTrackHist->GetAxis(6)->SetTitle("#eta");
 171   fRecMCTrackHist->GetAxis(7)->SetTitle("#phi (rad)");
 172   fRecMCTrackHist->GetAxis(8)->SetTitle("p_{T} (GeV/c)");
 173   fRecMCTrackHist->GetAxis(9)->SetTitle("hasStrangeMother");
 174   fRecMCTrackHist->GetAxis(10)->SetTitle("isFromMaterial");
 175   fRecMCTrackHist->GetAxis(11)->SetTitle("isPrim");
 176   fRecMCTrackHist->GetAxis(12)->SetTitle("charge");
 177   fRecMCTrackHist->Sumw2();
 178
 179   //nClust:chi2PerClust:nClust/nFindableClust:DCAy:DCAz:eta:phi:pt:kinkIdx:isPrim:polarity
 180   /*
 181   Int_t binsRecMCTrackHist[11]={160,80,80,100,100,90,90,ptNbins, 3, 2, 2};
 182   Double_t minRecMCTrackHist[11]={0., 0., 0., -1.,-1.,-1.5, 0., ptMin, -1., 0., 0.};
 183   Double_t maxRecMCTrackHist[11]={160.,10.,1.2, 1.,1.,1.5, 2.*TMath::Pi(), ptMax, 2., 2., 2.};
 184
 185   fRecMCTrackHist = new THnSparseF("fRecMCTrackHist","nClust:chi2PerClust:nClust/nFindableClust:DCAy:DCAz:eta:phi:pt:kinkIdx:isPrim:polarity",11,binsRecMCTrackHist,minRecMCTrackHist,maxRecMCTrackHist);
 186   fRecMCTrackHist->SetBinEdges(7,binsPt);
 187
 188   fRecMCTrackHist->GetAxis(0)->SetTitle("nClust");
 189   fRecMCTrackHist->GetAxis(1)->SetTitle("chi2PerClust");
 190   fRecMCTrackHist->GetAxis(2)->SetTitle("nClust/nFindableClust");
 191   fRecMCTrackHist->GetAxis(3)->SetTitle("DCAy (cm)");
 192   fRecMCTrackHist->GetAxis(4)->SetTitle("DCAz (cm)");
 193   fRecMCTrackHist->GetAxis(5)->SetTitle("#eta");
 194   fRecMCTrackHist->GetAxis(6)->SetTitle("#phi (rad)");
 195   fRecMCTrackHist->GetAxis(7)->SetTitle("p_{T} (GeV/c)");
 196   fRecMCTrackHist->GetAxis(8)->SetTitle("kinkIdx"); // 0 - no kink, -1 - kink mother, 1 - kink daugther
 197   fRecMCTrackHist->GetAxis(9)->SetTitle("isPrim");
 198   fRecMCTrackHist->GetAxis(10)->SetTitle("polarity");
 199   fRecMCTrackHist->Sumw2();
 200   */
 201
 202   // init output folder
 203   fAnalysisFolder = CreateFolder("folderdNdPt","Analysis dNdPt Folder");
 204
 205 }
 206
 207 //_____________________________________________________________________________
 208 void AlidNdPtCutAnalysis::Process(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent)
 209 {
 210   //
 211   // Process real and/or simulated events
 212   //
 213   if(!esdEvent) {
 214     AliDebug(AliLog::kError, "esdEvent not available");
 215     return;
 216   }
 217
 218   // get selection cuts
 219   AlidNdPtEventCuts *evtCuts = GetEventCuts();
 220   AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts();
 221   AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
 222
 223   if(!evtCuts || !accCuts  || !esdTrackCuts) {
 224     AliDebug(AliLog::kError, "cuts not available");
 225     return;
 226   }
 227
 228   // trigger selection
 229
 230   Bool_t isEventTriggered = kTRUE;
 231   AliPhysicsSelection *physicsSelection = NULL;
 232   AliTriggerAnalysis* triggerAnalysis = NULL;
 233
 234   //
 235   AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
 236   if (!inputHandler)
 237   {
 238     Printf("ERROR: Could not receive input handler");
 239     return;
 240   }
 241
 242   if(evtCuts->IsTriggerRequired())
 243   {
 244     // always MB
 245     isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
 246
 247     physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
 248     if(!physicsSelection) return;
 249     //SetPhysicsTriggerSelection(physicsSelection);
 250
 251     if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
 252       // set trigger (V0AND)
 253       triggerAnalysis = physicsSelection->GetTriggerAnalysis();
 254       if(!triggerAnalysis) return;
 255       isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
 256     }
 257   }
 258
 259   // use MC information
 260   AliHeader* header = 0;
 261   AliGenEventHeader* genHeader = 0;
 262   AliStack* stack = 0;
 263   TArrayF vtxMC(3);
 264   AliPWG0Helper::MCProcessType evtType = AliPWG0Helper::kInvalidProcess;
 265
 266   if(IsUseMCInfo())
 267   {
 268     if(!mcEvent) {
 269       AliDebug(AliLog::kError, "mcEvent not available");
 270       return;
 271     }
 272
 273     // get MC event header
 274     header = mcEvent->Header();
 275     if (!header) {
 276       AliDebug(AliLog::kError, "Header not available");
 277       return;
 278     }
 279
 280     // MC particle stack
 281     stack = mcEvent->Stack();
 282     if (!stack) {
 283       AliDebug(AliLog::kError, "Stack not available");
 284       return;
 285     }
 286
 287     // get event type (ND=0x1, DD=0x2, SD=0x4)
 288     evtType = AliPWG0Helper::GetEventProcessType(header);
 289     AliDebug(AliLog::kDebug+1, Form("Found process type %d", evtType));
 290
 291     // get MC vertex
 292     genHeader = header->GenEventHeader();
 293     if (!genHeader) {
 294       AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
 295       return;
 296     }
 297     genHeader->PrimaryVertex(vtxMC);
 298
 299     // Fill MC event histogram
 300     Double_t vMCEventHist[3]={vtxMC[0],vtxMC[1],vtxMC[2]};
 301     fMCEventHist->Fill(vMCEventHist);
 302
 303   } // end bUseMC
 304
 305   // get reconstructed vertex
 306   Bool_t bRedoTPCVertex = evtCuts->IsRedoTPCVertex();
 307   Bool_t bUseConstraints = evtCuts->IsUseBeamSpotConstraint();
 308   const AliESDVertex* vtxESD = AlidNdPtHelper::GetVertex(esdEvent,evtCuts,accCuts,esdTrackCuts,GetAnalysisMode(),kFALSE,bRedoTPCVertex,bUseConstraints);
 309   Bool_t isRecVertex = AlidNdPtHelper::TestRecVertex(vtxESD, esdEvent->GetPrimaryVertexSPD(), GetAnalysisMode(), kFALSE);
 310   Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD) && isRecVertex;
 311
 312   TObjArray *allChargedTracks=0;
 313   Int_t multAll=0;
 314
 315   //
 316   // event counter
 317   //
 318   //printf("isEventOK %d, isEventTriggered %d \n",isEventOK,isEventTriggered);
 319
 320   Bool_t isTrigAndVertex = isEventTriggered && isEventOK;
 321   Double_t vEventCount[2] = { isEventTriggered, isTrigAndVertex};
 322   fEventCount->Fill(vEventCount);
 323
 324   //
 325   // cosmic background and splitted tracks
 326   //
 327   if(GetParticleMode() == AlidNdPtHelper::kBackgroundTrack)
 328   {
 329     AlidNdPtBackgroundCuts *backCuts = GetBackgroundCuts();
 330     if(!backCuts) return;
 331
 332     for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
 333     {
 334       AliESDtrack *track1 = esdEvent->GetTrack(iTrack);
 335       if(!track1) continue;
 336       if(track1->Charge()==0) continue;
 337
 338       for (Int_t jTrack = iTrack+1; jTrack < esdEvent->GetNumberOfTracks(); jTrack++)
 339       {
 340         AliESDtrack *track2 = esdEvent->GetTrack(jTrack);
 341         if(!track2) continue;
 342         if(track2->Charge()==0) continue;
 343
 344         //printf("track2->Charge() %d\n",track2->Charge());
 345
 346         backCuts->IsBackgroundTrack(track1, track2);
 347       }
 348     }
 349   }
 350
 351   // check event cuts
 352   if(isEventOK && isEventTriggered)
 353   {
 354     // get all charged tracks
 355     allChargedTracks = AlidNdPtHelper::GetAllChargedTracks(esdEvent,GetAnalysisMode());
 356     if(!allChargedTracks) return;
 357
 358     Int_t entries = allChargedTracks->GetEntries();
 359     for(Int_t i=0; i<entries;++i)
 360     {
 361       AliESDtrack *track = (AliESDtrack*)allChargedTracks->At(i);
 362       if(!track) continue;
 363
 364       if(!esdTrackCuts->AcceptTrack(track)) continue;
 365
 366       //
 367       Bool_t isOK = kFALSE;
 368       Double_t x[3]; track->GetXYZ(x);
 369       Double_t b[3]; AliTracker::GetBxByBz(x,b);
 370
 371       //
 372       // if TPC-ITS hybrid tracking (kTPCITSHybrid)
 373       // replace track parameters with TPC-ony track parameters
 374       //
 375       if( GetAnalysisMode() == AlidNdPtHelper::kTPCITSHybrid || GetAnalysisMode() == AlidNdPtHelper::kTPCITSHybridTrackSPDvtx || GetAnalysisMode() == AlidNdPtHelper::kTPCITSHybridTrackSPDvtxDCArPt)
 376       {
 377         // Relate TPC-only tracks to SPD vertex
 378         isOK = track->RelateToVertexTPCBxByBz(vtxESD, b, kVeryBig);
 379         if(!isOK) continue;
 380
 381         // replace esd track parameters with TPCinner
 382         AliExternalTrackParam  *tpcTrack  = new AliExternalTrackParam(*(track->GetTPCInnerParam()));
 383         if (!tpcTrack) return;
 384         track->Set(tpcTrack->GetX(),tpcTrack->GetAlpha(),tpcTrack->GetParameter(),tpcTrack->GetCovariance());
 385
 386         if(tpcTrack) delete tpcTrack;
 387       }
 388
 389       //
 390       if (GetAnalysisMode()==AlidNdPtHelper::kTPCSPDvtxUpdate || GetAnalysisMode() == AlidNdPtHelper::kTPCTrackSPDvtxUpdate)
 391       {
 392         //
 393         // update track parameters
 394         //
 395         AliExternalTrackParam cParam;
 396         isOK = track->RelateToVertexTPCBxByBz(vtxESD, b, kVeryBig, &cParam);
 397         if(!isOK) continue;
 398
 399         track->Set(cParam.GetX(),cParam.GetAlpha(),cParam.GetParameter(),cParam.GetCovariance());
 400       }
 401
 402       FillHistograms(track, stack);
 403       multAll++;
 404     }
 405
 406     Double_t vRecEventHist[5] = {vtxESD->GetXv(),vtxESD->GetYv(),vtxESD->GetZv(),vtxESD->GetZRes(),multAll};
 407     fRecEventHist->Fill(vRecEventHist);
 408
 409     if(IsUseMCInfo()) {
 410       Double_t vRecMCEventHist[5] = {vtxESD->GetXv()-vtxMC[0],vtxESD->GetYv()-vtxMC[1],vtxESD->GetZv()-vtxMC[2],multAll};
 411       fRecMCEventHist->Fill(vRecMCEventHist);
 412     }
 413   }
 414
 415   if(allChargedTracks) delete allChargedTracks; allChargedTracks = 0;
 416
 417 }
 418
 419 //_____________________________________________________________________________
 420 void AlidNdPtCutAnalysis::FillHistograms(AliESDtrack *const esdTrack, AliStack *const stack) const
 421 {
 422   //
 423   // Fill ESD track and MC histograms
 424   //
 425   if(!esdTrack) return;
 426   if(esdTrack->Charge() == 0.) return;
 427
 428   Float_t pt = esdTrack->Pt();
 429   Float_t eta = esdTrack->Eta();
 430   Float_t phi = esdTrack->Phi();
 431
 432
 433   Int_t nClust = 0;
 434   if(GetAnalysisMode() == AlidNdPtHelper::kTPC) {
 435     nClust = esdTrack->GetTPCNclsIter1();
 436   } else {
 437     nClust = esdTrack->GetTPCclusters(0);
 438   }
 439
 440   Float_t chi2PerCluster = 0.;
 441   if(GetAnalysisMode() == AlidNdPtHelper::kTPC) {
 442     if(nClust>0.) chi2PerCluster = esdTrack->GetTPCchi2Iter1()/Float_t(nClust);
 443   } else {
 444     chi2PerCluster = esdTrack->GetTPCchi2()/Float_t(nClust);
 445   }
 446
 447   Int_t nFindableClust = esdTrack->GetTPCNclsF();
 448
 449
 450   Float_t clustPerFindClust = 0.;
 451   if(nFindableClust>0.) clustPerFindClust = Float_t(nClust)/nFindableClust;
 452
 453   Float_t b[2], bCov[3];
 454   esdTrack->GetImpactParameters(b,bCov);
 455
 456   //
 457   Float_t nCrossedRowsTPC = esdTrack->GetTPCClusterInfo(2,1);
 458
 459   Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
 460   if (esdTrack->GetTPCNclsF()>0) {
 461      ratioCrossedRowsOverFindableClustersTPC = esdTrack->GetTPCClusterInfo(2,1)/esdTrack->GetTPCNclsF();
 462   }
 463
 464   //
 465   Int_t nClustersTPCShared = esdTrack->GetTPCnclsS();
 466   Float_t fracClustersTPCShared = -1.;
 467   fracClustersTPCShared = Float_t(nClustersTPCShared)/Float_t(nClust);
 468
 469
 470   // kink idx
 471   Int_t kinkIdx = 0;
 472   //if(esdTrack->GetKinkIndex(0) > 0.)   isKink  = kTRUE;
 473   if(esdTrack->GetKinkIndex(0) > 0)      kinkIdx = 1;   // kink daughter
 474   else if(esdTrack->GetKinkIndex(0) < 0) kinkIdx = -1;  // kink mother
 475   else kinkIdx = 0; // not kink
 476
 477   //printf("esdTrack->GetKinkIndex(0) %d \n", esdTrack->GetKinkIndex(0));
 478   //printf("esdTrack->GetKinkIndex(1) %d \n", esdTrack->GetKinkIndex(1));
 479   //printf("esdTrack->GetKinkIndex(2) %d \n", esdTrack->GetKinkIndex(2));
 480   //printf("kinkIdx %d \n", kinkIdx);
 481
 482   //
 483   // Fill rec vs MC information
 484   //
 485   Bool_t isPrim = kTRUE;
 486   Bool_t hasStrangeMother = kFALSE;
 487   Bool_t isFromMaterial = kFALSE;
 488
 489   if(IsUseMCInfo()) {
 490     if(!stack) return;
 491     Int_t label = TMath::Abs(esdTrack->GetLabel());
 492     TParticle* particle = stack->Particle(label);
 493     if(!particle) return;
 494     if(particle->GetPDG() && particle->GetPDG()->Charge()==0.) return;
 495     isPrim = stack->IsPhysicalPrimary(label);
 496
 497     // check whether has stange mother
 498     //
 499     Int_t motherPdg = -1;
 500     TParticle* mother = 0;
 501
 502     Int_t motherLabel = particle->GetMother(0);
 503     if(motherLabel>0) mother = stack->Particle(motherLabel);
 504     if(mother) motherPdg = TMath::Abs(mother->GetPdgCode()); // take abs for visualisation only
 505     Int_t mech = particle->GetUniqueID(); // production mechanism
 506
 507     if( (motherPdg == 3122) || (motherPdg == -3122) || (motherPdg == 310)) // lambda, antilambda, k0s
 508     {
 509       if( (mech == 4) || (mech == 5) ) hasStrangeMother = kTRUE;
 510     }
 511     else {
 512       //if(isPrim==0 && mech == 13)
 513       //printf("mech %d \n", mech);
 514       if(!isPrim) isFromMaterial = kTRUE;
 515     }
 516
 517     //if(isPrim && pt > 1.5 && kinkIdx == -1) printf("nClust  %d \n", nClust);
 518   }
 519
 520   // fill histo
 521   Int_t charge = esdTrack->Charge();
 522
 523   //Double_t vRecMCTrackHist[11] = { nClust,chi2PerCluster,clustPerFindClust,b[0],b[1],eta,phi,pt,kinkIdx,isPrim, polarity };
 524   //fRecMCTrackHist->Fill(vRecMCTrackHist);
 525
 526   Double_t vRecMCTrackHist[13] = { nCrossedRowsTPC, chi2PerCluster, ratioCrossedRowsOverFindableClustersTPC, fracClustersTPCShared , b[0], b[1], eta, phi, pt, hasStrangeMother, isFromMaterial, isPrim, charge };
 527   fRecMCTrackHist->Fill(vRecMCTrackHist);
 528
 529 }
 530
 531 //_____________________________________________________________________________
 532 Long64_t AlidNdPtCutAnalysis::Merge(TCollection* const list)
 533 {
 534   // Merge list of objects (needed by PROOF)
 535
 536   if (!list)
 537   return 0;
 538
 539   if (list->IsEmpty())
 540   return 1;
 541
 542   TIterator* iter = list->MakeIterator();
 543   TObject* obj = 0;
 544
 545   //TList *collPhysSelection = new TList;
 546
 547   // collection of generated histograms
 548   Int_t count=0;
 549   while((obj = iter->Next()) != 0) {
 550     AlidNdPtCutAnalysis* entry = dynamic_cast<AlidNdPtCutAnalysis*>(obj);
 551     if (entry == 0) continue;
 552
 553     // event histo
 554     fEventCount->Add(entry->fEventCount);
 555     fRecEventHist->Add(entry->fRecEventHist);
 556     fRecMCEventHist->Add(entry->fRecMCEventHist);
 557     fMCEventHist->Add(entry->fMCEventHist);
 558
 559     // track histo
 560     fRecMCTrackHist->Add(entry->fRecMCTrackHist);
 561
 562     // physics selection
 563     //collPhysSelection->Add(entry->GetPhysicsTriggerSelection());
 564
 565   count++;
 566   }
 567
 568   //AliPhysicsSelection *trigSelection = GetPhysicsTriggerSelection();
 569   //trigSelection->Merge(collPhysSelection);
 570
 571   //if(collPhysSelection) delete collPhysSelection;
 572
 573 return count;
 574 }
 575
 576 //_____________________________________________________________________________
 577 void AlidNdPtCutAnalysis::Analyse()
 578 {
 579   //
 580   // Analyse histograms
 581   //
 582   TH1::AddDirectory(kFALSE);
 583   TObjArray *aFolderObj = new TObjArray;
 584   if(!aFolderObj) return;
 585
 586   TH1D *h1D = 0;
 587   TH2D *h2D = 0;
 588
 589
 590   //
 591   // get cuts
 592   //
 593   AlidNdPtEventCuts *evtCuts = GetEventCuts();
 594   AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts();
 595   AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
 596
 597   if(!evtCuts || !accCuts || !esdTrackCuts) {
 598     Error("AlidNdPtCutAnalysis::Analyse()", "cuts not available");
 599     return;
 600   }
 601
 602   //
 603   // set min and max values
 604   //
 605   Double_t minPt = accCuts->GetMinPt();
 606   Double_t maxPt = accCuts->GetMaxPt();
 607   Double_t minEta = accCuts->GetMinEta();
 608   Double_t maxEta = accCuts->GetMaxEta()-0.00001;
 609
 610   Double_t maxDCAr = accCuts->GetMaxDCAr();
 611
 612   //
 613   // Event counters
 614   //
 615   h2D = (TH2D*)fEventCount->Projection(0,1);
 616   if(!h2D) return;
 617   h2D->SetName("trig_vs_trigANDvertex");
 618   aFolderObj->Add(h2D);
 619
 620   fEventCount->GetAxis(0)->SetRange(2,2); // triggered
 621   h1D = (TH1D*)fEventCount->Projection(1);
 622   if(!h1D) return;
 623   h1D->SetTitle("rec. vertex for triggered events");
 624   h1D->SetName("trigANDvertex");
 625   aFolderObj->Add(h1D);
 626
 627   //
 628   // Create rec. event histograms
 629   //
 630   h1D = (TH1D *)fRecEventHist->Projection(0);
 631   if(!h1D) return;
 632   h1D->SetName("rec_xv");
 633   aFolderObj->Add(h1D);
 634
 635   h1D = (TH1D *)fRecEventHist->Projection(1);
 636   if(!h1D) return;
 637   h1D->SetName("rec_yv");
 638   aFolderObj->Add(h1D);
 639
 640   h1D = (TH1D *)fRecEventHist->Projection(2);
 641   if(!h1D) return;
 642   h1D->SetName("rec_zv");
 643   aFolderObj->Add(h1D);
 644
 645   h2D = (TH2D *)fRecEventHist->Projection(3,4);
 646   if(!h2D) return;
 647   h2D->SetName("rec_resZv_vs_Mult");
 648   aFolderObj->Add(h2D);
 649
 650   h2D = (TH2D *)fRecEventHist->Projection(0,1);
 651   if(!h2D) return;
 652   h2D->SetName("rec_xv_vs_yv");
 653   aFolderObj->Add(h2D);
 654
 655   h2D = (TH2D *)fRecEventHist->Projection(0,2);
 656   if(!h2D) return;
 657   h2D->SetName("rec_xv_vs_zv");
 658   aFolderObj->Add(h2D);
 659
 660   h2D = (TH2D *)fRecEventHist->Projection(3,4);
 661   if(!h2D) return;
 662   h2D->SetName("rec_resZv_vs_Mult");
 663   aFolderObj->Add(h2D);
 664
 665   //
 666   // MC available
 667   //
 668   if(IsUseMCInfo()) {
 669
 670   //
 671   // Create mc event histograms
 672   //
 673   h2D = (TH2D *)fMCEventHist->Projection(0,1);
 674   if(!h2D) return;
 675   h2D->SetName("mc_xv_vs_yv");
 676   aFolderObj->Add(h2D);
 677
 678   h2D = (TH2D *)fMCEventHist->Projection(0,2);
 679   if(!h2D) return;
 680   h2D->SetName("mc_xv_vs_zv");
 681   aFolderObj->Add(h2D);
 682
 683   //
 684   // Create rec-mc event histograms
 685   //
 686   h2D = (TH2D *)fRecMCEventHist->Projection(0,3);
 687   if(!h2D) return;
 688   h2D->SetName("rec_mc_deltaXv_vs_mult");
 689   aFolderObj->Add(h2D);
 690
 691   h2D = (TH2D *)fRecMCEventHist->Projection(1,3);
 692   if(!h2D) return;
 693   h2D->SetName("rec_mc_deltaYv_vs_mult");
 694   aFolderObj->Add(h2D);
 695
 696   h2D = (TH2D *)fRecMCEventHist->Projection(2,3);
 697   if(!h2D) return;
 698   h2D->SetName("rec_mc_deltaZv_vs_mult");
 699   aFolderObj->Add(h2D);
 700
 701   } // end use MC info
 702
 703
 704
 705   //
 706   // Create rec-mc track track histograms
 707   //
 708
 709   // DCA cuts
 710   fRecMCTrackHist->GetAxis(3)->SetRangeUser(-maxDCAr,maxDCAr);
 711   fRecMCTrackHist->GetAxis(4)->SetRangeUser(-maxDCAr,maxDCAr);
 712
 713   h2D = (TH2D *)fRecMCTrackHist->Projection(7,5);
 714   if(!h2D) return;
 715   h2D->SetName("pt_vs_eta");
 716   aFolderObj->Add(h2D);
 717
 718   fRecMCTrackHist->GetAxis(7)->SetRangeUser(minPt,maxPt);
 719
 720   h2D = (TH2D *)fRecMCTrackHist->Projection(0,5);
 721   if(!h2D) return;
 722   h2D->SetName("nClust_vs_eta");
 723   aFolderObj->Add(h2D);
 724
 725   h2D = (TH2D *)fRecMCTrackHist->Projection(1,5);
 726   if(!h2D) return;
 727   h2D->SetName("chi2PerClust_vs_eta");
 728   aFolderObj->Add(h2D);
 729
 730   h2D = (TH2D *)fRecMCTrackHist->Projection(2,5);
 731   if(!h2D) return;
 732   h2D->SetName("ratio_nClust_nFindableClust_vs_eta");
 733   aFolderObj->Add(h2D);
 734
 735   h2D = (TH2D *)fRecMCTrackHist->Projection(5,6);
 736   if(!h2D) return;
 737   h2D->SetName("eta_vs_phi");
 738   aFolderObj->Add(h2D);
 739
 740   //
 741   fRecMCTrackHist->GetAxis(5)->SetRangeUser(minEta,maxEta);
 742
 743   h2D = (TH2D *)fRecMCTrackHist->Projection(0,6);
 744   if(!h2D) return;
 745   h2D->SetName("nClust_vs_phi");
 746   aFolderObj->Add(h2D);
 747
 748   h2D = (TH2D *)fRecMCTrackHist->Projection(1,6);
 749   if(!h2D) return;
 750   h2D->SetName("chi2PerClust_vs_phi");
 751   aFolderObj->Add(h2D);
 752
 753   h2D = (TH2D *)fRecMCTrackHist->Projection(2,6);
 754   if(!h2D) return;
 755   h2D->SetName("ratio_nClust_nFindableClust_vs_phi");
 756   aFolderObj->Add(h2D);
 757
 758   //
 759   fRecMCTrackHist->GetAxis(7)->SetRangeUser(0.0,maxPt);
 760
 761   h2D = (TH2D *)fRecMCTrackHist->Projection(0,7);
 762   if(!h2D) return;
 763   h2D->SetName("nClust_vs_pt");
 764   aFolderObj->Add(h2D);
 765
 766   h2D = (TH2D *)fRecMCTrackHist->Projection(1,7);
 767   if(!h2D) return;
 768   h2D->SetName("chi2PerClust_vs_pt");
 769   aFolderObj->Add(h2D);
 770
 771   h2D = (TH2D *)fRecMCTrackHist->Projection(2,7);
 772   if(!h2D) return;
 773   h2D->SetName("ratio_nClust_nFindableClust_vs_pt");
 774   aFolderObj->Add(h2D);
 775
 776   h2D = (TH2D *)fRecMCTrackHist->Projection(6,7);
 777   if(!h2D) return;
 778   h2D->SetName("phi_vs_pt");
 779   aFolderObj->Add(h2D);
 780
 781
 782   // fiducial volume
 783   fRecMCTrackHist->GetAxis(5)->SetRangeUser(minEta,maxEta);
 784   fRecMCTrackHist->GetAxis(7)->SetRangeUser(minPt,maxPt);
 785
 786   // DCA cuts
 787   fRecMCTrackHist->GetAxis(3)->SetRangeUser(-maxDCAr,maxDCAr);
 788   fRecMCTrackHist->GetAxis(4)->SetRangeUser(-maxDCAr,maxDCAr);
 789
 790   h2D = (TH2D *)fRecMCTrackHist->Projection(0,1);
 791   if(!h2D) return;
 792   h2D->SetName("nClust_vs_chi2PerClust");
 793   aFolderObj->Add(h2D);
 794
 795   h2D = (TH2D *)fRecMCTrackHist->Projection(0,2);
 796   if(!h2D) return;
 797   h2D->SetName("nClust_vs_ratio_nClust_nFindableClust");
 798   aFolderObj->Add(h2D);
 799
 800   //
 801   // DCAy cuts
 802   //
 803   fRecMCTrackHist->GetAxis(0)->SetRange(50,160); // nClust/track > 50
 804   fRecMCTrackHist->GetAxis(1)->SetRangeUser(0.,3.9999); // chi2/cluster < 4.0
 805   fRecMCTrackHist->GetAxis(3)->SetRange(1,fRecMCTrackHist->GetAxis(3)->GetNbins());
 806   //fRecMCTrackHist->GetAxis(4)->SetRangeUser(-1.0,1.0);
 807   fRecMCTrackHist->GetAxis(4)->SetRange(1,fRecMCTrackHist->GetAxis(4)->GetNbins());
 808
 809   // sec
 810   fRecMCTrackHist->GetAxis(9)->SetRange(1,1);
 811   h1D = (TH1D *)fRecMCTrackHist->Projection(3);
 812   if(!h1D) return;
 813   h1D->SetName("dcay_sec");
 814   aFolderObj->Add(h1D);
 815
 816   // prim
 817   fRecMCTrackHist->GetAxis(9)->SetRange(2,2);
 818   h1D = (TH1D *)fRecMCTrackHist->Projection(3);
 819   if(!h1D) return;
 820   h1D->SetName("dcay_prim");
 821   aFolderObj->Add(h1D);
 822
 823   // DCAz cuts
 824   //fRecMCTrackHist->GetAxis(3)->SetRangeUser(-1.0,1.0);
 825   fRecMCTrackHist->GetAxis(4)->SetRange(1,fRecMCTrackHist->GetAxis(4)->GetNbins());
 826
 827   // sec
 828   fRecMCTrackHist->GetAxis(9)->SetRange(1,1);
 829   h1D = (TH1D *)fRecMCTrackHist->Projection(4);
 830   if(!h1D) return;
 831   h1D->SetName("dcaz_sec");
 832   aFolderObj->Add(h1D);
 833
 834   // prim
 835   fRecMCTrackHist->GetAxis(9)->SetRange(2,2);
 836   h1D = (TH1D *)fRecMCTrackHist->Projection(4);
 837   if(!h1D) return;
 838   h1D->SetName("dcaz_prim");
 839   aFolderObj->Add(h1D);
 840
 841
 842   // export objects to analysis folder
 843   fAnalysisFolder = ExportToFolder(aFolderObj);
 844   if(!fAnalysisFolder) {
 845       if(aFolderObj) delete aFolderObj;
 846       return;
 847   }
 848
 849   // delete only TObjArray
 850   if(aFolderObj) delete aFolderObj;
 851 }
 852
 853 //_____________________________________________________________________________
 854 TFolder* AlidNdPtCutAnalysis::ExportToFolder(TObjArray * const array)
 855 {
 856   // recreate folder avery time and export objects to new one
 857   //
 858   AlidNdPtCutAnalysis * comp=this;
 859   TFolder *folder = comp->GetAnalysisFolder();
 860
 861   TString name, title;
 862   TFolder *newFolder = 0;
 863   Int_t i = 0;
 864   Int_t size = array->GetSize();
 865
 866   if(folder) {
 867      // get name and title from old folder
 868      name = folder->GetName();
 869      title = folder->GetTitle();
 870
 871          // delete old one
 872      delete folder;
 873
 874          // create new one
 875      newFolder = CreateFolder(name.Data(),title.Data());
 876      newFolder->SetOwner();
 877
 878          // add objects to folder
 879      while(i < size) {
 880            newFolder->Add(array->At(i));
 881            i++;
 882          }
 883   }
 884
 885 return newFolder;
 886 }
 887
 888 //_____________________________________________________________________________
 889 TFolder* AlidNdPtCutAnalysis::CreateFolder(TString name,TString title) {
 890 // create folder for analysed histograms
 891 //
 892 TFolder *folder = 0;
 893   folder = new TFolder(name.Data(),title.Data());
 894
 895   return folder;
 896 }