PWGLF/SPECTRA/ChargedHadrons/multPbPb/AliAnalysisTaskMultPbTracks.cxx

   1 // AliAnalysisTaskMultPbTracks
   2
   3 // Author: Michele Floris, CERN
   4 // TODO:
   5 // - Add chi2/cluster plot for primary, secondaries and fakes
   6
   7 #include "AliAnalysisTaskMultPbTracks.h"
   8 #include "AliESDInputHandler.h"
   9 #include "AliAnalysisMultPbTrackHistoManager.h"
  10 #include "AliAnalysisManager.h"
  11 #include "AliESDtrackCuts.h"
  12 #include "AliMCEvent.h"
  13 #include "AliStack.h"
  14 #include "TH1I.h"
  15 #include "TH3D.h"
  16 #include "TH2D.h"
  17 #include "AliMCParticle.h"
  18 #include "AliGenEventHeader.h"
  19 #include "AliCentrality.h"
  20 #include "AliMultiplicity.h"
  21 #include "TFile.h"
  22 #include <iostream>
  23 #include "AliAnalysisMultPbCentralitySelector.h"
  24 #include "AliTriggerAnalysis.h"
  25 #include "AliPIDResponse.h"
  26
  27 using namespace std;
  28
  29 ClassImp(AliAnalysisTaskMultPbTracks)
  30
  31 AliAnalysisTaskMultPbTracks::AliAnalysisTaskMultPbTracks()
  32 : AliAnalysisTaskSE("TaskMultPbTracks"),
  33   fESD(0),fHistoManager(0),fCentrSelector(0),fTrackCuts(0),fTrackCutsNoDCA(0),fOfflineTrigger(0), fIsMC(0),fIsTPCOnly(0), fTriggerAnalysis(0),fPIDResponse(0),fRejectElectrons(0)
  34 {
  35   // constructor
  36
  37   DefineOutput(1, AliAnalysisMultPbTrackHistoManager::Class());
  38   DefineOutput(2, AliESDtrackCuts::Class());
  39   DefineOutput(3, AliAnalysisMultPbCentralitySelector::Class());
  40
  41   fHistoManager = new AliAnalysisMultPbTrackHistoManager("histoManager","Hitograms, Multiplicity, Track analysis");
  42   if(fIsMC) fHistoManager->SetSuffix("MC");
  43
  44 }
  45 AliAnalysisTaskMultPbTracks::AliAnalysisTaskMultPbTracks(const char * name)
  46   : AliAnalysisTaskSE(name),
  47     fESD(0),fHistoManager(0),fCentrSelector(0),fTrackCuts(0),fTrackCutsNoDCA(0),fOfflineTrigger(0),fIsMC(0),fIsTPCOnly(0), fTriggerAnalysis(0),fPIDResponse(0),fRejectElectrons(0)
  48 {
  49   //
  50   // Standard constructur which should be used
  51   //
  52
  53   DefineOutput(1, AliAnalysisMultPbTrackHistoManager::Class());
  54   DefineOutput(2, AliESDtrackCuts::Class());
  55   DefineOutput(3, AliAnalysisMultPbCentralitySelector::Class());
  56
  57   fHistoManager = new AliAnalysisMultPbTrackHistoManager("histoManager","Hitograms, Multiplicity, Track analysis");
  58   if(fIsMC) fHistoManager->SetSuffix("MC");
  59
  60 }
  61
  62 AliAnalysisTaskMultPbTracks::AliAnalysisTaskMultPbTracks(const AliAnalysisTaskMultPbTracks& obj) :
  63   AliAnalysisTaskSE(obj) ,fESD (0), fHistoManager(0), fCentrSelector(0), fTrackCuts(0),fTrackCutsNoDCA(0),fOfflineTrigger(0),fIsMC(0),fIsTPCOnly(0), fTriggerAnalysis(0),fPIDResponse(0),fRejectElectrons(0)
  64 {
  65   //copy ctor
  66   fESD = obj.fESD ;
  67   fHistoManager= obj.fHistoManager;
  68   fTrackCuts  = obj.fTrackCuts;
  69   fTrackCutsNoDCA  = obj.fTrackCutsNoDCA;
  70   fCentrSelector = obj.fCentrSelector;
  71   fOfflineTrigger = obj.fOfflineTrigger;
  72   fIsMC = obj.fIsMC;
  73   fIsTPCOnly = obj.fIsTPCOnly;
  74   fTriggerAnalysis = obj.fTriggerAnalysis;
  75
  76 }
  77
  78 AliAnalysisTaskMultPbTracks::~AliAnalysisTaskMultPbTracks(){
  79   // destructor
  80
  81   if(!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
  82     if(fHistoManager) {
  83       delete fHistoManager;
  84       fHistoManager = 0;
  85       delete fTriggerAnalysis;
  86       fTriggerAnalysis=0;
  87     }
  88   }
  89   // Histo list should not be destroyed: fListWrapper is owner!
  90
  91 }
  92 void AliAnalysisTaskMultPbTracks::UserCreateOutputObjects()
  93 {
  94   // Called once
  95
  96   // For the DCA distributions, we want to use exactly the same cuts
  97   // as for the other distributions, with the exception of the DCA cut
  98   fTrackCutsNoDCA = new AliESDtrackCuts(*fTrackCuts); // clone cuts
  99   // Reset all DCA cuts; FIXME: is this all?
 100   fTrackCutsNoDCA->SetMaxDCAToVertexXY();
 101   fTrackCutsNoDCA->SetMaxDCAToVertexZ ();
 102   fTrackCutsNoDCA->SetMaxDCAToVertexXYPtDep();
 103   fTrackCutsNoDCA->SetMaxDCAToVertexZPtDep();
 104
 105   fTriggerAnalysis = new AliTriggerAnalysis();
 106   fTriggerAnalysis->SetAnalyzeMC(fIsMC);
 107
 108   //The common PID object can then be retrieved from the input handler. This can naturally be done in the UserCreateOutputObjects:
 109   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 110   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 111   fPIDResponse = inputHandler->GetPIDResponse();
 112
 113   PostData(1,fHistoManager);
 114   PostData(2,fTrackCuts);
 115   PostData(3,fCentrSelector);
 116
 117
 118 }
 119
 120
 121 void AliAnalysisTaskMultPbTracks::UserExec(Option_t *)
 122 {
 123   // User code
 124
 125
 126   /* PostData(0) is taken care of by AliAnalysisTaskSE */
 127   PostData(1,fHistoManager);
 128   PostData(2,fTrackCuts);
 129   PostData(3,fCentrSelector);
 130
 131   // Cache (some) histogram pointers to avoid loop in the histo manager
 132   static TH3D * hTracks  [AliAnalysisMultPbTrackHistoManager::kNHistos];
 133   static TH2D * hDCA     [AliAnalysisMultPbTrackHistoManager::kNHistos];
 134   static TH1D * hNTracks [AliAnalysisMultPbTrackHistoManager::kNHistos];
 135   hTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen]        = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoGen       );
 136   hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRec]        = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec       );
 137   hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecPrim]    = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim   );
 138   hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecFake]    = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecFake   );
 139   hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat]  = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat );
 140   hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak] = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak);
 141
 142   hDCA[AliAnalysisMultPbTrackHistoManager::kHistoGen]        = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoGen       );
 143   hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRec]        = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec       );
 144   hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecPrim]    = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim   );
 145   hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecFake]    = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecFake   );
 146   hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat]  = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat );
 147   hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak] = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak);
 148
 149   hNTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen]        = fHistoManager->GetHistoMult(AliAnalysisMultPbTrackHistoManager::kHistoGen );
 150   hNTracks[AliAnalysisMultPbTrackHistoManager::kHistoRec]        = fHistoManager->GetHistoMult(AliAnalysisMultPbTrackHistoManager::kHistoRec );
 151
 152
 153   fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec)->Reset();
 154
 155   fESD = dynamic_cast<AliESDEvent*>(fInputEvent);
 156   if (strcmp(fESD->ClassName(),"AliESDEvent")) {
 157     AliFatal("Not processing ESDs");
 158   }
 159
 160   fHistoManager->GetHistoStats()->Fill(AliAnalysisMultPbTrackHistoManager::kStatAll);
 161
 162
 163   // Centrality selection
 164   Bool_t isCentralitySelected = fCentrSelector->IsCentralityBinSelected(fESD,fTrackCuts);
 165   if(!isCentralitySelected) return;
 166
 167   // AliESDCentrality *centrality = fESD->GetCentrality();
 168   // if(!centrality) {
 169   //   AliFatal("Centrality object not available");
 170   // }
 171   // else {
 172   //   Int_t centrBin = centrality->GetCentralityClass5(fCentralityEstimator.Data()) ;
 173
 174   //   if (centrBin != fCentrBin && fCentrBin != -1) return;
 175   // }
 176
 177   fHistoManager->GetHistoStats()->Fill(AliAnalysisMultPbTrackHistoManager::kStatCentr);
 178
 179   Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & fOfflineTrigger);
 180
 181   if(!isSelected) return;
 182   fHistoManager->GetHistoStats()->Fill(AliAnalysisMultPbTrackHistoManager::kStatPhysSel);
 183
 184
 185   if (fIsMC) {
 186     // Get MC vertex
 187     //FIXME: which vertex do I take for MC?
 188     TArrayF   vertex;
 189     fMCEvent->GenEventHeader()->PrimaryVertex(vertex);
 190     Float_t zvGen = vertex[2];
 191
 192     if (!fMCEvent) {
 193       AliError("No MC info found");
 194     } else {
 195
 196       //loop on the MC event, only over primaries, which are always
 197       //      the first in stack.
 198       // WRONG: D&B decays are produced in the transport... Need To Loop over all particles
 199       //      Int_t nMCTracks = fMCEvent->GetNumberOfPrimaries();
 200       Int_t nMCTracks = fMCEvent->GetNumberOfTracks();
 201       Int_t nPhysicalPrimaries = 0;
 202       for (Int_t ipart=0; ipart<nMCTracks; ipart++) {
 203
 204         AliMCParticle *mcPart  = (AliMCParticle*)fMCEvent->GetTrack(ipart);
 205
 206         // We don't care about neutrals and non-physical primaries
 207         if(mcPart->Charge() == 0) continue;
 208
 209         // If we are cutting away electrons, also exclude them here:
 210         if(fRejectElectrons)
 211           if (TMath::Abs(mcPart->PdgCode())==11) continue;
 212
 213         //check if current particle is a physical primary and fill map of physical primaries
 214         if(!IsPhysicalPrimaryAndTransportBit(ipart)) continue;
 215         fHistoManager->FillSpeciesMap(mcPart->PdgCode());
 216
 217         if(TMath::Abs(mcPart->Zv()-zvGen)>1e-6) {
 218           // This cures a bug in Hijing
 219           // A little hack here: I put those in the underflow bin of the process type to keep track of them
 220           fHistoManager->GetHistoProcess(AliAnalysisMultPbTrackHistoManager::kHistoGen)->Fill(-1);
 221           continue;
 222         }
 223
 224         nPhysicalPrimaries++;
 225         // Fill species histo and particle species
 226         fHistoManager->GetHistoProcess(AliAnalysisMultPbTrackHistoManager::kHistoGen)->Fill(mcPart->Particle()->GetUniqueID());
 227         fHistoManager->FillParticleID(AliAnalysisMultPbTrackHistoManager::kHistoGen, mcPart);
 228
 229         //      Float_t zv = vtxESD->GetZ();
 230         // Fill generated histo
 231         hTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen]->Fill(mcPart->Pt(),mcPart->Eta(),zvGen);
 232         Int_t partCode = fHistoManager->GetLocalParticleID(mcPart);
 233         //      cout << "Part " << partCode << endl;
 234         if (partCode == AliAnalysisMultPbTrackHistoManager::kPartPiPlus  ||
 235             partCode == AliAnalysisMultPbTrackHistoManager::kPartPiMinus ||
 236             partCode == AliAnalysisMultPbTrackHistoManager::kPartKPlus   ||
 237             partCode == AliAnalysisMultPbTrackHistoManager::kPartKMinus  ||
 238             partCode == AliAnalysisMultPbTrackHistoManager::kPartP       ||
 239             partCode == AliAnalysisMultPbTrackHistoManager::kPartPBar
 240             ){
 241           //      cout << "Found " << partCode << endl;
 242
 243           fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoGen, partCode);
 244         }
 245       }
 246       hNTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen]->Fill(nPhysicalPrimaries);
 247       fHistoManager->GetHistoVzEvent(AliAnalysisMultPbTrackHistoManager::kHistoGen)->Fill(zvGen);
 248
 249
 250     }
 251   }
 252
 253
 254   const AliESDVertex* vtxESD = fESD->GetPrimaryVertex();
 255   if(!vtxESD) return;
 256   // TODO: leave the cuts here or find a better place?)
 257   // Quality cut on vertexer Z, as suggested by Francesco Prino
 258   if(vtxESD->IsFromVertexerZ()) {
 259     if (vtxESD->GetNContributors() <= 0) return;
 260     if (vtxESD->GetDispersion() > 0.04) return;
 261     if (vtxESD->GetZRes() > 0.25) return;
 262   }
 263   // "Beam gas" vertex cut
 264   const AliESDVertex * vtxESDTPC= fESD->GetPrimaryVertexTPC();
 265   if(vtxESDTPC->GetNContributors()<1) return;
 266   if (vtxESDTPC->GetNContributors()<(-10.+0.25*fESD->GetMultiplicity()->GetNumberOfITSClusters(0)))     return;
 267
 268   // Fill  statistics
 269   fHistoManager->GetHistoStats()->Fill(AliAnalysisMultPbTrackHistoManager::kStatVtx);
 270
 271   // ZDC cut, only ZNs
 272   Bool_t zdcA   = fTriggerAnalysis->ZDCTDCTrigger(fESD, AliTriggerAnalysis::kASide, kTRUE, kFALSE) ;
 273   Bool_t zdcC   = fTriggerAnalysis->ZDCTDCTrigger(fESD, AliTriggerAnalysis::kCSide, kTRUE, kFALSE) ;
 274
 275   if (!(zdcA && zdcC) && (!fIsMC)) return;
 276   fHistoManager->GetHistoStats()->Fill(AliAnalysisMultPbTrackHistoManager::kStatZDCCut);
 277
 278
 279   if(TMath::Abs(vtxESD->GetZ()) > 10) return;
 280   fHistoManager->GetHistoStats()->Fill(AliAnalysisMultPbTrackHistoManager::kStatVtxRangeCut);
 281
 282   // FIXME
 283   Float_t ntracksOk = fTrackCuts->CountAcceptedTracks(fESD);
 284   const AliMultiplicity * mult = fESD->GetMultiplicity();
 285   if (ntracksOk < 1000 &&  mult->GetNumberOfITSClusters(1) > 4500) {
 286     Float_t dummy;
 287     Printf("IEV: %d, Orbit: %x, Period: %d, BC: %d\n",fESD->GetEventNumberInFile(), fESD->GetOrbitNumber(),fESD->GetPeriodNumber(),fESD->GetBunchCrossNumber());
 288     printf("%s ----> Processing event # %lld\n", CurrentFileName(), Entry());
 289     cout << "File " << AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName() << endl;
 290     cout << "Nt " << ntracksOk << ", "<< fESD->GetNumberOfTracks() << " V0 " <<  fCentrSelector->GetCorrV0(fESD,dummy)  << " SPD1 " << mult->GetNumberOfITSClusters(1) << endl;
 291   }
 292
 293
 294
 295   // Fill Vertex
 296
 297   fHistoManager->GetHistoVzEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec)->Fill(vtxESD->GetZ());
 298
 299   // loop on tracks
 300   Int_t ntracks = fESD->GetNumberOfTracks();
 301   Int_t acceptedTracks = 0;
 302
 303   for (Int_t itrack = 0; itrack<ntracks; itrack++) {
 304     AliESDtrack * esdTrack = fIsTPCOnly ? fTrackCuts->GetTPCOnlyTrack(fESD,itrack) : fESD->GetTrack(itrack);// FIXMEL TrackCuts or TrackCutsNoDCA for the TPC?
 305     if (!esdTrack) continue;
 306
 307     // Fill DCA distibutions, without the DCA cut, Fill the other stuff, with all cuts!
 308     Bool_t accepted = fTrackCuts->AcceptTrack(esdTrack);
 309     Bool_t acceptedNoDCA = fTrackCutsNoDCA->AcceptTrack(esdTrack);
 310
 311     //fHistoManager->GetHistoElectronCutQA()->Fill(esdTrack->Pt(), esdTrack->GetTPCsignal()); FIXME: temporary hack
 312
 313     if(fRejectElectrons && acceptedNoDCA) { // don't check this for tracks to be anyways rejected
 314       Bool_t isEl = IsElectron(esdTrack);
 315       if(isEl) {
 316         //      cout << " --> FOUND!" << endl;
 317         fHistoManager->GetHistoElectronCutQA()->Fill(esdTrack->Pt(), esdTrack->GetTPCsignal());
 318       }
 319       accepted = accepted && !isEl;
 320       acceptedNoDCA = acceptedNoDCA && !isEl;
 321     }
 322
 323     if(accepted) acceptedTracks++;
 324
 325     // Compute weighted offset
 326     // TODO: other choiches for the DCA?
 327     Double_t weightedDCA = 10;
 328
 329
 330 #if defined WEIGHTED_DCA
 331     Double_t b = fESD->GetMagneticField();
 332     Double_t dca[2];
 333     Double_t cov[3];
 334     if (esdTrack->PropagateToDCA(vtxESD, b,3., dca, cov)) {
 335       Double_t det = cov[0]*cov[2]-cov[1]*cov[1];
 336       if (det<=0) {
 337         AliError("DCA Covariance matrix is not positive definite!");
 338       }
 339       else {
 340         weightedDCA = (dca[0]*dca[0]*cov[2]+dca[1]*dca[1]*cov[0]-2*dca[0]*dca[1]*cov[1])/det;
 341         weightedDCA = weightedDCA>0 ? TMath::Sqrt(weightedDCA) : 10;
 342       }
 343       //      printf("dR:%e dZ%e  sigR2:%e sigRZ:%e sigZ2:%e\n",dca[0],dca[1],cov[0],cov[1],cov[2]);
 344     }
 345
 346 #elif defined TRANSVERSE_DCA
 347     Float_t xz[2];
 348     esdTrack->GetDZ(vtxESD->GetX(),vtxESD->GetY(),vtxESD->GetZ(), fESD->GetMagneticField(), xz);
 349     weightedDCA = xz[0];
 350 #endif
 351
 352
 353     // Alternative: p*DCA
 354     // Float_t xz[2];
 355     // esdTrack->GetDZ(vtxESD->GetX(),vtxESD->GetY(),vtxESD->GetZ(), fESD->GetMagneticField(), xz);
 356     // Float_t dca = xz[0]*esdTrack->P();
 357
 358   // This cures a bug in Hijing (displaced primaries)
 359   if (fIsMC) {
 360     Int_t label = TMath::Abs(esdTrack->GetLabel()); // no fakes!!!
 361     if (IsPhysicalPrimaryAndTransportBit(label)) {
 362       AliMCParticle *mcPart  =  (AliMCParticle*)fMCEvent->GetTrack(label);
 363       if(!mcPart) AliFatal("No particle");
 364       TArrayF   vertex;
 365       fMCEvent->GenEventHeader()->PrimaryVertex(vertex);
 366       Float_t zvGen = vertex[2];
 367       if(TMath::Abs(mcPart->Zv()-zvGen)>1e-6) continue;
 368     }
 369   }
 370   // for each track
 371   if(accepted) {
 372     hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRec]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
 373     if (TMath::Abs(esdTrack->Eta())<0.5 && TMath::Abs(vtxESD->GetZ())<7)
 374       fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec)->Fill(esdTrack->Pt());
 375   }
 376   if(acceptedNoDCA)
 377     hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRec]->Fill(weightedDCA,esdTrack->Pt());
 378
 379     // Fill separately primaries and secondaries
 380     // FIXME: fakes? Is this the correct way to account for them?
 381     // Get label and corresponding mcPart;
 382     if (fIsMC) {
 383       Int_t label = TMath::Abs(esdTrack->GetLabel()); // no fakes!!!
 384       AliMCParticle *mcPart  = label < 0 ? 0 : (AliMCParticle*)fMCEvent->GetTrack(label);
 385       if (!mcPart)  {
 386         if(accepted)
 387           hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecFake]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
 388         if(acceptedNoDCA)
 389           hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecFake]->Fill(weightedDCA,esdTrack->Pt());
 390       }
 391       else {
 392         if(IsPhysicalPrimaryAndTransportBit(label)) {
 393           // Fill species histo
 394           fHistoManager->GetHistoProcess(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->Fill(mcPart->Particle()->GetUniqueID());
 395           if(accepted) {
 396             hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecPrim]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
 397             Int_t partCode = fHistoManager->GetLocalParticleID(mcPart);
 398             if (partCode == AliAnalysisMultPbTrackHistoManager::kPartPiPlus  ||
 399                 partCode == AliAnalysisMultPbTrackHistoManager::kPartPiMinus ||
 400                 partCode == AliAnalysisMultPbTrackHistoManager::kPartKPlus   ||
 401                 partCode == AliAnalysisMultPbTrackHistoManager::kPartKMinus  ||
 402                 partCode == AliAnalysisMultPbTrackHistoManager::kPartP       ||
 403                 partCode == AliAnalysisMultPbTrackHistoManager::kPartPBar    ||
 404                 partCode == AliAnalysisMultPbTrackHistoManager::kPartLPlus   ||
 405                 partCode == AliAnalysisMultPbTrackHistoManager::kPartKMinus
 406                 )
 407               fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim, partCode);
 408           }
 409           if(acceptedNoDCA)
 410             hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecPrim]->Fill(weightedDCA,esdTrack->Pt());
 411         }
 412         else {
 413           Int_t mfl=0;
 414           Int_t indexMoth=mcPart->Particle()->GetFirstMother();
 415           if(indexMoth>=0){
 416             TParticle* moth = fMCEvent->Stack()->Particle(indexMoth);
 417             Float_t codemoth = TMath::Abs(moth->GetPdgCode());
 418             mfl = Int_t (codemoth/ TMath::Power(10, Int_t(TMath::Log10(codemoth))));
 419           }
 420           if(mfl==3){ // strangeness
 421             fHistoManager->GetHistoProcess(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak)->Fill(mcPart->Particle()->GetUniqueID());
 422             if(accepted)
 423               hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
 424             if(acceptedNoDCA)
 425               hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak]->Fill(weightedDCA,esdTrack->Pt());
 426           }else{ // material
 427             fHistoManager->GetHistoProcess(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat)->Fill(mcPart->Particle()->GetUniqueID());
 428             if(accepted)
 429               hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
 430             if(acceptedNoDCA)
 431               hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat]->Fill(weightedDCA,esdTrack->Pt());
 432
 433           }
 434
 435
 436         }
 437       }
 438     }
 439
 440
 441   }
 442   //  cout << acceptedTracks << endl;
 443
 444   hNTracks[AliAnalysisMultPbTrackHistoManager::kHistoRec]  ->Fill(acceptedTracks);
 445
 446   // FIXME: this used to be filled with rescaled V0. I'm using raw V0 to test some newer production here.
 447   AliESDVZERO* esdV0 = fESD->GetVZEROData();
 448   Float_t multV0A=esdV0->GetMTotV0A();
 449   Float_t multV0C=esdV0->GetMTotV0C();
 450   Float_t multV0 = multV0A+multV0C;
 451   fHistoManager->GetHistoV0vsNtracks(AliAnalysisMultPbTrackHistoManager::kHistoRec)->Fill(acceptedTracks, multV0);
 452   // FIXME: the old version:
 453   // Float_t v0;
 454   // fHistoManager->GetHistoV0vsNtracks(AliAnalysisMultPbTrackHistoManager::kHistoRec)->Fill(acceptedTracks, fCentrSelector->GetCorrV0(fESD,v0));
 455
 456   // Fill <pt> analysis histos
 457   fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec)->Scale(1,"width");// correct bin width
 458   if (fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec)->GetEntries()>0)
 459     fHistoManager->GetHistoMeanPt(AliAnalysisMultPbTrackHistoManager::kHistoRec)->Fill(fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec)->GetMean());
 460   if (fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec)->GetMean() >
 461       fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRecHighestMeanPt)->GetMean()) {
 462     // Found a new highest pt: first reset than add it to the container for the highest
 463     fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRecHighestMeanPt)->Reset();
 464     fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRecHighestMeanPt)
 465       ->Add(fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec));
 466   }
 467   if ((fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec)->GetMean() <
 468        fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRecLowestMeanPt)->GetMean() &&
 469        fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec)->GetEntries()>0) ||
 470       !(fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRecLowestMeanPt)->GetEntries()>0)) {
 471     // Found a new lowest pt: first reset than add it to the container for the lowest
 472     fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRecLowestMeanPt)->Reset();
 473     fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRecLowestMeanPt)
 474       ->Add(fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec));
 475   }
 476 }
 477
 478 void   AliAnalysisTaskMultPbTracks::Terminate(Option_t *){
 479   // terminate
 480
 481 }
 482
 483
 484 Bool_t AliAnalysisTaskMultPbTracks::IsPhysicalPrimaryAndTransportBit(Int_t ipart) {
 485   // Check if it is a primary and if it was transported
 486   Bool_t physprim=fMCEvent->IsPhysicalPrimary(ipart);
 487   if (!physprim) return kFALSE;
 488   Bool_t transported = ((AliMCParticle*)fMCEvent->GetTrack(ipart))->Particle()->TestBit(kTransportBit);
 489   if(!transported) return kFALSE;
 490
 491   return kTRUE;
 492
 493 }
 494
 495 Bool_t AliAnalysisTaskMultPbTracks::IsElectron(AliESDtrack * esdTrack) {
 496   Bool_t isElectron = kFALSE;
 497     // ((fPIDResponse->NumberOfSigmasTPC(esdTrack,AliPID::kElectron) < 2 ) &&
 498     //  (fPIDResponse->NumberOfSigmasTPC(esdTrack,AliPID::kPion    ) > 1 ) &&
 499     //  (fPIDResponse->NumberOfSigmasTPC(esdTrack,AliPID::kProton  ) > 1 ) &&
 500     //  (fPIDResponse->NumberOfSigmasTPC(esdTrack,AliPID::kKaon    ) > 1 )
 501     //  );//FIXME SKIP ELECTRONS below p = 1.2 gev, make configurable. Keep particles if they are in the crossing
 502
 503   // cout << " - TPC: " << fPIDResponse->NumberOfSigmasTPC(esdTrack,AliPID::kElectron) << ", TOF: " << fPIDResponse->NumberOfSigmasTOF(esdTrack,AliPID::kElectron);
 504   // printf("[%f,%f,%d]\n", esdTrack->Pt(), esdTrack->Eta(), esdTrack->GetTPCclusters(0));
 505   // if (fPIDResponse->NumberOfSigmasTPC(esdTrack,AliPID::kElectron) <0){
 506
 507
 508   // }
 509   isElectron = (TMath::Abs(fPIDResponse->NumberOfSigmasTPC(esdTrack,AliPID::kElectron)) < 2) && (TMath::Abs(fPIDResponse->NumberOfSigmasTOF(esdTrack,AliPID::kElectron) < 3));
 510
 511   return isElectron;
 512 }
 513
 514 // void AliAnalysisTaskEvil::PrintProcInfo()
 515 // {
 516 //   ProcInfo_t info;
 517 //   gSystem->GetProcInfo(&info);
 518 //   AliInfo(Form("fMemResident=%10ld kB  fMemVirtual=%10ld kB",info.fMemResident,info.fMemVirtual));
 519 // }