PWGHF/hfe/AliAnalysisTaskHFEemcQA.cxx

   1 #include "TChain.h"
   2 #include "TTree.h"
   3 #include "TH1F.h"
   4 #include "TCanvas.h"
   5
   6 #include "AliAnalysisTask.h"
   7 #include "AliAnalysisManager.h"
   8
   9 #include "AliESDEvent.h"
  10 #include "AliESDInputHandler.h"
  11 #include "AliESDtrackCuts.h"
  12
  13 #include "AliAODEvent.h"
  14 #include "AliAODHandler.h"
  15
  16 #include "AliAnalysisTaskHFEemcQA.h"
  17
  18 //QA task for EMCAL electron analysis
  19
  20 ClassImp(AliAnalysisTaskHFEemcQA)
  21   //________________________________________________________________________
  22   AliAnalysisTaskHFEemcQA::AliAnalysisTaskHFEemcQA(const char *name)
  23 : AliAnalysisTaskSE(name),
  24   fVevent(0),
  25   fESD(0),
  26   fAOD(0),
  27   fOutputList(0),
  28   fVtxZ(0),
  29   fVtxX(0),
  30   fVtxY(0),
  31   fTrigMulti(0),
  32   fHistClustE(0),
  33   fEMCClsEtaPhi(0),
  34   fNegTrkIDPt(0),
  35   fTrkPt(0),
  36   fTrketa(0),
  37   fTrkphi(0),
  38   fdEdx(0),
  39   fTPCNpts(0),
  40   fHistPtMatch(0),
  41   fEMCTrkMatch(0),
  42   fEMCTrkPt(0),
  43   fEMCTrketa(0),
  44   fEMCTrkphi(0),
  45   fEMCdEdx(0),
  46   fEMCTPCNpts(0),
  47   fHistdEdxEop(0),
  48   fHistEop(0),
  49   fEleCanTPCNpts(0),
  50   fEleCanTPCNCls(0),
  51   fEleCanITSNCls(0),
  52   fEleCanITShit(0),
  53   fEleCanSPD1(0),
  54   fEleCanSPD2(0),
  55   fEleCanSPDBoth(0),
  56   fEleCanSPDOr(0)
  57 {
  58   // Constructor
  59
  60   // Define input and output slots here
  61   // Input slot #0 works with a TChain
  62   DefineInput(0, TChain::Class());
  63   // Output slot #0 id reserved by the base class for AOD
  64   // Output slot #1 writes into a TH1 container
  65   DefineOutput(1, TList::Class());
  66 }
  67
  68 //________________________________________________________________________
  69 void AliAnalysisTaskHFEemcQA::UserCreateOutputObjects()
  70 {
  71   // Create histograms
  72   // Called once
  73
  74   fOutputList = new TList();
  75   fOutputList->SetOwner();
  76
  77   fVtxZ = new TH1F("fVtxZ","Z vertex position;Vtx_{z};counts",1000,-50,50);
  78   fOutputList->Add(fVtxZ);
  79
  80   fVtxY = new TH1F("fVtxY","Y vertex position;Vtx_{y};counts",1000,-50,50);
  81   fOutputList->Add(fVtxY);
  82
  83   fVtxX = new TH1F("fVtxX","X vertex position;Vtx_{x};counts",1000,-50,50);
  84   fOutputList->Add(fVtxX);
  85
  86   fTrigMulti = new TH2F("fTrigMulti","Multiplicity distribution for different triggers; Trigger type; multiplicity",8,-1,7,2000,0,2000);
  87   fOutputList->Add(fTrigMulti);
  88
  89   fHistClustE = new TH1F("fHistClustE", "EMCAL cluster energy distribution; Cluster E;counts", 500, 0.0, 50.0);
  90   fOutputList->Add(fHistClustE);
  91
  92   fEMCClsEtaPhi = new TH2F("fEMCClsEtaPhi","EMCAL cluster #eta and #phi distribution;#eta;#phi",100,-0.9,0.9,200,0,6.3);
  93   fOutputList->Add(fEMCClsEtaPhi);
  94
  95   fNegTrkIDPt = new TH1F("fNegTrkIDPt", "p_{T} distribution of tracks with negative track id;p_{T} (GeV/c);counts", 500, 0.0, 50.0);
  96   fOutputList->Add(fNegTrkIDPt);
  97
  98   fTrkPt = new TH1F("fTrkPt","p_{T} distribution of all tracks;p_{T} (GeV/c);counts",1000,0,100);
  99   fOutputList->Add(fTrkPt);
 100
 101   fTrketa = new TH1F("fTrketa","All Track #eta distribution;#eta;counts",100,-1.5,1.5);
 102   fOutputList->Add(fTrketa);
 103
 104   fTrkphi = new TH1F("fTrkphi","All Track #phi distribution;#phi;counts",100,0,6.3);
 105   fOutputList->Add(fTrkphi);
 106
 107   fdEdx = new TH2F("fdEdx","All Track dE/dx distribution;p (GeV/c);dE/dx",200,0,20,500,0,100);
 108   fOutputList->Add(fdEdx);
 109
 110   fTPCNpts = new TH2F("fTPCNpts","All track TPC Npoints used for dE/dx calculation;p (GeV/c);N points",200,0,20,200,0.,200.);
 111   fOutputList->Add(fTPCNpts);
 112
 113   fHistPtMatch = new TH1F("fHistPtMatch", "p_{T} distribution of tracks matched to EMCAL;p_{T} (GeV/c);counts",1000, 0.0, 100.0);
 114   fOutputList->Add(fHistPtMatch);
 115
 116   fEMCTrkMatch = new TH2F("fEMCTrkMatch","Distance of EMCAL cluster to its closest track;#phi;z",100,-0.3,0.3,100,-0.3,0.3);
 117   fOutputList->Add(fEMCTrkMatch);
 118
 119   fEMCTrkPt = new TH1F("fEMCTrkPt","p_{T} distribution of tracks with EMCAL cluster;p_{T} (GeV/c);counts",1000,0,100);
 120   fOutputList->Add(fEMCTrkPt);
 121
 122   fEMCTrketa = new TH1F("fEMCTrketa","#eta distribution of tracks matched to EMCAL;#eta;counts",100,-1.5,1.5);
 123   fOutputList->Add(fEMCTrketa);
 124
 125   fEMCTrkphi = new TH1F("fEMCTrkphi","#phi distribution of tracks matched to EMCAL;#phi;counts",100,0,6.3);
 126   fOutputList->Add(fEMCTrkphi);
 127
 128   fEMCdEdx = new TH2F("fEMCdEdx","dE/dx distribution of tracks matched to EMCAL;p (GeV/c);dE/dx",200,0,20,500,0,100);
 129   fOutputList->Add(fEMCdEdx);
 130
 131   fEMCTPCNpts = new TH2F("fEMCTPCNpts","TPC Npoints used for dE/dx for tracks matched to EMCAL;p (GeV/c);N points",200,0,20,200,0.,200.);
 132   fOutputList->Add(fEMCTPCNpts);
 133
 134   fHistEop = new TH2F("fHistEop", "E/p distribution;p_{T} (GeV/c);E/p", 200,0,20,60, 0.0, 3.0);
 135   fOutputList->Add(fHistEop);
 136
 137   fHistdEdxEop = new TH2F("fHistdEdxEop", "E/p vs. dE/dx;E/p;dE/dx", 60, 0.0, 3.0, 500,0,100);
 138   fOutputList->Add(fHistdEdxEop);
 139
 140   fEleCanTPCNpts = new TH2F("fEleCanTPCNpts","TPC Npoints used for dE/dx for electron candidates;p_{T} (GeV/c);N points",200,0,20,200,0,200);
 141   fOutputList->Add(fEleCanTPCNpts);
 142
 143   fEleCanTPCNCls = new TH2F("fEleCanTPCNCls","TPC N clusters for electron candidates;p_{T} (GeV/c);N TPC clusters",200,0,20,171,-0.5,170.5);
 144   fOutputList->Add(fEleCanTPCNCls);
 145
 146   fEleCanITSNCls = new TH2F("fEleCanITSNCls","ITS N clusters for electron candidates;p_{T} (GeV/c);N ITS clusters",200,0,20,8,-0.5,7.5);
 147   fOutputList->Add(fEleCanITSNCls);
 148
 149   fEleCanITShit = new TH1F("fEleCanITShit","ITS hit map;ITS layer;counts",7,-0.5,6.5);
 150   fOutputList->Add(fEleCanITShit);
 151
 152   fEleCanSPD1 = new TH2F("fEleCanSPD1","Hit on SPD layer 1;p_{T} (GeV/c);Hit",200,0,20,1,0,1);
 153   fOutputList->Add(fEleCanSPD1);
 154
 155   fEleCanSPD2 = new TH2F("fEleCanSPD2","Hit on SPD layer 2;p_{T} (GeV/c);Hit",200,0,20,1,0,1);
 156   fOutputList->Add(fEleCanSPD2);
 157
 158   fEleCanSPDBoth = new TH2F("fEleCanSPDBoth","Tracks with hits on both SPD layer;p_{T} (GeV/c);Hit",200,0,20,1,0,1);
 159   fOutputList->Add(fEleCanSPDBoth);
 160
 161   fEleCanSPDOr = new TH2F("fEleCanSPDOr","Tracks with hits on both SPD layer;p_{T} (GeV/c);Hit",200,0,20,1,0,1);
 162   fOutputList->Add(fEleCanSPDOr);
 163
 164   PostData(1,fOutputList);
 165 }
 166
 167 //________________________________________________________________________
 168 void AliAnalysisTaskHFEemcQA::UserExec(Option_t *)
 169 {
 170   // Main loop
 171   // Called for each event
 172   // Post output data.
 173
 174   UInt_t evSelMask=((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 175
 176   fVevent = dynamic_cast<AliVEvent*>(InputEvent());
 177   if (!fVevent) {
 178     printf("ERROR: fVEvent not available\n");
 179     return;
 180   }
 181
 182   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 183   if (fESD) {
 184     printf("fESD available\n");
 185     //return;
 186   }
 187
 188   AliESDtrackCuts* esdTrackCutsH = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kFALSE);
 189   esdTrackCutsH->SetMaxDCAToVertexXY(2.4);
 190   esdTrackCutsH->SetMaxDCAToVertexZ(3.2);
 191   esdTrackCutsH->SetDCAToVertex2D(kTRUE);
 192
 193   fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 194   if (fAOD) {
 195     printf("fAOD available\n");
 196     //return;
 197   }
 198
 199   int ntracks;
 200   ntracks = fVevent->GetNumberOfTracks();
 201   printf("There are %d tracks in this event\n",ntracks);
 202   //  if(fESD)printf("There are %d tracks in this event\n", fESD->GetNumberOfTracks());
 203   //  if(fAOD)printf("There are %d tracks in this event\n", fAOD->GetNumberOfTracks());
 204
 205   double Zvertex = -100, Xvertex = -100, Yvertex = -100;
 206   const AliVVertex *pVtx = fVevent->GetPrimaryVertex();
 207   Zvertex = pVtx->GetZ();
 208   Yvertex = pVtx->GetY();
 209   Xvertex = pVtx->GetX();
 210   fVtxZ->Fill(Zvertex);
 211   fVtxX->Fill(Xvertex);
 212   fVtxY->Fill(Yvertex);
 213   /*
 214      if(fESD)
 215      {
 216      const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
 217      Zvertex = pVtx->GetZ();
 218      }
 219      if(fAOD)
 220      {
 221      const AliAODVertex *pVtx = fAOD->GetPrimaryVertex();
 222      Zvertex = pVtx->GetZ();
 223      }
 224    */
 225   // trigger check
 226
 227   TString firedTrigger;
 228   TString TriggerEG1("EG1");
 229   TString TriggerEG2("EG2");
 230   fVevent->GetFiredTriggerClasses();
 231
 232   //if(fAOD) firedTrigger = fAOD->GetFiredTriggerClasses();
 233   //else if(fESD) firedTrigger = fESD->GetFiredTriggerClasses();
 234
 235   Bool_t EG1tr = kFALSE;
 236   Bool_t EG2tr = kFALSE;
 237
 238   if(firedTrigger.Contains(TriggerEG1))EG1tr = kTRUE;
 239   if(firedTrigger.Contains(TriggerEG2))EG2tr = kTRUE;
 240
 241
 242   if (fAOD){
 243     Double_t multiplicity=fAOD->GetHeader()->GetRefMultiplicity();
 244     fTrigMulti->Fill(-0.5, multiplicity);
 245     if(evSelMask & AliVEvent::kAny) fTrigMulti->Fill(0.5, multiplicity);
 246     if(evSelMask & AliVEvent::kMB) fTrigMulti->Fill(1.5, multiplicity);
 247     if(evSelMask & AliVEvent::kINT7) fTrigMulti->Fill(1.5, multiplicity);
 248     if(evSelMask & AliVEvent::kINT8) fTrigMulti->Fill(1.5, multiplicity);
 249     if(evSelMask & AliVEvent::kEMC1) fTrigMulti->Fill(2.5, multiplicity);
 250     if(evSelMask & AliVEvent::kEMC8) fTrigMulti->Fill(3.5, multiplicity);
 251     if(evSelMask & AliVEvent::kEMCEJE) fTrigMulti->Fill(4.5, multiplicity);
 252     if(evSelMask & AliVEvent::kEMCEGA) fTrigMulti->Fill(5.5, multiplicity);
 253     if(evSelMask & AliVEvent::kEMCEGA & EG2tr) fTrigMulti->Fill(5.5, multiplicity);
 254     if(evSelMask & AliVEvent::kEMCEGA & EG2tr) fTrigMulti->Fill(5.5, multiplicity);
 255   }
 256
 257   // event selection
 258   if(fabs(Zvertex>10.0))return;
 259
 260   //EMCAL cluster information
 261   int Nclust = 0;
 262   Nclust = fVevent->GetNumberOfCaloClusters();
 263   //  if(fESD)Nclust = fESD->GetNumberOfCaloClusters();
 264   //  if(fAOD)Nclust = fAOD->GetNumberOfCaloClusters();
 265   for(Int_t icl=0; icl<Nclust; icl++)
 266   {
 267     AliVCluster *clust = 0x0;
 268     clust = fVevent->GetCaloCluster(icl);
 269     //    if(fESD)clust = (AliVCluster*) fESD->GetCaloCluster(icl);
 270     //    if(fAOD)clust = (AliVCluster*) fAOD->GetCaloCluster(icl);
 271     if(clust && clust->IsEMCAL())
 272     {
 273       double clustE = clust->E();
 274       float  emcx[3]; // cluster pos
 275       clust->GetPosition(emcx);
 276       TVector3 clustpos(emcx[0],emcx[1],emcx[2]);
 277       double emcphi = clustpos.Phi();
 278       double emceta = clustpos.Eta();
 279       fHistClustE->Fill(clustE);
 280       fEMCClsEtaPhi->Fill(emceta,emcphi);
 281     }
 282   }
 283
 284   // Track loop
 285   for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
 286
 287     //---------combine both esd and aod tracks -------
 288     AliVParticle* Vtrack = fVevent->GetTrack(iTracks);
 289     if (!Vtrack) {
 290       printf("ERROR: Could not receive track %d\n", iTracks);
 291       continue;
 292     }
 293     AliVTrack *track = dynamic_cast<AliVTrack*>(Vtrack);
 294     AliESDtrack *etrack = dynamic_cast<AliESDtrack*>(Vtrack);
 295     AliAODTrack *atrack = dynamic_cast<AliAODTrack*>(Vtrack);
 296
 297     if(fAOD)
 298       if(!atrack->TestFilterMask(AliAODTrack::kTrkGlobalNoDCA)) continue;
 299
 300     if(fESD)
 301       if(!esdTrackCutsH->AcceptTrack(etrack))continue;
 302
 303     //track properties
 304     double dEdx = track->GetTPCsignal();
 305     if(track->GetID()<0) fNegTrkIDPt->Fill(track->Pt());
 306     fTrkPt->Fill(track->Pt());
 307     fTrketa->Fill(track->Eta());
 308     fTrkphi->Fill(track->Phi());
 309     fdEdx->Fill(track->P(),dEdx);
 310     fTPCNpts->Fill(track->P(),track->GetTPCsignalN());
 311
 312     //track matching to EMCAL
 313     int EMCalIndex = -1;
 314     EMCalIndex = track->GetEMCALcluster();
 315     if(EMCalIndex < 0) continue;
 316     fHistPtMatch->Fill(track->Pt());
 317
 318     AliVCluster *clustMatch = (AliVCluster*)fAOD->GetCaloCluster(EMCalIndex);
 319     if(clustMatch && clustMatch->IsEMCAL()) //isnt clustMatch always EMCal cluster?
 320     {
 321       fEMCTrkMatch->Fill(clustMatch->GetTrackDx(),clustMatch->GetTrackDz());
 322       //properties of tracks matched to the EMCAL
 323       fEMCTrkPt->Fill(track->Pt());
 324       fEMCTrketa->Fill(track->Eta());
 325       fEMCTrkphi->Fill(track->Phi());
 326       fEMCdEdx->Fill(track->P(),dEdx);
 327       fEMCTPCNpts->Fill(track->P(),track->GetTPCsignalN());
 328
 329       //E/p distribution
 330       double clustMatchE = clustMatch->E();
 331       double eop = -1.0;
 332       if(track->P()>0)eop = clustMatchE/track->P();
 333
 334       if(track->Pt()>2.0)fHistdEdxEop->Fill(eop,dEdx);
 335       fHistEop->Fill(track->Pt(),eop);
 336
 337       //track properties of EMCAL electron cadidates
 338       if(eop>0.8 and eop<1.2){
 339         fEleCanTPCNpts->Fill(track->Pt(),track->GetTPCsignalN());
 340         fEleCanTPCNCls->Fill(track->Pt(),track->GetTPCNcls());
 341
 342         Int_t fITSncls=0;
 343         for(Int_t l=0;l<6;l++) {
 344           if(TESTBIT(track->GetITSClusterMap(),l)) {
 345             fEleCanITShit->Fill(l);
 346             if(l==0) fEleCanSPD1->Fill(track->Pt(),0.5);
 347             if(l==1) fEleCanSPD2->Fill(track->Pt(),0.5);
 348             if(l==0 && l==1) fEleCanSPDBoth->Fill(track->Pt(),0.5);
 349             if(l==0 || l==1) fEleCanSPDOr->Fill(track->Pt(),0.5);
 350             fITSncls++;
 351           }
 352         }
 353         fEleCanITSNCls->Fill(track->Pt(),fITSncls++);
 354       }
 355     }
 356   } //track loop
 357
 358   PostData(1, fOutputList);
 359 }
 360 //________________________________________________________________________
 361 void AliAnalysisTaskHFEemcQA::Terminate(Option_t *)
 362 {
 363   // Draw result to the screen
 364   // Called once at the end of the query
 365
 366   fOutputList = dynamic_cast<TList*> (GetOutputData(1));
 367   if (!fOutputList) {
 368     printf("ERROR: Output list not available\n");
 369     return;
 370   }
 371
 372 }