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