PWG0/FirstPhysics/AliAnalysisTaskFirstPhysics.cxx

   1
   2 #include "AliAnalysisTaskFirstPhysics.h"
   3
   4 #include "TH1D.h"
   5 #include "TH2D.h"
   6 #include "TList.h"
   7 #include "TString.h"
   8
   9 #include "AliAnalysisTaskSE.h"
  10 #include "AliAnalysisManager.h"
  11 #include "AliESDtrackCuts.h"
  12 #include "AliESDEvent.h"
  13 #include "AliESDHeader.h"
  14 #include "AliMCEvent.h"
  15 #include "AliMultiplicity.h"
  16
  17 #include "AliGenEventHeader.h"
  18 #include "AliGenPythiaEventHeader.h"
  19 #include "AliGenDPMjetEventHeader.h"
  20 #include "AliStack.h"
  21
  22 ClassImp(AliAnalysisTaskFirstPhysics)
  23
  24 //________________________________________________________________________
  25 AliAnalysisTaskFirstPhysics::AliAnalysisTaskFirstPhysics(const char *name) :
  26   AliAnalysisTaskSE(name),
  27   fESD(0),
  28   fMCEvent(0),
  29   fOutput(0),
  30   fbReadMC(0),
  31   fMCProcessType(kProcND),
  32   fTrigger(0),
  33   fCutTrackPtMin(0.15),
  34   fCutTrackPtMax(100),
  35   fCutEta(0.8),
  36   fCutVertexZ(10)
  37 {
  38   DefineOutput(1, TList::Class());
  39   for (Int_t i = 0; i < knTrackCuts; i ++) {
  40     fTrackCuts[i] = 0;
  41   }
  42 }
  43
  44 //________________________________________________________________________
  45 AliAnalysisTaskFirstPhysics::~AliAnalysisTaskFirstPhysics()
  46 {
  47   // Destructor. Clean-up the output list, but not the histograms that are put inside
  48   // (the list is owner and will clean-up these histograms). Protect in PROOF case.
  49   if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
  50     delete fOutput;
  51   }
  52   for (Int_t i = 0; i < knTrackCuts; i ++) {
  53     delete fTrackCuts[i];
  54     fTrackCuts[i] = 0;
  55   }
  56   if (fTrigger) {
  57     delete fTrigger;
  58   }
  59 }
  60
  61 //________________________________________________________________________
  62 void AliAnalysisTaskFirstPhysics::PrepareOutputList()
  63 {
  64   if (fOutput) {
  65     AliError("fOutput already initialised.");
  66     return;
  67   }
  68   fOutput = new TList();
  69   fOutput->SetOwner();
  70   TH1::SetDefaultSumw2(kTRUE);
  71   if (fTrigger) {
  72     AliError("fTrigger is already initialised.");
  73     return;
  74   }
  75   fTrigger = new AliTriggerAnalysis;
  76 }
  77
  78 //________________________________________________________________________
  79 void AliAnalysisTaskFirstPhysics::PrepareDefaultTrackCuts()
  80 {
  81   // quality cut on ITS+TPC tracks
  82   fTrackCuts[kTrackCutQGlo] = new AliESDtrackCuts();
  83   // TPC
  84   fTrackCuts[kTrackCutQGlo]->SetMinNClustersTPC(70);
  85   fTrackCuts[kTrackCutQGlo]->SetMaxChi2PerClusterTPC(4);
  86   fTrackCuts[kTrackCutQGlo]->SetAcceptKinkDaughters(kFALSE);
  87   fTrackCuts[kTrackCutQGlo]->SetRequireTPCRefit(kTRUE);
  88   // ITS
  89   fTrackCuts[kTrackCutQGlo]->SetRequireITSRefit(kTRUE);
  90   fTrackCuts[kTrackCutQGlo]->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff);
  91   fTrackCuts[kTrackCutQGlo]->SetEtaRange(-fCutEta, fCutEta);
  92   fTrackCuts[kTrackCutQGlo]->SetPtRange(fCutTrackPtMin, fCutTrackPtMax);
  93   // quality cut on ITS_SA tracks (complementary to ITS+TPC)
  94   fTrackCuts[kTrackCutQITS] = new AliESDtrackCuts();
  95   fTrackCuts[kTrackCutQITS]->SetRequireITSRefit(kTRUE);
  96   fTrackCuts[kTrackCutQITS]->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff);
  97   fTrackCuts[kTrackCutQITS]->SetEtaRange(-fCutEta, fCutEta);
  98   // primary selection for tracks with SPD hits
  99   fTrackCuts[kTrackCutDCAwSPD] = new AliESDtrackCuts();
 100   fTrackCuts[kTrackCutDCAwSPD]->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
 101   fTrackCuts[kTrackCutDCAwSPD]->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
 102   fTrackCuts[kTrackCutDCAwSPD]->SetMaxDCAToVertexZ(0.5);
 103   fTrackCuts[kTrackCutDCAwSPD]->SetEtaRange(-fCutEta, fCutEta);
 104   // primary selection for tracks w/o SPD hits
 105   fTrackCuts[kTrackCutDCAwoSPD] = new AliESDtrackCuts();
 106   fTrackCuts[kTrackCutDCAwoSPD]->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kNone);
 107   fTrackCuts[kTrackCutDCAwoSPD]->SetMaxDCAToVertexXYPtDep("1.5*(0.0182+0.0350/pt^1.01)");
 108   fTrackCuts[kTrackCutDCAwoSPD]->SetMaxDCAToVertexZ(0.5);
 109   fTrackCuts[kTrackCutDCAwoSPD]->SetEtaRange(-fCutEta, fCutEta);
 110
 111   // tracks without SPD hits
 112   fTrackCuts[kTrackCutNoSPD] = new AliESDtrackCuts();
 113   fTrackCuts[kTrackCutNoSPD]->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kNone);
 114   // tracks from TPC SA
 115   fTrackCuts[kTrackCutTPConly] = new AliESDtrackCuts();
 116   fTrackCuts[kTrackCutTPConly]->SetRequireTPCStandAlone(kTRUE);
 117 }
 118
 119 //________________________________________________________________________
 120 void AliAnalysisTaskFirstPhysics::UserCreateOutputObjects()
 121 {
 122   PrepareOutputList();
 123   // Define cuts
 124   PrepareDefaultTrackCuts();
 125
 126   // create more histograms here using UserHisto1d and UserHisto2d
 127
 128   PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
 129 }
 130
 131 //________________________________________________________________________
 132 bool AliAnalysisTaskFirstPhysics::GetESDEvent()
 133 {
 134   fESD = 0;
 135   // Get a pointer to the reconstructed event
 136   AliVEvent *event = InputEvent();
 137   if (!event) {
 138     AliError("ERROR: Could not retrieve event.");
 139     return false;
 140   }
 141   // try to access the ESD information
 142   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
 143   if (!esd) {
 144     AliError("InputEvent() is not AliESDEvent.");
 145     return false;
 146   }
 147   fESD = esd;
 148   return true;
 149 }
 150
 151 //________________________________________________________________________
 152 bool AliAnalysisTaskFirstPhysics::CheckVertex()
 153 {
 154   // check for good reconstructed vertex
 155   if (!(fESD->GetPrimaryVertex()->GetStatus()) ||
 156       !(fESD->GetPrimaryVertexSPD()->GetStatus())) {
 157     return false;
 158   }
 159   // if vertex is from spd vertexZ, require more stringent cut
 160   if (fESD->GetPrimaryVertex()->IsFromVertexerZ() &&
 161       (fESD->GetPrimaryVertex()->GetDispersion() > 0.02 ||
 162        fESD->GetPrimaryVertex()->GetZRes() > 0.25 )) {
 163     return false;
 164   }
 165   if (TMath::Abs(fESD->GetPrimaryVertex()->GetZ()) > fCutVertexZ) {
 166     return false;
 167   }
 168   return true;
 169 }
 170
 171 //________________________________________________________________________
 172 bool AliAnalysisTaskFirstPhysics::CheckVertexMC() {
 173   if (!fbReadMC) {
 174     return false;
 175   }
 176   if (TMath::Abs(fMCEvent->GetPrimaryVertex()->GetZ()) > fCutVertexZ) {
 177     return false;
 178   }
 179   return true;
 180 }
 181
 182 //________________________________________________________________________
 183 bool AliAnalysisTaskFirstPhysics::PrepareMCInfo()
 184 {
 185   fbReadMC = true;
 186   fMCEvent = MCEvent();
 187   if (!fMCEvent) {
 188     fbReadMC = false;
 189   }
 190   if (fbReadMC) {
 191     AliDebug(4, TString::Format("MC particles: %d", fMCEvent->GetNumberOfTracks()));
 192     AliDebug(4, TString::Format("number of particles in the stack: %d", fMCEvent->Stack()->GetNtrack()));
 193     AliGenPythiaEventHeader* pythiaGenHeader = 0;
 194     AliGenDPMjetEventHeader* dpmHeader = 0;
 195     pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(fMCEvent->GenEventHeader());
 196     dpmHeader = dynamic_cast<AliGenDPMjetEventHeader*>(fMCEvent->GenEventHeader());
 197     // MC event classification
 198     Int_t iProcessType = 0;
 199     Int_t iSD1flag = 0, iSD2flag = 0, iDDflag = 0,
 200       iCDflag = 0, iELflag = 0, iNDflag = 0;
 201     if (pythiaGenHeader) {
 202       iProcessType = pythiaGenHeader->ProcessType();
 203       iSD1flag = 92;
 204       iSD2flag = 93;
 205       iDDflag = 94;
 206       iELflag = 91;
 207       iCDflag = -2;
 208       iNDflag = -1;
 209       AliDebug(5, TString::Format("Pythia event, %d type", iProcessType));
 210     } else if (dpmHeader) {
 211       iProcessType = dpmHeader->ProcessType();
 212       iSD1flag = 5;
 213       iSD2flag = 6;
 214       iDDflag = 7;
 215       iELflag = 2;
 216       iCDflag = 4;
 217       iNDflag = 1;
 218       AliDebug(5, TString::Format("Phojet event, %d type", iProcessType));
 219     } else {
 220       AliError("neither pythia nor phojet event");
 221     }
 222     if (iProcessType == iSD1flag) {
 223       fMCProcessType = kProcSD1;
 224     } else if (iProcessType == iSD2flag) {
 225       fMCProcessType = kProcSD2;
 226     } else if (iProcessType == iDDflag) {
 227       fMCProcessType = kProcDD;
 228     } else if (iProcessType == iELflag) {
 229       fMCProcessType = kProcEL;
 230     } else if (iProcessType == iCDflag) {
 231       fMCProcessType = kProcCD;
 232     } else if (iProcessType == iNDflag) {
 233       fMCProcessType = kProcND;
 234     } else {
 235       fMCProcessType = kProcIndef;
 236     }
 237   }
 238   return fbReadMC;
 239 }
 240
 241 //________________________________________________________________________
 242 void AliAnalysisTaskFirstPhysics::UserExec(Option_t *)
 243 {
 244   if (!GetESDEvent()) {
 245     AliError("Unable to read the ESD");
 246     return;
 247   }
 248   if (!CheckVertex()) {
 249     return;
 250   }
 251
 252   // custom code comes here
 253
 254   PostData(1, fOutput);
 255 }
 256
 257 //________________________________________________________________________
 258 TH1D* AliAnalysisTaskFirstPhysics::UserHisto1d(const char *name, const char *title, const char *xlabel, Int_t nbinsx, Double_t xlow, Double_t xup)
 259 {
 260   TH1D *hist =  new TH1D(name, title, nbinsx, xlow, xup);
 261   hist->GetXaxis()->SetTitle(xlabel);
 262   fOutput->Add(hist);
 263   return hist;
 264 }
 265
 266 //________________________________________________________________________
 267 TH2D* AliAnalysisTaskFirstPhysics::UserHisto2d(const char *name, const char *title, const char *xlabel, Int_t nbinsx, Double_t xlow, Double_t xup, const char *ylabel, Int_t nbinsy, Double_t ylow, Double_t yup)
 268 {
 269   TH2D *hist =  new TH2D(name, title, nbinsx, xlow, xup, nbinsy, ylow, yup);
 270   hist->GetXaxis()->SetTitle(xlabel);
 271   hist->GetYaxis()->SetTitle(ylabel);
 272   fOutput->Add(hist);
 273   return hist;
 274 }
 275
 276 //________________________________________________________________________
 277 bool AliAnalysisTaskFirstPhysics::GetHisto1FromOutput(const char *name, TH1D *&h) const
 278 {
 279   h = dynamic_cast<TH1D*> (fOutput->FindObject(name));
 280   if (!h) {
 281     AliError(TString::Format("Unable to load histogram from output: %s", name));
 282     return false;
 283   } else {
 284     return true;
 285   }
 286 }
 287
 288 //________________________________________________________________________
 289 bool AliAnalysisTaskFirstPhysics::GetHisto2FromOutput(const char *name, TH2D *&h) const
 290 {
 291   h = dynamic_cast<TH2D*> (fOutput->FindObject(name));
 292   if (!h) {
 293     AliError(TString::Format("Unable to load histogram from output: %s", name));
 294     return false;
 295   } else {
 296     return true;
 297   }
 298 }
 299
 300 //________________________________________________________________________
 301 void AliAnalysisTaskFirstPhysics::Terminate(Option_t *)
 302 {
 303   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 304   if(!fOutput) {
 305     AliError("Could not retrieve TList fOutput.");
 306     return;
 307   }
 308 }