| 3722909b |
1 | /************************************************************************** |
| 2 | * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. * |
| 3 | * * |
| 4 | * Author: The ALICE Off-line Project. * |
| 5 | * Contributors are mentioned in the code where appropriate. * |
| 6 | * * |
| 7 | * Permission to use, copy, modify and distribute this software and its * |
| 8 | * documentation strictly for non-commercial purposes is hereby granted * |
| 9 | * without fee, provided that the above copyright notice appears in all * |
| 10 | * copies and that both the copyright notice and this permission notice * |
| 11 | * appear in the supporting documentation. The authors make no claims * |
| 12 | * about the suitability of this software for any purpose. It is * |
| 13 | * provided "as is" without express or implied warranty. * |
| 14 | **************************************************************************/ |
| 15 | |
| 16 | //----------------------------------------------------------------- |
| 17 | // AliAnalysisTaskPhiEffMc class |
| 18 | //----------------------------------------------------------------- |
| 19 | |
| 20 | #include "TChain.h" |
| 21 | #include "TTree.h" |
| 22 | #include "TLegend.h" |
| 23 | #include "TH1F.h" |
| 24 | #include "TH2F.h" |
| 25 | #include "TH3F.h" |
| 26 | #include "THnSparse.h" |
| 27 | #include "TCanvas.h" |
| 28 | #include "TObjArray.h" |
| 29 | #include "AliAnalysisTask.h" |
| 30 | #include "AliAnalysisManager.h" |
| 31 | #include "AliVTrack.h" |
| 32 | #include "AliAODMCParticle.h" |
| 33 | #include "AliVParticle.h" |
| 34 | #include "AliAODEvent.h" |
| 35 | #include "AliAODInputHandler.h" |
| 36 | #include "AliAnalysisTaskPhiEffMc.h" |
| 37 | #include "AliAnalysisTaskESDfilter.h" |
| 38 | #include "AliAnalysisDataContainer.h" |
| 39 | #include "AliCentrality.h" |
| 40 | #include "TProof.h" |
| 41 | #include "AliPID.h" |
| 42 | #include "AliVEvent.h" |
| 43 | #include "AliPIDResponse.h" |
| 44 | #include "AliSpectraAODTrackCuts.h" |
| 45 | #include "AliSpectraAODEventCuts.h" |
| 46 | #include "AliStack.h" |
| 47 | #include <TMCProcess.h> |
| 48 | #include "AliAnalyseLeadingTrackUE.h" |
| 49 | |
| 50 | #include <iostream> |
| 51 | |
| 52 | using namespace std; |
| 53 | |
| 54 | ClassImp(AliAnalysisTaskPhiEffMc) |
| 55 | |
| 56 | //________________________________________________________________________ |
| 57 | AliAnalysisTaskPhiEffMc::AliAnalysisTaskPhiEffMc(const char *name) : AliAnalysisTaskSE(name), |
| 58 | fAOD(0x0), |
| 59 | fIsMC(0), |
| 60 | fOutput(0x0), |
| 61 | fHelperPID(0x0), |
| 62 | fTrackCuts(0x0), |
| 63 | fEventCuts(0x0) |
| 64 | { |
| 65 | // Default constructor |
| 66 | |
| 67 | |
| 68 | DefineInput(0, TChain::Class()); |
| 69 | DefineOutput(1, TList::Class()); |
| 70 | DefineOutput(2, AliHelperPID::Class()); |
| 71 | DefineOutput(3, AliSpectraAODEventCuts::Class()); |
| 72 | DefineOutput(4, AliSpectraAODTrackCuts::Class()); |
| 73 | } |
| 74 | //________________________________________________________________________ |
| 75 | //________________________________________________________________________ |
| 76 | void AliAnalysisTaskPhiEffMc::UserCreateOutputObjects() |
| 77 | { |
| 78 | Printf("\n\n\n\n\n\n In CreateOutput Object:"); |
| 79 | |
| 80 | if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro"); |
| 81 | if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro"); |
| 82 | if (!fHelperPID) AliFatal("HelperPID should be set in the steering macro"); |
| 83 | |
| 84 | fOutput = new TList(); |
| 85 | fOutput->SetOwner(); |
| 86 | fOutput->SetName("BlaBla"); |
| 87 | |
| 88 | // event histogram |
| 89 | const Int_t nEventDim = 4; |
| 90 | // cent vz MCmult RecoMult |
| 91 | Int_t nEventBins[4] = { 10, 41, 50, 50}; |
| 92 | Double_t nEventMin[4] = { 0, -10, 0, 0}; |
| 93 | Double_t nEventMax[4] = {100, 10, 100, 100}; |
| 94 | |
| 95 | THnSparseF* hEvent = new THnSparseF("hEvent", "hEvent", nEventDim, nEventBins, nEventMin, nEventMax); |
| 96 | hEvent->GetAxis(0)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data())); |
| 97 | hEvent->GetAxis(1)->SetTitle("v_{z} (cm)"); |
| 98 | hEvent->GetAxis(2)->SetTitle("# MC particles"); |
| 99 | hEvent->GetAxis(3)->SetTitle("# reconstructed tracks"); |
| 100 | fOutput->Add(hEvent); |
| 101 | |
| 102 | TH3F* hKcorr = new TH3F("hKcorr","hKcorr",100,0,10,100,0,10,100,0,10); |
| 103 | hKcorr->GetXaxis()->SetTitle("#phi p_{T} (GeV/c)"); |
| 104 | hKcorr->GetYaxis()->SetTitle("k+ p_{T} (GeV/c)"); |
| 105 | hKcorr->GetZaxis()->SetTitle("k- p_{T} (GeV/c)"); |
| 106 | fOutput->Add(hKcorr); |
| 107 | |
| 108 | // single charged particles -- bins for THnSparse histograms |
| 109 | const Int_t nTrackDim = 6; |
| 110 | // PID ID, pt, y, cent, eta, phi |
| 111 | Int_t nTrackBins[nTrackDim] = { 7, 50, 20, 10, 20, 30}; |
| 112 | Double_t nTrackMin[nTrackDim] = {-3.5, 0, -1, 0, -1, 0}; |
| 113 | Double_t nTrackMax[nTrackDim] = { 3.5, 5, 1, 100, 1, 2*TMath::Pi()}; |
| 114 | |
| 115 | // single charged particles -- Monte Carlo particles |
| 116 | THnSparseF* hTrackMc = new THnSparseF("hTrackMc", "hTrackMc", nTrackDim, nTrackBins, nTrackMin, nTrackMax); |
| 117 | hTrackMc->GetAxis(0)->SetTitle("PID ID * Charge"); |
| 118 | hTrackMc->GetAxis(1)->SetTitle("p_{T} (GeV/c)"); |
| 119 | hTrackMc->GetAxis(2)->SetTitle("y"); |
| 120 | hTrackMc->GetAxis(3)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data())); |
| 121 | hTrackMc->GetAxis(4)->SetTitle("#eta"); |
| 122 | hTrackMc->GetAxis(5)->SetTitle("#phi (rad)"); |
| 123 | fOutput->Add(hTrackMc); |
| 124 | |
| 125 | // single charged particles -- reconstructed tracks |
| 126 | THnSparseF* hTrackReco = new THnSparseF("hTrackReco", "hTrackReco", nTrackDim, nTrackBins, nTrackMin, nTrackMax); |
| 127 | hTrackReco->GetAxis(0)->SetTitle("PID ID * Charge"); |
| 128 | hTrackReco->GetAxis(1)->SetTitle("p_{T} (GeV/c)"); |
| 129 | hTrackReco->GetAxis(2)->SetTitle("y"); |
| 130 | hTrackReco->GetAxis(3)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data())); |
| 131 | hTrackReco->GetAxis(4)->SetTitle("#eta"); |
| 132 | hTrackReco->GetAxis(5)->SetTitle("#phi (rad)"); |
| 133 | fOutput->Add(hTrackReco); |
| 134 | |
| 135 | // single charged particles -- Monte Carlo particles matched to reconstructed tracks |
| 136 | THnSparseF* hTrackMatch = new THnSparseF("hTrackMatch", "hTrackMatch", nTrackDim, nTrackBins, nTrackMin, nTrackMax); |
| 137 | hTrackMatch->GetAxis(0)->SetTitle("PID ID * Charge"); |
| 138 | hTrackMatch->GetAxis(1)->SetTitle("p_{T} (GeV/c)"); |
| 139 | hTrackMatch->GetAxis(2)->SetTitle("y"); |
| 140 | hTrackMatch->GetAxis(3)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data())); |
| 141 | hTrackMatch->GetAxis(4)->SetTitle("#eta"); |
| 142 | hTrackMatch->GetAxis(5)->SetTitle("#phi (rad)"); |
| 143 | fOutput->Add(hTrackMatch); |
| 144 | |
| 145 | // kaon pairs -- bins for THnSparse histograms |
| 146 | const Int_t nPairDim = 7; |
| 147 | // InvMass, pt, y, cent, eta, phi, pair ID |
| 148 | Int_t nPairBins[nPairDim] = { 100, 50, 20, 10, 20, 30, 3}; |
| 149 | Double_t nPairMin[nPairDim] = {0.98, 0, -1, 0, -1, 0, -1.5}; |
| 150 | Double_t nPairMax[nPairDim] = { 1.1, 5, 1, 100, 1, 2*TMath::Pi(), 1.5}; |
| 151 | // pair ID = 0 unlike sign k+k- |
| 152 | // 1 like sign pos k+k+ |
| 153 | // -1 like sign neg k-k- |
| 154 | |
| 155 | // kaon pairs -- Monte Carlo particles -- real phi's |
| 156 | THnSparseF* hPairMcPhi = new THnSparseF("hPairMcPhi", "hPairMcPhi", nPairDim, nPairBins, nPairMin, nPairMax); |
| 157 | hPairMcPhi->GetAxis(0)->SetTitle("Invariant Mass (GeV/c/c)"); |
| 158 | hPairMcPhi->GetAxis(1)->SetTitle("p_{T} (GeV/c)"); |
| 159 | hPairMcPhi->GetAxis(2)->SetTitle("y"); |
| 160 | hPairMcPhi->GetAxis(3)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data())); |
| 161 | hPairMcPhi->GetAxis(4)->SetTitle("#eta"); |
| 162 | hPairMcPhi->GetAxis(5)->SetTitle("#phi (rad)"); |
| 163 | hPairMcPhi->GetAxis(6)->SetTitle("pair ID"); |
| 164 | fOutput->Add(hPairMcPhi); |
| 165 | |
| 166 | // kaon pairs -- Monte Carlo particles -- real phi's with kaon daughters that pass pt and eta cuts |
| 167 | THnSparseF* hPairMcPhiCuts = new THnSparseF("hPairMcPhiCuts", "hPairMcPhiCuts", nPairDim, nPairBins, nPairMin, nPairMax); |
| 168 | hPairMcPhiCuts->GetAxis(0)->SetTitle("Invariant Mass (GeV/c/c)"); |
| 169 | hPairMcPhiCuts->GetAxis(1)->SetTitle("p_{T} (GeV/c)"); |
| 170 | hPairMcPhiCuts->GetAxis(2)->SetTitle("y"); |
| 171 | hPairMcPhiCuts->GetAxis(3)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data())); |
| 172 | hPairMcPhiCuts->GetAxis(4)->SetTitle("#eta"); |
| 173 | hPairMcPhiCuts->GetAxis(5)->SetTitle("#phi (rad)"); |
| 174 | hPairMcPhiCuts->GetAxis(6)->SetTitle("pair ID"); |
| 175 | fOutput->Add(hPairMcPhiCuts); |
| 176 | |
| 177 | // kaon pairs -- Monte Carlo particles -- phi's reconstructed from MC kaons |
| 178 | THnSparseF* hPairMc = new THnSparseF("hPairMc", "hPairMc", nPairDim, nPairBins, nPairMin, nPairMax); |
| 179 | hPairMc->GetAxis(0)->SetTitle("Invariant Mass (GeV/c/c)"); |
| 180 | hPairMc->GetAxis(1)->SetTitle("p_{T} (GeV/c)"); |
| 181 | hPairMc->GetAxis(2)->SetTitle("y"); |
| 182 | hPairMc->GetAxis(3)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data())); |
| 183 | hPairMc->GetAxis(4)->SetTitle("#eta"); |
| 184 | hPairMc->GetAxis(5)->SetTitle("#phi (rad)"); |
| 185 | hPairMc->GetAxis(6)->SetTitle("pair ID"); |
| 186 | fOutput->Add(hPairMc); |
| 187 | |
| 188 | // kaon pairs -- reconstructed tracks |
| 189 | THnSparseF* hPairReco = new THnSparseF("hPairReco", "hPairReco", nPairDim, nPairBins, nPairMin, nPairMax); |
| 190 | hPairReco->GetAxis(0)->SetTitle("Invariant Mass (GeV/c/c)"); |
| 191 | hPairReco->GetAxis(1)->SetTitle("p_{T} (GeV/c)"); |
| 192 | hPairReco->GetAxis(2)->SetTitle("y"); |
| 193 | hPairReco->GetAxis(3)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data())); |
| 194 | hPairReco->GetAxis(4)->SetTitle("#eta"); |
| 195 | hPairReco->GetAxis(5)->SetTitle("#phi (rad)"); |
| 196 | hPairReco->GetAxis(6)->SetTitle("pair ID"); |
| 197 | fOutput->Add(hPairReco); |
| 198 | |
| 199 | // kaon pairs -- Monte Carlo particles matched to reconstructed tracks |
| 200 | THnSparseF* hPairMatch = new THnSparseF("hPairMatch", "hPairMatch", nPairDim, nPairBins, nPairMin, nPairMax); |
| 201 | hPairMatch->GetAxis(0)->SetTitle("Invariant Mass (GeV/c/c)"); |
| 202 | hPairMatch->GetAxis(1)->SetTitle("p_{T} (GeV/c)"); |
| 203 | hPairMatch->GetAxis(2)->SetTitle("y"); |
| 204 | hPairMatch->GetAxis(3)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data())); |
| 205 | hPairMatch->GetAxis(4)->SetTitle("#eta"); |
| 206 | hPairMatch->GetAxis(5)->SetTitle("#phi (rad)"); |
| 207 | hPairMatch->GetAxis(6)->SetTitle("pair ID"); |
| 208 | fOutput->Add(hPairMatch); |
| 209 | |
| 210 | PostData(1, fOutput); |
| 211 | PostData(2, fHelperPID); |
| 212 | PostData(3, fEventCuts); |
| 213 | PostData(4, fTrackCuts); |
| 214 | } |
| 215 | |
| 216 | //________________________________________________________________________ |
| 217 | void AliAnalysisTaskPhiEffMc::UserExec(Option_t *) |
| 218 | { |
| 219 | // main event loop |
| 220 | fAOD = dynamic_cast<AliAODEvent*>(fInputEvent); |
| 221 | if (!fAOD) { |
| 222 | AliWarning("ERROR: AliAODEvent not available \n"); |
| 223 | return; |
| 224 | } |
| 225 | |
| 226 | if (strcmp(fAOD->ClassName(), "AliAODEvent")) |
| 227 | { |
| 228 | AliFatal("Not processing AODs"); |
| 229 | } |
| 230 | |
| 231 | if(!fEventCuts->IsSelected(fAOD,fTrackCuts))return;//event selection |
| 232 | |
| 233 | Double_t cent = fEventCuts->GetCent(); |
| 234 | |
| 235 | THnSparseF* hEvent = (THnSparseF*)fOutput->FindObject("hEvent"); |
| 236 | THnSparseF* hTrackMc = (THnSparseF*)fOutput->FindObject("hTrackMc"); |
| 237 | THnSparseF* hTrackReco = (THnSparseF*)fOutput->FindObject("hTrackReco"); |
| 238 | THnSparseF* hTrackMatch = (THnSparseF*)fOutput->FindObject("hTrackMatch"); |
| 239 | THnSparseF* hPairMcPhi = (THnSparseF*)fOutput->FindObject("hPairMcPhi"); |
| 240 | THnSparseF* hPairMcPhiCuts = (THnSparseF*)fOutput->FindObject("hPairMcPhiCuts"); |
| 241 | THnSparseF* hPairMc = (THnSparseF*)fOutput->FindObject("hPairMc"); |
| 242 | THnSparseF* hPairReco = (THnSparseF*)fOutput->FindObject("hPairReco"); |
| 243 | THnSparseF* hPairMatch = (THnSparseF*)fOutput->FindObject("hPairMatch"); |
| 244 | TH3F* hKcorr = (TH3F*)fOutput->FindObject("hKcorr"); |
| 245 | |
| 246 | TObjArray* kaonsPosMc = new TObjArray(); |
| 247 | TObjArray* kaonsNegMc = new TObjArray(); |
| 248 | TObjArray* kaonsPosData = new TObjArray(); |
| 249 | TObjArray* kaonsNegData = new TObjArray(); |
| 250 | TObjArray* kaonsPosGen = new TObjArray(); |
| 251 | TObjArray* kaonsNegGen = new TObjArray(); |
| 252 | |
| 253 | Int_t nMc = 0; // number of charged MC particles which satisfy cuts (except PID) |
| 254 | Int_t nReco = 0; // number of reco tracks which satisfy cuts (except PID) |
| 255 | |
| 256 | //MC Loop |
| 257 | TClonesArray *arrayMC = 0; |
| 258 | if (fIsMC) |
| 259 | { |
| 260 | arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName()); |
| 261 | if (!arrayMC) { |
| 262 | AliFatal("Error: MC particles branch not found!\n"); |
| 263 | } |
| 264 | |
| 265 | |
| 266 | Int_t nMC = arrayMC->GetEntries(); |
| 267 | for (Int_t iMC = 0; iMC < nMC; iMC++) |
| 268 | { |
| 269 | AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC); |
| 270 | |
| 271 | if(partMC->Pt()<0.15) continue; |
| 272 | if(TMath::Abs(partMC->Eta())>0.8) continue; |
| 273 | |
| 274 | // PDG codes: pion(211), kaon(+/-321), proton(2212), phi(333) |
| 275 | |
| 276 | // make pt spectrum of only MC phis that decay to k+k- |
| 277 | if(partMC->GetPdgCode()==333) |
| 278 | { |
| 279 | if(partMC->GetNDaughters()==2) |
| 280 | { |
| 281 | if(TMath::Abs(((AliAODMCParticle*) arrayMC->At(partMC->GetDaughter(0)))->GetPdgCode()) == 321 && TMath::Abs(((AliAODMCParticle*) arrayMC->At(partMC->GetDaughter(1)))->GetPdgCode()) == 321) |
| 282 | { |
| 283 | AliAODMCParticle *d1 = (AliAODMCParticle*) arrayMC->At(partMC->GetDaughter(0)); |
| 284 | AliAODMCParticle *d2 = (AliAODMCParticle*) arrayMC->At(partMC->GetDaughter(1)); |
| 285 | if(d1->Charge() > 0) hKcorr->Fill(partMC->Pt(),d1->Pt(),d2->Pt()); |
| 286 | else hKcorr->Fill(partMC->Pt(),d2->Pt(),d1->Pt()); |
| 287 | TLorentzVector *d3 = (TLorentzVector*)makePhi(d1,d2); |
| 288 | // InvMass, pt, y, eta, phi, pair ID |
| 289 | Double_t invMass = (d3->M2() > 0 ? sqrt(d3->M2()) : 0); |
| 290 | if(invMass < 1.11) |
| 291 | { |
| 292 | Double_t varfill1[7] = {invMass, d3->Pt(), d3->Rapidity(), cent, d3->Eta(), (d3->Phi() > 0 ? d3->Phi() : d3->Phi()+2*TMath::Pi()), 0}; |
| 293 | hPairMcPhi->Fill(varfill1); |
| 294 | if(d1->Pt() > 0.15 && d2->Pt() > 0.15 && TMath::Abs(d1->Eta()) < 0.8 && TMath::Abs(d2->Eta()) < 0.8) |
| 295 | hPairMcPhiCuts->Fill(varfill1); |
| 296 | } |
| 297 | delete d3; |
| 298 | } |
| 299 | } |
| 300 | } |
| 301 | |
| 302 | if(partMC->Charge()==0) continue; // remove neutrals |
| 303 | if(!partMC->IsPhysicalPrimary()) continue; // need to remove all non-final state particles |
| 304 | |
| 305 | nMc++; |
| 306 | |
| 307 | Int_t mcID = fHelperPID->GetParticleSpecies(partMC); |
| 308 | if(mcID>3 || mcID < -3) continue; |
| 309 | // PID ID, pt, y, eta, phi |
| 310 | Double_t varfill2[6] = {(partMC->Charge() > 0 ? mcID+1 : -1*(mcID+1)), partMC->Pt(), partMC->Y(), cent, partMC->Eta(), partMC->Phi()}; |
| 311 | hTrackMc->Fill(varfill2); |
| 312 | |
| 313 | if(mcID == 1) |
| 314 | { |
| 315 | if(partMC->Charge() > 0) kaonsPosMc->Add(partMC); |
| 316 | else if(partMC->Charge() < 0) kaonsNegMc->Add(partMC); |
| 317 | } |
| 318 | |
| 319 | } // end MC track loop |
| 320 | |
| 321 | UnlikeSign(kaonsPosMc,kaonsNegMc,hPairMc,cent); |
| 322 | LikeSign(kaonsPosMc,kaonsNegMc,hPairMc,cent); |
| 323 | } |
| 324 | |
| 325 | |
| 326 | //track loop |
| 327 | for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) { |
| 328 | AliAODTrack* track = fAOD->GetTrack(iTracks); |
| 329 | if (!fTrackCuts->IsSelected(track,kTRUE)) continue; //track selection (rapidity selection NOT in the standard cuts?) |
| 330 | if(track->Charge()==0) continue; |
| 331 | if(track->Pt()<0.15) continue; |
| 332 | if(TMath::Abs(track->Eta())>0.8) continue; |
| 333 | |
| 334 | nReco++; |
| 335 | |
| 336 | Int_t dataID=fHelperPID->GetParticleSpecies(track,kTRUE); |
| 337 | if(dataID>3 || dataID < -3) continue; |
| 338 | |
| 339 | // PID ID, pt, y, eta, phi |
| 340 | Double_t varfill3[6] = {(track->Charge() > 0 ? dataID+1 : -1*(dataID+1)), track->Pt(), track->Y(), cent, track->Eta(), track->Phi()}; |
| 341 | hTrackReco->Fill(varfill3); |
| 342 | |
| 343 | if(dataID==1) |
| 344 | { |
| 345 | if(track->Charge() > 0) kaonsPosData->Add(track); |
| 346 | else if(track->Charge() < 0) kaonsNegData->Add(track); |
| 347 | } |
| 348 | |
| 349 | if (fIsMC) |
| 350 | { |
| 351 | Int_t mcLabel = track->GetLabel(); |
| 352 | AliAODMCParticle* tempMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(mcLabel)); |
| 353 | if(!tempMC) |
| 354 | { |
| 355 | AliError("Cannot get MC particle"); |
| 356 | continue; |
| 357 | } |
| 358 | |
| 359 | Int_t genID = fHelperPID->GetParticleSpecies(tempMC); |
| 360 | if(genID>3 || genID < -3) continue; |
| 361 | |
| 362 | // PID ID, pt, y, eta, phi |
| 363 | Double_t varfill4[6] = {(tempMC->Charge() > 0 ? genID+1 : -1*(genID+1)), tempMC->Pt(), tempMC->Y(), cent, tempMC->Eta(), tempMC->Phi()}; |
| 364 | hTrackMatch->Fill(varfill4); |
| 365 | |
| 366 | if(genID==1) |
| 367 | { |
| 368 | |
| 369 | if(tempMC->GetMother() >= 0) |
| 370 | { |
| 371 | if(((AliAODMCParticle*)arrayMC->At(tempMC->GetMother()))->GetPdgCode() == 333) |
| 372 | { |
| 373 | if(tempMC->Charge() > 0) |
| 374 | { |
| 375 | for(Int_t k = 0; k < kaonsNegGen->GetEntries(); k++) |
| 376 | { |
| 377 | if(tempMC->GetMother() != ((AliAODMCParticle*)arrayMC->At(((AliVParticle*)kaonsNegGen->UncheckedAt(k))->GetLabel()))->GetMother()) continue; |
| 378 | TLorentzVector*c = (TLorentzVector*)makePhi(track,(AliVParticle*)kaonsNegGen->UncheckedAt(k)); |
| 379 | Double_t invMass = (c->M2() > 0 ? sqrt(c->M2()) : 0); |
| 380 | if(invMass < 1.11) |
| 381 | { |
| 382 | Double_t fillgen[7] = {invMass, c->Pt(), c->Rapidity(), cent, c->Eta(), (c->Phi() > 0 ? c->Phi() : c->Phi()+2*TMath::Pi()), 0}; |
| 383 | hPairMatch->Fill(fillgen); |
| 384 | } |
| 385 | delete c; |
| 386 | } |
| 387 | } |
| 388 | else if(tempMC->Charge() < 0) |
| 389 | { |
| 390 | for(Int_t k = 0; k < kaonsPosGen->GetEntries(); k++) |
| 391 | { |
| 392 | if(tempMC->GetMother() != ((AliAODMCParticle*)arrayMC->At(((AliVParticle*)kaonsPosGen->UncheckedAt(k))->GetLabel()))->GetMother()) continue; |
| 393 | TLorentzVector*c = (TLorentzVector*)makePhi(track,(AliVParticle*)kaonsPosGen->UncheckedAt(k)); |
| 394 | Double_t invMass = (c->M2() > 0 ? sqrt(c->M2()) : 0); |
| 395 | if(invMass < 1.11) |
| 396 | { |
| 397 | Double_t fillgen2[7] = {invMass, c->Pt(), c->Rapidity(), cent, c->Eta(), (c->Phi() > 0 ? c->Phi() : c->Phi()+2*TMath::Pi()), 0}; |
| 398 | hPairMatch->Fill(fillgen2); |
| 399 | } |
| 400 | delete c; |
| 401 | } |
| 402 | } |
| 403 | } |
| 404 | } |
| 405 | |
| 406 | if(track->Charge() > 0) kaonsPosGen->Add(track); |
| 407 | else if(track->Charge() < 0) kaonsNegGen->Add(track); |
| 408 | } |
| 409 | } |
| 410 | |
| 411 | } // end loop on tracks |
| 412 | |
| 413 | UnlikeSign(kaonsPosData,kaonsNegData,hPairReco,cent); |
| 414 | LikeSign(kaonsPosData,kaonsNegData,hPairReco,cent); |
| 415 | |
| 416 | Double_t varfill5[4] = {cent, fAOD->GetPrimaryVertex()->GetZ(), nMc, nReco}; |
| 417 | hEvent->Fill(varfill5); |
| 418 | |
| 419 | kaonsPosMc->Clear(); |
| 420 | kaonsNegMc->Clear(); |
| 421 | kaonsPosData->Clear(); |
| 422 | kaonsNegData->Clear(); |
| 423 | kaonsPosGen->Clear(); |
| 424 | kaonsNegGen->Clear(); |
| 425 | |
| 426 | if(kaonsPosMc) delete kaonsPosMc; |
| 427 | if(kaonsNegMc) delete kaonsNegMc; |
| 428 | if(kaonsPosData) delete kaonsPosData; |
| 429 | if(kaonsNegData) delete kaonsNegData; |
| 430 | if(kaonsPosGen) delete kaonsPosGen; |
| 431 | if(kaonsNegGen) delete kaonsNegGen; |
| 432 | |
| 433 | PostData(1,fOutput); |
| 434 | PostData(2, fHelperPID); |
| 435 | PostData(3, fEventCuts); |
| 436 | PostData(4, fTrackCuts); |
| 437 | //Printf("............. end of Exec"); |
| 438 | |
| 439 | } |
| 440 | |
| 441 | //_________________________________________________________________ |
| 442 | void AliAnalysisTaskPhiEffMc::Terminate(Option_t *) |
| 443 | { |
| 444 | // Terminate analysis |
| 445 | // |
| 446 | fOutput = dynamic_cast<TList*>(GetOutputData(1)); |
| 447 | if (!fOutput) { |
| 448 | printf("ERROR: fOutput not available\n"); |
| 449 | return; |
| 450 | } |
| 451 | |
| 452 | printf("AliAnalysisTaskPhiEffMc: Terminate() \n"); |
| 453 | |
| 454 | } |
| 455 | |
| 456 | //_________________________________________________________________ |
| 457 | void AliAnalysisTaskPhiEffMc::UnlikeSign(TObjArray* kaonsPos, TObjArray* kaonsNeg, THnSparseF* h, Double_t cent) |
| 458 | { |
| 459 | Int_t countPos = kaonsPos->GetEntriesFast(); |
| 460 | Int_t countNeg = kaonsNeg->GetEntriesFast(); |
| 461 | |
| 462 | if(countPos+countNeg < 2) return; |
| 463 | |
| 464 | for(Int_t cp = 0; cp < countPos; cp++) // form k+k- pairs and compute invariant mass |
| 465 | { |
| 466 | AliVParticle* temp1 = (AliVParticle*)kaonsPos->UncheckedAt(cp); |
| 467 | for(Int_t cn = 0; cn < countNeg; cn++) |
| 468 | { |
| 469 | AliVParticle* temp2 = (AliVParticle*)kaonsNeg->UncheckedAt(cn); |
| 470 | TLorentzVector*c = (TLorentzVector*)makePhi(temp1,temp2); |
| 471 | // InvMass, pt, y, eta, phi, pair ID |
| 472 | Double_t invMass = (c->M2() > 0 ? sqrt(c->M2()) : 0); |
| 473 | if(invMass < 1.11) |
| 474 | { |
| 475 | Double_t varfill[7] = {invMass, c->Pt(), c->Rapidity(), cent, c->Eta(), (c->Phi() > 0 ? c->Phi() : c->Phi()+2*TMath::Pi()), 0}; |
| 476 | h->Fill(varfill); |
| 477 | } |
| 478 | delete c; |
| 479 | } |
| 480 | } |
| 481 | } |
| 482 | |
| 483 | void AliAnalysisTaskPhiEffMc::LikeSign(TObjArray* kaonsPos, TObjArray* kaonsNeg, THnSparseF* h, Double_t cent) |
| 484 | { |
| 485 | Int_t countPos = kaonsPos->GetEntriesFast(); |
| 486 | Int_t countNeg = kaonsNeg->GetEntriesFast(); |
| 487 | |
| 488 | if(countPos < 2 && countNeg < 2) return; |
| 489 | |
| 490 | for(Int_t cp = 0; cp < countPos; cp++) // form k+k+ pairs and compute invariant mass |
| 491 | { |
| 492 | AliVParticle* temp1 = (AliVParticle*)kaonsPos->UncheckedAt(cp); |
| 493 | for(Int_t cn = cp+1; cn < countPos; cn++) |
| 494 | { |
| 495 | AliVParticle* temp2 = (AliVParticle*)kaonsPos->UncheckedAt(cn); |
| 496 | TLorentzVector*c = (TLorentzVector*)makePhi(temp1,temp2); |
| 497 | // InvMass, pt, y, eta, phi, pair ID |
| 498 | Double_t invMass = (c->M2() > 0 ? sqrt(c->M2()) : 0); |
| 499 | if(invMass < 1.11) |
| 500 | { |
| 501 | //cout << "y = " << 0.5*log((c->E()+c->Pz())/(c->E()-c->Pz())) << " " << c->Rapidity() << endl; |
| 502 | //cout << "eta = " << -1.*log(tan(c->Theta()/2.)) << " " << c->Eta() << endl; |
| 503 | Double_t varfill1[7] = {invMass, c->Pt(), c->Rapidity(), cent, c->Eta(), (c->Phi() > 0 ? c->Phi() : c->Phi()+2*TMath::Pi()), 1}; |
| 504 | h->Fill(varfill1); |
| 505 | } |
| 506 | delete c; |
| 507 | } |
| 508 | } |
| 509 | for(Int_t cp = 0; cp < countNeg; cp++) // form k-k- pairs and compute invariant mass |
| 510 | { |
| 511 | AliVParticle* temp1 = (AliVParticle*)kaonsNeg->UncheckedAt(cp); |
| 512 | for(Int_t cn = cp+1; cn < countNeg; cn++) |
| 513 | { |
| 514 | AliVParticle* temp2 = (AliVParticle*)kaonsNeg->UncheckedAt(cn); |
| 515 | TLorentzVector*c = (TLorentzVector*)makePhi(temp1,temp2); |
| 516 | Double_t invMass = (c->M2() > 0 ? sqrt(c->M2()) : 0); |
| 517 | if(invMass < 1.11) |
| 518 | { |
| 519 | Double_t varfill2[7] = {invMass, c->Pt(), c->Rapidity(), cent, c->Eta(), (c->Phi() > 0 ? c->Phi() : c->Phi()+2*TMath::Pi()), -1}; |
| 520 | h->Fill(varfill2); |
| 521 | } |
| 522 | delete c; |
| 523 | } |
| 524 | } |
| 525 | } |
| 526 | |
| 527 | TLorentzVector* AliAnalysisTaskPhiEffMc::makePhi(AliVParticle* p1, AliVParticle* p2) |
| 528 | { |
| 529 | TLorentzVector* a = new TLorentzVector(p1->Px(),p1->Py(),p1->Pz(),sqrt(pow(p1->P(),2)+pow(4.93676999999999977e-01, 2))); |
| 530 | TLorentzVector* b = new TLorentzVector(p2->Px(),p2->Py(),p2->Pz(),sqrt(pow(p2->P(),2)+pow(4.93676999999999977e-01, 2))); |
| 531 | TLorentzVector* c = new TLorentzVector((*a)+(*b)); |
| 532 | delete a; |
| 533 | delete b; |
| 534 | return c; |
| 535 | } |
| 536 | |
git://git.uio.no / u / mrichter / AliRoot.git / blame