1 /**************************************************************************
\r
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
\r
4 * Author: The ALICE Off-line Project. *
\r
5 * Contributors are mentioned in the code where appropriate. *
\r
7 * Permission to use, copy, modify and distribute this software and its *
\r
8 * documentation strictly for non-commercial purposes is hereby granted *
\r
9 * without fee, provided that the above copyright notice appears in all *
\r
10 * copies and that both the copyright notice and this permission notice *
\r
11 * appear in the supporting documentation. The authors make no claims *
\r
12 * about the suitability of this software for any purpose. It is *
\r
13 * provided "as is" without express or implied warranty. *
\r
14 **************************************************************************/
\r
18 #include <TPDGCode.h>
\r
19 #include <TDatabasePDG.h>
\r
22 #include "TTreeStream.h"
\r
25 #include "TCanvas.h"
\r
28 #include "TMatrixD.h"
\r
29 #include "TRandom3.h"
\r
31 #include "AliHeader.h"
\r
32 #include "AliGenEventHeader.h"
\r
33 #include "AliInputEventHandler.h"
\r
34 #include "AliStack.h"
\r
35 #include "AliTrackReference.h"
\r
37 #include "AliPhysicsSelection.h"
\r
38 #include "AliAnalysisTask.h"
\r
39 #include "AliAnalysisManager.h"
\r
40 #include "AliESDEvent.h"
\r
41 #include "AliESDfriend.h"
\r
42 #include "AliMCEvent.h"
\r
43 #include "AliESDInputHandler.h"
\r
44 #include "AliESDVertex.h"
\r
45 #include "AliTracker.h"
\r
46 #include "AliGeomManager.h"
\r
48 #include "AliCentrality.h"
\r
49 #include "AliESDVZERO.h"
\r
50 #include "AliMultiplicity.h"
\r
52 #include "AliESDtrackCuts.h"
\r
53 #include "AliMCEventHandler.h"
\r
54 #include "AlidNdPt.h"
\r
55 #include "AlidNdPtEventCuts.h"
\r
56 #include "AlidNdPtAcceptanceCuts.h"
\r
58 #include "AlidNdPtTrackDumpTask.h"
\r
59 #include "AliKFParticle.h"
\r
60 #include "AliESDv0.h"
\r
62 using namespace std;
\r
64 ClassImp(AlidNdPtTrackDumpTask)
\r
66 //_____________________________________________________________________________
\r
67 AlidNdPtTrackDumpTask::AlidNdPtTrackDumpTask(const char *name)
\r
68 : AliAnalysisTaskSE(name)
\r
74 , fUseMCInfo(kFALSE)
\r
75 , fdNdPtEventCuts(0)
\r
76 , fdNdPtAcceptanceCuts(0)
\r
77 , fdNdPtRecAcceptanceCuts(0)
\r
79 , fTrigger(AliTriggerAnalysis::kMB1)
\r
80 , fAnalysisMode(AlidNdPtHelper::kTPC)
\r
82 , fTreeSRedirector(0)
\r
83 , fCentralityEstimator(0)
\r
84 , fLowPtTrackDownscaligF(0)
\r
85 , fLowPtV0DownscaligF(0)
\r
89 // Define input and output slots here
\r
90 DefineOutput(0, TTree::Class());
\r
91 //DefineOutput(1, TList::Class());
\r
95 //_____________________________________________________________________________
\r
96 AlidNdPtTrackDumpTask::~AlidNdPtTrackDumpTask()
\r
98 if(fOutput) delete fOutput; fOutput =0;
\r
99 //if(fOutputSummary) delete fOutputSummary; fOutputSummary =0;
\r
100 if(fTreeSRedirector) delete fTreeSRedirector; fTreeSRedirector =0;
\r
102 if(fdNdPtEventCuts) delete fdNdPtEventCuts; fdNdPtEventCuts=NULL;
\r
103 if(fdNdPtAcceptanceCuts) delete fdNdPtAcceptanceCuts; fdNdPtAcceptanceCuts=NULL;
\r
104 if(fdNdPtRecAcceptanceCuts) delete fdNdPtRecAcceptanceCuts; fdNdPtRecAcceptanceCuts=NULL;
\r
105 if(fEsdTrackCuts) delete fEsdTrackCuts; fEsdTrackCuts=NULL;
\r
108 //____________________________________________________________________________
\r
109 Bool_t AlidNdPtTrackDumpTask::Notify()
\r
111 static Int_t count = 0;
\r
113 TTree *chain = (TChain*)GetInputData(0);
\r
116 Printf("Processing %d. file: %s", count, chain->GetCurrentFile()->GetName());
\r
122 //_____________________________________________________________________________
\r
123 void AlidNdPtTrackDumpTask::UserCreateOutputObjects()
\r
125 // Create histograms
\r
127 fOutput = new TList;
\r
128 fOutput->SetOwner();
\r
131 // create output tree
\r
133 fTreeSRedirector = new TTreeSRedirector("jotwinow_HighPt_TrackAndV0_Trees.root");
\r
135 PostData(0, fOutputSummary);
\r
136 //PostData(1, fOutput);
\r
139 //_____________________________________________________________________________
\r
140 void AlidNdPtTrackDumpTask::UserExec(Option_t *)
\r
143 // Called for each event
\r
147 fESD = (AliESDEvent*) (InputEvent());
\r
149 Printf("ERROR: ESD event not available");
\r
157 Printf("ERROR: MC event not available");
\r
162 fESDfriend = static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
\r
164 Printf("ERROR: ESD friends not available");
\r
168 Process(fESD,fMC,fESDfriend);
\r
169 ProcessV0(fESD,fMC,fESDfriend);
\r
171 // Post output data.
\r
172 PostData(0, fOutputSummary);
\r
173 //PostData(1, fOutput);
\r
176 //_____________________________________________________________________________
\r
177 void AlidNdPtTrackDumpTask::Process(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const esdFriend)
\r
180 // Process real and/or simulated events
\r
183 AliDebug(AliLog::kError, "esdEvent not available");
\r
187 // get selection cuts
\r
188 AlidNdPtEventCuts *evtCuts = GetEventCuts();
\r
189 AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts();
\r
190 AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
\r
192 if(!evtCuts || !accCuts || !esdTrackCuts) {
\r
193 AliDebug(AliLog::kError, "cuts not available");
\r
197 // trigger selection
\r
198 Bool_t isEventTriggered = kTRUE;
\r
199 AliPhysicsSelection *physicsSelection = NULL;
\r
200 AliTriggerAnalysis* triggerAnalysis = NULL;
\r
203 AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
\r
206 Printf("ERROR: Could not receive input handler");
\r
211 TTree *chain = (TChain*)GetInputData(0);
\r
213 Printf("ERROR: Could not receive input chain");
\r
216 TObjString fileName(chain->GetCurrentFile()->GetName());
\r
219 if(evtCuts->IsTriggerRequired())
\r
222 isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
\r
224 physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
\r
225 if(!physicsSelection) return;
\r
226 //SetPhysicsTriggerSelection(physicsSelection);
\r
228 if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
\r
229 // set trigger (V0AND)
\r
230 triggerAnalysis = physicsSelection->GetTriggerAnalysis();
\r
231 if(!triggerAnalysis) return;
\r
232 isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
\r
236 // centrality determination
\r
237 Float_t centralityF = -1;
\r
238 AliCentrality *esdCentrality = esdEvent->GetCentrality();
\r
239 centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());
\r
241 // use MC information
\r
242 AliHeader* header = 0;
\r
243 AliGenEventHeader* genHeader = 0;
\r
244 AliStack* stack = 0;
\r
247 Int_t multMCTrueTracks = 0;
\r
252 AliDebug(AliLog::kError, "mcEvent not available");
\r
255 // get MC event header
\r
256 header = mcEvent->Header();
\r
258 AliDebug(AliLog::kError, "Header not available");
\r
261 // MC particle stack
\r
262 stack = mcEvent->Stack();
\r
264 AliDebug(AliLog::kError, "Stack not available");
\r
269 genHeader = header->GenEventHeader();
\r
271 AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
\r
274 genHeader->PrimaryVertex(vtxMC);
\r
276 // multipliticy of all MC primary tracks
\r
277 // in Zv, pt and eta ranges)
\r
278 multMCTrueTracks = AlidNdPtHelper::GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);
\r
282 // get reconstructed vertex
\r
283 //const AliESDVertex* vtxESD = 0;
\r
284 const AliESDVertex* vtxESD = 0;
\r
285 if(GetAnalysisMode() == AlidNdPtHelper::kTPC) {
\r
286 vtxESD = esdEvent->GetPrimaryVertexTPC();
\r
288 else if(GetAnalysisMode() == AlidNdPtHelper::kTPCITS) {
\r
289 vtxESD = esdEvent->GetPrimaryVertexTracks();
\r
295 if(!vtxESD) return;
\r
297 Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
\r
298 //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
\r
299 //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
\r
301 // check event cuts
\r
302 if(isEventOK && isEventTriggered)
\r
306 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
\r
308 AliESDtrack *track = esdEvent->GetTrack(iTrack);
\r
309 if(!track) continue;
\r
310 if(track->Charge()==0) continue;
\r
311 if(!esdTrackCuts->AcceptTrack(track)) continue;
\r
312 if(!accCuts->AcceptTrack(track)) continue;
\r
314 // downscale low-pT tracks
\r
315 Double_t scalempt= TMath::Min(track->Pt(),10.);
\r
316 if(TMath::Exp(2*scalempt)<fLowPtTrackDownscaligF*random.Rndm()) continue;
\r
318 // Dump to the tree
\r
320 // TPC constrained tracks
\r
321 // InnerParams constrained tracks
\r
323 // ITSout-InnerParams tracks
\r
324 // chi2 distance between TPC constrained and TPC-ITS tracks
\r
325 // chi2 distance between TPC refitted constrained and TPC-ITS tracks
\r
326 // chi2 distance between ITSout and InnerParams tracks
\r
329 Double_t x[3]; track->GetXYZ(x);
\r
330 Double_t b[3]; AliTracker::GetBxByBz(x,b);
\r
331 Bool_t isOK = kFALSE;
\r
334 // Constrain TPC-only tracks (TPCinner) to vertex
\r
336 AliExternalTrackParam * tpcInner = (AliExternalTrackParam *)(track->GetTPCInnerParam());
\r
337 if (!tpcInner) continue;
\r
338 // transform to the track reference frame
\r
339 isOK = tpcInner->Rotate(track->GetAlpha());
\r
340 isOK = tpcInner->PropagateTo(track->GetX(),esdEvent->GetMagneticField());
\r
341 if(!isOK) continue;
\r
343 // clone TPCinner has to be deleted
\r
344 AliExternalTrackParam * tpcInnerC = new AliExternalTrackParam(*(track->GetTPCInnerParam()));
\r
345 if (!tpcInnerC) continue;
\r
347 // constrain TPCinner
\r
348 //isOK = ConstrainTPCInner(tpcInnerC,vtxESD,esdEvent->GetMagneticField());
\r
349 isOK = ConstrainTPCInner(tpcInnerC,vtxESD,b);
\r
351 // transform to the track reference frame
\r
352 isOK = tpcInnerC->Rotate(track->GetAlpha());
\r
353 isOK = tpcInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField());
\r
356 if(tpcInnerC) delete tpcInnerC;
\r
362 // Constrain TPC refitted tracks at inner TPC wall (InnerParams)
\r
365 // clone track InnerParams has to be deleted
\r
366 AliExternalTrackParam * trackInnerC = new AliExternalTrackParam(*(track->GetInnerParam()));
\r
367 if (!trackInnerC) continue;
\r
369 // constrain track InnerParams
\r
370 isOK = ConstrainTrackInner(trackInnerC,vtxESD,track->GetMass(),b);
\r
372 // transform to the track reference frame
\r
373 isOK = trackInnerC->Rotate(track->GetAlpha());
\r
374 isOK = trackInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField());
\r
377 if(trackInnerC) delete trackInnerC;
\r
382 // calculate chi2 between vi and vj vectors
\r
383 // with covi and covj covariance matrices
\r
384 // chi2ij = (vi-vj)^(T)*(covi+covj)^(-1)*(vi-vj)
\r
386 TMatrixD deltaT(5,1), deltaTtrackC(5,1);
\r
387 TMatrixD delta(1,5), deltatrackC(1,5);
\r
388 TMatrixD covarM(5,5), covarMtrackC(5,5);
\r
390 for (Int_t ipar=0; ipar<5; ipar++) {
\r
391 deltaT(ipar,0)=tpcInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];
\r
392 delta(0,ipar)=tpcInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];
\r
394 deltaTtrackC(ipar,0)=trackInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];
\r
395 deltatrackC(0,ipar)=trackInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];
\r
397 for (Int_t jpar=0; jpar<5; jpar++) {
\r
398 Int_t index=track->GetIndex(ipar,jpar);
\r
399 covarM(ipar,jpar)=track->GetCovariance()[index]+tpcInnerC->GetCovariance()[index];
\r
400 covarMtrackC(ipar,jpar)=track->GetCovariance()[index]+trackInnerC->GetCovariance()[index];
\r
403 // chi2 distance TPC constrained and TPC+ITS
\r
404 TMatrixD covarMInv = covarM.Invert();
\r
405 TMatrixD mat2 = covarMInv*deltaT;
\r
406 TMatrixD chi2 = delta*mat2;
\r
409 // chi2 distance TPC refitted constrained and TPC+ITS
\r
410 TMatrixD covarMInvtrackC = covarMtrackC.Invert();
\r
411 TMatrixD mat2trackC = covarMInvtrackC*deltaTtrackC;
\r
412 TMatrixD chi2trackC = deltatrackC*mat2trackC;
\r
413 //chi2trackC.Print();
\r
417 // Propagate ITSout to TPC inner wall
\r
418 // and calculate chi2 distance to track (InnerParams)
\r
420 const Double_t kTPCRadius=85;
\r
421 const Double_t kStep=3;
\r
423 // clone track InnerParams has to be deleted
\r
424 AliExternalTrackParam *trackInnerC2 = new AliExternalTrackParam(*(track->GetInnerParam()));
\r
425 if (!trackInnerC2) continue;
\r
426 if(!AliTracker::PropagateTrackToBxByBz(trackInnerC2,kTPCRadius,track->GetMass(),kStep,kFALSE))
\r
428 if(trackInnerC2) delete trackInnerC2;
\r
432 AliExternalTrackParam *outerITSc = new AliExternalTrackParam();
\r
433 if(!outerITSc) continue;
\r
435 TMatrixD chi2OuterITS(1,1);
\r
436 chi2OuterITS(0,0) = 0;
\r
438 if(esdFriend && esdFriend->TestSkipBit()==kFALSE)
\r
440 // propagate ITSout to TPC inner wall
\r
441 AliESDfriendTrack *friendTrack = esdFriend->GetTrack(iTrack);
\r
445 if( (outerITSc = new AliExternalTrackParam(*friendTrack->GetITSOut())) )
\r
447 if(AliTracker::PropagateTrackToBxByBz(outerITSc,kTPCRadius,track->GetMass(),kStep,kFALSE))
\r
449 // transform ITSout to the track InnerParams reference frame
\r
451 isOK = outerITSc->Rotate(trackInnerC2->GetAlpha());
\r
452 isOK = outerITSc->PropagateTo(trackInnerC2->GetX(),esdEvent->GetMagneticField());
\r
455 if(outerITSc) delete outerITSc;
\r
456 if(trackInnerC2) delete trackInnerC2;
\r
461 // calculate chi2 between outerITS and innerParams
\r
462 // cov without z-coordinate at the moment
\r
464 TMatrixD deltaTouterITS(4,1);
\r
465 TMatrixD deltaouterITS(1,4);
\r
466 TMatrixD covarMouterITS(4,4);
\r
470 for (Int_t ipar=0; ipar<5; ipar++) {
\r
472 deltaTouterITS(kipar,0)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar];
\r
473 deltaouterITS(0,kipar)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar];
\r
477 for (Int_t jpar=0; jpar<5; jpar++) {
\r
478 Int_t index=outerITSc->GetIndex(ipar,jpar);
\r
479 if(ipar !=1 || jpar!=1) {
\r
480 covarMouterITS(kipar,kjpar)=outerITSc->GetCovariance()[index]+trackInnerC2->GetCovariance()[index];
\r
482 if(jpar!=1) kjpar++;
\r
484 if(ipar!=1) kipar++;
\r
487 // chi2 distance ITSout and InnerParams
\r
488 TMatrixD covarMInvouterITS = covarMouterITS.Invert();
\r
489 TMatrixD mat2outerITS = covarMInvouterITS*deltaTouterITS;
\r
490 chi2OuterITS = deltaouterITS*mat2outerITS;
\r
491 //chi2OuterITS.Print();
\r
500 TParticle *particle=NULL, *particleTPC=NULL, *particleITS=NULL;
\r
501 TParticle *particleMother=NULL, *particleMotherTPC=NULL, *particleMotherITS=NULL;
\r
502 Int_t mech=-1, mechTPC=-1, mechITS=-1;
\r
503 Bool_t isPrim=kFALSE, isPrimTPC=kFALSE, isPrimITS=kFALSE;
\r
504 Bool_t isFromStrangess=kFALSE, isFromStrangessTPC=kFALSE, isFromStrangessITS=kFALSE;
\r
505 Bool_t isFromConversion=kFALSE, isFromConversionTPC=kFALSE, isFromConversionITS=kFALSE;
\r
506 Bool_t isFromMaterial=kFALSE, isFromMaterialTPC=kFALSE, isFromMaterialITS=kFALSE;
\r
508 AliTrackReference *refTPCIn = NULL;
\r
509 AliTrackReference *refITS = NULL;
\r
510 AliExternalTrackParam *trackInnerC3 = new AliExternalTrackParam(*(track->GetInnerParam()));
\r
511 if(!trackInnerC3) continue;
\r
520 Int_t label = TMath::Abs(track->GetLabel());
\r
521 particle = stack->Particle(label);
\r
522 if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0.)
\r
524 particleMother = GetMother(particle,stack);
\r
525 mech = particle->GetUniqueID();
\r
526 isPrim = stack->IsPhysicalPrimary(label);
\r
527 isFromStrangess = IsFromStrangeness(label,stack);
\r
528 isFromConversion = IsFromConversion(label,stack);
\r
529 isFromMaterial = IsFromMaterial(label,stack);
\r
535 Int_t labelTPC = TMath::Abs(track->GetTPCLabel());
\r
536 particleTPC = stack->Particle(labelTPC);
\r
537 if(particleTPC && particleTPC->GetPDG() && particleTPC->GetPDG()->Charge()!=0.)
\r
539 particleMotherTPC = GetMother(particleTPC,stack);
\r
540 mechTPC = particleTPC->GetUniqueID();
\r
541 isPrimTPC = stack->IsPhysicalPrimary(labelTPC);
\r
542 isFromStrangessTPC = IsFromStrangeness(labelTPC,stack);
\r
543 isFromConversionTPC = IsFromConversion(labelTPC,stack);
\r
544 isFromMaterialTPC = IsFromMaterial(labelTPC,stack);
\r
548 // store first track reference
\r
552 TClonesArray *trefs=0;
\r
553 Int_t status = mcEvent->GetParticleAndTR(track->GetTPCLabel(), part, trefs);
\r
555 if(status>0 && part && trefs && part->GetPDG() && part->GetPDG()->Charge()!=0.)
\r
557 Int_t nTrackRef = trefs->GetEntries();
\r
558 //printf("nTrackRef %d \n",nTrackRef);
\r
560 Int_t countITS = 0;
\r
561 for (Int_t iref = 0; iref < nTrackRef; iref++)
\r
563 AliTrackReference *ref = (AliTrackReference *)trefs->At(iref);
\r
564 //printf("ref %p \n",ref);
\r
565 //if(ref) printf("ref->DetectorId() %d \n",ref->DetectorId());
\r
566 //printf("AliTrackReference::kTPC %d \n",AliTrackReference::kTPC);
\r
569 // Ref. in the middle ITS
\r
570 if(ref && ref->DetectorId()==AliTrackReference::kITS)
\r
572 if(!refITS && countITS==2) {
\r
574 //printf("refITS %p \n",refITS);
\r
580 if(ref && ref->DetectorId()==AliTrackReference::kTPC)
\r
584 //printf("refTPCIn %p \n",refTPCIn);
\r
590 // transform inner params to TrackRef
\r
592 if(refTPCIn && trackInnerC3)
\r
594 Double_t kRefPhi = TMath::ATan2(refTPCIn->Y(),refTPCIn->X());
\r
596 isOK = trackInnerC3->Rotate(kRefPhi);
\r
597 isOK = AliTracker::PropagateTrackToBxByBz(trackInnerC3,refTPCIn->R(),track->GetMass(),kStep,kFALSE);
\r
600 if(trackInnerC3) delete trackInnerC3;
\r
601 if(refTPCIn) delete refTPCIn;
\r
609 Int_t labelITS = TMath::Abs(track->GetITSLabel());
\r
610 particleITS = stack->Particle(labelITS);
\r
611 if(particleITS && particleITS->GetPDG() && particleITS->GetPDG()->Charge()!=0.)
\r
613 particleMotherITS = GetMother(particleITS,stack);
\r
614 mechITS = particleITS->GetUniqueID();
\r
615 isPrimITS = stack->IsPhysicalPrimary(labelITS);
\r
616 isFromStrangessITS = IsFromStrangeness(labelITS,stack);
\r
617 isFromConversionITS = IsFromConversion(labelITS,stack);
\r
618 isFromMaterialITS = IsFromMaterial(labelITS,stack);
\r
623 Double_t vert[3] = {0};
\r
624 vert[0] = vtxESD->GetXv();
\r
625 vert[1] = vtxESD->GetYv();
\r
626 vert[2] = vtxESD->GetZv();
\r
627 Int_t mult = vtxESD->GetNContributors();
\r
628 Double_t bz = esdEvent->GetMagneticField();
\r
629 Double_t runNumber = esdEvent->GetRunNumber();
\r
630 Double_t evtTimeStamp = esdEvent->GetTimeStamp();
\r
631 Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
\r
635 if(!fTreeSRedirector) return;
\r
636 (*fTreeSRedirector)<<"dNdPtTree"<<
\r
637 "fileName.="<<&fileName<<
\r
638 "runNumber="<<runNumber<<
\r
639 "evtTimeStamp="<<evtTimeStamp<<
\r
640 "evtNumberInFile="<<evtNumberInFile<<
\r
642 "vertX="<<vert[0]<<
\r
643 "vertY="<<vert[1]<<
\r
644 "vertZ="<<vert[2]<<
\r
646 "esdTrack.="<<track<<
\r
647 "extTPCInnerC.="<<tpcInnerC<<
\r
648 "extInnerParamC.="<<trackInnerC<<
\r
649 "extInnerParam.="<<trackInnerC2<<
\r
650 "extOuterITS.="<<outerITSc<<
\r
651 "extInnerParamRef.="<<trackInnerC3<<
\r
652 "refTPCIn.="<<refTPCIn<<
\r
653 "refITS.="<<refITS<<
\r
654 "chi2TPCInnerC="<<chi2(0,0)<<
\r
655 "chi2InnerC="<<chi2trackC(0,0)<<
\r
656 "chi2OuterITS="<<chi2OuterITS(0,0)<<
\r
657 "centralityF="<<centralityF<<
\r
658 "particle.="<<particle<<
\r
659 "particleMother.="<<particleMother<<
\r
661 "isPrim="<<isPrim<<
\r
662 "isFromStrangess="<<isFromStrangess<<
\r
663 "isFromConversion="<<isFromConversion<<
\r
664 "isFromMaterial="<<isFromMaterial<<
\r
665 "particleTPC.="<<particleTPC<<
\r
666 "particleMotherTPC.="<<particleMotherTPC<<
\r
667 "mechTPC="<<mechTPC<<
\r
668 "isPrimTPC="<<isPrimTPC<<
\r
669 "isFromStrangessTPC="<<isFromStrangessTPC<<
\r
670 "isFromConversionTPC="<<isFromConversionTPC<<
\r
671 "isFromMaterialTPC="<<isFromMaterialTPC<<
\r
672 "particleITS.="<<particleITS<<
\r
673 "particleMotherITS.="<<particleMotherITS<<
\r
674 "mechITS="<<mechITS<<
\r
675 "isPrimITS="<<isPrimITS<<
\r
676 "isFromStrangessITS="<<isFromStrangessITS<<
\r
677 "isFromConversionITS="<<isFromConversionITS<<
\r
678 "isFromMaterialITS="<<isFromMaterialITS<<
\r
681 if(tpcInnerC) delete tpcInnerC;
\r
682 if(trackInnerC) delete trackInnerC;
\r
683 if(trackInnerC2) delete trackInnerC2;
\r
684 if(outerITSc) delete outerITSc;
\r
686 if(trackInnerC3) delete trackInnerC3;
\r
690 PostData(0, fOutputSummary);
\r
691 //PostData(1, fOutput);
\r
694 //_____________________________________________________________________________
\r
695 void AlidNdPtTrackDumpTask::ProcessV0(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)
\r
698 // Process real and/or simulated events
\r
701 AliDebug(AliLog::kError, "esdEvent not available");
\r
705 // get selection cuts
\r
706 AlidNdPtEventCuts *evtCuts = GetEventCuts();
\r
707 AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts();
\r
708 AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
\r
710 if(!evtCuts || !accCuts || !esdTrackCuts) {
\r
711 AliDebug(AliLog::kError, "cuts not available");
\r
718 // trigger selection
\r
719 Bool_t isEventTriggered = kTRUE;
\r
720 AliPhysicsSelection *physicsSelection = NULL;
\r
721 AliTriggerAnalysis* triggerAnalysis = NULL;
\r
724 AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
\r
727 Printf("ERROR: Could not receive input handler");
\r
732 TTree *chain = (TChain*)GetInputData(0);
\r
734 Printf("ERROR: Could not receive input chain");
\r
737 TObjString fileName(chain->GetCurrentFile()->GetName());
\r
740 if(evtCuts->IsTriggerRequired())
\r
743 isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
\r
745 physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
\r
746 if(!physicsSelection) return;
\r
747 //SetPhysicsTriggerSelection(physicsSelection);
\r
749 if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
\r
750 // set trigger (V0AND)
\r
751 triggerAnalysis = physicsSelection->GetTriggerAnalysis();
\r
752 if(!triggerAnalysis) return;
\r
753 isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
\r
757 // centrality determination
\r
758 Float_t centralityF = -1;
\r
759 AliCentrality *esdCentrality = esdEvent->GetCentrality();
\r
760 centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());
\r
763 // get reconstructed vertex
\r
764 //const AliESDVertex* vtxESD = 0;
\r
765 const AliESDVertex* vtxESD = 0;
\r
766 if(GetAnalysisMode() == AlidNdPtHelper::kTPC) {
\r
767 vtxESD = esdEvent->GetPrimaryVertexTPC();
\r
769 else if(GetAnalysisMode() == AlidNdPtHelper::kTPCITS) {
\r
770 vtxESD = esdEvent->GetPrimaryVertexTracks();
\r
776 if(!vtxESD) return;
\r
778 Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
\r
779 //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
\r
780 //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
\r
782 // check event cuts
\r
783 if(isEventOK && isEventTriggered) {
\r
785 // Dump the pt downscaled V0 into the tree
\r
788 Int_t ntracks = esdEvent->GetNumberOfTracks();
\r
789 Int_t nV0s = esdEvent->GetNumberOfV0s();
\r
790 Int_t run = esdEvent->GetRunNumber();
\r
791 Int_t time= esdEvent->GetTimeStamp();
\r
792 Int_t evNr=esdEvent->GetEventNumberInFile();
\r
794 for (Int_t iv0=0; iv0<nV0s; iv0++){
\r
795 AliESDv0 * v0 = esdEvent->GetV0(iv0);
\r
797 AliESDtrack * track0 = esdEvent->GetTrack(v0->GetIndex(0));
\r
798 AliESDtrack * track1 = esdEvent->GetTrack(v0->GetIndex(1));
\r
799 if (!track0) continue;
\r
800 if (!track1) continue;
\r
801 if (track0->GetSign()<0) {
\r
802 track1 = esdEvent->GetTrack(v0->GetIndex(0));
\r
803 track0 = esdEvent->GetTrack(v0->GetIndex(1));
\r
806 Bool_t isDownscaled = IsV0Downscaled(v0);
\r
807 if (isDownscaled) continue;
\r
808 AliKFParticle kfparticle; //
\r
809 Int_t type=GetKFParticle(v0,esdEvent,kfparticle);
\r
810 if (type==0) continue;
\r
812 if(!fTreeSRedirector) return;
\r
813 (*fTreeSRedirector)<<"V0s"<<
\r
814 "isDownscaled="<<isDownscaled<<
\r
816 "fname="<<&fileName<<
\r
820 "ntracks="<<ntracks<<
\r
822 "kf.="<<&kfparticle<<
\r
823 "track0.="<<track0<<
\r
824 "track1.="<<track1<<
\r
825 "centralityF="<<centralityF<<
\r
829 PostData(0, fOutputSummary);
\r
832 //_____________________________________________________________________________
\r
833 Int_t AlidNdPtTrackDumpTask::GetKFParticle(AliESDv0 *const v0, AliESDEvent * const event, AliKFParticle & kfparticle)
\r
836 // Create KF particle in case the V0 fullfill selection criteria
\r
838 // Selection criteria
\r
839 // 0. algorithm cut
\r
842 // 4. rough mass cut
\r
843 // 5. Normalized pointing angle cut
\r
845 const Double_t cutMass=0.2;
\r
846 const Double_t kSigmaDCACut=3;
\r
848 // 0.) algo cut - accept only on the fly
\r
850 if (v0->GetOnFlyStatus() ==kFALSE) return 0;
\r
854 AliESDtrack * track0 = event->GetTrack(v0->GetIndex(0));
\r
855 AliESDtrack * track1 = event->GetTrack(v0->GetIndex(1));
\r
857 TCut cutD="abs(track0.fD/sqrt(track0.fCdd))>2&&abs(track1.fD/sqrt(track1.fCdd))>2";
\r
858 TCut cutTheta="abs(track0.fP[3])<1&&abs(track1.fP[3])<1";
\r
859 TCut cutNcl="track0.GetTPCClusterInfo(2,1)>100&&track1.GetTPCClusterInfo(2,1)>100";
\r
861 if (TMath::Abs(track0->GetTgl())>1) return 0;
\r
862 if (TMath::Abs(track1->GetTgl())>1) return 0;
\r
863 if ((track0->GetTPCClusterInfo(2,1))<100) return 0;
\r
864 if ((track1->GetTPCClusterInfo(2,1))<100) return 0;
\r
865 //if ((track0->GetITSclusters(0))<2) return 0;
\r
866 //if ((track1->GetITSclusters(0))<2) return 0;
\r
867 Float_t pos0[2]={0}, cov0[3]={0};
\r
868 Float_t pos1[2]={0}, cov1[3]={0};
\r
869 track0->GetImpactParameters(pos0,cov0);
\r
870 track0->GetImpactParameters(pos1,cov1);
\r
872 if (TMath::Abs(pos0[0])<kSigmaDCACut*TMath::Sqrt(cov0[0])) return 0;
\r
873 if (TMath::Abs(pos1[0])<kSigmaDCACut*TMath::Sqrt(cov1[0])) return 0;
\r
878 Double_t chi2KF = v0->GetKFInfo(2,2,2);
\r
879 if (chi2KF>25) return 0;
\r
881 // 4.) Rough mass cut - 0.200 GeV
\r
883 static Double_t masses[2]={-1};
\r
885 masses[0] = TDatabasePDG::Instance()->GetParticle("K_S0")->Mass();
\r
886 masses[1] = TDatabasePDG::Instance()->GetParticle("Lambda0")->Mass();
\r
888 Double_t mass00= v0->GetEffMass(0,0);
\r
889 Double_t mass22= v0->GetEffMass(2,2);
\r
890 Double_t mass42= v0->GetEffMass(4,2);
\r
891 Double_t mass24= v0->GetEffMass(2,4);
\r
892 Bool_t massOK=kFALSE;
\r
897 if (TMath::Abs(mass00-0)<cutMass) {
\r
898 massOK=kTRUE; type+=1;
\r
899 if (TMath::Abs(mass00-0)<dmass) {
\r
901 dmass=TMath::Abs(mass00-0);
\r
905 if (TMath::Abs(mass24-masses[1])<cutMass) {
\r
906 massOK=kTRUE; type+=2;
\r
907 if (TMath::Abs(mass24-masses[1])<dmass){
\r
908 dmass = TMath::Abs(mass24-masses[1]);
\r
913 if (TMath::Abs(mass42-masses[1])<cutMass) {
\r
914 massOK=kTRUE; type+=4;
\r
915 if (TMath::Abs(mass42-masses[1])<dmass){
\r
916 dmass = TMath::Abs(mass42-masses[1]);
\r
921 if (TMath::Abs(mass22-masses[0])<cutMass) {
\r
922 massOK=kTRUE; type+=8;
\r
923 if (TMath::Abs(mass22-masses[0])<dmass){
\r
924 dmass = TMath::Abs(mass22-masses[0]);
\r
929 if (type==0) return 0;
\r
931 const Int_t spdg[5]={kPositron,kMuonPlus,kPiPlus, kKPlus, kProton};
\r
932 const AliExternalTrackParam *paramP = v0->GetParamP();
\r
933 const AliExternalTrackParam *paramN = v0->GetParamN();
\r
934 if (paramP->GetSign()<0){
\r
935 paramP=v0->GetParamP();
\r
936 paramN=v0->GetParamN();
\r
938 //Double_t *pparam1 = (Double_t*)paramP->GetParameter();
\r
939 //Double_t *pparam2 = (Double_t*)paramN->GetParameter();
\r
941 AliKFParticle kfp1( *paramP, spdg[p1] );
\r
942 AliKFParticle kfp2( *paramN, -1 *spdg[p2] );
\r
943 AliKFParticle V0KF;
\r
950 Double_t errPhi = V0KF.GetErrPhi();
\r
951 Double_t pointAngle= TMath::ACos(v0->GetV0CosineOfPointingAngle());
\r
952 if (pointAngle/errPhi>10) return 0;
\r
957 //_____________________________________________________________________________
\r
958 Bool_t AlidNdPtTrackDumpTask::IsV0Downscaled(AliESDv0 *const v0)
\r
961 // Downscale randomly low pt V0
\r
964 Double_t maxPt= TMath::Max(v0->GetParamP()->Pt(), v0->GetParamN()->Pt());
\r
965 Double_t scalempt= TMath::Min(maxPt,10.);
\r
966 if (TMath::Exp(2*scalempt)<fLowPtV0DownscaligF*gRandom->Rndm()) return kTRUE;
\r
970 TH1F his1("his1","his1",100,0,10);
\r
971 TH1F his2("his2","his2",100,0,10);
\r
972 {for (Int_t i=0; i<10000; i++){
\r
973 Double_t rnd=gRandom->Exp(1);
\r
974 Bool_t isDownscaled =TMath::Exp(rnd)<100*gRandom->Rndm();
\r
976 if (!isDownscaled) his2->Fill(rnd);
\r
991 //_____________________________________________________________________________
\r
992 Bool_t AlidNdPtTrackDumpTask::ConstrainTPCInner(AliExternalTrackParam *const tpcInnerC, const AliESDVertex* vtx, Double_t b[3])
\r
994 // Constrain TPC inner params constrained
\r
996 if(!tpcInnerC) return kFALSE;
\r
997 if(!vtx) return kFALSE;
\r
999 Double_t dz[2],cov[3];
\r
1000 //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();
\r
1001 //if(!tpcInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE;
\r
1002 //if(!tpcInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE;
\r
1003 if(!tpcInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE;
\r
1006 Double_t covar[6]; vtx->GetCovMatrix(covar);
\r
1007 Double_t p[2]={tpcInnerC->GetParameter()[0]-dz[0],tpcInnerC->GetParameter()[1]-dz[1]};
\r
1008 Double_t c[3]={covar[2],0.,covar[5]};
\r
1009 Double_t chi2C=tpcInnerC->GetPredictedChi2(p,c);
\r
1010 if (chi2C>kVeryBig) return kFALSE;
\r
1012 if(!tpcInnerC->Update(p,c)) return kFALSE;
\r
1017 //_____________________________________________________________________________
\r
1018 Bool_t AlidNdPtTrackDumpTask::ConstrainTrackInner(AliExternalTrackParam *const trackInnerC, const AliESDVertex* vtx, Double_t mass, Double_t b[3])
\r
1020 // Constrain TPC inner params constrained
\r
1022 if(!trackInnerC) return kFALSE;
\r
1023 if(!vtx) return kFALSE;
\r
1025 const Double_t kRadius = 2.8;
\r
1026 const Double_t kMaxStep = 1.0;
\r
1028 Double_t dz[2],cov[3];
\r
1030 //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();
\r
1031 //if(!trackInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE;
\r
1032 //if(!trackInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE;
\r
1034 if(!AliTracker::PropagateTrackToBxByBz(trackInnerC,kRadius,mass,kMaxStep,kFALSE)) return kFALSE;
\r
1035 if(!trackInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE;
\r
1038 Double_t covar[6]; vtx->GetCovMatrix(covar);
\r
1039 Double_t p[2]={trackInnerC->GetParameter()[0]-dz[0],trackInnerC->GetParameter()[1]-dz[1]};
\r
1040 Double_t c[3]={covar[2],0.,covar[5]};
\r
1041 Double_t chi2C=trackInnerC->GetPredictedChi2(p,c);
\r
1042 if (chi2C>kVeryBig) return kFALSE;
\r
1044 if(!trackInnerC->Update(p,c)) return kFALSE;
\r
1050 //_____________________________________________________________________________
\r
1051 TParticle *AlidNdPtTrackDumpTask::GetMother(TParticle *const particle, AliStack *const stack)
\r
1053 if(!particle) return NULL;
\r
1054 if(!stack) return NULL;
\r
1056 Int_t motherLabel = TMath::Abs(particle->GetMother(0));
\r
1057 TParticle* mother = NULL;
\r
1058 mother = stack->Particle(motherLabel);
\r
1063 //_____________________________________________________________________________
\r
1064 Bool_t AlidNdPtTrackDumpTask::IsFromConversion(const Int_t label, AliStack *const stack)
\r
1066 Bool_t isFromConversion = kFALSE;
\r
1069 TParticle* particle = stack->Particle(label);
\r
1071 if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0)
\r
1073 Int_t mech = particle->GetUniqueID(); // production mechanism
\r
1074 Bool_t isPrim = stack->IsPhysicalPrimary(label);
\r
1076 Int_t motherLabel = TMath::Abs(particle->GetMother(0));
\r
1077 TParticle* mother = stack->Particle(motherLabel);
\r
1079 Int_t motherPdg = mother->GetPdgCode();
\r
1081 if(!isPrim && mech==5 && motherPdg==kGamma) {
\r
1082 isFromConversion=kTRUE;
\r
1088 return isFromConversion;
\r
1091 //_____________________________________________________________________________
\r
1092 Bool_t AlidNdPtTrackDumpTask::IsFromMaterial(const Int_t label, AliStack *const stack)
\r
1094 Bool_t isFromMaterial = kFALSE;
\r
1097 TParticle* particle = stack->Particle(label);
\r
1099 if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0)
\r
1101 Int_t mech = particle->GetUniqueID(); // production mechanism
\r
1102 Bool_t isPrim = stack->IsPhysicalPrimary(label);
\r
1104 Int_t motherLabel = TMath::Abs(particle->GetMother(0));
\r
1105 TParticle* mother = stack->Particle(motherLabel);
\r
1107 if(!isPrim && mech==13) {
\r
1108 isFromMaterial=kTRUE;
\r
1114 return isFromMaterial;
\r
1117 //_____________________________________________________________________________
\r
1118 Bool_t AlidNdPtTrackDumpTask::IsFromStrangeness(const Int_t label, AliStack *const stack)
\r
1120 Bool_t isFromStrangeness = kFALSE;
\r
1123 TParticle* particle = stack->Particle(label);
\r
1125 if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0)
\r
1127 Int_t mech = particle->GetUniqueID(); // production mechanism
\r
1128 Bool_t isPrim = stack->IsPhysicalPrimary(label);
\r
1130 Int_t motherLabel = TMath::Abs(particle->GetMother(0));
\r
1131 TParticle* mother = stack->Particle(motherLabel);
\r
1133 Int_t motherPdg = mother->GetPdgCode();
\r
1135 // K+-, lambda, antilambda, K0s decays
\r
1136 if(!isPrim && mech==4 &&
\r
1137 (TMath::Abs(motherPdg)==kKPlus || TMath::Abs(motherPdg)==kLambda0 || motherPdg==kK0Short))
\r
1139 isFromStrangeness = kTRUE;
\r
1145 return isFromStrangeness;
\r
1149 //_____________________________________________________________________________
\r
1150 void AlidNdPtTrackDumpTask::FinishTaskOutput()
\r
1153 // Called one at the end
\r
1154 // locally on working node
\r
1157 // Post output data.
\r
1158 PostData(1, fOutput);
\r
1159 //PostData(0, fOutputSummary);
\r
1162 //_____________________________________________________________________________
\r
1163 void AlidNdPtTrackDumpTask::Terminate(Option_t *)
\r
1165 // Called one at the end
\r
1166 if(fTreeSRedirector) delete fTreeSRedirector; fTreeSRedirector=0;
\r
1167 fOutputSummary = dynamic_cast<TTree*> (GetOutputData(0));
\r
1168 if(fOutputSummary) delete fOutputSummary; fOutputSummary=0;
\r
1170 TChain* chain = new TChain("dNdPtTree");
\r
1171 if(!chain) return;
\r
1172 chain->Add("jotwinow_HighPt_TrackAndV0_Trees.root");
\r
1173 TTree *tree = chain->CopyTree("1");
\r
1174 if (chain) { delete chain; chain=0; }
\r
1177 fOutputSummary = tree;
\r
1179 if (!fOutputSummary) {
\r
1180 Printf("ERROR: AlidNdPtTrackDumpTask::Terminate(): Output data not avaiable GetOutputData(0)==0x0 ..." );
\r
1184 PostData(0, fOutputSummary);
\r
1185 //PostData(1, fOutput);
\r