PWGGA/EMCALJetTasks/AliHadCorrTask.cxx

   1 // $Id$
   2
   3 #include <TClonesArray.h>
   4 #include <TParticle.h>
   5 #include <TList.h>
   6 #include <TH1F.h>
   7 #include <TH2F.h>
   8 #include <TChain.h>
   9 #include <TLorentzVector.h>
  10
  11 #include "AliCentrality.h"
  12 #include "AliEmcalJet.h"
  13 #include "AliAnalysisManager.h"
  14 #include "AliESDtrack.h"
  15 #include "AliFJWrapper.h"
  16 #include "AliESDCaloCluster.h"
  17 #include "AliAODCaloCluster.h"
  18 #include "AliPicoTrack.h"
  19 #include "AliVEventHandler.h"
  20
  21 #include "AliHadCorrTask.h"
  22
  23 ClassImp(AliHadCorrTask)
  24
  25 //________________________________________________________________________
  26 AliHadCorrTask::AliHadCorrTask(const char *name) :
  27   AliAnalysisTaskSE("AliHadCorrTask"),
  28   fTracksName("Tracks"),
  29   fCaloName("CaloClusters"),
  30   fHadCorr(0),
  31   fJets(0),
  32   fMinPt(0.15),
  33   fOutCaloName("CaloClustersOut"),
  34   fOutClusters(0),
  35   fOutputList(0),
  36   fHistNclusvsCent(0),
  37   fHistNclusMatchvsCent(0),
  38   fHistEbefore(0),
  39   fHistEafter(0)
  40 {
  41
  42   for(Int_t i=0; i<4; i++){
  43     fHistMatchEtaPhi[i]=0x0;
  44     fHistMatchEvsP[i]=0x0;
  45     fHistMatchdRvsEP[i]=0x0;
  46   }
  47
  48   // Standard constructor.
  49   cout << "Constructor for HadCorrTask " << name <<endl;
  50   if (!name)
  51     return;
  52
  53   SetName(name);
  54   fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
  55
  56   DefineInput(0,TChain::Class());
  57   DefineOutput(1,TList::Class());
  58 }
  59
  60 AliHadCorrTask::AliHadCorrTask() :
  61   AliAnalysisTaskSE("AliHadCorrTask"),
  62   fTracksName("Tracks"),
  63   fCaloName("CaloClusters"),
  64   fHadCorr(0),
  65   fJets(0),
  66   fMinPt(0.15),
  67   fOutCaloName("CaloClustersOut"),
  68   fOutClusters(0),
  69   fOutputList(0),
  70   fHistNclusvsCent(0),
  71   fHistNclusMatchvsCent(0),
  72   fHistEbefore(0),
  73   fHistEafter(0)
  74 {
  75   // Standard constructor.
  76
  77   fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
  78 }
  79
  80
  81 //________________________________________________________________________
  82 AliHadCorrTask::~AliHadCorrTask()
  83 {
  84   // Destructor
  85 }
  86
  87 //________________________________________________________________________
  88 void AliHadCorrTask::UserCreateOutputObjects()
  89 {
  90
  91   AliVEventHandler* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
  92
  93   if (!handler) {
  94     AliError("Input handler not available!");
  95     return;
  96   }
  97
  98   if (handler->InheritsFrom("AliESDInputHandler")) {
  99     fOutClusters = new TClonesArray("AliESDCaloCluster");
 100   }
 101   else if (handler->InheritsFrom("AliAODInputHandler")) {
 102     fOutClusters = new TClonesArray("AliAODCaloCluster");
 103   }
 104   else {
 105     AliError("Input handler not recognized!");
 106     return;
 107   }
 108
 109   fOutClusters->SetName(fOutCaloName);
 110
 111   OpenFile(1);
 112   fOutputList = new TList();
 113   fOutputList->SetOwner();
 114
 115   char name[200];
 116
 117   for(int icent=0; icent<4; icent++){
 118     sprintf(name,"fHistMatchEtaPhi_%i",icent);
 119     fHistMatchEtaPhi[icent] = new TH2F(name,name,400,-0.2,0.2,1600,-0.8,0.8);
 120     sprintf(name,"fHistMatchEvsP_%i",icent);
 121     fHistMatchEvsP[icent]=new TH2F(name,name,400,0.,200.,1000,0.,10.);
 122     sprintf(name,"fHistMatchdRvsEP_%i",icent);
 123     fHistMatchdRvsEP[icent]=new TH2F(name,name,1000,0.,1.,1000,0.,10.);
 124
 125
 126   fOutputList->Add(fHistMatchEtaPhi[icent]);
 127   fOutputList->Add(fHistMatchEvsP[icent]);
 128   fOutputList->Add(fHistMatchdRvsEP[icent]);
 129
 130   }
 131   fHistNclusvsCent=new TH1F("Nclusvscent","NclusVsCent",100,0,100);
 132   fHistNclusMatchvsCent=new TH1F("NclusMatchvscent","NclusMatchVsCent",100,0,100);
 133   fHistEbefore=new TH1F("Ebefore","Ebefore",100,0,100);
 134   fHistEafter=new TH1F("Eafter","Eafter",100,0,100);
 135
 136   fOutputList->Add(fHistNclusMatchvsCent);
 137   fOutputList->Add(fHistNclusvsCent);
 138   fOutputList->Add(fHistEbefore);
 139   fOutputList->Add(fHistEafter);
 140
 141   PostData(1, fOutputList);
 142 }
 143
 144 //________________________________________________________________________
 145 void AliHadCorrTask::UserExec(Option_t *)
 146 {
 147
 148   if (!(InputEvent()->FindListObject(fOutCaloName)))
 149     InputEvent()->AddObject(fOutClusters);
 150   else fOutClusters->Delete();
 151
 152
 153   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
 154   TClonesArray *tracks = 0;
 155   TClonesArray *clus   = 0;
 156   TList *l = InputEvent()->GetList();
 157
 158   AliCentrality *centrality = 0;
 159   centrality = dynamic_cast<AliCentrality*>(l->FindObject("Centrality"));
 160   Float_t fCent=-1;
 161   if(centrality)
 162     fCent = centrality->GetCentralityPercentile("V0M");
 163   else
 164     fCent=99;//probably pp data
 165
 166   if(fCent<0) return;
 167
 168   if (fTracksName == "Tracks")
 169     am->LoadBranch("Tracks");
 170   tracks = dynamic_cast<TClonesArray*>(l->FindObject(fTracksName));
 171   if (!tracks) {
 172     AliError(Form("Pointer to tracks %s == 0", fTracksName.Data() ));
 173     return;
 174   }
 175   const Int_t Ntrks = tracks->GetEntries();
 176
 177   if (fCaloName == "CaloClusters")
 178     am->LoadBranch("CaloClusters");
 179   clus = dynamic_cast<TClonesArray*>(l->FindObject(fCaloName));
 180   if (!clus) {
 181     AliError(Form("Pointer to clus %s == 0", fCaloName.Data() ));
 182     return;
 183   }
 184
 185   Int_t centbin=-1;
 186   if(fCent>=0 && fCent<10) centbin=0;
 187   else if(fCent>=10 && fCent<30) centbin=1;
 188   else if(fCent>=30 && fCent<50) centbin=2;
 189   else if(fCent>=50 && fCent<100) centbin=3;
 190
 191   if (clus) {
 192     Double_t vertex[3] = {0, 0, 0};
 193     InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
 194     const Int_t Nclus = clus->GetEntries();
 195     for (Int_t iClus = 0, clusCount=0; iClus < Nclus; ++iClus) {
 196       AliVCluster *c = dynamic_cast<AliVCluster*>(clus->At(iClus));
 197       if (!c)
 198         continue;
 199       if (!c->IsEMCAL())
 200         continue;
 201       TLorentzVector nPart;
 202
 203       c->SetEmcCpvDistance(-1);
 204       c->SetTrackDistance(999,999);
 205       Double_t dEtaMin  = 1e9;
 206       Double_t dPhiMin  = 1e9;
 207       Int_t    imin     = -1;
 208       for(Int_t t = 0; t<Ntrks; ++t) {
 209         AliVTrack *track = dynamic_cast<AliVTrack*>(tracks->At(t));
 210         Double_t etadiff=999;
 211         Double_t phidiff=999;
 212         AliPicoTrack::GetEtaPhiDiff(track,c,phidiff,etadiff);
 213         Double_t dR = TMath::Sqrt(etadiff*etadiff+phidiff*phidiff);
 214         Double_t dRmin = TMath::Sqrt(dEtaMin*dEtaMin+dPhiMin*dPhiMin);
 215         if(dR > 25)
 216           continue;
 217         if(dR<dRmin){
 218           dEtaMin = etadiff;
 219           dPhiMin = phidiff;
 220           imin = t;
 221         }
 222       }
 223       c->SetEmcCpvDistance(imin);
 224       c->SetTrackDistance(dPhiMin, dEtaMin);
 225       c->GetMomentum(nPart, vertex);
 226       Double_t energy = nPart.P();
 227       if(energy<fMinPt) continue;
 228       if (imin>=0) {
 229         dPhiMin = c->GetTrackDx();
 230         dEtaMin = c->GetTrackDz();
 231         fHistMatchEtaPhi[centbin]->Fill(dEtaMin,dPhiMin);
 232       }
 233
 234       fHistNclusvsCent->Fill(fCent);
 235
 236       if (fHadCorr>0) {
 237         if (imin>=0) {
 238           dPhiMin = c->GetTrackDx();
 239           dEtaMin = c->GetTrackDz();
 240           Double_t dR=TMath::Sqrt(dEtaMin*dEtaMin+dPhiMin*dPhiMin);
 241
 242           AliVTrack *t = dynamic_cast<AliVTrack*>(tracks->At(imin));
 243           if (t) {
 244             if (t->Pt()<fMinPt)
 245               continue;
 246             if (t->P()>0) fHistMatchEvsP[centbin]->Fill(energy,energy/t->P());
 247             if (t->P()>0) fHistMatchdRvsEP[centbin]->Fill(dR,energy/t->P());
 248             fHistEbefore->Fill(fCent,energy);
 249             if (dPhiMin<0.05 && dEtaMin<0.025) { // pp cuts!!!
 250               energy -= fHadCorr*t->P();
 251               fHistNclusMatchvsCent->Fill(fCent);
 252             }
 253             if (energy<0)
 254               continue;
 255           }
 256         }
 257         fHistEafter->Fill(fCent,energy);
 258
 259       }//end had correction if
 260       if (energy>0){//Output new corrected clusters
 261         AliVCluster *oc;
 262         if (c->InheritsFrom("AliESDCaloCluster")) {
 263           oc = new ((*fOutClusters)[clusCount]) AliESDCaloCluster(*(dynamic_cast<AliESDCaloCluster*>(c)));
 264         }
 265         else if (c->InheritsFrom("AliAODCaloCluster")) {
 266           oc = new ((*fOutClusters)[clusCount]) AliAODCaloCluster(*(dynamic_cast<AliAODCaloCluster*>(c)));
 267         }
 268         else {
 269           AliError("Cluster type not recognized (nor ESD neither AOD)!");
 270           continue;
 271         }
 272         oc->SetE(energy);
 273         clusCount++;
 274       }
 275     }
 276
 277
 278   }
 279 }
 280
 281
 282 //________________________________________________________________________
 283 void AliHadCorrTask::Terminate(Option_t *)
 284 {
 285
 286 }
 287
 288 //________________________________________________________________________