PWGCF/EBYE/LRC/AliAnalysisTaskLRC.cxx

   1 /**************************************************************************
   2  * Author: Andrey Ivanov.                                           *
   3  * Contributors are mentioned in the code where appropriate.              *
   4  *                                                                        *
   5  * Permission to use, copy, modify and distribute this software and its   *
   6  * documentation strictly for non-commercial purposes is hereby granted   *
   7  * without fee, provided that the above copyright notice appears in all   *
   8  * copies and that both the copyright notice and this permission notice   *
   9  * appear in the supporting documentation. The authors make no claims     *
  10  * about the suitability of this software for any purpose. It is          *
  11  * provided "as is" without express or implied warranty.                  *
  12  **************************************************************************/
  13 // Analysis task for Long Range Correlation (LRC) analysis using TPC data
  14 // This includes a TList of AliLRCBase objects that are processing LRC analysis
  15 // for a given Eta window
  16 // This task is worcking with ESD data only
  17
  18 // Authors : Andrey Ivanov , Igor Altsybeev, St.Peterburg State University
  19 // Email: Andrey.Ivanov@cern.ch, Igor.Altsybeev@cern.ch
  20
  21 #include <AliAnalysisManager.h>
  22 #include <AliESDInputHandler.h>
  23 #include "AliAnalysisTaskLRC.h"
  24 #include <AliLRCBase.h>
  25 //#include <AliLRCProcess.h>
  26 #include <AliVEvent.h>
  27 #include <AliMCEvent.h>
  28 #include <AliESDEvent.h>
  29 #include "AliAODEvent.h"
  30 #include "AliAODTrack.h"
  31 #include "AliAODInputHandler.h"
  32
  33 #include "AliGenEventHeader.h"
  34 #include "AliGenHijingEventHeader.h"
  35
  36 #include <AliPID.h>
  37 #include <AliPIDCombined.h>
  38 #include <AliPIDResponse.h>
  39 #include <AliMultiplicity.h>
  40 #include <AliESDv0.h>
  41
  42 #include <AliRunLoader.h>
  43 #include <AliRun.h>
  44 #include <AliStack.h>
  45 #include <AliESDtrackCuts.h>
  46 #include <AliKalmanTrack.h>
  47 #include <AliPhysicsSelection.h>
  48 #include <AliCentrality.h>
  49 #include <AliESDZDC.h>
  50 #include <AliESDFMD.h>
  51
  52 #include "AliEventplane.h"
  53
  54
  55 #include <TParticle.h>
  56 #include <TH1.h>
  57 #include <TH2.h>
  58 #include "TH1I.h"
  59 #include "TRandom3.h"
  60 #include <TComplex.h>
  61
  62 #include <TF1.h>
  63
  64 //#include "AliSimpleEvent.h"
  65
  66 //flag use my Event Tree
  67 //#define fSetIncludeEventTreeInOutput 1
  68
  69
  70 #include "TStopwatch.h"
  71
  72 using std::endl;
  73 using std::cout;
  74
  75 ClassImp(AliAnalysisTaskLRC)
  76
  77 //________________________________________________________________________
  78 AliAnalysisTaskLRC::AliAnalysisTaskLRC( const char *name, Bool_t runKine)
  79     : AliAnalysisTaskSE(name)
  80     ,fNumberOfPhiSectors(1)
  81     ,fFlagSuppressAddingSomeHistos(kTRUE)
  82     ,fAnalysisLevel("ESD")
  83     ,fEsdTrackCuts(0)
  84     ,fAODtrackCutBit(128)
  85     ,fArrTrackCuts(0x0)
  86     ,fHistCutsNamesBins(0)
  87     //    ,fNumberOfCutsToRemember(0)
  88     ,fSwitchToListingCuts(0)
  89
  90     //    ,fAnalysisType(en_AnalysisType_ESD)
  91     ,fMinNumberOfSPDtracklets(0)
  92     ,fMaxPtLimit(100.0)
  93     ,fMinPtLimit(0.0)
  94     ,fKineLowPtCut(0.0)
  95     ,fMinAcceptedTracksCut(0)
  96     ,fMaxAcceptedTracksCut(0)
  97     ,fCheckForkVtx(kTRUE)
  98     ,fCheckForVtxPosition(kFALSE)
  99     ,fVxMax(0)
 100     ,fVyMax(0)
 101     ,fVzMax(0)
 102     ,fLRCproc(0)
 103     ,fOutList(0)
 104     ,fRunKine(0)
 105     ,fAnalyseMCESD(0)
 106     ,fShowEventStats(kFALSE)
 107     ,fShowPerTrackStats(kFALSE)
 108     ,fEtaRegionForTests(0.8)
 109     ,fMultCutInEtaRegion(0)
 110     ,fHistEventCutStats(0)
 111     ,fHistTrackCutStats(0)
 112     ,fHistAODTrackStats(0)
 113     ,fHistVx(0)
 114     ,fHistVy(0)
 115     ,fHistVz(0)
 116     ,fHistVxMCrecoDiff(0)
 117     ,fHistVyMCrecoDiff(0)
 118     ,fHistVzMCrecoDiff(0)
 119     ,fHistVertexNconributors(0)
 120     ,fHistNumberOfPileupVerticesTracks(0)
 121     ,fHistNumberOfPileupVerticesSPD(0)
 122     ,fHistEventPlane(0)
 123     ,fHistPt(0)
 124     ,fHistEta(0)
 125     ,fHistEtaAODpure(0)
 126     ,fHistPhi(0)
 127     ,fHistEtaPhi(0)
 128     ,fHistPhiLRCrotationsCheck(0)
 129     ,fHistPhiArtificialProfilingCheck(0)
 130     ,fHistPhiArtificialProfilingCheckWrtEvPlane(0)
 131     ,fHistPhiArtificialEvPlane(0)
 132     ,fHistEtaVsZvCoverage(0)
 133     ,fHistEtaVsZvCoverageAccepted(0)
 134     ,fHistMultBeforeCuts(0)
 135     ,fHistAcceptedMult(0)
 136     ,fHistAcceptedTracks(0)
 137     ,fHistMultiplicityInEtaRegion(0)
 138     ,fHistMultiplicityInEtaRegionAfterPtCuts(0)
 139     ,fHist2DMultiplicityMCESDInEtaRegion(0)
 140     ,fHistAcceptedTracksAfterPtCuts(0)
 141     ,fHistAcceptedTPCtracks(0)
 142     ,fHistClustersTPC(0)
 143     ,fHistClustersTPCafterCuts(0)
 144     ,fHistCrossedRowsTPC(0)
 145     ,fHistCrossedRowsTPCafterCuts(0)
 146     ,fHistClustersITS(0)
 147     ,fHistTrackletsITS(0)
 148     ,fHist2DClustersTPCvsPt(0)
 149     ,fHist2DClustersTPCvsEta(0)
 150     ,fHist2DAcceptedTracksPtvsEta(0)
 151     ,fHistMClabels(0)
 152     ,fHistRejectedTracksCharge(0)
 153     ,fHistTracksCharge(0)
 154     ,fHistPtPlus(0)
 155     ,fHistPtMinus(0)
 156     ,fHistNetChargeVsPt(0)
 157     ,fHistChargePlusVsPtTmp(0)
 158     ,fHistChargeMinusVsPtTmp(0)
 159     ,fHist2DNetChargeVsPt(0)
 160     ,fHist2DNetChargeVsPtCorrectedOnEventMean(0)
 161     ,fHist2DNetChargeVsPtCorrectedOnEventMeanNormOnNch(0)
 162     ,fHistProbabilitiesPID(0)
 163     ,fHistESDtrackMass(0)
 164     ,fHistProbabilityPion(0)
 165     ,fHistProbabilityKaon(0)
 166     ,fHistProbabilityProton(0)
 167     ,fHistParticlesDistr(0)
 168     ,fHistParticlesDistrBeforeCuts(0)
 169     //    ,fNumberOfSectors(1)
 170     ,fHistCentralityPercentile(0)
 171     ,fHistCentralityClass10(0)
 172     ,fHistCentralityClass5(0)
 173     //,fHistZDCenergy(0x0)
 174     ,fHistZDCparticipants(0)
 175     ,fHistV0multiplicity(0)
 176     ,fHistV0Amultiplicity(0)
 177     ,fHistV0Cmultiplicity(0)
 178     ,fHist2DV0ACmultiplicity(0)
 179     ,fHist2DTracksAcceptedVsV0multiplicity(0)
 180     //,fHistV0spectra(0)
 181
 182     ,fHistV0cells    (0)
 183     ,fHistV0Acells   (0)
 184     ,fHistV0Ccells   (0)
 185     ,fHist2DV0ACcells(0)
 186
 187     ,fThresholdOnV0mult(0)
 188
 189
 190     ,fMinCentralityClass(-0.01)
 191     ,fMaxCentralityClass(98)
 192     ,fIsIonsAnalysis(kFALSE)
 193     ,fEtInsteadOfPt(kFALSE)
 194     ,fUsePhiShufflingByHand(kFALSE)
 195     ,fUseToyEvents(kFALSE)
 196     ,fNumberOfToyEvents(1)
 197     ,fTmpCounter(0)
 198     ,fPIDResponse(0x0)
 199     ,fPIDCombined(0x0)
 200
 201     ,fHistPidMaxProbability(0)
 202     ,fHistPidPureMaxProbability(0)
 203     ,fStrPIDforFwd("any")
 204     ,fStrPIDforBwd("any")
 205
 206     ,fPIDsensingFlag(kFALSE)
 207     ,fArtificialInefficiency(-1.)
 208     ,fHistNumberOfDroppedByHandTracks(0)
 209     ,fRand(0)
 210
 211
 212     ,fPhiArtificialGapBegin(0)
 213     ,fPhiArtificialGapEnd(0)
 214
 215     ,fFlagWatchZDC(0)
 216     ,fFlagWatchV0 (0)
 217     ,fFlagWatchFMD(0)
 218
 219     ,fAnalysisTimer(0)
 220
 221     //MC qa and studies
 222     //    ,fHistMCvertexRdeltaFromParent(0)
 223     //    ,fHistMCparentsStat(0)
 224     //,fHistMCparentsEta(0)
 225     //,fHistMCchildsEta(0)
 226     //,fHistMCdeltaEtaChildParent(0)
 227     //,fProbabilitiesPID(0)
 228
 229
 230     //    ,fSimpleEvent(0)
 231     //    ,fNsimpleEvents(0)
 232     //    ,fEventTree(0)
 233     //    ,fSetIncludeEventTreeInOutput(0)
 234 {
 235     fAnalysisTimer = new TStopwatch;
 236     //Init
 237     fRunKine = runKine;
 238     for ( int i = 0; i < 5; i++ )
 239     {
 240         fHistZDCenergy[i] = 0x0;
 241     }
 242     for ( int i = 0; i < 4; i++ )
 243     {
 244         fHistV0AmultiplicityRing[i] = 0x0;
 245         fHistV0CmultiplicityRing[i] = 0x0;
 246         fHist2DV0ACmultiplicityRing[i] = 0x0;
 247         fHist2DTracksAcceptedVsV0AmultiplicityRing[i] = 0x0;
 248         fHist2DTracksAcceptedVsV0CmultiplicityRing[i] = 0x0;
 249     }
 250     //fProbabilitiesPID = new Double_t[AliPID::kSPECIES];
 251     // Output slot #1 writes into a TList container for common task data and QA
 252     DefineOutput(1, TList::Class());
 253
 254     //Defining output slots for each LRC processor (required to avoid TList of TLists on merging)
 255     for(Int_t i=0; i < fLRCproc.GetEntries(); i++)
 256     {
 257         DefineOutput(Proc(i)->GetOutputSlotNumber(),TList::Class());
 258     }
 259 }
 260
 261
 262 // ---------------------------------------  Setters ------------------
 263
 264 void AliAnalysisTaskLRC::SetMaxPtLimit(Double_t MaxPtLimit)
 265 {
 266     //Sets  Max Pt filter
 267     fMaxPtLimit = MaxPtLimit;
 268 }
 269 void AliAnalysisTaskLRC::SetMinPtLimit(Double_t MinPtLimit)
 270 {
 271     //Sets  Min Pt filter
 272     fMinPtLimit = MinPtLimit;
 273 }
 274 void AliAnalysisTaskLRC::SetKineLowPtCut(Double_t kineLowPtCut)
 275 {
 276     //Sets  Min Pt filter for Kine tracks
 277     fKineLowPtCut = kineLowPtCut;
 278 }
 279 AliLRCBase*  AliAnalysisTaskLRC::Proc(Int_t index)
 280 {
 281     // Get Processor i
 282     return (dynamic_cast<AliLRCBase*> (fLRCproc.At(index)));
 283 }
 284
 285 void AliAnalysisTaskLRC::SetMinNumberOfSPDtracklets( Int_t MinSPDtracklets )
 286 {
 287     //Sets  Min SPD tracklets number
 288     fMinNumberOfSPDtracklets = MinSPDtracklets;
 289 }
 290
 291 //________________________________________________________________________
 292 void AliAnalysisTaskLRC::UserCreateOutputObjects()
 293 {
 294     // --------- Output list
 295     fOutList = new TList();
 296
 297     //### added 13.12.2011
 298     fOutList->SetOwner();  // IMPORTANT!
 299
 300     fAnalysisTimer->Start();
 301
 302
 303     //##### PID stuff
 304     // ------- setup PIDCombined
 305     fPIDCombined=new AliPIDCombined;
 306     fPIDCombined->SetDefaultTPCPriors();
 307     fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC+AliPIDResponse::kDetTOF);
 308
 309     // no light nuclei
 310     fPIDCombined->SetSelectedSpecies(AliPID::kSPECIES);
 311     //Int_t lNumberOfPidSpecies = (Int_t)AliPID::kSPECIES;
 312     for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec)
 313     {
 314         fProbTPCTOF[ispec]=new TH2D(Form("prob%s_mom_TPCTOF",AliPID::ParticleName(ispec)),
 315                                     Form("%s probability vs. momentum;momentum;probability",AliPID::ParticleName(ispec)),
 316                                     100,0.,20.,50,0.,1.);
 317         fProbAllDets[ispec]=new TH2D(Form("prob%s_mom_AllDets",AliPID::ParticleName(ispec)),
 318                                      Form("%s probability vs. momentum;momentum;probability",AliPID::ParticleName(ispec)),
 319                                      100,0.,20.,50,0.,1.);
 320         if ( !fFlagSuppressAddingSomeHistos )
 321         {
 322             fOutList->Add(fProbTPCTOF[ispec]);
 323             fOutList->Add(fProbAllDets[ispec]);
 324         }
 325
 326
 327         // basic priors
 328         fPriors[ispec] = new TH1F(Form("%s_priors",AliPID::ParticleName(ispec)),
 329                                   Form("%s priors vs momentum",AliPID::ParticleName(ispec)),
 330                                   100,0.,20.);
 331         if ( !fFlagSuppressAddingSomeHistos )
 332             fOutList->Add(fPriors[ispec]);
 333         switch (ispec) {
 334         case AliPID::kElectron:
 335             for (Int_t ich=1;ich<=100;ich++) fPriors[ispec]->SetBinContent(ich,0.02);
 336             break;
 337         case AliPID::kMuon:
 338             for (Int_t ich=1;ich<=100;ich++) fPriors[ispec]->SetBinContent(ich,0.02);
 339             break;
 340         case AliPID::kPion:
 341             for (Int_t ich=1;ich<=100;ich++) fPriors[ispec]->SetBinContent(ich,0.56);
 342             break;
 343         case AliPID::kKaon:
 344             for (Int_t ich=1;ich<=100;ich++) fPriors[ispec]->SetBinContent(ich,0.20);
 345             break;
 346         case AliPID::kProton:
 347             for (Int_t ich=1;ich<=100;ich++) fPriors[ispec]->SetBinContent(ich,0.20);
 348             break;
 349         default:
 350             break;
 351         }
 352         fPIDCombined->SetPriorDistribution((AliPID::EParticleType)ispec,fPriors[ispec]);
 353
 354         // priors used
 355         fPriorsUsed[ispec] = new TH2D(Form("%s_priorsUsed",AliPID::ParticleName(ispec)),
 356                                       Form("%s priors vs transverse momentum;p_{t} (GeV/c);priors",AliPID::ParticleName(ispec)),
 357                                       100,0.,20.,101,0,1.01);
 358         if ( !fFlagSuppressAddingSomeHistos )
 359             fOutList->Add(fPriorsUsed[ispec]);
 360     }
 361
 362     fHistPidMaxProbability = new TH1D("fHistPidMaxProbability","HistPidMaxProbability;Probability;Entries",100,0,1);
 363     fOutList->Add(fHistPidMaxProbability);
 364
 365     fHistPidPureMaxProbability = new TH1D("fHistPidPureMaxProbability","HistPidMaxProbability;Probability;Entries",100,0,1);
 366     fOutList->Add(fHistPidPureMaxProbability);
 367
 368     // ##### end PID part
 369     Printf("UserCreateOutputObjects.........");
 370     //Disabling "replacing existing TH2D ...etc" warning
 371
 372     Bool_t lTH1oldStatus = TH1::AddDirectoryStatus();
 373     TH1::AddDirectory(kFALSE);
 374
 375     Bool_t lTH2oldStatus = TH2::AddDirectoryStatus();
 376     TH2::AddDirectory(kFALSE);
 377
 378
 379
 380     //LRC processors init
 381     Int_t lLrcNum = fLRCproc.GetEntries();
 382     for(Int_t i = 0; i < lLrcNum; i++)
 383     {
 384         AliLRCBase *lrcBase = dynamic_cast<AliLRCBase*> (fLRCproc.At(i));
 385         if(lrcBase)
 386             lrcBase->InitDataMembers();
 387         else continue;
 388         //remember pointer to be used in the analysis
 389         //fLRCprocArrayPointers[i] = lrcBase;
 390     }
 391
 392
 393     //fOutList->Add( fEsdTrackCuts );
 394
 395
 396     //create here array with bin names
 397     //15.03.2013 - list with cuts to apply on tracks and remember decisions
 398     if ( fSwitchToListingCuts )
 399     {
 400         int lNumberOfCutsToRemember = fArrTrackCuts.GetEntries();
 401         fHistCutsNamesBins = new TH1I("fHistCutsNamesBins","Tracks in different cut conditions;;N_{tracks}", lNumberOfCutsToRemember,-0.5,lNumberOfCutsToRemember-0.5);
 402         for(Int_t i = 1; i <= lNumberOfCutsToRemember; i++)
 403         {
 404             //            TString cutsName = (TString*)fArrCutsNames[i-1];
 405             //            fHistCutsNamesBins->GetXaxis()->SetBinLabel(i,((TString*)fArrCutsNames[i-1])->Data());
 406             fHistCutsNamesBins->GetXaxis()->SetBinLabel(i,fArrCutsNames[i-1].Data());
 407
 408         }
 409         fOutList->Add(fHistCutsNamesBins);
 410     }
 411
 412
 413     //Event level QA
 414     const Int_t nEventStatBins = 10;
 415     fHistEventCutStats = new TH1I("fHistEventCutStats","Event statistics;;N_{events}", nEventStatBins,-0.5,nEventStatBins-0.5);
 416     TString gEventCutBinNames[nEventStatBins] = {"Total","No trigger","Wrong centrality","No reco vertex","Bad vertex params", "Bad vertex position","Few SPD tracklets","HighMult cut","LowMult cut","Analyzed"};
 417     for(Int_t i = 1; i <= nEventStatBins; i++)fHistEventCutStats->GetXaxis()->SetBinLabel(i,gEventCutBinNames[i-1].Data());
 418     fOutList->Add(fHistEventCutStats);
 419
 420     //Track level QA
 421     const Int_t nTracksStatBins = 6;
 422     fHistTrackCutStats = new TH1I("fHistTrackCutStats","Track statistics;;N_{tracks}", nTracksStatBins,-0.5,nTracksStatBins-0.5);
 423     TString gTrackCutBinNames[nTracksStatBins] = {"Total","QA track cuts","HighPt cut","LowPt cut","Good", "Not phys primary"};
 424     for(Int_t i = 1; i <= nTracksStatBins; i++)fHistTrackCutStats->GetXaxis()->SetBinLabel(i,gTrackCutBinNames[i-1].Data());
 425     fOutList->Add(fHistTrackCutStats);
 426
 427     //AOD track cut bits stats
 428     const int nAODtrackStats = 16;
 429     fHistAODTrackStats = new TH1I("fHistAODTrackStats","AOD tracks statistics;TrackFilterBit;N_{tracks}",nAODtrackStats,-0.5,nAODtrackStats-0.5);
 430     if ( fAnalysisLevel == "AOD" )
 431         fOutList->Add(fHistAODTrackStats);
 432
 433
 434     //Vertex distributions
 435     fHistVx = new TH1D("fHistVx","Primary vertex distribution - x coordinate;V_{x} (cm);Entries",200,-0.4,0.4);
 436     fOutList->Add(fHistVx);
 437     fHistVy = new TH1D("fHistVy","Primary vertex distribution - y coordinate;V_{y} (cm);Entries",200,-0.4,0.4);
 438     fOutList->Add(fHistVy);
 439     fHistVz = new TH1D("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Entries",200,-20.,20.);
 440     fOutList->Add(fHistVz);
 441
 442     if ( fRunKine )
 443     {
 444         fHistVxMCrecoDiff = new TH1D("fHistVxMCrecoDiff","Primary vertex mc-reco diff - x coordinate;V_{x} (cm);Entries",200,-0.2,0.2);
 445         fOutList->Add(fHistVxMCrecoDiff);
 446         fHistVyMCrecoDiff = new TH1D("fHistVyMCrecoDiff","Primary vertex mc-reco diff - y coordinate;V_{y} (cm);Entries",200,-0.2,0.2);
 447         fOutList->Add(fHistVyMCrecoDiff);
 448         fHistVzMCrecoDiff = new TH1D("fHistVMCrecoDiffz","Primary vertex mc-reco diff - z coordinate;V_{z} (cm);Entries",200,-10.,10.);
 449         fOutList->Add(fHistVzMCrecoDiff);
 450     }
 451
 452     fHistVertexNconributors = new TH1I("fHistVertexNconributors","Primary vertex n contributors;N contributors;Entries",101,-0.5,100.5);
 453     fOutList->Add(fHistVertexNconributors);
 454
 455     fHistNumberOfPileupVerticesTracks = new TH1I("fHistNumberOfPileupVerticesTracks","Number of pilup verteces (by tracks);N verteces;Entries",11,-0.5,10.5);
 456     fOutList->Add(fHistNumberOfPileupVerticesTracks);
 457
 458     fHistNumberOfPileupVerticesSPD = new TH1I("fHistNumberOfPileupVerticesSPD","Number of pilup verteces (by SPD);N verteces;Entries",11,-0.5,10.5);
 459     fOutList->Add(fHistNumberOfPileupVerticesSPD);
 460
 461
 462     if ( fIsIonsAnalysis )
 463     {
 464         //Event plane
 465         fHistEventPlane = new TH2F("fHistEventPlane",";#Psi, rad.;Centrality percentile;Counts",80,0,TMath::Pi()/*0,360.*/,101,-0.5,100.5);
 466         fOutList->Add(fHistEventPlane);
 467     }
 468     //pt, eta, phi checkplots
 469     fHistPt = new TH1F("fHistPt", "p_{T} distribution", 120, 0.0, 6.0);
 470     fHistPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 471     fHistPt->GetYaxis()->SetTitle("dN/dp_{T} (c/GeV)");
 472     fHistPt->SetMarkerStyle(kFullCircle);
 473     fOutList->Add(fHistPt);
 474
 475     fHistEta = new TH1F("fHistEta", "#eta distribution", 200, -4, 4);
 476     fHistEta->GetXaxis()->SetTitle("#eta");
 477     fHistEta->GetYaxis()->SetTitle("dN/#eta");
 478     fHistEta->SetMarkerStyle(kFullCircle);
 479     fOutList->Add(fHistEta);
 480
 481     fHistEtaAODpure = new TH1F("fHistEtaAODpure", "#eta distribution AOD before cuts;#eta;dN/d#eta", 200, -4, 4);
 482     fOutList->Add(fHistEtaAODpure);
 483
 484     fHistPhi = new TH1F("fHistPhi", "#phi distribution", 200, 0, 2*TMath::Pi() );
 485     fHistPhi->GetXaxis()->SetTitle("#phi");
 486     fHistPhi->GetYaxis()->SetTitle("dN/#phi");
 487     fHistPhi->SetMarkerStyle(kFullCircle);
 488     fOutList->Add(fHistPhi);
 489
 490     fHistEtaPhi = new TH2D("fHistEtaPhi","N tracks in (#eta, #phi);#eta;#phi",25,-0.8,0.8,25,0,2*TMath::Pi());
 491     fOutList->Add(fHistEtaPhi);
 492
 493     fHistPhiLRCrotationsCheck = new TH1F("fHistPhiLRCrotationsCheck", "#phi distribution in LRC rotations;#phi;dN/#phi", 200, -4*TMath::Pi(), 4*TMath::Pi() );
 494     fOutList->Add(fHistPhiLRCrotationsCheck);
 495
 496     //some histos for tracks manipulations by hand
 497     fHistPhiArtificialProfilingCheck = new TH1F("fHistPhiArtificialProfilingCheck", "#phi distribution tracks distr by hand;#phi;dN/#phi", 200, -4*TMath::Pi(), 4*TMath::Pi() );
 498     fHistPhiArtificialProfilingCheckWrtEvPlane = new TH1F("fHistPhiArtificialProfilingCheckWrtEvPlane", "#phi distribution by hand wrt event plane;#phi;dN/#phi", 200, -4*TMath::Pi(), 4*TMath::Pi() );
 499     fHistPhiArtificialEvPlane = new TH1F("fHistPhiArtificialEvPlane", "#phi distribution of event plane;#phi;dN/#phi", 200, -4*TMath::Pi(), 4*TMath::Pi() );
 500     if ( fUsePhiShufflingByHand )
 501     {
 502         fOutList->Add(fHistPhiArtificialProfilingCheck);
 503         fOutList->Add(fHistPhiArtificialProfilingCheckWrtEvPlane);
 504         fOutList->Add(fHistPhiArtificialEvPlane);
 505     }
 506
 507     //Eta vs Zv
 508     fHistEtaVsZvCoverage = new TH2D("fHistEtaVsZvCoverage","TPC tracks #eta vs Zv;V_{z} (cm);#eta",100,-20,20,50,-2,2);
 509     fHistEtaVsZvCoverageAccepted = new TH2D("fHistEtaVsZvCoverageAccepted","Accepted TPC tracks #eta vs Zv;V_{z} (cm);#eta",100,-20,20,50,-2,2);
 510     fOutList->Add(fHistEtaVsZvCoverage);
 511     fOutList->Add(fHistEtaVsZvCoverageAccepted);
 512
 513     int lMaxAcceptedTracksInHist = fIsIonsAnalysis ? 4001 : 101;
 514     fHistMultBeforeCuts = new TH1D("fHistMultBeforeCuts","N_{ch} - tracks;N_{ch} tracks;Entries"
 515                                    ,lMaxAcceptedTracksInHist,-0.5,lMaxAcceptedTracksInHist-0.5);
 516     fOutList->Add(fHistMultBeforeCuts);
 517
 518     fHistAcceptedMult = new TH1D("fHistAcceptedMult","N_{ch} - accepted tracks;N_{ch} accepted;Entries"
 519                                  ,lMaxAcceptedTracksInHist,-0.5,lMaxAcceptedTracksInHist-0.5);
 520     fOutList->Add(fHistAcceptedMult);
 521
 522     fHistAcceptedTracks = new TH1D("fHistAcceptedTracks","N_{ch} - accepted tracks for LRC;N_{ch} accepted;Entries"
 523                                    ,lMaxAcceptedTracksInHist,-0.5,lMaxAcceptedTracksInHist-0.5);
 524     fOutList->Add(fHistAcceptedTracks);
 525
 526     fHistMultiplicityInEtaRegion = new TH1D("fHistMultiplicityInEtaRegion","N_{ch} in |#eta| region;N_{ch};Entries"
 527                                             ,lMaxAcceptedTracksInHist,-0.5,lMaxAcceptedTracksInHist-0.5);
 528     fOutList->Add(fHistMultiplicityInEtaRegion);
 529
 530     fHistMultiplicityInEtaRegionAfterPtCuts = new TH1D("fHistMultiplicityInEtaRegionAfterPtCuts","N_{ch} in |#eta| region after p_{T} cuts;N_{ch};Entries"
 531                                                        ,lMaxAcceptedTracksInHist,-0.5,lMaxAcceptedTracksInHist-0.5);
 532     fOutList->Add(fHistMultiplicityInEtaRegionAfterPtCuts);
 533
 534     if ( fAnalyseMCESD )
 535     {
 536         fHist2DMultiplicityMCESDInEtaRegion = new TH2D("fHist2DMultiplicityMCESDInEtaRegion","N_{ch} in |#eta| region;N_{ch} MC;N_{ch} ESD;Entries"
 537                                                        ,101,-0.5,100.5,101,-0.5,100.5);
 538         fOutList->Add(fHist2DMultiplicityMCESDInEtaRegion);
 539     }
 540
 541
 542     fHistAcceptedTracksAfterPtCuts = new TH1D("fHistAcceptedTracksAfterPtCuts","N_{ch} - accepted tracks for LRC after Pt cuts;N_{ch} accepted;Entries"
 543                                               ,lMaxAcceptedTracksInHist,-0.5,lMaxAcceptedTracksInHist-0.5);
 544     fOutList->Add(fHistAcceptedTracksAfterPtCuts);
 545
 546     fHistAcceptedTPCtracks = new TH1D("fHistAcceptedTPCtracks","N_{ch} - accepted tracks with TPC inner param;N_{ch} accepted;Entries"
 547                                       ,lMaxAcceptedTracksInHist,-0.5,lMaxAcceptedTracksInHist-0.5);
 548     fOutList->Add(fHistAcceptedTPCtracks);
 549
 550     fHistClustersTPC = new TH1D("fHistClustersTPC","N Clusters TPC;N_{TPC clusters};Entries",161,-0.5,160.5);
 551     fOutList->Add(fHistClustersTPC);
 552
 553     fHistClustersTPCafterCuts = new TH1D("fHistClustersTPCafterCuts","N Clusters TPC after cuts;N_TPC_clusters_after_cuts;Entries",161,-0.5,160.5);
 554     fOutList->Add(fHistClustersTPCafterCuts);
 555
 556
 557     fHistCrossedRowsTPC = new TH1D("fHistCrossedRowsTPC","N Crossed Rows TPC;N_{TPC CrossedRows};Entries",161,-0.5,160.5);
 558     fOutList->Add(fHistCrossedRowsTPC);
 559
 560     fHistCrossedRowsTPCafterCuts = new TH1D("fHistCrossedRowsTPCafterCuts","N CrossedRows TPC after cuts;N_TPC_clusters_after_cuts;Entries",161,-0.5,160.5);
 561     fOutList->Add(fHistCrossedRowsTPCafterCuts);
 562
 563
 564
 565     fHistTrackletsITS = new TH1D("fHistTrackletsITS","N Tracklets ITS;N_ITS_tracklets;Entries",101,-0.5,100.5);
 566     fOutList->Add(fHistTrackletsITS);
 567
 568
 569     fHistClustersITS = new TH1D("fHistClustersITS","N Clusters ITS;N_ITS_clusters;Entries",11,-0.5,10.5);
 570     fOutList->Add(fHistClustersITS);
 571
 572
 573     fHist2DClustersTPCvsPt = new TH2D("fHist2DClustersTPCvsPt","Num TPC clusters vs Pt;P_t;N_clusters",50,0,2,161,-0.5,160.5);
 574     fOutList->Add(fHist2DClustersTPCvsPt);
 575
 576     fHist2DClustersTPCvsEta = new TH2D("fHist2DClustersTPCvsEta","Num TPC clusters vs Eta;#eta;N_clusters",50,-1.5,1.5,161,-0.5,160.5);
 577     fOutList->Add(fHist2DClustersTPCvsEta);
 578
 579     fHist2DAcceptedTracksPtvsEta = new TH2D("fHist2DAcceptedTracksPtvsEta","Accepted Tracks Pt vs Eta;pt;#eta",50,0,2,50,-1.5,1.5);
 580     fOutList->Add(fHist2DAcceptedTracksPtvsEta);
 581
 582     fHistProbabilitiesPID = new TH1D("fHistProbabilitiesPID","PID Probabilities;pid;Entries",10,-0.5,9.5 );//AliPID::kSPECIES,-0.5,AliPID::kSPECIES-0.5);
 583     fOutList->Add(fHistProbabilitiesPID);
 584
 585     fHistESDtrackMass = new TH1D("fHistESDtrackMass","ESD Mass;Mass;Entries",100,0,1);
 586     fOutList->Add(fHistESDtrackMass);
 587
 588     fHistProbabilityPion = new TH1D("fHistProbabilityPion","Probability Pion;Probability;Entries",100,0,1);
 589     fOutList->Add(fHistProbabilityPion);
 590
 591     fHistProbabilityKaon = new TH1D("fHistProbabilityKaon","Probability Kaon;Probability;Entries",100,0,1);
 592     fOutList->Add(fHistProbabilityKaon);
 593
 594     fHistProbabilityProton = new TH1D("fHistProbabilityProton","Probability Proton;Probability;Entries",100,0,1);
 595     fOutList->Add(fHistProbabilityProton);
 596
 597
 598
 599     fHistParticlesDistrBeforeCuts = new TH1D("fHistParticlesDistrBeforeCuts","Particles Distr;Particle;Entries",5,-0.5,4.5);
 600     TString gBinParticleNames[5] = {/*"Other",*/"Electron","Muon","Pion","Kaon", "Proton"};
 601     for(Int_t i = 1; i <= 5; i++)fHistParticlesDistrBeforeCuts->GetXaxis()->SetBinLabel(i,gBinParticleNames[i-1].Data());
 602     fOutList->Add(fHistParticlesDistrBeforeCuts);
 603
 604     fHistParticlesDistr = new TH1D("fHistParticlesDistr","Particles Distr;Particle;Entries",5,-0.5,4.5);
 605     //TString gBinParticleNames[6] = {"Undefined","Electron","Muon","Pion","Kaon", "Proton"};
 606     for(Int_t i = 1; i <= 5; i++)fHistParticlesDistr->GetXaxis()->SetBinLabel(i,gBinParticleNames[i-1].Data());
 607     fOutList->Add(fHistParticlesDistr);
 608
 609     fHistCentralityPercentile = new TH1D("fHistCentralityPercentile","Centrality Percentile;Centrality;Entries",101,-1.01,100.01);
 610     if ( fIsIonsAnalysis ) fOutList->Add(fHistCentralityPercentile);
 611
 612     fHistCentralityClass10 = new TH1D("fHistCentralityClass10","Centrality Class 10;Centrality;Entries",10,-0.5,9.5);
 613     if ( fIsIonsAnalysis ) fOutList->Add(fHistCentralityClass10);
 614
 615     fHistCentralityClass5 = new TH1D("fHistCentralityClass5","Centrality Class 5;Centrality;Entries",20,-0.5,19.5);
 616     if ( fIsIonsAnalysis ) fOutList->Add(fHistCentralityClass5);
 617
 618
 619     fRand = new TRandom3;
 620     //ZDC additional study
 621
 622     fMultForZDCstudy[0] = 0;
 623     fMultForZDCstudy[1] = 200;
 624     fMultForZDCstudy[2] = 220;
 625     fMultForZDCstudy[3] = 240;
 626     fMultForZDCstudy[4] = 250;
 627
 628     TString strZDCenergyName;
 629     TString strZDCenergyTitle;
 630     for ( int i = 0; i < 5; i++ )
 631     {
 632         strZDCenergyName = Form( "fHistZDCenergy%d", fMultForZDCstudy[i] );
 633         strZDCenergyTitle = Form(  "ZDC Energy, multThreshold %d;Energy;Entries", fMultForZDCstudy[i] );
 634         fHistZDCenergy[i] = new TH1D(strZDCenergyName, strZDCenergyTitle, 200, 0, 10000);
 635         if ( fFlagWatchZDC && fIsIonsAnalysis )
 636             fOutList->Add(fHistZDCenergy[i]);
 637     }
 638
 639
 640     fHistZDCparticipants = new TH1D("fHistZDCparticipants","ZDC Participants;Participants;Entries",800,0,800);
 641     if ( fFlagWatchZDC && fIsIonsAnalysis )
 642         fOutList->Add(fHistZDCparticipants);
 643
 644     if ( fFlagWatchV0 )//fIsIonsAnalysis )
 645     {
 646         const int nBinsV0multForNonIons = 700;
 647         const int nBinsV0multForIons = 220;
 648         const int nMaxV0multForIons = 22000;
 649         fHistV0multiplicity = new TH1D("fHistV0multiplicity","V0 Multiplicity;Multiplicity;Entries",
 650                                        fIsIonsAnalysis ? nBinsV0multForIons : nBinsV0multForNonIons, 0, fIsIonsAnalysis ? nMaxV0multForIons : nBinsV0multForNonIons);
 651         fHistV0Amultiplicity = new TH1D("fHistV0Amultiplicity","V0-A Multiplicity;Multiplicity;Entries",
 652                                         fIsIonsAnalysis ? nBinsV0multForIons : nBinsV0multForNonIons, 0, fIsIonsAnalysis ? nMaxV0multForIons : nBinsV0multForNonIons);
 653         fHistV0Cmultiplicity = new TH1D("fHistV0Cmultiplicity","V0-C Multiplicity;Multiplicity;Entries",
 654                                         fIsIonsAnalysis ? nBinsV0multForIons : nBinsV0multForNonIons, 0, fIsIonsAnalysis ? nMaxV0multForIons : nBinsV0multForNonIons);
 655         fHist2DV0ACmultiplicity = new TH2D("fHist2DV0ACmultiplicity","V0 A-C Multiplicity;Multiplicity A;Multiplicity C;Entries"
 656                                            ,50, 0, fIsIonsAnalysis ? nMaxV0multForIons/2 : nBinsV0multForNonIons/2
 657                                                                      ,50, 0, fIsIonsAnalysis ? nMaxV0multForIons/2 : nBinsV0multForNonIons/2);
 658
 659         fHist2DTracksAcceptedVsV0multiplicity = new TH2D("fHist2DTracksAcceptedVsV0multiplicity","V0 Multiplicity vs Accepted;N Accepted tracks;Multiplicity V0;Entries"
 660                                                          ,70, 0, fIsIonsAnalysis ? 3000 : 70
 661                                                                                    ,50, 0, fIsIonsAnalysis ? nMaxV0multForIons/2 : 100 );
 662
 663         //number of fired cells
 664         const int nV0cells = 32;
 665         //fHistV0cells = new TH1D("fHistV0cells","V0 cells;N cells;Entries", nV0cells+1, -0.5, nV0cells + 0.5 );
 666         fHistV0Acells = new TH1D("fHistV0Acells","V0-A cells;N cells;Entries", nV0cells+1, -0.5, nV0cells + 0.5 );
 667         fHistV0Ccells = new TH1D("fHistV0Ccells","V0-C cells;N cells;Entries", nV0cells+1, -0.5, nV0cells + 0.5 );
 668         fHist2DV0ACcells = new TH2D("fHist2DV0ACcells","V0 A-C cells;N cells A;N cells C;Entries", nV0cells+1, -0.5, nV0cells + 0.5, nV0cells+1, -0.5, nV0cells + 0.5 );
 669
 670
 671         fOutList->Add(fHistV0multiplicity);
 672         fOutList->Add(fHistV0Amultiplicity);
 673         fOutList->Add(fHistV0Cmultiplicity);
 674         fOutList->Add(fHist2DV0ACmultiplicity);
 675         fOutList->Add(fHist2DTracksAcceptedVsV0multiplicity);
 676         //cells
 677         //fOutList->Add(fHistV0cells    );
 678         fOutList->Add(fHistV0Acells   );
 679         fOutList->Add(fHistV0Ccells   );
 680         fOutList->Add(fHist2DV0ACcells);
 681
 682         //mult in rings
 683         //        char strV0ringName[200];
 684         //        char strV0ringTitle[200];
 685         TString strV0ringName;
 686         TString strV0ringTitle;
 687         for ( int i = 0; i < 4; i++ )
 688         {
 689             strV0ringName =  Form( "fHistV0AmultiplicityRing%d", i );
 690             strV0ringTitle = Form( strV0ringTitle, "V0-A Multiplicity Ring %d;Multiplicity;Entries", i );
 691             fHistV0AmultiplicityRing[i] = new TH1D( strV0ringName, strV0ringTitle,
 692                                                     fIsIonsAnalysis ? nBinsV0multForIons : nBinsV0multForNonIons, 0, fIsIonsAnalysis ? nMaxV0multForIons : nBinsV0multForNonIons);
 693             strV0ringName =  Form(  "fHistV0CmultiplicityRing%d", i );
 694             strV0ringTitle = Form( strV0ringTitle, "V0-C Multiplicity Ring %d;Multiplicity;Entries", i );
 695             fHistV0CmultiplicityRing[i] = new TH1D( strV0ringName, strV0ringTitle,
 696                                                     fIsIonsAnalysis ? nBinsV0multForIons : nBinsV0multForNonIons, 0, fIsIonsAnalysis ? nMaxV0multForIons : nBinsV0multForNonIons);
 697             strV0ringName =  Form( "fHist2DV0ACmultiplicityRing%d", i );
 698             strV0ringTitle = Form(  "V0-AC Multiplicity Ring %d;Multiplicity A;Multiplicity C;Entries", i );
 699             fHist2DV0ACmultiplicityRing[i] = new TH2D( strV0ringName, strV0ringTitle
 700                                                        , 100, 0, fIsIonsAnalysis ? nMaxV0multForIons/2 : 100
 701                                                                                    , 100, 0, fIsIonsAnalysis ? nMaxV0multForIons/2 : 100 );
 702             fOutList->Add(fHistV0AmultiplicityRing[i]);
 703             fOutList->Add(fHistV0CmultiplicityRing[i]);
 704             fOutList->Add(fHist2DV0ACmultiplicityRing[i]);
 705
 706             //mult in barrel vs V0 rings
 707             strV0ringName =  Form( "fHist2DTracksAcceptedVsV0AmultiplicityRing%d", i );
 708             strV0ringTitle = Form(  "Accepted tracks vs V0-A Multiplicity in Ring %d;N Accepted tracks;Multiplicity V0A;Entries", i );
 709             fHist2DTracksAcceptedVsV0AmultiplicityRing[i] = new TH2D( strV0ringName, strV0ringTitle
 710                                                                       , 100, 0, fIsIonsAnalysis ? nMaxV0multForIons/2 : 100
 711                                                                                                   , 100, 0, fIsIonsAnalysis ? nMaxV0multForIons/2 : 100 );
 712             strV0ringName =  Form( "fHist2DTracksAcceptedVsV0CmultiplicityRing%d", i );
 713             strV0ringTitle = Form(  "Accepted tracks vs V0-C Multiplicity in Ring %d;N Accepted tracks;Multiplicity V0C;Entries", i );
 714             fHist2DTracksAcceptedVsV0CmultiplicityRing[i] = new TH2D( strV0ringName, strV0ringTitle
 715                                                                       , 100, 0, fIsIonsAnalysis ? nMaxV0multForIons/2 : 100
 716                                                                                                   , 100, 0, fIsIonsAnalysis ? nMaxV0multForIons/2 : 100 );
 717
 718             fOutList->Add(fHist2DTracksAcceptedVsV0AmultiplicityRing[i]);
 719             fOutList->Add(fHist2DTracksAcceptedVsV0CmultiplicityRing[i]);
 720         }
 721
 722     }
 723
 724
 725
 726
 727     //    fHistV0spectra = new TH1D("fHistV0spectra","V0 spectra;Mass, GeV;Entries",500,0,5000);
 728     //    fOutList->Add(fHistV0spectra);
 729
 730
 731     fHistMClabels = new TH1D("fHistMClabels","MC label;label;Entries",102,-100.5,100.5);
 732     fOutList->Add(fHistMClabels);
 733
 734     fHistRejectedTracksCharge = new TH1D("fHistRejectedTracksCharge","Rejected tracks charge;charge;Entries",3,-1.5,1.5);
 735     fOutList->Add(fHistRejectedTracksCharge);
 736
 737     fHistTracksCharge = new TH1D("fHistTracksCharge","Accepted tracks charge;charge;Entries",3,-1.5,1.5);
 738     fOutList->Add(fHistTracksCharge);
 739
 740     // ##### net charge vs pt study
 741     const int kPtNetChargePtBins = 200;
 742     const double kPtNetChargePtMin = 0.1;
 743     const double kPtNetChargePtMax = 2.1;
 744
 745     //pt plus
 746     fHistPtPlus = new TH1D("fHistPtPlus","p_{T} +;p_{T};dN/dpT"
 747                                   , kPtNetChargePtBins,  kPtNetChargePtMin, kPtNetChargePtMax
 748                                   );
 749     fOutList->Add(fHistPtPlus);
 750
 751     //pt minus
 752     fHistPtMinus = new TH1D("fHistPtMinus","p_{T} -;p_{T};dN/dpT"
 753                                   , kPtNetChargePtBins,  kPtNetChargePtMin, kPtNetChargePtMax
 754                                   );
 755     fOutList->Add(fHistPtMinus);
 756
 757     //net charge vs pT
 758     fHistNetChargeVsPt = new TH1D("fHistNetChargeVsPt","charge vs p_{T};p_{T};Q"
 759                                   , kPtNetChargePtBins,  kPtNetChargePtMin, kPtNetChargePtMax
 760                                   );
 761     fOutList->Add(fHistNetChargeVsPt);
 762
 763     //plus
 764     fHistChargePlusVsPtTmp = new TH1D( "fHistChargePlusVsPtTmp","charge plus vs p_{T};p_{T};q plus"
 765                                        , kPtNetChargePtBins,  kPtNetChargePtMin, kPtNetChargePtMax
 766                                        );
 767     //    fOutList->Add(fHistChargePlusVsPtTmp);
 768
 769     //minus
 770     fHistChargeMinusVsPtTmp = new TH1D( "fHistChargeMinusVsPtTmp","charge minus vs p_{T};p_{T};q minus"
 771                                         , kPtNetChargePtBins,  kPtNetChargePtMin, kPtNetChargePtMax
 772                                         );
 773     //    fOutList->Add(fHistChargeMinusVsPtTmp);
 774
 775     //NetChargeVsPt
 776     fHist2DNetChargeVsPt = new TH2D( "fHist2DNetChargeVsPt","Net charge vs p_{T};p_{T};Q"
 777                                      , kPtNetChargePtBins,  kPtNetChargePtMin, kPtNetChargePtMax
 778                                      , 40,  -20, 20
 779                                      );
 780     fOutList->Add(fHist2DNetChargeVsPt);
 781
 782     //NetChargeVsPt CorrectedOnEventMean
 783     fHist2DNetChargeVsPtCorrectedOnEventMean = new TH2D( "fHist2DNetChargeVsPtCorrectedOnEventMean","Net charge corrected on mean e-by-e vs p_{T};p_{T};Q"
 784                                                          , kPtNetChargePtBins,  kPtNetChargePtMin, kPtNetChargePtMax
 785                                                          , 40,  -20, 20
 786                                                          );
 787     fOutList->Add(fHist2DNetChargeVsPtCorrectedOnEventMean);
 788
 789     //NetChargeVsPt CorrectedOnEventMean normalized on Nch
 790     fHist2DNetChargeVsPtCorrectedOnEventMeanNormOnNch = new TH2D( "fHist2DNetChargeVsPtCorrectedOnEventMeanNormOnNch","Net charge vs p_{T} corrected on mean e-by-e normalized on nCh;p_{T};Q"
 791                                                                   , kPtNetChargePtBins,  kPtNetChargePtMin, kPtNetChargePtMax
 792                                                                   , 40,  -20, 20
 793                                                                   );
 794     fOutList->Add(fHist2DNetChargeVsPtCorrectedOnEventMeanNormOnNch);
 795
 796
 797
 798     if ( fArtificialInefficiency >= 0 ) //i.e. have this kind of analysis
 799     {
 800         fHistNumberOfDroppedByHandTracks = new TH2D("fHistNumberOfDroppedByHandTracks","Accepted tracks vs Dropped artificially;N_{ch} accepted;N_{ch} dropped;Entries", 71, -0.5, 70.5,   71, -0.5, 70.5);
 801         fOutList->Add(fHistNumberOfDroppedByHandTracks);
 802     }
 803
 804
 805     // Returning TH1,TH2 AddDirectory status
 806     TH1::AddDirectory(lTH1oldStatus);
 807     TH2::AddDirectory(lTH2oldStatus);
 808
 809     //fHistProbabilitiesPID->Fill(1,fStrPIDforFwd);
 810     //fHistProbabilitiesPID->Fill(2,fStrPIDforBwd);
 811     Printf("UserCreateOutputObjects done.");
 812
 813     // NEW HISTO added to fOutput here
 814     PostData(1, fOutList); // Post data for ALL output slots >0 here, to get at least an empty histogram
 815
 816     for( Int_t i = 0; i < fLRCproc.GetEntries(); i++)
 817     {
 818         PostData( Proc(i)->GetOutputSlotNumber(),Proc(i)->CreateOutput() );
 819     }
 820
 821
 822     //if ( fSetIncludeEventTreeInOutput )
 823     //    PostData(2, fEventTree);
 824     //int a;
 825     //cin >> a;
 826
 827
 828 }
 829
 830 //________________________________________________________________________
 831 //void AliAnalysisTaskLRC::UserExec(Option_t *)
 832 //{
 833
 834 //    //cout << "TEST in Task" << endl;
 835 //    // ###### tuning phi
 836
 837 //    //rotate phi when N_phi_sectors > 1 and call UserExecLoop:
 838 //    double phiStep = 2 * TMath::Pi();
 839 //    if ( fNumberOfSectors > 1 )
 840 //        phiStep /= fNumberOfSectors;
 841
 842 //    double lPhiRotatedExtra = 0; //additional phi rotation of all tracks
 843 //    for ( Int_t sectorId = 0; sectorId < fNumberOfSectors; sectorId++ )
 844 //    {
 845 //        //cout << "loop " << sectorId << endl;
 846 //        UserExecLoop( lPhiRotatedExtra );
 847 //        lPhiRotatedExtra += phiStep; //rotate track
 848 //    }
 849
 850
 851 //}
 852
 853 //________________________________________________________________________
 854 void AliAnalysisTaskLRC::UserExec(Option_t *)   //UserExecLoop( Double_t phiAdditional )//Option_t *)
 855 {
 856     // ########### if use toy events
 857     if ( fUseToyEvents )
 858     {
 859         //generate all events here and fill LRC processors
 860         UseToyEvents();
 861
 862         //just return
 863         return;
 864     }
 865
 866     // Main loop
 867     // Called for each event
 868     //printf( "starting UserExec...\n" );
 869     //Pid setting
 870     if ( fPIDsensingFlag )
 871     {
 872         SetParticleTypeToProcessors( 0 /*fwd*/, fStrPIDforFwd );
 873         SetParticleTypeToProcessors( 1 /*backward*/, fStrPIDforBwd );
 874         fPIDsensingFlag = kFALSE;
 875     }
 876
 877     //printf( "starting event...\n" );
 878     //The common PID object can then be retrieved from the input handler:
 879     AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
 880     AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 881     fPIDResponse = inputHandler->GetPIDResponse();
 882     //if (!fPIDResponse) !!!!!!!!!!!!
 883     //  AliFatal("This Task needs the PID response attached to the inputHandler");
 884
 885
 886     AliVEvent *event = InputEvent();
 887     AliESDEvent *lESD = 0x0;
 888     AliAODEvent *lAOD = 0x0;
 889
 890     AliStack *stack = 0x0;
 891     AliMCEvent *eventMC = 0x0;
 892
 893
 894     if ( fAnalysisLevel == "ESD" ) //fAnalysisType == en_AnalysisType_AOD )
 895     {
 896         lESD = dynamic_cast<AliESDEvent*> (event) ;
 897         //printf( "cast to ESD is Ok\n." );
 898     }
 899     else if ( fAnalysisLevel == "AOD" ) //fAnalysisType == en_AnalysisType_AOD )
 900     {
 901         //cout << "test AOD analysis... " << endl;
 902         lAOD = dynamic_cast<AliAODEvent*>(event); // from TaskSE
 903         //Int_t bMagneticFieldSign = (lAOD->GetMagneticField() > 0) ? 1 : -1;
 904         //printf( "Number of AOD tracks = %d, magnetic field = %f\n", lAOD->GetNumberOfTracks(), lAOD->GetMagneticField() );
 905         //return;
 906
 907     }
 908
 909     if( fRunKine ) //|| fAnalyseMCESD )
 910     {
 911         eventMC = MCEvent();
 912         stack = eventMC->Stack();
 913         //Printf("Number of primaries: %d",stack->GetNprimary());
 914     }
 915
 916     if (!event)
 917     {
 918         //Printf("ERROR: Could not retrieve event");
 919         return;
 920     }
 921     //cout << "test start Event!" << endl;
 922
 923     if( (lESD ) && ( !fEsdTrackCuts ) )
 924     {
 925         AliDebug(AliLog::kError, "No ESD track cuts avalible");
 926         return;
 927     }
 928
 929     //make combined flag for esd selection during MC kine analysis
 930     Bool_t lRunKine = ( fRunKine && !fAnalyseMCESD );
 931
 932     // Processing event selection
 933     fHistEventCutStats->Fill( "Total", 1 );
 934
 935
 936
 937     //cout << "test: nOfTracks = " << event->GetNumberOfTracks() << endl;
 938
 939     //Trigger
 940     Bool_t lTrigger = kTRUE;
 941     if( lESD && !lRunKine ) //just for tests; usually it is done before UserExecLoop in PhysSel task
 942     {
 943         if(fShowEventStats)
 944             Printf("Trigger classes: %s:", lESD->GetFiredTriggerClasses().Data());
 945
 946         lTrigger = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 947
 948         //8 TeV runs: TEST triggers! (29.12.2012)
 949         //        Bool_t lTrigger1 = (lESD->IsTriggerClassFired("CINT7-B-NOPF-ALLNOTRD")) ? 1 : 0;
 950         //        Bool_t lTrigger2 = (lESD->IsTriggerClassFired("CINT7WU-B-NOPF-ALL")) ? 1 : 0;
 951         //        lTrigger = lTrigger1 && lTrigger2;
 952
 953         if( !lTrigger )
 954         {
 955             if( fShowEventStats )
 956                 Printf("Rejected!");
 957             fHistEventCutStats->Fill( "No trigger", 1 );
 958             PostData(1, fOutList);
 959             return;
 960         }
 961     }
 962     //Centrality      //https://twiki.cern.ch/twiki/bin/viewauth/ALICE/CentStudies
 963     Double_t lCentralityPercentile      = 0.;
 964     Int_t lCentralityClass10            = 0;
 965     Int_t lCentralityClass5             = 0;
 966     Float_t lReactionPlane       = -1.;
 967
 968
 969     if( !lRunKine && fIsIonsAnalysis )
 970     {
 971         AliCentrality *centrality = 0x0;
 972
 973         if ( lESD )
 974             centrality = lESD->GetCentrality();
 975         else if ( lAOD )
 976             centrality = lAOD->GetCentrality();
 977
 978         if ( centrality )
 979         {
 980             lCentralityPercentile       = centrality->GetCentralityPercentile("V0M");   // returns the centrality percentile, a float from 0 to 100 (or to the trigger efficiency)
 981             lCentralityClass10          = centrality->GetCentralityClass10("V0M");//"FMD");             // returns centrality class for 10% width (a integer from 0 to 10)
 982             lCentralityClass5           = centrality->GetCentralityClass5("V0M");//"TKL");              // returns centrality class for 5% width (a integer from 0 to 20)
 983             Bool_t lCentralityInClass = centrality->IsEventInCentralityClass(fMinCentralityClass,fMaxCentralityClass,"V0M"); // returns kTRUE if the centrality of the event is between a and b, otherwise kFALSE
 984             //cout << lCentralityPercentile << " "
 985             //<< fMinCentralityClass << " "
 986             //<< fMaxCentralityClass << " "
 987             //<< lCentralityInClass << endl;
 988
 989             if ( !lCentralityInClass )
 990             {
 991                 //cout << "outside of centrality class!" << endl;
 992                 fHistEventCutStats->Fill("Wrong centrality", 1);
 993                 PostData(1, fOutList);
 994                 return;
 995             }
 996
 997             fHistCentralityPercentile->Fill( lCentralityPercentile );
 998             fHistCentralityClass10->Fill( lCentralityClass10 );
 999             fHistCentralityClass5->Fill( lCentralityClass5 );
1000
1001             // get the reaction plane
1002             lReactionPlane = GetEventPlane( event );
1003             fHistEventPlane->Fill( lReactionPlane, lCentralityPercentile );
1004         }
1005     }
1006
1007     //number of verteces in ESD (pileup)
1008     if ( lESD )
1009     {
1010         int lNumberOfPileUpVerteces = 0;
1011         TClonesArray *lPileupVertecesTracks = lESD->GetPileupVerticesTracks();
1012         lNumberOfPileUpVerteces = lPileupVertecesTracks->GetSize();
1013         fHistNumberOfPileupVerticesTracks->Fill( lNumberOfPileUpVerteces );
1014
1015         TClonesArray *lPileupVertecesSPD = lESD->GetPileupVerticesSPD();
1016         lNumberOfPileUpVerteces = lPileupVertecesSPD->GetSize();
1017         fHistNumberOfPileupVerticesSPD->Fill( lNumberOfPileUpVerteces );
1018     }
1019     else if ( lAOD ) //number of verteces in AOD (pileup)
1020     {
1021         int lNumberOfPileUpVerteces = 0;
1022         lNumberOfPileUpVerteces = lAOD->GetNumberOfPileupVerticesTracks();
1023         fHistNumberOfPileupVerticesTracks->Fill( lNumberOfPileUpVerteces );
1024         lNumberOfPileUpVerteces = lAOD->GetNumberOfPileupVerticesSPD();
1025         fHistNumberOfPileupVerticesSPD->Fill( lNumberOfPileUpVerteces );
1026     }
1027
1028
1029     // Vertex present
1030     //    const AliESDVertex *vertex = lESD->GetPrimaryVertex();
1031     Double_t lVertexX = -1000;
1032     Double_t lVertexY = -1000;
1033     Double_t lVertexZ = -1000;
1034     Double_t lMCVertexX = -1000;
1035     Double_t lMCVertexY = -1000;
1036     Double_t lMCVertexZ = -1000;
1037
1038     Bool_t lReconstructedVertexPresent = kFALSE;
1039     Bool_t lVertexAcceptable = kFALSE;
1040
1041     if ( fAnalysisLevel == "ESD" ||  fAnalysisLevel == "AOD" )
1042     {
1043         const AliVVertex *vertex = event->GetPrimaryVertex();
1044         fHistVertexNconributors->Fill( vertex->GetNContributors() );
1045
1046         Double32_t lCov[6];
1047         vertex->GetCovarianceMatrix(lCov);
1048
1049         //check nContributors and z*z>0
1050         lReconstructedVertexPresent = ( ( vertex->GetNContributors() > 0 ) && ( lCov[5] != 0 ) );
1051         //        if ( (vertex) && (vertex->GetNContributors() > 0 ) && ( vertex->GetZRes() == 0 ) )
1052         //        {
1053         //            int aa;
1054         //            cout << "stop with nContributors = " << vertex->GetNContributors() << endl;
1055         //            cin >> aa;
1056         //        }
1057         lVertexX = vertex->GetX();
1058         lVertexY = vertex->GetY();
1059         lVertexZ = vertex->GetZ();
1060     }
1061
1062     if ( fRunKine ) //take MC vertex params and rewrite position
1063     {
1064         AliGenEventHeader *header = eventMC->GenEventHeader();
1065         if(!header) return;
1066
1067         TArrayF gVertexArray;
1068         header->PrimaryVertex(gVertexArray);
1069         lMCVertexX = gVertexArray.At(0);
1070         lMCVertexY = gVertexArray.At(1);
1071         lMCVertexZ = gVertexArray.At(2);
1072
1073         //fill mc-reco diff
1074         if ( lReconstructedVertexPresent )
1075         {
1076             fHistVxMCrecoDiff->Fill( lMCVertexX - lVertexX );
1077             fHistVyMCrecoDiff->Fill( lMCVertexY - lVertexY );
1078             fHistVzMCrecoDiff->Fill( lMCVertexZ - lVertexZ );
1079         }
1080
1081     }
1082
1083     //cut on reco vertex params
1084     if ( !lRunKine )
1085     {
1086         if( ( !lReconstructedVertexPresent ) && fCheckForkVtx )
1087         {
1088             if( fShowEventStats )
1089                 Printf("Bad vertex params");
1090             fHistEventCutStats->Fill("Bad vertex params",1);
1091             PostData(1, fOutList);
1092             return;
1093         }
1094     }
1095     else
1096     {
1097         //rewrite vert x,y,z into mc values
1098         lVertexX = lMCVertexX;
1099         lVertexY = lMCVertexY;
1100         lVertexZ = lMCVertexZ;
1101     }
1102
1103
1104     // Vertex in range
1105     lVertexAcceptable = (TMath::Abs(lVertexX) < fVxMax) && (TMath::Abs(lVertexY) < fVyMax);
1106     if( lVertexAcceptable )
1107     {
1108         if( fVzMax > 0 )   //   fVzMax < 0 -> select Zv outside selected range
1109             lVertexAcceptable = (TMath::Abs(lVertexZ) < fVzMax);
1110         else
1111             lVertexAcceptable = (TMath::Abs(lVertexZ) > -fVzMax);
1112     }
1113
1114     if( (!lVertexAcceptable) && fCheckForVtxPosition && !lRunKine )
1115     {
1116         if(fShowEventStats)
1117             Printf("Vertex out of range");
1118         fHistEventCutStats->Fill("Bad vertex position",1);
1119         PostData(1, fOutList);
1120         return;
1121     }
1122
1123     //fill bin if no reco vertex in MC analysis
1124     if ( lRunKine && !lReconstructedVertexPresent )
1125         fHistEventCutStats->Fill( "No reco vertex", 1 );
1126
1127     fHistVx->Fill(lVertexX);
1128     fHistVy->Fill(lVertexY);
1129     fHistVz->Fill(lVertexZ);
1130     //end with vertex
1131
1132
1133
1134
1135     //cut on number of SPD tracklets (25.03.2012)
1136     if( lESD && !lRunKine )
1137     {
1138         //How to get tracklets
1139         const AliMultiplicity *tracklets = lESD->GetMultiplicity();
1140         Int_t multSPD = tracklets->GetNumberOfTracklets();
1141         //Int_t nITStracklets = kalmanTrack->GetNumberOfTracklets();
1142         fHistTrackletsITS->Fill( multSPD );
1143         if ( multSPD < fMinNumberOfSPDtracklets )
1144         {
1145             if(fShowEventStats)
1146                 Printf("Too few SPD tracklets");
1147             fHistEventCutStats->Fill( "Few SPD tracklets", 1 );
1148             PostData(1, fOutList);
1149             return;
1150         }
1151     }
1152
1153
1154     //tmp eta ranges (for ZDC respond study)
1155     double etaBmin = -0.8;
1156     double etaBmax = -0.6;
1157     double etaFmin = 0.6;
1158     double etaFmax = 0.8;
1159     int countTracksEtaB = 0;
1160     int countTracksEtaF = 0;
1161
1162     //n accepted tracks
1163     int lNchTrigger = 0;
1164     int lTPCtracks  = 0; // added 25.03.2012
1165     fHistMultBeforeCuts->Fill( event->GetNumberOfTracks() );
1166
1167     // Pre event loop
1168     if( !lRunKine )
1169     {
1170         for ( Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++)
1171         {
1172             if ( fAnalysisLevel == "ESD" )
1173             {
1174                 AliESDtrack *track = lESD->GetTrack(iTracks);
1175                 if ( fEsdTrackCuts->AcceptTrack(track) )
1176                 {
1177                     lNchTrigger++;
1178                     //todo: try to implement pt cuts here?....
1179                     //How to get TPC standalone tracks
1180                     AliExternalTrackParam *tpcTrack
1181                             = (AliExternalTrackParam *)track->GetTPCInnerParam();
1182                     if ( tpcTrack )
1183                         lTPCtracks++;
1184                 }
1185                 else
1186                 {
1187                     fHistRejectedTracksCharge->Fill( track->Charge() );
1188                 }
1189             }
1190             else if ( fAnalysisLevel == "AOD" )
1191             {
1192                 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));
1193                 if( aodTrack->TestFilterBit( fAODtrackCutBit ) )
1194                     lNchTrigger++;
1195             }
1196         }
1197     }
1198     fHistAcceptedTPCtracks->Fill( lTPCtracks );
1199
1200
1201
1202
1203     //if(lRunKine)lNchTrigger=eventMC->GetNumberOfPrimaries();
1204
1205     // Nch bins cut
1206     if( ( lNchTrigger > fMaxAcceptedTracksCut ) && ( fMaxAcceptedTracksCut != 0 ) )
1207     {
1208         fHistEventCutStats->Fill( "HighMult cut", 1 );
1209         PostData(1, fOutList);
1210         return;
1211     }
1212
1213     if( lNchTrigger < fMinAcceptedTracksCut )
1214     {
1215         fHistEventCutStats->Fill( "LowMult cut", 1 );
1216         PostData(1, fOutList);
1217         return;
1218     }
1219     fHistAcceptedMult->Fill(lNchTrigger);
1220     fHistEventCutStats->Fill("Analyzed",1);
1221
1222     if ( lESD && fFlagWatchV0 ) // cut on V0 multiplicity "radius" in 2D-hist for both A and C sides
1223     {
1224         const AliESDVZERO* vzrData = lESD->GetVZEROData(); //aod the same
1225
1226         //        const int lThresholdMultV0RingId = 3; //which ring is considered
1227         //        Double_t lThisEventV0MultSumRing3 = vzrData->GetMRingV0A(2) + vzrData->GetMRingV0C(2);
1228         //        Double_t lThisEventV0MultSumRing4 = vzrData->GetMRingV0A(3) + vzrData->GetMRingV0C(3);
1229         //        Double_t lThisEventV0MultSum = lThisEventV0MultSumRing3 + lThisEventV0MultSumRing4;
1230         double sumV0Amult = vzrData->GetMTotV0A();
1231         double sumV0Cmult = vzrData->GetMTotV0C();
1232
1233         Double_t lThisEventV0MultSum = sumV0Amult + sumV0Cmult;
1234
1235         //sqrt ( pow ( vzrData->GetMRingV0A(lThresholdMultV0RingId), 2 ) + pow ( vzrData->GetMRingV0C(lThresholdMultV0RingId), 2 ) );
1236         if ( lThisEventV0MultSum < fThresholdOnV0mult
1237              || sumV0Amult < fThresholdOnV0mult / 5 //additionally, to "exclude" diffraction
1238              || sumV0Cmult < fThresholdOnV0mult / 5 //additionally, to "exclude" diffraction
1239              )
1240         {
1241             PostData(1, fOutList);
1242             return;
1243         }
1244     }
1245
1246     //    if ( fAnalysisLevel == "AOD" )
1247     //        return;
1248
1249
1250     //reset tmp histos
1251     fHistChargePlusVsPtTmp->Reset();
1252     fHistChargeMinusVsPtTmp->Reset();
1253
1254     //Track selection counters
1255     int lNaccept=0;
1256     int lNacceptEtaInRegion=0;
1257     int lNacceptAfterPtCuts=0;
1258     int lNacceptEtaInRegionAfterPtCuts=0;
1259     //    int lNacceptPlusEtaInRegion=0;
1260     //    int lNacceptMinusEtaInRegion=0;
1261     int lPtOver=0;
1262     int lPtUnder=0;
1263     int lNoCharge=0;
1264     int lCutfail=0;
1265     int lDecay=0;
1266     //AliLRCBase::LRCparticleType lParticleType = AliLRCBase::kLRCany;
1267
1268     int lNumberOfDropByHandTracks = 0; //number of artificially dropped tracks
1269
1270     //Double_t probTPC[AliPID::kSPECIES]={0.};
1271     Double_t probTOF[AliPID::kSPECIES]={0.};
1272     Double_t probTPCTOF[AliPID::kSPECIES]={0.};
1273
1274     Int_t nTracksInEvent = ( !lRunKine )
1275             ? event->GetNumberOfTracks()
1276             : eventMC->GetNumberOfPrimaries();//GetNumberOfTracks();
1277
1278     //        if ( fSetIncludeEventTreeInOutput )
1279     //        {
1280     //            fSimpleEvent->SetHeader( fNsimpleEvents, -1, -1, lCentralityPercentile, lReactionPlane );
1281     //            fSimpleEvent->SetVertexPos( vertex->GetX(), vertex->GetY(), vertex->GetZ() );
1282     //            fNsimpleEvents++;
1283     //            fSimpleEvent->StartEventFilling();
1284     //        }
1285     //if ( !lRunKine )
1286     int lNumberOfAcceptedTracksForLRC = 0;
1287     // Track loop -----------------------------------------------------------------
1288     for (Int_t iTracks = 0; iTracks < nTracksInEvent/*event->GetNumberOfTracks()*/; iTracks++)
1289     {
1290
1291         //Track variables
1292         double lPt = 0;   // Temp Pt
1293         double lEta = -100;       // Temp ETA
1294         double lPhi = -100;    // Temp Phi
1295         double lMass = 0;    // Temp Mass
1296         double lEt = 0;
1297         Short_t lCharge = 0;
1298
1299         AliVParticle* track = ( !lRunKine )
1300                 ? event->GetTrack(iTracks)
1301                 : eventMC->GetTrack(iTracks);
1302
1303
1304         if ( !track ) {
1305             Printf("ERROR: Could not receive track %d", iTracks);
1306             continue;
1307         }
1308
1309         fHistTrackCutStats->Fill("Total",1);
1310         fHistEtaVsZvCoverage->Fill(lVertexZ,track->Eta());
1311
1312         lPt = track->Pt();
1313         //if ( iTracks == 0 ) cout << "pt of 1st track is " << lPt << endl;
1314         lEta = track->Eta();
1315         lPhi = track->Phi();
1316
1317         if ( !lRunKine && fAnalysisLevel == "AOD" )
1318         {
1319             AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));
1320             if ( aodTrack->TestFilterBit( fAODtrackCutBit ) )
1321                 fHistEtaAODpure->Fill( lEta );
1322         }
1323
1324         //1.10.2012 - "inefficiency" in phi-acceptance "by hand"
1325         if ( lPhi > fPhiArtificialGapBegin &&
1326              lPhi < fPhiArtificialGapEnd )
1327             continue; // drop track in this "ineffective" area
1328
1329
1330         //lPhi += phiAdditional; // add here extra phi angle! (when applying tracks rotation by hand)
1331
1332
1333         //        if ( lPhi > 2 * TMath::Pi() )
1334         //            lPhi -= 2 * TMath::Pi();
1335         //cout << "track pt = " <<  lPt << endl;
1336         lCharge = track->Charge();
1337
1338         Int_t lTriggerMask = -1; //for manual event tree extraction
1339
1340         // ESD or AOD track cuts
1341         if( !lRunKine )
1342         {
1343             if ( fAnalysisLevel == "ESD" )
1344             {
1345                 AliESDtrack* lESDtrack = lESD->GetTrack(iTracks);
1346                 if( fShowPerTrackStats )
1347                     Printf("ESD Track N%d , Eta=%f, Pt=%f , TPC Nclusters =%d Sigma=%f ",  iTracks , lEta , lPt, lESDtrack-> GetTPCNcls(),fEsdTrackCuts->GetSigmaToVertex( lESDtrack) );
1348                 //cluster histograms (for watching)
1349                 fHistClustersTPC->Fill( lESDtrack->GetTPCNcls() );
1350                 fHistCrossedRowsTPC->Fill( lESDtrack->GetTPCCrossedRows() );
1351                 Int_t nITSclusters = lESDtrack->GetITSclusters(0);
1352                 fHistClustersITS->Fill( nITSclusters );//kalmanTrack->GetNumberOfClusters()  );
1353                 //Printf( " after \n" );
1354
1355                 fHist2DClustersTPCvsPt->Fill( lPt, lESDtrack->GetTPCNcls() );
1356                 fHist2DClustersTPCvsEta->Fill( lEta, lESDtrack->GetTPCNcls() );
1357
1358                 //now check track cuts: either look at fEsdTrackCuts or take cuts from array and look
1359                 if ( !fSwitchToListingCuts )
1360                 {
1361                     if( !fEsdTrackCuts->AcceptTrack(lESDtrack) )
1362                     {
1363                         lCutfail++;
1364                         if( fShowPerTrackStats )
1365                             Printf("Rejected by cuts");
1366                         fHistTrackCutStats->Fill( "QA track cuts", 1 );
1367                         continue;
1368                     }
1369                     else
1370                     {
1371                         if(fShowPerTrackStats)
1372                             Printf("OK");
1373                     }
1374                 }
1375                 else
1376                 {
1377                     lTriggerMask = 0;
1378                     //                    cout << "cuts: " << endl;
1379                     for( Int_t cutsId = 0; cutsId < fArrTrackCuts.GetEntries()/*fNumberOfCutsToRemember*/; cutsId++ )
1380                     {
1381                         //                        if ( !fArrTrackCuts[cutsId] )
1382                         //                            continue;
1383                         int cutsPassed = ( ((AliESDtrackCuts*)fArrTrackCuts[cutsId])->AcceptTrack( lESDtrack ) );
1384                         if ( cutsPassed )
1385                         {
1386                             fHistCutsNamesBins->Fill( cutsId );
1387                             Int_t cutMask = 1;
1388                             for( int iPow = 0; iPow < cutsId; iPow++ )
1389                                 cutMask *= 2;
1390                             lTriggerMask ^=  cutMask;//(cutsId+1);
1391                         }
1392                         //                        cout << cutsPassed;
1393                     }
1394                 }
1395
1396                 //fHist2DRejectedTracksPtvsEta->Fill( lPt, lEta );
1397                 fHistClustersTPCafterCuts->Fill( lESDtrack->GetTPCNcls() );
1398                 fHistCrossedRowsTPCafterCuts->Fill( lESDtrack->GetTPCCrossedRows() );
1399             }
1400             else if ( fAnalysisLevel == "AOD" )
1401             {
1402                 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));
1403                 if (!aodTrack)
1404                 {
1405                     AliError(Form("Could not receive track %d", iTracks));
1406                     continue;
1407                 }
1408
1409                 for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){
1410                     fHistAODTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));
1411                     //                    cout << aodTrack->TestFilterBit(1<<iTrackBit) << " ";
1412                 }
1413                 //                cout << endl;
1414                 if( !aodTrack->TestFilterBit( fAODtrackCutBit ) )
1415                 {
1416                     //                    cout << "fAODtrackCutBit=" << fAODtrackCutBit << endl;
1417                     fHistTrackCutStats->Fill( "QA track cuts", 1 );
1418                     continue;
1419                 }
1420                 //                int a;
1421                 //                cin >> a;
1422             }
1423
1424             fHistTracksCharge->Fill(lCharge);
1425
1426         }
1427         //        lEta = fRand->Uniform(-1,1); //!!!!!!!!!!!!!!!!!!!!!!! TEST!!!!!!
1428         // end of ESD or AOD track cuts
1429
1430
1431         //MC truth
1432         //Int_t label = -1000;
1433         TParticle * part = 0x0;
1434         if( lRunKine )   // Dropping undetectable particles in Kine
1435         {
1436             //            if ( fAnalyseMCESD ) //run MC+ESD
1437             //            {
1438             //                label = track->GetLabel();//TMath::Abs(track->GetLabel()); //abs!!!!
1439             //                fHistMClabels->Fill( label );
1440             //                if ( label < 0 ) //by Peter Hristov - it's ghost tracks
1441             //                    continue;
1442             //                part = stack->Particle(label);
1443             //            }
1444             //            else
1445             {
1446                 AliMCParticle* trackMC = dynamic_cast<AliMCParticle *>(eventMC->GetTrack(iTracks));
1447                 if(!trackMC)
1448                     continue;
1449                 part = trackMC->Particle();
1450                 if( !part ) continue;
1451             }
1452
1453             Int_t lNoD = part->GetNDaughters();
1454             if(fShowPerTrackStats)
1455             {
1456                 printf("Track %d, Mother %d, track Pt=%f, MC Pt=%f, PDG=%d  Nof Dothers=%d  ETA=%f   ## ",  iTracks , part->GetFirstMother() ,lPt,      part->Pt() ,part->GetPdgCode(),lNoD,lEta);
1457                 printf("%s", part->GetName());
1458                 printf(" ");
1459                 printf("%s", part->GetPDG()->ParticleClass());
1460             }
1461
1462             if( !stack->IsPhysicalPrimary(iTracks) )//label))
1463             {
1464                 fHistTrackCutStats->Fill( "Not phys primary", 1 );
1465                 continue;
1466             }
1467
1468             //charge in TParticlePDG is in units of |e|/3
1469             fHistTracksCharge->Fill( part->GetPDG()->Charge() / 3 ); // 0-charge bin fills only for MC truth events
1470             //cout << " charge = " << part->GetPDG()->Charge();
1471             if ( part->GetPDG()->Charge() == 0 )
1472             {
1473                 if(fShowPerTrackStats)
1474                     Printf(" ChargeReject");
1475                 lNoCharge++;
1476                 continue;
1477             }
1478
1479             //cut on low-pt kine particles
1480             if ( lPt < fKineLowPtCut )
1481                 continue;
1482         }    //End of  Kine particle filter
1483
1484         //now decide that we have a good track
1485         lNaccept++;
1486
1487         if( fabs(lEta) < fEtaRegionForTests ) //look at eta region
1488         {
1489             lNacceptEtaInRegion++;
1490             //net charge fillings
1491             fHistNetChargeVsPt->Fill( lPt, lCharge );
1492             if ( lCharge > 0 )
1493             {
1494                 fHistChargePlusVsPtTmp->Fill( lPt );
1495                 fHistPtPlus->Fill( lPt );
1496             }
1497             else if ( lCharge < 0 )
1498             {
1499                 fHistChargeMinusVsPtTmp->Fill( lPt );
1500                 fHistPtMinus->Fill( lPt );
1501             }
1502         }
1503
1504         if( lPt > fMaxPtLimit )
1505         {
1506             lPtOver++;
1507             fHistTrackCutStats->Fill(2);
1508             continue;
1509         } // Dropping tracks with hi Pt
1510
1511         if( lPt < fMinPtLimit )
1512         {
1513             lPtUnder++;
1514             fHistTrackCutStats->Fill(3);
1515             continue;
1516         } // Dropping tracks with low Pt
1517         lNacceptAfterPtCuts++;
1518         fHistTrackCutStats->Fill(4);
1519         fHistPt->Fill( lPt );
1520         fHistEta->Fill( lEta );
1521         fHistPhi->Fill( lPhi );
1522
1523         if( fabs(lEta) < fEtaRegionForTests ) //look at eta region
1524             lNacceptEtaInRegionAfterPtCuts++;
1525
1526         //fill counters for ZDC-comparing with a 'shelve'
1527         if ( lEta > etaBmin && lEta < etaBmax )
1528             countTracksEtaB++;
1529         if ( lEta > etaFmin && lEta < etaFmax )
1530             countTracksEtaF++;
1531
1532         if ( fabs( lEta ) < 0.8 ) //eta-phi for tracks in "ALICE barrel acceptance"
1533             fHistEtaPhi->Fill(lEta, lPhi);
1534
1535
1536         // ###### New PID
1537         Int_t lMostProbablePIDPure = -1;
1538         Int_t lMostProbablePIDdirty = -1;
1539         if ( (lESD || lAOD) && !lRunKine && fPIDResponse )
1540         {
1541             AliVTrack *trackV = (AliVTrack*)event->GetTrack(iTracks);
1542             //##### from Pid Combined Task
1543             // compute priors for TPC+TOF, even if we ask just TOF for PID
1544             fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
1545             UInt_t detUsed = fPIDCombined->ComputeProbabilities(trackV, fPIDResponse, probTOF);
1546             Double_t priors[5];         // check priors used for TOF
1547             fPIDCombined->GetPriors(trackV,priors,fPIDResponse,detUsed);
1548             for(Int_t ispec=0;ispec<5;ispec++) fPriorsUsed[ispec]->Fill(TMath::Abs(trackV->Pt()),priors[ispec]);
1549
1550             fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);
1551             detUsed = fPIDCombined->ComputeProbabilities(trackV, fPIDResponse, probTPCTOF);
1552             fHistProbabilityPion->Fill( probTPCTOF[AliPID::kPion] );////probDensity[iParticle] );
1553             fHistProbabilityKaon->Fill( probTPCTOF[AliPID::kKaon] );
1554             fHistProbabilityProton->Fill( probTPCTOF[AliPID::kProton] );////probDensity[iParticle] );
1555             //cout << "chto-to s detId: " << detUsed << endl;
1556 //            for ( Int_t ispec=0; ispec < (int)(AliPID::kSPECIES); ++ispec ) {
1557 //                //cout << probTPCTOF[ispec] << "  " << endl;
1558 //            }
1559             //cout << endl;
1560             if ( detUsed == (UInt_t)fPIDCombined->GetDetectorMask() )
1561             {
1562                 Double_t lMaxProb = 0;
1563                 for ( Int_t ispec=0; ispec < (int)(AliPID::kSPECIES); ++ispec )
1564                 {
1565                     //Double_t nSigmaTPC = fPIDResponse->NumberOfSigmasTPC(trackV,(AliPID::EParticleType)ispec);
1566                     fProbTPCTOF[ispec]->Fill( lPt, probTPCTOF[ispec] );
1567                     //fProbTPCTOFnSigmaTPC[ispec]->Fill(nSigmaTPC,probTPCTOF[ispec]);
1568                     //fProbTPCTOFnSigTPCMom[ibin][ispec]->Fill(nSigmaTPC,probTPCTOF[ispec]);
1569
1570                     if ( probTPCTOF[ispec] > lMaxProb )
1571                     {
1572                         lMaxProb = probTPCTOF[ispec];
1573                         lMostProbablePIDPure = ispec;
1574                     }
1575                 }
1576                 //fill for watching at max prob-s in pid species arrays
1577                 fHistPidPureMaxProbability->Fill( lMaxProb );
1578             }
1579
1580             const bool useDirtyPID = true;
1581             if ( useDirtyPID )
1582             {
1583                 //Int_t lMostProbablePIDdirty = -1;
1584
1585                 Double_t lMaxProb = 0;
1586                 for ( Int_t ispec=0; ispec < (int)(AliPID::kSPECIES); ++ispec )
1587                 {
1588                     fProbAllDets[ispec]->Fill( lPt, probTPCTOF[ispec] );
1589                     if ( probTPCTOF[ispec] > lMaxProb )
1590                     {
1591                         lMaxProb = probTPCTOF[ispec];
1592                         lMostProbablePIDdirty = ispec;  // define most probable particle!!!
1593                     }
1594                 }
1595                 //fill for watching at max prob-s in pid species arrays
1596                 fHistPidMaxProbability->Fill( lMaxProb );
1597                 fHistParticlesDistr->Fill( lMostProbablePIDdirty );
1598             }
1599
1600
1601         }
1602         // end of new PID filling
1603
1604         fHist2DAcceptedTracksPtvsEta->Fill( lPt, lEta );
1605
1606         //get particle mass and Et
1607         if ( lMostProbablePIDPure >= 0 )
1608         {
1609             lMass = AliPID::ParticleMass( lMostProbablePIDPure );
1610             fHistESDtrackMass->Fill( lMass );
1611             lEt = sqrt( lPt*lPt + lMass*lMass );
1612             //cout << lEt << " " << lMass << " " << lPt << endl;
1613         }
1614
1615         //eta vs vertex z coverage
1616         fHistEtaVsZvCoverageAccepted->Fill(lVertexZ,track->Eta());
1617
1618         //artificial inefficiency (27.09.12) - check, what if we drop some fraction of tracks randomly
1619         Bool_t lDropDecision = ( ( 1.0 - fArtificialInefficiency ) < fRand->Uniform(0,1) );
1620
1621         if ( lDropDecision )
1622         {
1623             lNumberOfDropByHandTracks++;
1624             continue;
1625         }
1626
1627         if ( lNumberOfAcceptedTracksForLRC >= kMaxParticlesNumber )
1628         {
1629             AliDebug(AliLog::kError, "lNumberOfAcceptedTracksForLRC is too large...");
1630             break;
1631         }
1632
1633         //fill arrays with track data for LRC
1634         fArrayTracksPt[lNumberOfAcceptedTracksForLRC]       = (fEtInsteadOfPt ? lEt : lPt);
1635         fArrayTracksEta[lNumberOfAcceptedTracksForLRC]      = lEta;
1636         fArrayTracksPhi[lNumberOfAcceptedTracksForLRC]      = lPhi;
1637         fArrayTracksCharge[lNumberOfAcceptedTracksForLRC]   = lCharge;
1638         fArrayTracksPID[lNumberOfAcceptedTracksForLRC]      = lMostProbablePIDPure;
1639
1640         lNumberOfAcceptedTracksForLRC++;
1641
1642     } //end of track loop
1643
1644     //net charge vs pt
1645 //    if ( lNacceptEtaInRegion > 20 )
1646 //    {
1647 //        fHistChargePlusVsPtTmp->DrawClone();
1648 //        fHistChargeMinusVsPtTmp->SetLineColor( kBlue );
1649 //        fHistChargeMinusVsPtTmp->DrawClone("same");
1650 //    }
1651
1652
1653     //final actions for net-charge
1654     double integralPlus     = fHistChargePlusVsPtTmp->ComputeIntegral();
1655     double integralMinus    = fHistChargeMinusVsPtTmp->ComputeIntegral();
1656     double xAxisLowEdge     = fHistChargePlusVsPtTmp->GetBinLowEdge(1);
1657     double xAxisHighEdge    = fHistChargePlusVsPtTmp->GetBinLowEdge(fHistChargePlusVsPtTmp->GetNbinsX()) + fHistChargePlusVsPtTmp->GetBinWidth(fHistChargePlusVsPtTmp->GetNbinsX());
1658     double ptRange =  xAxisHighEdge - xAxisLowEdge;
1659 //    cout << "integralPlus=" << integralPlus << ", integralMinus=" << integralMinus << ", ptRange=" << ptRange << endl;
1660
1661     for ( int bin = 0; bin < fHistChargePlusVsPtTmp->GetNbinsX(); bin++)
1662     {
1663         double binCenter = fHistChargePlusVsPtTmp->GetBinCenter(bin+1);
1664         double binContentPlus = fHistChargePlusVsPtTmp->GetBinContent(bin+1);
1665         double binContentMinus = fHistChargeMinusVsPtTmp->GetBinContent(bin+1);
1666
1667
1668
1669         double lNetCharge = binContentPlus-binContentMinus;
1670         fHist2DNetChargeVsPt->Fill( binCenter, lNetCharge );
1671         if ( ptRange>0 )
1672             fHist2DNetChargeVsPtCorrectedOnEventMean->Fill( binCenter, lNetCharge - (integralPlus-integralMinus)/ptRange );
1673         //        fHist2DNetChargeVsPtCorrectedOnEventMeanNormOnNch->Fill( );
1674
1675     }
1676
1677     //############ MC ESD QA
1678     if ( fRunKine && fAnalyseMCESD ) //run MC+ESD
1679     {
1680         Double_t lEta;
1681         Double_t lPt;
1682 //        Double_t lY;
1683         Int_t lMCtracksAcceptedEtaInRegionAfterPtCuts = 0;
1684         for (Int_t iMCtracks = 0; iMCtracks < eventMC->GetNumberOfPrimaries(); iMCtracks++)
1685         {
1686             AliMCParticle* track = dynamic_cast<AliMCParticle *>(eventMC->GetTrack(iMCtracks));
1687             if (!track) {
1688                 Printf("ERROR: Could not receive particle %d", iMCtracks);
1689                 continue;
1690             }
1691
1692             //exclude non stable particles
1693             if(!eventMC->IsPhysicalPrimary(iMCtracks)) continue;
1694
1695             if ( track->Particle()->GetPDG()->Charge() == 0 )
1696                 continue;
1697
1698             lEta    = track->Eta();
1699             lPt     = track->Pt();
1700 //            lY      = track->Y();
1701
1702             if( fabs(lEta) > fEtaRegionForTests ) //look at eta region
1703                 continue;
1704
1705             if( lPt > fMaxPtLimit || lPt < fMinPtLimit )
1706                 continue;
1707
1708
1709             lMCtracksAcceptedEtaInRegionAfterPtCuts++;
1710         }
1711         fHist2DMultiplicityMCESDInEtaRegion->Fill( lMCtracksAcceptedEtaInRegionAfterPtCuts, lNacceptEtaInRegionAfterPtCuts );
1712     }
1713
1714
1715
1716
1717     //profiling tracks in phi BY HAND if requested
1718     //    if ( fUsePhiShufflingByHand )
1719     //    {
1720     //        double lFictionEventPlanePhi = fRand->Uniform(0,2*TMath::Pi());
1721     //        fHistPhiArtificialEvPlane->Fill( lFictionEventPlanePhi );
1722     //        TF1 lFictionEventPlaneFunc( "fictionEventPlane", "0.75 + 0.25*TMath::Cos(x)", 0, 2*TMath::Pi() );
1723     //        for ( Int_t trackId = 0; trackId < lNumberOfAcceptedTracksForLRC; trackId++ )
1724     //        {
1725     //            //check phi wrt event plane
1726     //            Double_t lProfiledPhiWrtEvPlane = lFictionEventPlaneFunc.GetRandom();
1727     //            //FixAngleInTwoPi( lProfiledPhiWrtEvPlane );
1728     //            fHistPhiArtificialProfilingCheckWrtEvPlane->Fill( lProfiledPhiWrtEvPlane );
1729
1730     //            //  double lOppositePhi = ( ( fRand->Uniform(0,1) > 0.5 ) ? -TMath::Pi() : 0 );
1731     //            //  double lRandomConeAngleEta = fRand->Gaus(0,0.5);
1732     //            //  double lRandomConeAnglePhi = fRand->Gaus(0,TMath::PiOver2());
1733
1734     //            //double lFictionPhi = fRand->Uniform(0,2*TMath::Pi());
1735     //            //double lPhiSign = ( ( fRand->Uniform(0,1) > 0.5 ) ? TMath::Pi() : 0 );
1736
1737     //            //add event plane phi angle
1738     //            Double_t lProfiledPhi = lProfiledPhiWrtEvPlane + lFictionEventPlanePhi;
1739     //            FixAngleInTwoPi( lProfiledPhi );
1740     //            fHistPhiArtificialProfilingCheck->Fill( lProfiledPhi );
1741
1742     //            fArrayTracksPhi[trackId] = lProfiledPhi; //lPhi + lRandomConeAnglePhi; //lPhi;
1743
1744     //            //        fArrayTracksPhi[lNumberOfAcceptedTracksForLRC]      = lPhi - lRandomConeAnglePhi; //lPhi;
1745     //            //        if ( fArrayTracksPhi[lNumberOfAcceptedTracksForLRC] > 2 * TMath::Pi() )
1746     //            //            fArrayTracksPhi[lNumberOfAcceptedTracksForLRC] -= 2 * TMath::Pi();
1747     //        }
1748     //    }
1749
1750
1751     // ######### fill some QA plots
1752     fHistAcceptedTracks->Fill(lNaccept);
1753     fHistMultiplicityInEtaRegion->Fill( lNacceptEtaInRegion );
1754     fHistMultiplicityInEtaRegionAfterPtCuts->Fill( lNacceptEtaInRegionAfterPtCuts );
1755     fHistAcceptedTracksAfterPtCuts->Fill(lNacceptAfterPtCuts);
1756     if ( fArtificialInefficiency >= 0 ) //i.e. have ineff analysis ON
1757         fHistNumberOfDroppedByHandTracks->Fill( lNacceptAfterPtCuts, lNumberOfDropByHandTracks );
1758
1759     // ########### ProfilePhiByHand
1760     if ( fUsePhiShufflingByHand )
1761         ProfilePhiByHand( lNumberOfAcceptedTracksForLRC);
1762
1763     // ########### LRC processors filling
1764     Bool_t lFlagAcceptEventForLRC = ( lNacceptEtaInRegion >= fMultCutInEtaRegion );
1765     if ( lFlagAcceptEventForLRC )
1766         FillLRCProcessors( lNumberOfAcceptedTracksForLRC, lCentralityPercentile );
1767
1768     // ##### ZDC, VZERO, etc
1769     if ( fFlagWatchZDC && fIsIonsAnalysis )
1770     {
1771         AliESDZDC *zdcData = lESD->GetESDZDC();
1772         //fHistZDCenergy->Fill( zdcData->GetZDCEMEnergy(0) );
1773         fHistZDCparticipants->Fill( zdcData->GetZDCParticipants ());
1774         for ( int i = 0; i < 5; i++ )
1775         {
1776             //if multiplicity in eta windows above threshold - fill ZDC energy hist
1777             if ( countTracksEtaB >= fMultForZDCstudy[i] && countTracksEtaF >= fMultForZDCstudy[i] )
1778                 fHistZDCenergy[i]->Fill( zdcData->GetZDCEMEnergy(0) );
1779         }
1780     }
1781     if ( fFlagWatchV0 )
1782     {
1783         const AliESDVZERO* vzrData = lESD->GetVZEROData(); //aod the same
1784
1785
1786         //V0 multiplicity
1787         float lV0mult[64];
1788         //float lV0Amult[64];
1789         //float lV0Cmult[64];
1790         double sumV0mult = 0;
1791         double sumV0Amult = 0;
1792         double sumV0Cmult = 0;
1793
1794         for (int i = 0; i < 64; i++ )
1795         {
1796             lV0mult[i] = (float)(vzrData->GetMultiplicity(i));
1797             sumV0mult += lV0mult[i];
1798         }
1799         fHistV0multiplicity->Fill( sumV0mult );
1800
1801         //            for (int i=0; i<32; i++)
1802         //            {
1803         //                lV0Cmult[i] = (float)(vzrData->GetMultiplicityV0C(i));
1804         //                sumV0Cmult += lV0Cmult[i];
1805         //                lV0Amult[i] = (float)(vzrData->GetMultiplicityV0A(i));
1806         //                sumV0Amult += lV0Amult[i];
1807         //            }
1808         sumV0Amult = vzrData->GetMTotV0A();
1809         sumV0Cmult = vzrData->GetMTotV0C();
1810         fHistV0Amultiplicity->Fill( sumV0Amult );
1811         fHistV0Cmultiplicity->Fill( sumV0Cmult );
1812         fHist2DV0ACmultiplicity->Fill( sumV0Amult, sumV0Cmult );
1813
1814         fHist2DTracksAcceptedVsV0multiplicity->Fill( lNaccept, sumV0Amult );
1815
1816         //cells
1817         //fHistV0cells->Fill(  );
1818         fHistV0Acells->Fill( vzrData->GetNbPMV0A() );
1819         fHistV0Ccells->Fill( vzrData->GetNbPMV0C() );
1820         fHist2DV0ACcells->Fill( vzrData->GetNbPMV0A(), vzrData->GetNbPMV0C() );
1821
1822         //rings
1823         for ( int i = 0; i < 4; i++ )
1824         {
1825             int lMultRingV0A = vzrData->GetMRingV0A(i);
1826             int lMultRingV0C = vzrData->GetMRingV0C(i);
1827             fHistV0AmultiplicityRing[i]->Fill( lMultRingV0A );
1828             fHistV0CmultiplicityRing[i]->Fill( lMultRingV0C );
1829             fHist2DV0ACmultiplicityRing[i]->Fill( lMultRingV0A, lMultRingV0C );
1830
1831             fHist2DTracksAcceptedVsV0AmultiplicityRing[i]->Fill( lNaccept, lMultRingV0A );
1832             fHist2DTracksAcceptedVsV0CmultiplicityRing[i]->Fill( lNaccept, lMultRingV0C );
1833         }
1834
1835     }
1836     if ( fFlagWatchFMD )
1837     {
1838         // multiplicity
1839         //const AliESDFMD *lFMDdata = lESD->GetFMDData();
1840
1841
1842     }
1843
1844     //        //look at V0s
1845     //        AliESDv0 *lV0 = 0x0;
1846     //        //cout << ">>>> showing v0's:" << endl;
1847     //        for ( int iV0 = 0; iV0 < lESD->GetNumberOfV0s(); iV0++ )
1848     //        {
1849     //            lV0 = lESD->GetV0( iV0 );
1850
1851     //            int lV0pdgCode = lV0->GetPdgCode();
1852     //            double lMassV0 = lV0->GetEffMass();//ChangeMassHypothesis( lV0pdgCode );  //GetEffMass();
1853     //            //if ( abs(lV0pdgCode) > 10 )//== 3122 )
1854     //            {
1855     //                //cout << ">>>>>>>>>>>>>>>>>>>>>>>>> " << lV0pdgCode << ": mass = " << lMassV0 << endl;
1856     //                fHistV0spectra->Fill(lV0pdgCode);//lMassV0);
1857     //            }
1858
1859     //            if ( lMassV0 > 1. )// abs(lPdgCode) == 3122 || abs(lPdgCode) == 421 ) //abs(lPdgCode) % 1000 == 122 )
1860     //            {
1861     //                //cout << ">>>>>>>>>>>>>>>>>>>>>>>>> " << lV0pdgCode << ": mass = " << lMassV0 << endl;
1862     //            }
1863     //        }
1864
1865
1866
1867
1868
1869
1870     //Debuging output of track filter
1871     if( fShowPerTrackStats )
1872         Printf("NofTracks= %d , accepted %d , LowPt %d, HighPt %d, LowCharge %d,  ESD cut fail %d , Decay Filer %d", event->GetNumberOfTracks(), lNaccept, lPtUnder, lPtOver ,lNoCharge , lCutfail  , lDecay);
1873
1874     // Post output data
1875     PostData(1, fOutList);
1876
1877     for( Int_t i = 0; i < fLRCproc.GetEntries(); i++)
1878     {
1879         PostData( Proc(i)->GetOutputSlotNumber(),Proc(i)->CreateOutput() );
1880     }
1881
1882 }
1883
1884
1885 void AliAnalysisTaskLRC::SetParticleTypeToProcessors( int windowId, TString strPid ) //char* strPid )//AliLRCBase::LRCparticleType whichParticleToFill )
1886 {
1887     //dummy actions
1888     windowId+=1;
1889     if(0) printf("%s", strPid.Data());
1890
1891     //Set pid types for LRC processors (windowId = 0 - fwd window, =1 - backward window)
1892     //    Int_t lLrcNum=fLRCproc.GetEntries();
1893
1894     //    for(Int_t i=0; i < lLrcNum; i++)
1895     //    {
1896     //        AliLRCBase *lrcProc = dynamic_cast<AliLRCBase*> (fLRCproc.At(i));
1897     //        if(lrcProc)
1898     //        {
1899     //            //tmp lines!
1900     //            windowId++;
1901     //            if (0)
1902     //                printf("%s",strPid);
1903     //            //commented for this commit
1904     //            //            if ( windowId == 0 )
1905     //            //                lrcProc->SetParticleType( "fwd", strPid );
1906     //            //            else if ( windowId == 1 )
1907     //            //                lrcProc->SetParticleType( "bkwd", strPid );
1908
1909     //            //fTmpCounter++;
1910     //        }
1911     //        else continue;
1912     //        //fTmpCounter++;
1913     //    }
1914
1915 }
1916
1917 void AliAnalysisTaskLRC::SetParticleTypeForTask( TString strF, TString strB ) //char* strF, char* strB )
1918 {
1919     // Set PID sensitivity to 'true' and write pids for windows
1920     fPIDsensingFlag = kTRUE;
1921     fStrPIDforFwd = strF;
1922     fStrPIDforBwd = strB;
1923
1924     //    sprintf ( fStrPIDforFwd, "%s", strF );
1925     //    sprintf ( fStrPIDforBwd, "%s", strB );
1926     //cout << "we just set fwd win to " << fStrPIDforFwd << " and bwd win to " << fStrPIDforBwd << endl;
1927 }
1928
1929
1930
1931 //________________________________________________________________________
1932 void AliAnalysisTaskLRC::Terminate(Option_t *)
1933 {
1934     // Draw result to the screen
1935     // Called once at the end of the query
1936     //fOutList = dynamic_cast<TList*> (GetOutputData(0));
1937
1938     //lESDTrackCuts->DrawHistograms();
1939
1940
1941     fAnalysisTimer->Stop();
1942     Double_t rtime = fAnalysisTimer->RealTime();
1943     Double_t ctime = fAnalysisTimer->CpuTime();
1944
1945     printf("RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime);
1946
1947 }
1948
1949
1950 void AliAnalysisTaskLRC::AddLRCProcess(AliLRCBase *newProc)
1951 {
1952     // Add new AliLRCBase (Main LRC processor per ETA window) to the processing list
1953     // Used to add new ETA window to AnalysisTask
1954     if(!newProc)
1955     {
1956         Printf("ERROR:No AliLRCBase object -  NULL pointer!");
1957         return;
1958     }
1959
1960     fLRCproc.Add(newProc);
1961     newProc->SetOutputSlotNumber(fLRCproc.GetEntries() + 1 );//( fSetIncludeEventTreeInOutput ? 2 : 1 ) );
1962     DefineOutput(newProc->GetOutputSlotNumber(),TList::Class());
1963     return ;
1964 }
1965
1966
1967 Double_t AliAnalysisTaskLRC::GetEventPlane(AliVEvent *event)
1968 {
1969     // Get the event plane
1970
1971     if ( !event )
1972         return -1;
1973     //TString gAnalysisLevel = fBalance->GetAnalysisLevel();
1974
1975     Float_t lVZEROEventPlane    = -10.;
1976     Float_t lReactionPlane      = -10.;
1977     Double_t qxTot = 0.0, qyTot = 0.0;
1978
1979     //MC: from reaction plane
1980     if( fRunKine )//gAnalysisLevel == "MC"){
1981     {
1982         AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(dynamic_cast<AliMCEvent*>(event)->GenEventHeader());
1983         if (headerH)
1984         {
1985             lReactionPlane = headerH->ReactionPlaneAngle();
1986             //lReactionPlane *= TMath::RadToDeg();
1987         }
1988     }//MC
1989
1990     // AOD,ESD,ESDMC: from VZERO Event Plane
1991     else
1992     {
1993         AliEventplane *ep = event->GetEventplane();
1994         if( ep )
1995         {
1996             lVZEROEventPlane = ep->CalculateVZEROEventPlane( event, 10, 2, qxTot, qyTot);
1997             if( lVZEROEventPlane < 0. )
1998                 lVZEROEventPlane += TMath::Pi();
1999             //lReactionPlane = lVZEROEventPlane*TMath::RadToDeg();
2000             lReactionPlane = lVZEROEventPlane;
2001             //            cout << "lReactionPlane: " << lReactionPlane << endl;
2002             //            int aa;
2003             //            cin >> aa;
2004         }
2005     }//AOD,ESD,ESDMC
2006
2007     return lReactionPlane;
2008 }
2009
2010
2011 void AliAnalysisTaskLRC::AddTrackCutForBits(AliESDtrackCuts*  const cuts, TString cutsName )
2012 {
2013     //TString *strPtrCutsName = new TString(cutsName.Data());
2014 //    cout << "ADDING CUTS " << cutsName << endl;
2015     //    fArrTrackCuts[fNumberOfCutsToRemember] = cuts;
2016     int lNumberOfCutsToRemember = fArrTrackCuts.GetEntries();
2017
2018     fArrCutsNames[lNumberOfCutsToRemember] = cutsName;
2019     fArrTrackCuts.Add( cuts );// = cuts;
2020     //    fArrCutsNames.Add( strPtrCutsName );
2021     //    fNumberOfCutsToRemember++;
2022 }
2023
2024 void AliAnalysisTaskLRC::UseToyEvents()
2025 {
2026     //generate all events here and fill LRC processors
2027     for ( Int_t toyEventId = 0; toyEventId < fNumberOfToyEvents; toyEventId++ )
2028     {
2029         if ( toyEventId % 10000 == 0 )
2030             printf("Processing %d event...\n", toyEventId );
2031
2032         //        cout << "toyEventId=" << toyEventId << endl;
2033         int nToyTracks =  0;
2034         const double kExpParam = 0.42;
2035         const double kEtaGausSigma = 0.3;
2036         const double kPhiGausSigma = TMath::PiOver4();
2037
2038         //try to do "jets" with some particles
2039         const int nJets = 2;
2040
2041         //jets eta angles
2042         double lJets_Eta[nJets];
2043         lJets_Eta[0] = fRand->Uniform(0, 0.8);
2044         lJets_Eta[1] = -lJets_Eta[0];
2045
2046         //jets phi angles
2047         double lJets_Phi[nJets];
2048         lJets_Phi[0] = fRand->Uniform(0, 2*TMath::Pi());
2049         lJets_Phi[1] = lJets_Phi[0] + TMath::Pi();
2050         FixAngleInTwoPi( lJets_Phi[1] );
2051
2052
2053         //jets loop
2054         for ( Int_t jetId = 0; jetId < nJets; jetId++ )
2055         {
2056             int nPartInJet = fRand->Gaus(4,1);
2057             //jet tracks params
2058             for ( Int_t toyTrackId = 0; toyTrackId < nPartInJet; toyTrackId++ )
2059             {
2060                 fArrayTracksPt[nToyTracks] = fRand->Exp( kExpParam );
2061                 fArrayTracksEta[nToyTracks] = fRand->Gaus( lJets_Eta[jetId], kEtaGausSigma );
2062                 fArrayTracksPhi[nToyTracks] = fRand->Gaus( lJets_Phi[jetId], kPhiGausSigma );
2063                 FixAngleInTwoPi( fArrayTracksPhi[nToyTracks] );
2064
2065                 nToyTracks++;
2066             }
2067         }
2068
2069         FillLRCProcessors( nToyTracks, 0 );
2070
2071         //        for ( Int_t toyTrackId = 0; toyTrackId < nToyTracks; toyTrackId++ )
2072         //        {
2073         //            // ########### ProfilePhiByHand
2074         //            ProfilePhiByHand( nToyTracks);
2075         //        }
2076
2077     }
2078 }
2079
2080
2081 void AliAnalysisTaskLRC::FillLRCProcessors( int numberOfAcceptedTracksForLRC, Double_t eventCentrality )
2082 {
2083     //pass signal to LRC-based analysis to start event
2084     Int_t lLrcNum = fLRCproc.GetEntries(); // Number of processors attached
2085
2086     //prepare phi rotation step
2087     Double_t phiStep = 2 * TMath::Pi();
2088     if ( fNumberOfPhiSectors > 1 )
2089         phiStep /= fNumberOfPhiSectors;
2090     Double_t lPhiRotatedExtra = 0; //additional phi rotation of all tracks in case
2091
2092
2093     //start LRC filling
2094     for ( Int_t sectorId = 0; sectorId < fNumberOfPhiSectors; sectorId++ )
2095     {
2096         if ( lLrcNum > kMaxParticlesNumber )
2097         {
2098             AliDebug(AliLog::kError, "lLrcNum is too large... Mistake???");
2099             break;
2100         }
2101
2102         //pass signal to LRC-based analysis to start event
2103         for( Int_t lrcProcessorId = 0; lrcProcessorId < lLrcNum; lrcProcessorId++ )
2104         {
2105             AliLRCBase *lrcBase = dynamic_cast<AliLRCBase*> (fLRCproc.At(lrcProcessorId));
2106             if (!lrcBase)
2107             {
2108                 AliDebug(AliLog::kError, "can't cast lrcBase to AliLRCBase!" );
2109                 break;
2110             }
2111             lrcBase->StartEvent();
2112             //pass the centrality
2113             lrcBase->SetEventCentrality( eventCentrality );
2114             //        }
2115
2116             //pass track data to LRC-based analysis
2117             for ( Int_t trackId = 0; trackId < numberOfAcceptedTracksForLRC; trackId++ )
2118             {
2119                 if ( fEtInsteadOfPt && (fArrayTracksPt[trackId]  == 0 ) ) //in Et-mode we didn't get Et! exit
2120                 {
2121                     AliDebug(AliLog::kError, "pT=0 Mistake???" );
2122                     break;
2123                 }               //if ( fEtInsteadOfPt && lMostProbablePIDPure == -1 ) //don't have pure PID
2124                 //      continue;
2125
2126                 //rotate track phi
2127                 Double_t lPhi = fArrayTracksPhi[trackId] + lPhiRotatedExtra;
2128                 FixAngleInTwoPi( lPhi );
2129
2130                 fHistPhiLRCrotationsCheck->Fill( lPhi );
2131
2132                 //            for(Int_t i = 0; i < lLrcNum; i++ )
2133                 //            {
2134                 lrcBase->AddTrackPtEta(
2135                             fArrayTracksPt[trackId]
2136                             , fArrayTracksEta[trackId]
2137                             , lPhi
2138                             , fArrayTracksCharge[trackId]
2139                             , fArrayTracksPID[trackId]
2140                             );
2141                 //            }
2142
2143             }
2144
2145             //take event only if at least 1 track in this event fulfill the requirements! //21.11.11
2146             Bool_t lDontTakeEventDecision =  kFALSE; //lNaccept > 0 ? kFALSE : kTRUE;
2147             //pass signal to LRC-based analysis to finish the event
2148             //        for( Int_t i = 0; i < lLrcNum; i++ )
2149             //        {
2150             lrcBase->FinishEvent( lDontTakeEventDecision );
2151         }
2152
2153         lPhiRotatedExtra += phiStep; //increase phi step to rotate tracks
2154     }
2155 }
2156
2157 void AliAnalysisTaskLRC::ProfilePhiByHand( int numberOfAcceptedTracksForLRC )
2158 {
2159     //profiling tracks in phi BY HAND if requested
2160
2161     double lFictionEventPlanePhi = fRand->Uniform(0,2*TMath::Pi());
2162     fHistPhiArtificialEvPlane->Fill( lFictionEventPlanePhi );
2163     TF1 lFictionEventPlaneFunc( "fictionEventPlane", "0.75 + 0.25*TMath::Cos(x)", 0, 2*TMath::Pi() );
2164     for ( Int_t trackId = 0; trackId < numberOfAcceptedTracksForLRC; trackId++ )
2165     {
2166         //check phi wrt event plane
2167         Double_t lProfiledPhiWrtEvPlane = lFictionEventPlaneFunc.GetRandom();
2168         //FixAngleInTwoPi( lProfiledPhiWrtEvPlane );
2169         fHistPhiArtificialProfilingCheckWrtEvPlane->Fill( lProfiledPhiWrtEvPlane );
2170
2171         //  double lOppositePhi = ( ( fRand->Uniform(0,1) > 0.5 ) ? -TMath::Pi() : 0 );
2172         //  double lRandomConeAngleEta = fRand->Gaus(0,0.5);
2173         //  double lRandomConeAnglePhi = fRand->Gaus(0,TMath::PiOver2());
2174
2175         //double lFictionPhi = fRand->Uniform(0,2*TMath::Pi());
2176         //double lPhiSign = ( ( fRand->Uniform(0,1) > 0.5 ) ? TMath::Pi() : 0 );
2177
2178         //add event plane phi angle
2179         Double_t lProfiledPhi = lProfiledPhiWrtEvPlane + lFictionEventPlanePhi;
2180         FixAngleInTwoPi( lProfiledPhi );
2181         fHistPhiArtificialProfilingCheck->Fill( lProfiledPhi );
2182
2183         fArrayTracksPhi[trackId] = lProfiledPhi; //lPhi + lRandomConeAnglePhi; //lPhi;
2184
2185         //        fArrayTracksPhi[lNumberOfAcceptedTracksForLRC]      = lPhi - lRandomConeAnglePhi; //lPhi;
2186         //        if ( fArrayTracksPhi[lNumberOfAcceptedTracksForLRC] > 2 * TMath::Pi() )
2187         //            fArrayTracksPhi[lNumberOfAcceptedTracksForLRC] -= 2 * TMath::Pi();
2188     }
2189 }