PWGGA/EMCALJetTasks/AliAnalysisTaskSOH.cxx

   1 // $Id$\r
   2 //\r
   3 // Simulation EMCal task.\r
   4 //\r
   5 // Author: Saehanseul Oh\r
   6 \r
   7 #include "AliAnalysisTaskSOH.h"\r
   8 \r
   9 #include <TH1F.h>\r
  10 #include <TH2F.h>\r
  11 \r
  12 #include "AliAnalysisManager.h"\r
  13 #include "AliAnalysisTask.h"\r
  14 #include "AliEMCALRecoUtils.h"\r
  15 #include "AliEMCALTrack.h"\r
  16 #include "AliESDCaloCluster.h"\r
  17 #include "AliESDEvent.h"\r
  18 #include "AliESDInputHandler.h"\r
  19 #include "AliESDtrack.h"\r
  20 #include "AliESDtrackCuts.h"\r
  21 #include "AliExternalTrackParam.h"\r
  22 #include "AliLog.h"\r
  23 #include "AliMCEvent.h"\r
  24 \r
  25 ClassImp(AliAnalysisTaskSOH)\r
  26 \r
  27 //________________________________________________________________________\r
  28 AliAnalysisTaskSOH::AliAnalysisTaskSOH() :\r
  29   AliAnalysisTaskSE(), \r
  30   fESD(0), \r
  31   fMC(0), \r
  32   fEsdTrackCuts(0x0),\r
  33   fHybridTrackCuts1(0x0),\r
  34   fHybridTrackCuts2(0x0),\r
  35   fTrackIndices(0x0),\r
  36   fOutputList(0x0),        \r
  37   fHEventStat(0),      \r
  38   fHTrkEffParGenPt(0), \r
  39   fHTrkEffDetGenPt(0), \r
  40   fHTrkEffDetRecPt(0), \r
  41   fHEOverPVsPt(0x0),\r
  42   fHEMCalResponsePion(0x0), \r
  43   fHEMCalResponseElec(0x0),\r
  44   fHEMCalResponseProton(0x0), \r
  45   fHEMCalRecPhiEtaClus(0x0),\r
  46   fHEMCalRecPhiEtaTrk(0x0), \r
  47   fHClsPhiEta(0x0),\r
  48   fHEMCalRecdPhidEta(0x0),\r
  49   fHEMCalRecdPhidEtaP(0x0),\r
  50   fHEMCalRecdPhidEtaM(0x0),\r
  51   fHEMCalRecdPhidEta_Truth(0x0),\r
  52   fHEMCalRecdPhidEtaP_Truth(0x0),\r
  53   fHEMCalRecdPhidEtaM_Truth(0x0),\r
  54   fHEMCalRecdPhidEtaposEta(0x0),\r
  55   fHEMCalRecdPhidEtanegEta(0x0)\r
  56 {\r
  57   // Constructor\r
  58 \r
  59   // Output slot #1 writes into a TH1 container\r
  60   DefineOutput(1, TList::Class());\r
  61 }\r
  62 \r
  63 \r
  64 //________________________________________________________________________\r
  65 AliAnalysisTaskSOH::AliAnalysisTaskSOH(const char *name) :\r
  66   AliAnalysisTaskSE(name), \r
  67   fESD(0), \r
  68   fMC(0), \r
  69   fEsdTrackCuts(0x0),\r
  70   fHybridTrackCuts1(0x0),\r
  71   fHybridTrackCuts2(0x0),\r
  72   fTrackIndices(0x0),\r
  73   fOutputList(0x0),        \r
  74   fHEventStat(0),      \r
  75   fHTrkEffParGenPt(0), \r
  76   fHTrkEffDetGenPt(0), \r
  77   fHTrkEffDetRecPt(0), \r
  78   fHEOverPVsPt(0x0),\r
  79   fHEMCalResponsePion(0x0), \r
  80   fHEMCalResponseElec(0x0),\r
  81   fHEMCalResponseProton(0x0), \r
  82   fHEMCalRecPhiEtaClus(0x0),\r
  83   fHEMCalRecPhiEtaTrk(0x0), \r
  84   fHClsPhiEta(0x0),\r
  85   fHEMCalRecdPhidEta(0x0),\r
  86   fHEMCalRecdPhidEtaP(0x0),\r
  87   fHEMCalRecdPhidEtaM(0x0),\r
  88   fHEMCalRecdPhidEta_Truth(0x0),\r
  89   fHEMCalRecdPhidEtaP_Truth(0x0),\r
  90   fHEMCalRecdPhidEtaM_Truth(0x0),\r
  91   fHEMCalRecdPhidEtaposEta(0x0),\r
  92   fHEMCalRecdPhidEtanegEta(0x0)\r
  93 {\r
  94   // Constructor\r
  95 \r
  96   // Output slot #1 writes into a TH1 container\r
  97   DefineOutput(1, TList::Class());\r
  98 }\r
  99 \r
 100 \r
 101 //________________________________________________________________________\r
 102 AliAnalysisTaskSOH::~AliAnalysisTaskSOH()\r
 103 {\r
 104   // Destructor.\r
 105 \r
 106   delete fEsdTrackCuts;\r
 107   delete fHybridTrackCuts1;\r
 108   delete fHybridTrackCuts2;\r
 109   delete fTrackIndices;\r
 110 }\r
 111 \r
 112 \r
 113 //________________________________________________________________________\r
 114 void AliAnalysisTaskSOH::UserCreateOutputObjects()\r
 115 {\r
 116   // Create histograms, called once.\r
 117 \r
 118   OpenFile(1);\r
 119   \r
 120   fOutputList = new TList();\r
 121   fOutputList->SetOwner(1);\r
 122 \r
 123   fHEventStat = new TH1F("fHEventStat","Event statistics for analysis",1,0,1);\r
 124   fHEventStat->GetXaxis()->SetBinLabel(1,"Event");\r
 125   fOutputList->Add(fHEventStat);\r
 126 \r
 127   fHTrkEffParGenPt = new TH1F("fHTrkEffParGenPt", "Particle level truth p_{T} distribution of generated primary charged pions;p_{T}^{gen} (GeV/c)",15,0.15,3.1);\r
 128   fOutputList->Add(fHTrkEffParGenPt);\r
 129 \r
 130   fHTrkEffDetGenPt = new TH1F("fHTrkEffDetGenPt", "Detector level truth p_{T} distribution of primary charged pions;p_{T}^{gen} (GeV/c)",15,0.15,3.1);\r
 131   fOutputList->Add(fHTrkEffDetGenPt);\r
 132 \r
 133   fHTrkEffDetRecPt = new TH1F("fHTrkEffDetRecPt", "Reconstructed track p_{T} distribution of primary charged pions;p_{T}^{rec} (GeV/c)",15,0.15,3.1);\r
 134   fOutputList->Add(fHTrkEffDetRecPt);\r
 135 \r
 136   fHEOverPVsPt = new TH2F("fHEOverPVsPt", "E/P vs track p_{T}; p_{T} (GeV/c); E/P", 20 , 0, 4, 40, 0, 3.2);\r
 137   fOutputList->Add(fHEOverPVsPt);\r
 138 \r
 139   fHEMCalResponsePion = new TH2F("fHEMCalResponsePion", "Pion E/P vs track p_{T}; p_{T} (GeV/c); E/P", 100 , 0, 4, 100, 0, 3.2);\r
 140   fOutputList->Add(fHEMCalResponsePion);\r
 141 \r
 142   fHEMCalResponseElec = new TH2F("fHEMCalResponseElec", "Electron E/P vs track p_{T};  p_{T} (GeV/c); E/P", 100 , 0, 4, 100, 0, 3.2);\r
 143   fOutputList->Add(fHEMCalResponseElec);\r
 144 \r
 145   fHEMCalResponseProton = new TH2F("fHEMCalResponseProton", "Proton E/P vs track p_{T};  p_{T} (GeV/c); E/P", 100 , 0, 4, 100, 0, 3.2);\r
 146   fOutputList->Add(fHEMCalResponseElec);\r
 147 \r
 148   fHEMCalRecPhiEtaClus = new TH2F("fHEMCalRecPhiEtaClus","EMCAL Cluster#phi-#eta; #eta; #phi",1000,-3,3,1000,-4,4);\r
 149   fOutputList->Add(fHEMCalRecPhiEtaClus);\r
 150 \r
 151   fHEMCalRecPhiEtaTrk = new TH2F("fHEMCalRecPhiEtaTrk","EMCAL Track #phi-#eta; #eta; #phi",1000,-3,3,1000,-4,4);\r
 152   fOutputList->Add(fHEMCalRecPhiEtaTrk);\r
 153 \r
 154   fHClsPhiEta = new TH2F("fHClsPhiEta", "cluster #phi-#eta; #eta; #phi",1000,-3,3,1000,-4,4);\r
 155   fOutputList->Add(fHClsPhiEta);\r
 156 \r
 157   fHEMCalRecdPhidEta = new TH2F("fHEMCalRecdPhidEta","EMCAL Cluster-Track #Delta#phi-#Delta#eta; #Delta#eta; #Delta#phi",1000,-0.1,0.1,1000,-0.5,0.5);\r
 158   fOutputList->Add(fHEMCalRecdPhidEta);\r
 159 \r
 160   fHEMCalRecdPhidEtaP = new TH2F("fHEMCalRecdPhidEtaP","EMCAL Charge+ Cluster-Track #Delta#phi-#Delta#eta; #Delta#eta; #Delta#phi",1000,-0.1,0.1,1000,-0.5,0.5);\r
 161   fOutputList->Add(fHEMCalRecdPhidEtaP);\r
 162 \r
 163   fHEMCalRecdPhidEtaM = new TH2F("fHEMCalRecdPhidEtaM","EMCAL Charge- Cluster-Track #Delta#phi-#Delta#eta; #Delta#eta; #Delta#phi",1000,-0.1,0.1,1000,-0.5,0.5);\r
 164   fOutputList->Add(fHEMCalRecdPhidEtaM);\r
 165 \r
 166   fHEMCalRecdPhidEta_Truth = new TH2F("fHEMCalRecdPhidEta_Truth","EMCAL Cluster-Track(Truth matched) #Delta#phi-#Delta#eta; #Delta#eta; #Delta#phi",1000,-0.1,0.1,1000,-0.5,0.5);\r
 167   fOutputList->Add(fHEMCalRecdPhidEta_Truth);\r
 168 \r
 169   fHEMCalRecdPhidEtaP_Truth = new TH2F("fHEMCalRecdPhidEtaP_Truth","EMCAL Charge+ Cluster-Track(Truth matched) #Delta#phi-#Delta#eta; #Delta#eta; #Delta#phi",1000,-0.1,0.1,1000,-0.5,0.5);\r
 170   fOutputList->Add(fHEMCalRecdPhidEtaP_Truth);\r
 171 \r
 172   fHEMCalRecdPhidEtaM_Truth = new TH2F("fHEMCalRecdPhidEtaM_Truth","EMCAL Charge- Cluster-Track(Truth matched) #Delta#phi-#Delta#eta; #Delta#eta; #Delta#phi",1000,-0.1,0.1,1000,-0.5,0.5);\r
 173   fOutputList->Add(fHEMCalRecdPhidEtaM_Truth);\r
 174 \r
 175   fHEMCalRecdPhidEtaposEta = new TH2F("fHEMCalRecdPhidEtaposEta","(+eta track) EMCAL Cluster-Track #Delta#phi-#Delta#eta; #Delta#eta; #Delta#phi",1000,-0.1,0.1,1000,-0.5,0.5);\r
 176   fOutputList->Add(fHEMCalRecdPhidEtaposEta);\r
 177 \r
 178   fHEMCalRecdPhidEtanegEta = new TH2F("fHEMCalRecdPhidEtanegEta","(-eta track) EMCAL Cluster-Track #Delta#phi-#Delta#eta; #Delta#eta; #Delta#phi",1000,-0.1,0.1,1000,-0.5,0.5);\r
 179   fOutputList->Add(fHEMCalRecdPhidEtanegEta);\r
 180   fTrackIndices = new TArrayI();\r
 181 \r
 182   PostData(1, fOutputList);\r
 183 }\r
 184 \r
 185 //________________________________________________________________________\r
 186 void AliAnalysisTaskSOH::UserExec(Option_t *) \r
 187 {\r
 188   // Main loop, called for each event.\r
 189 \r
 190   // Post output data.\r
 191   fESD = dynamic_cast<AliESDEvent*>(InputEvent());\r
 192   if (!fESD) {\r
 193     printf("ERROR: fESD not available\n");\r
 194     return;\r
 195   }\r
 196 \r
 197   fMC = MCEvent();\r
 198   if (!fMC) {\r
 199     printf("ERROR: fMC not available\n");\r
 200     return;\r
 201   }\r
 202 \r
 203   fHEventStat->Fill(0.5);\r
 204 \r
 205   if(fTrackIndices) \r
 206     fTrackIndices->Reset();\r
 207 \r
 208   ProcessTrack();\r
 209   ProcessMc();\r
 210   \r
 211   PostData(1, fOutputList);\r
 212 }   \r
 213    \r
 214 //________________________________________________________________________\r
 215 void  AliAnalysisTaskSOH::ProcessTrack()\r
 216 {\r
 217   // Process track.\r
 218 \r
 219   fTrackIndices->Set(fESD->GetNumberOfTracks());\r
 220   AliDebug(3,Form("%s:%d Selecting tracks",(char*)__FILE__,__LINE__));\r
 221 \r
 222   Int_t isMth = 0;\r
 223   Int_t nTracks = 0;\r
 224 \r
 225   Float_t ClsPos[3] = {-999,-999,-999};\r
 226   Double_t emcTrkpos[3] = {-999,-999,-999};\r
 227 \r
 228   for(Int_t itr=0; itr<fESD->GetNumberOfCaloClusters(); itr++)\r
 229   {\r
 230     AliESDCaloCluster *cluster = fESD->GetCaloCluster(itr);\r
 231     if(!cluster) continue;\r
 232     cluster->GetPosition(ClsPos);\r
 233     TVector3 VClsPos(ClsPos[0], ClsPos[1], ClsPos[2]);\r
 234     fHClsPhiEta->Fill(VClsPos.Eta(), VClsPos.Phi());\r
 235   }\r
 236 \r
 237   for(Int_t itr=0; itr<fESD->GetNumberOfTracks(); itr++)\r
 238   {\r
 239     AliESDtrack *esdtrack = fESD->GetTrack(itr);\r
 240     if(!esdtrack)continue;\r
 241     AliESDtrack *newTrack = GetAcceptTrack(esdtrack);\r
 242     if(!newTrack) continue;\r
 243   \r
 244     Double_t clsE = -1;\r
 245     Int_t clsIndex = newTrack->GetEMCALcluster();\r
 246     if(newTrack->GetEMCALcluster()>-1)\r
 247     {\r
 248       AliESDCaloCluster *cluster = fESD->GetCaloCluster(clsIndex);\r
 249       if(IsGoodCluster(cluster))\r
 250       {\r
 251         isMth=1;\r
 252               \r
 253         cluster->GetPosition(ClsPos);\r
 254         TVector3 VClsPos(ClsPos[0], ClsPos[1], ClsPos[2]);\r
 255               \r
 256         AliEMCALTrack EMCTrk(*newTrack);\r
 257         if(!EMCTrk.PropagateToGlobal(ClsPos[0], ClsPos[1], ClsPos[2], 0.0, 0.0)) {\r
 258           continue;\r
 259         }\r
 260         EMCTrk.GetXYZ(emcTrkpos);\r
 261         TVector3 VemcTrkPos(emcTrkpos[0],emcTrkpos[1],emcTrkpos[2]);\r
 262               \r
 263         fHEMCalRecPhiEtaClus->Fill(VClsPos.Eta(), VClsPos.Phi());\r
 264         fHEMCalRecPhiEtaTrk->Fill(VemcTrkPos.Eta(), VemcTrkPos.Phi());\r
 265               \r
 266         Double_t dPhi = VClsPos.Phi() - VemcTrkPos.Phi();\r
 267         if (dPhi < -1*TMath::Pi()) dPhi += (2*TMath::Pi());\r
 268         else if (dPhi > TMath::Pi()) dPhi -= (2*TMath::Pi());\r
 269               \r
 270         Double_t dEta = VClsPos.Eta() - VemcTrkPos.Eta();\r
 271               \r
 272         fHEMCalRecdPhidEta->Fill(dEta, dPhi);   \r
 273 \r
 274         if((newTrack->GetLabel())>-1 && (newTrack->GetLabel()) < fMC->GetNumberOfTracks())\r
 275         {\r
 276           AliVParticle *vParticle = fMC->GetTrack(newTrack->GetLabel());\r
 277           if(IsGoodMcParticle(vParticle, newTrack->GetLabel()))\r
 278           {\r
 279             fHEMCalRecdPhidEta_Truth->Fill(dEta, dPhi);\r
 280             if(vParticle->Charge() > 0) fHEMCalRecdPhidEtaP_Truth->Fill(dEta, dPhi);\r
 281             if(vParticle->Charge() < 0) fHEMCalRecdPhidEtaM_Truth->Fill(dEta, dPhi);\r
 282           }\r
 283         }\r
 284 \r
 285         if(esdtrack->Charge() > 0) {fHEMCalRecdPhidEtaP->Fill(dEta, dPhi);}\r
 286         if(esdtrack->Charge() < 0) {fHEMCalRecdPhidEtaM->Fill(dEta, dPhi);}\r
 287 \r
 288         if(VemcTrkPos.Eta() > 0) fHEMCalRecdPhidEtaposEta->Fill(dEta, dPhi);\r
 289         if(VemcTrkPos.Eta() < 0) fHEMCalRecdPhidEtanegEta->Fill(dEta, dPhi);\r
 290               \r
 291         clsE = cluster->E();\r
 292         if(newTrack->P()>0) fHEOverPVsPt->Fill(newTrack->Pt(),clsE/newTrack->P());\r
 293       }\r
 294           \r
 295       Int_t ipart = newTrack->GetLabel();\r
 296       if(ipart>-1 && ipart<fMC->GetNumberOfTracks())\r
 297       {\r
 298         AliVParticle* vParticle = fMC->GetTrack(ipart);\r
 299         if(isMth && vParticle)\r
 300         {\r
 301           if(TMath::Abs(vParticle->PdgCode())==211)\r
 302           {\r
 303             fHEMCalResponsePion->Fill(newTrack->Pt(),clsE/newTrack->P());\r
 304           }\r
 305           if(TMath::Abs(vParticle->PdgCode())==11)\r
 306           {\r
 307             fHEMCalResponseElec->Fill(newTrack->Pt(),clsE/newTrack->P());\r
 308           }\r
 309           if(TMath::Abs(vParticle->PdgCode())==2212)\r
 310           {\r
 311             fHEMCalResponseProton->Fill(newTrack->Pt(),clsE/newTrack->P());\r
 312           }\r
 313         }\r
 314       }\r
 315     }\r
 316     fTrackIndices->AddAt(itr,nTracks);\r
 317     nTracks++;\r
 318   }\r
 319   fTrackIndices->Set(nTracks);\r
 320 }\r
 321 \r
 322 //________________________________________________________________________\r
 323 void AliAnalysisTaskSOH::ProcessMc()\r
 324 {\r
 325   // Process MC.\r
 326 \r
 327   for(Int_t ipart=0; ipart<fMC->GetNumberOfTracks(); ipart++)\r
 328   {\r
 329     AliVParticle* vParticle = fMC->GetTrack(ipart);\r
 330     if(!IsGoodMcParticle(vParticle, ipart)) continue;\r
 331     Int_t pdgCode = vParticle->PdgCode();\r
 332       \r
 333     //tracking effciency\r
 334     if(TMath::Abs(vParticle->Eta())<0.9 && TMath::Abs(pdgCode==211))\r
 335     {\r
 336       fHTrkEffParGenPt->Fill(vParticle->Pt());\r
 337       for(Int_t j=0; j<fTrackIndices->GetSize(); j++)\r
 338       {\r
 339         AliESDtrack *esdtrack = fESD->GetTrack(fTrackIndices->At(j));\r
 340         if(esdtrack && esdtrack->GetLabel()==ipart)\r
 341         {\r
 342           fHTrkEffDetGenPt->Fill(vParticle->Pt());\r
 343           fHTrkEffDetRecPt->Fill(esdtrack->Pt());\r
 344           break;\r
 345         }\r
 346       }\r
 347     }\r
 348   } \r
 349 }\r
 350 \r
 351 //________________________________________________________________________\r
 352 AliESDtrack *AliAnalysisTaskSOH::GetAcceptTrack(AliESDtrack *esdtrack)\r
 353 {\r
 354   // Get accepted track.\r
 355 \r
 356   static AliESDtrack newTrack;\r
 357   if(fEsdTrackCuts->AcceptTrack(esdtrack))\r
 358   {\r
 359     esdtrack->Copy(newTrack);\r
 360     newTrack.SetTRDQuality(0);\r
 361   }\r
 362   else if(fHybridTrackCuts1->AcceptTrack(esdtrack))\r
 363   {\r
 364     if(esdtrack->GetConstrainedParam())\r
 365     {\r
 366       esdtrack->Copy(newTrack);\r
 367       const AliExternalTrackParam* constrainParam = esdtrack->GetConstrainedParam();\r
 368       newTrack.Set(constrainParam->GetX(),constrainParam->GetAlpha(),constrainParam->GetParameter(),constrainParam->GetCovariance());\r
 369       newTrack.SetTRDQuality(1);                \r
 370     }\r
 371     else \r
 372       return 0x0;\r
 373   }\r
 374   else if(fHybridTrackCuts2->AcceptTrack(esdtrack))\r
 375   {\r
 376     if(esdtrack->GetConstrainedParam())\r
 377     {\r
 378       esdtrack->Copy(newTrack);\r
 379       const AliExternalTrackParam* constrainParam = esdtrack->GetConstrainedParam();\r
 380       newTrack.Set(constrainParam->GetX(),constrainParam->GetAlpha(),constrainParam->GetParameter(),constrainParam->GetCovariance());\r
 381       newTrack.SetTRDQuality(2);                \r
 382     }\r
 383     else \r
 384       return 0x0;\r
 385   }\r
 386   else\r
 387   {\r
 388     return 0x0;\r
 389   }\r
 390 \r
 391   return &newTrack;\r
 392 }\r
 393 \r
 394 //________________________________________________________________________\r
 395 Bool_t AliAnalysisTaskSOH::IsGoodMcParticle(AliVParticle* vParticle, Int_t ipart)\r
 396 {\r
 397   // Return true if good MC particle.\r
 398 \r
 399   if(!vParticle) return kFALSE;\r
 400   if(!fMC->IsPhysicalPrimary(ipart)) return kFALSE;\r
 401   if (TMath::Abs(vParticle->Eta())>2) return kFALSE;\r
 402   if(vParticle->Pt()<0.15) return kFALSE;\r
 403   return kTRUE;\r
 404 }\r
 405 \r
 406 \r
 407 //________________________________________________________________________\r
 408 Bool_t AliAnalysisTaskSOH::IsGoodCluster(AliESDCaloCluster *cluster)\r
 409 {\r
 410   // Return true if good cluster.\r
 411 \r
 412   if(!cluster) return kFALSE;\r
 413   if (!cluster->IsEMCAL()) return kFALSE;\r
 414   return kTRUE;\r
 415 }\r
 416 \r
 417 \r
 418 //________________________________________________________________________\r
 419 void AliAnalysisTaskSOH::Terminate(Option_t *) \r
 420 {\r
 421   // Terminate analysis.\r
 422 }\r