]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG0/multPbPb/AliAnalysisTaskMultPbTracks.cxx
Glauber files for real data (Alberica)
[u/mrichter/AliRoot.git] / PWG0 / multPbPb / AliAnalysisTaskMultPbTracks.cxx
index bd33b8e235cbeac92a6ee364c1c027b74027c8a6..8076d0cca7061146e50940c3dd9c4578ea346b38 100644 (file)
@@ -3,51 +3,54 @@
 // Author: Michele Floris, CERN
 // TODO:
 // - Add chi2/cluster plot for primary, secondaries and fakes
-
+// FIXME:
+// - re-implement centrality estimator to cut on tracks and multiplicity
 
 #include "AliAnalysisTaskMultPbTracks.h"
 #include "AliESDInputHandler.h"
 #include "AliAnalysisMultPbTrackHistoManager.h"
 #include "AliAnalysisManager.h"
 #include "AliESDtrackCuts.h"
-#include "AliAnalysisMultPbCentralitySelector.h"
 #include "AliMCEvent.h"
 #include "AliStack.h"
 #include "TH1I.h"
 #include "TH3D.h"
 #include "AliMCParticle.h"
 #include "AliGenEventHeader.h"
-
+#include "AliESDCentrality.h"
+#include "AliMultiplicity.h"
 #include <iostream>
+#include "AliAnalysisMultPbCentralitySelector.h"
 
 using namespace std;
 
 ClassImp(AliAnalysisTaskMultPbTracks)
 
 AliAnalysisTaskMultPbTracks::AliAnalysisTaskMultPbTracks()
-: AliAnalysisTaskSE("TaskMultPbTracks"),fESD(0),fHistoManager(0),fCentrSelector(0),fTrackCuts(0),fTrackCutsNoDCA(0),fIsMC(0)
+: AliAnalysisTaskSE("TaskMultPbTracks"),
+  fESD(0),fHistoManager(0),fCentrSelector(0),fTrackCuts(0),fTrackCutsNoDCA(0),fIsMC(0)
 {
   // constructor
 
   DefineOutput(1, AliAnalysisMultPbTrackHistoManager::Class());
-  DefineOutput(2, AliAnalysisMultPbCentralitySelector::Class());
-  DefineOutput(3, AliESDtrackCuts::Class());
-  //  DefineOutput(2, TH1I::Class());
+  DefineOutput(2, AliESDtrackCuts::Class());
+  DefineOutput(3, AliAnalysisMultPbCentralitySelector::Class());
 
   fHistoManager = new AliAnalysisMultPbTrackHistoManager("histoManager","Hitograms, Multiplicity, Track analysis");
   if(fIsMC) fHistoManager->SetSuffix("MC");
 
 }
 AliAnalysisTaskMultPbTracks::AliAnalysisTaskMultPbTracks(const char * name)
-  : AliAnalysisTaskSE(name),fESD(0),fHistoManager(0),fCentrSelector(0),fTrackCuts(0),fTrackCutsNoDCA(0),fIsMC(0)
+  : AliAnalysisTaskSE(name),
+    fESD(0),fHistoManager(0),fCentrSelector(0),fTrackCuts(0),fTrackCutsNoDCA(0),fIsMC(0)
 {
   //
   // Standard constructur which should be used
   //
 
   DefineOutput(1, AliAnalysisMultPbTrackHistoManager::Class());
-  DefineOutput(2, AliAnalysisMultPbCentralitySelector::Class());
-  DefineOutput(3, AliESDtrackCuts::Class());
+  DefineOutput(2, AliESDtrackCuts::Class());
+  DefineOutput(3, AliAnalysisMultPbCentralitySelector::Class());
 
   fHistoManager = new AliAnalysisMultPbTrackHistoManager("histoManager","Hitograms, Multiplicity, Track analysis");
   if(fIsMC) fHistoManager->SetSuffix("MC");
@@ -60,9 +63,9 @@ AliAnalysisTaskMultPbTracks::AliAnalysisTaskMultPbTracks(const AliAnalysisTaskMu
   //copy ctor
   fESD = obj.fESD ;
   fHistoManager= obj.fHistoManager;
-  fCentrSelector = obj.fCentrSelector;
   fTrackCuts  = obj.fTrackCuts;
   fTrackCutsNoDCA  = obj.fTrackCutsNoDCA;
+  fCentrSelector = obj.fCentrSelector;
 
 }
 
@@ -100,11 +103,13 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *)
 
   /* PostData(0) is taken care of by AliAnalysisTaskSE */
   PostData(1,fHistoManager);
-  PostData(2,fCentrSelector);
-  PostData(3,fTrackCuts);
+  PostData(2,fTrackCuts);
+  PostData(3,fCentrSelector);
 
-  TH3D * hTracks[AliAnalysisMultPbTrackHistoManager::kNHistos];
-  TH1D * hDCA   [AliAnalysisMultPbTrackHistoManager::kNHistos];
+  // Cache histogram pointers
+  static TH3D * hTracks  [AliAnalysisMultPbTrackHistoManager::kNHistos];
+  static TH1D * hDCA     [AliAnalysisMultPbTrackHistoManager::kNHistos];
+  static TH1D * hNTracks [AliAnalysisMultPbTrackHistoManager::kNHistos];
   hTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen]        = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoGen       );
   hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRec]        = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec       );
   hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecPrim]    = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim   );
@@ -119,6 +124,8 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *)
   hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat]  = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat );
   hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak] = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak);
 
+  hNTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen]        = fHistoManager->GetHistoMult(AliAnalysisMultPbTrackHistoManager::kHistoGen );
+  hNTracks[AliAnalysisMultPbTrackHistoManager::kHistoRec]        = fHistoManager->GetHistoMult(AliAnalysisMultPbTrackHistoManager::kHistoRec );
 
   fESD = dynamic_cast<AliESDEvent*>(fInputEvent);
   if (strcmp(fESD->ClassName(),"AliESDEvent")) {
@@ -128,9 +135,26 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *)
   // FIXME: use physics selection here to keep track of events lost?
   fHistoManager->GetHistoStats()->Fill(AliAnalysisMultPbTrackHistoManager::kStatAll);
 
-  Bool_t isCentralitySelected = fCentrSelector->IsSelected(fESD);  
+  Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & fOfflineTrigger);
+
+  if(!isSelected) return;
+  fHistoManager->GetHistoStats()->Fill(AliAnalysisMultPbTrackHistoManager::kStatPhysSel);
+
+
+  // Centrality selection
+  Bool_t isCentralitySelected = fCentrSelector->IsCentralityBinSelected(fESD,fTrackCuts);  
   if(!isCentralitySelected) return;
 
+  // AliESDCentrality *centrality = fESD->GetCentrality();
+  // if(!centrality) {
+  //   AliFatal("Centrality object not available"); 
+  // }
+  // else {
+  //   Int_t centrBin = centrality->GetCentralityClass5(fCentralityEstimator.Data()) ;
+    
+  //   if (centrBin != fCentrBin && fCentrBin != -1) return;
+  // }
+
   fHistoManager->GetHistoStats()->Fill(AliAnalysisMultPbTrackHistoManager::kStatCentr);
 
   if (fIsMC) {
@@ -140,8 +164,10 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *)
       AliError("No MC info found");
     } else {
       
-      //loop on the MC event
-      Int_t nMCTracks = fMCEvent->GetNumberOfTracks();
+      //loop on the MC event, only over primaries, which are always
+      //      the first in stack.
+      Int_t nMCTracks = fMCEvent->GetNumberOfPrimaries();
+      Int_t nPhysicalPrimaries = 0;
       for (Int_t ipart=0; ipart<nMCTracks; ipart++) { 
        
        AliMCParticle *mcPart  = (AliMCParticle*)fMCEvent->GetTrack(ipart);
@@ -149,8 +175,13 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *)
        // We don't care about neutrals and non-physical primaries
        if(mcPart->Charge() == 0) continue;
 
-       // FIXME: add kTransportBit
-       if(!fMCEvent->Stack()->IsPhysicalPrimary(ipart)) continue;
+       //check if current particle is a physical primary
+       if(!IsPhysicalPrimaryAndTransportBit(ipart)) continue;
+       nPhysicalPrimaries++;
+       // Fill species histo
+       fHistoManager->GetHistoSpecies(AliAnalysisMultPbTrackHistoManager::kHistoGen)->Fill(mcPart->Particle()->GetUniqueID());
+
        
        // Get MC vertex
        //FIXME: which vertex do I take for MC?
@@ -162,9 +193,12 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *)
        hTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen]->Fill(mcPart->Pt(),mcPart->Eta(),zv);
        
       }
+      hNTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen]->Fill(nPhysicalPrimaries);
+
     }
   }
   
+
   // FIXME: shall I take the primary vertex?
   const AliESDVertex* vtxESD = fESD->GetPrimaryVertex();
   if(!vtxESD) return;
@@ -173,7 +207,7 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *)
 
   // loop on tracks
   Int_t ntracks = fESD->GetNumberOfTracks();
-
+  Int_t acceptedTracks = 0;
 
   for (Int_t itrack = 0; itrack<ntracks; itrack++) {    
     AliESDtrack * esdTrack = fESD->GetTrack(itrack);
@@ -182,6 +216,8 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *)
     Bool_t accepted = fTrackCuts->AcceptTrack(esdTrack);
     Bool_t acceptedNoDCA = fTrackCutsNoDCA->AcceptTrack(esdTrack);
 
+    if(accepted) acceptedTracks++;
+
     // Compute weighted offset
     // TODO: other choiches for the DCA?
     Double_t b = fESD->GetMagneticField();
@@ -223,7 +259,7 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *)
     if (fIsMC) {
       //      Int_t label = TMath::Abs(esdTrack->GetLabel()); // no fakes!!!
       Int_t label = esdTrack->GetLabel(); // 
-      AliMCParticle *mcPart  = (AliMCParticle*)fMCEvent->GetTrack(label);
+      AliMCParticle *mcPart  = label < 0 ? 0 : (AliMCParticle*)fMCEvent->GetTrack(label);
       if (!mcPart)  {
        if(accepted)
          hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecFake]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
@@ -231,7 +267,9 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *)
          hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecFake]->Fill(weightedDCA);
       }
       else {
-       if(fMCEvent->Stack()->IsPhysicalPrimary(label)) {
+       if(IsPhysicalPrimaryAndTransportBit(label)) {
+         // Fill species histo
+         fHistoManager->GetHistoSpecies(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->Fill(mcPart->Particle()->GetUniqueID());
          if(accepted)
            hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecPrim]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
          if(acceptedNoDCA)
@@ -246,11 +284,13 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *)
            mfl = Int_t (codemoth/ TMath::Power(10, Int_t(TMath::Log10(codemoth))));
          }
          if(mfl==3){ // strangeness
+           fHistoManager->GetHistoSpecies(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak)->Fill(mcPart->Particle()->GetUniqueID());
            if(accepted)
              hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
            if(acceptedNoDCA)
              hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak]->Fill(weightedDCA);      
          }else{ // material
+           fHistoManager->GetHistoSpecies(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat)->Fill(mcPart->Particle()->GetUniqueID());
            if(accepted)
              hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
            if(acceptedNoDCA)
@@ -265,6 +305,11 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *)
 
 
   }
+  //  cout << acceptedTracks << endl;
+  
+  hNTracks[AliAnalysisMultPbTrackHistoManager::kHistoRec]  ->Fill(acceptedTracks);
+  // FIXME
+  //  hNTracks[AliAnalysisMultPbTrackHistoManager::kHistoRec]  ->Fill(fESD->GetMultiplicity()->GetNumberOfTracklets());
 
 
 }
@@ -275,3 +320,13 @@ void   AliAnalysisTaskMultPbTracks::Terminate(Option_t *){
 }
 
 
+Bool_t AliAnalysisTaskMultPbTracks::IsPhysicalPrimaryAndTransportBit(Int_t ipart) {
+
+  Bool_t physprim=fMCEvent->IsPhysicalPrimary(ipart);
+  if (!physprim) return kFALSE;
+  Bool_t transported = ((AliMCParticle*)fMCEvent->GetTrack(ipart))->Particle()->TestBit(kTransportBit);
+  if(!transported) return kFALSE;
+
+  return kTRUE;
+
+}