X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PWGPP%2FAliAnalysisTaskFilteredTree.cxx;h=5bcb187a2721d61e7b16724fe71e61028a9f5898;hb=35b193b453a046e5a5dfb0acf982699153966551;hp=adaddf34485d2eeff0524ab9ea4ef77e25718d85;hpb=8f63d5f478ed05171ac02d8ea4ae91ce298fba0f;p=u%2Fmrichter%2FAliRoot.git diff --git a/PWGPP/AliAnalysisTaskFilteredTree.cxx b/PWGPP/AliAnalysisTaskFilteredTree.cxx index adaddf34485..5bcb187a272 100644 --- a/PWGPP/AliAnalysisTaskFilteredTree.cxx +++ b/PWGPP/AliAnalysisTaskFilteredTree.cxx @@ -13,9 +13,29 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ +/* + Class to process/filter reconstruction information from ESD, ESD friends, MC and provide them for later reprocessing + Filtering schema - low pt part is downscaled - to have flat pt specatra of selected topologies (tracks and V0s) + Downscaling schema is controlled by downscaling factors + Usage: + 1.) Filtering on Lego train + 2.) expert QA for tracking (resolution efficnecy) + 3.) pt reoslution studies using V0s + 4.) dEdx calibration using V0s + 5.) pt resolution and dEdx studies using cosmic + + + 6.) Info used for later raw data OFFLINE triggering (highPt, V0, laser, cosmic, high dEdx) + + Exported trees (with full objects and dereived variables): + 1.) "highPt" - filtered trees with esd tracks, derived variables(propagated tracks), optional MC info +optional space points + 2.) "V0s" - - filtered trees with selected V0s (rough KF chi2 cut), KF particle and corresponding esd tracks + optionla space points + 3.) "Laser" - dump laser tracks with space points if exests + 4.) "CosmicTree" - cosmic track candidate (random or triggered) + esdTracks(up/down)+ optional points + 5.) "dEdx" - tree with high dEdx tpc tracks +*/ + #include "iostream" #include "TSystem.h" - #include #include @@ -36,6 +56,8 @@ #include "AliInputEventHandler.h" #include "AliStack.h" #include "AliTrackReference.h" +#include "AliTrackPointArray.h" +#include "AliSysInfo.h" #include "AliPhysicsSelection.h" #include "AliAnalysisTask.h" @@ -89,8 +111,10 @@ ClassImp(AliAnalysisTaskFilteredTree) , fCentralityEstimator(0) , fLowPtTrackDownscaligF(0) , fLowPtV0DownscaligF(0) + , fFriendDownscaling(-3.) , fProcessAll(kFALSE) , fProcessCosmics(kFALSE) + , fProcessITSTPCmatchOut(kFALSE) // swittch to process ITS/TPC standalone tracks , fHighPtTree(0) , fV0Tree(0) , fdEdxTree(0) @@ -105,7 +129,9 @@ ClassImp(AliAnalysisTaskFilteredTree) , fPtResEtaPtTPCITS(0) , fPtResCentPtTPC(0) , fPtResCentPtTPCc(0) - , fPtResCentPtTPCITS(0) + , fPtResCentPtTPCITS(0) + , fCurrentFileName("") + , fDummyTrack(0) { // Constructor @@ -122,16 +148,9 @@ ClassImp(AliAnalysisTaskFilteredTree) //_____________________________________________________________________________ AliAnalysisTaskFilteredTree::~AliAnalysisTaskFilteredTree() { - //the output trees not to be deleted in case of proof analysis - //Bool_t deleteTrees=kTRUE; - //if ((AliAnalysisManager::GetAnalysisManager())) - //{ - // if (AliAnalysisManager::GetAnalysisManager()->GetAnalysisType() == - // AliAnalysisManager::kProofAnalysis) - // deleteTrees=kFALSE; - //} - //if (deleteTrees) delete fTreeSRedirector; - + // + // Destructor + // delete fFilteredTreeEventCuts; delete fFilteredTreeAcceptanceCuts; delete fFilteredTreeRecAcceptanceCuts; @@ -141,14 +160,15 @@ AliAnalysisTaskFilteredTree::~AliAnalysisTaskFilteredTree() //____________________________________________________________________________ Bool_t AliAnalysisTaskFilteredTree::Notify() { + // + // static Int_t count = 0; count++; TTree *chain = (TChain*)GetInputData(0); - if(chain) - { + if(chain){ Printf("Processing %d. file: %s", count, chain->GetCurrentFile()->GetName()); - } - + } + fCurrentFileName=chain->GetCurrentFile()->GetName(); return kTRUE; } @@ -172,8 +192,9 @@ void AliAnalysisTaskFilteredTree::UserCreateOutputObjects() fMCEffTree = ((*fTreeSRedirector)<<"MCEffTree").GetTree(); fCosmicPairsTree = ((*fTreeSRedirector)<<"CosmicPairs").GetTree(); - - + if (!fDummyTrack) { + fDummyTrack=new AliESDtrack(); + } // histogram booking @@ -278,45 +299,42 @@ void AliAnalysisTaskFilteredTree::UserExec(Option_t *) Printf("ERROR: ESD event not available"); return; } - - //// MC event - //if(fUseMCInfo) { - // fMC = MCEvent(); - // if (!fMC) { - // Printf("ERROR: MC event not available"); - // return; - // } - //} - //if MC info available - use it. fMC = MCEvent(); - if(fUseESDfriends) { //fESDfriend = dynamic_cast(fESD->FindListObject("AliESDfriend")); fESDfriend = ESDfriend(); - if(!fESDfriend) { Printf("ERROR: ESD friends not available"); } } + AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler(); + if (!inputHandler){ + return; + } //if set, use the environment variables to set the downscaling factors //AliAnalysisTaskFilteredTree_fLowPtTrackDownscaligF //AliAnalysisTaskFilteredTree_fLowPtV0DownscaligF + //AliAnalysisTaskFilteredTree_fFriendDownscaling TString env; env = gSystem->Getenv("AliAnalysisTaskFilteredTree_fLowPtTrackDownscaligF"); - if (!env.IsNull()) - { + if (!env.IsNull()){ fLowPtTrackDownscaligF=env.Atof(); AliInfo(Form("fLowPtTrackDownscaligF=%f",fLowPtTrackDownscaligF)); } env = gSystem->Getenv("AliAnalysisTaskFilteredTree_fLowPtV0DownscaligF"); - if (!env.IsNull()) - { + if (!env.IsNull()){ fLowPtV0DownscaligF=env.Atof(); AliInfo(Form("fLowPtV0DownscaligF=%f",fLowPtTrackDownscaligF)); } - + env = gSystem->Getenv("AliAnalysisTaskFilteredTree_fFriendDownscaling"); + if (!env.IsNull()){ + fFriendDownscaling=env.Atof(); + AliInfo(Form(" fFriendDownscaling=%f",fFriendDownscaling)); + } + // + // // if(fProcessAll) { ProcessAll(fESD,fMC,fESDfriend); // all track stages and MC @@ -324,47 +342,21 @@ void AliAnalysisTaskFilteredTree::UserExec(Option_t *) else { Process(fESD,fMC,fESDfriend); // only global and TPC tracks } - // ProcessV0(fESD,fMC,fESDfriend); ProcessLaser(fESD,fMC,fESDfriend); ProcessdEdx(fESD,fMC,fESDfriend); - if (fProcessCosmics) { ProcessCosmics(fESD,fESDfriend); } if(fMC) { ProcessMCEff(fESD,fMC,fESDfriend);} + if (fProcessITSTPCmatchOut) ProcessITSTPCmatchOut(fESD, fESDfriend); + printf("processed event %d\n", Int_t(Entry())); } //_____________________________________________________________________________ void AliAnalysisTaskFilteredTree::ProcessCosmics(AliESDEvent *const event, AliESDfriend* esdFriend) { // - // Select real events with high-pT tracks - // - if(!event) { - AliDebug(AliLog::kError, "event not available"); - return; - } - - // - AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler(); - if (!inputHandler) - { - Printf("ERROR: Could not receive input handler"); - return; - } - - // get file name - TTree *chain = (TChain*)GetInputData(0); - if(!chain) { - Printf("ERROR: Could not receive input chain"); - return; - } - TObjString fileName(chain->GetCurrentFile()->GetName()); - - - // check for cosmic pairs - // - // find cosmic pairs trigger by random trigger + // Find cosmic pairs (triggered or random) // // AliESDVertex *vertexSPD = (AliESDVertex *)event->GetPrimaryVertexSPD(); @@ -374,12 +366,9 @@ void AliAnalysisTaskFilteredTree::ProcessCosmics(AliESDEvent *const event, AliES const Double_t kMinNcl=50; const Double_t kMaxDelta[5]={2,600,0.02,0.02,0.1}; Int_t ntracks=event->GetNumberOfTracks(); - // Float_t dcaTPC[2]={0,0}; - // Float_t covTPC[3]={0,0,0}; - UInt_t specie = event->GetEventSpecie(); // skip laser events if (specie==AliRecoParam::kCalib) return; - + Int_t ntracksFriend = esdFriend ? esdFriend->GetNumberOfTracks() : 0; for (Int_t itrack0=0;itrack0GetInnerParam(); - Bool_t newFriendTrack0=kFALSE; AliESDfriendTrack* friendTrack0=NULL; - if (esdFriend) {if (!esdFriend->TestSkipBit()) friendTrack0 = esdFriend->GetTrack(itrack0);} //this guy can be NULL - if (!friendTrack0) {friendTrack0=new AliESDfriendTrack(); newFriendTrack0=kTRUE;} + if (esdFriend &&!esdFriend->TestSkipBit()){ + if (itrack0GetTrack(itrack0); + } //this guy can be NULL + } for (Int_t itrack1=itrack0+1;itrack1GetTrack(itrack1); @@ -429,20 +420,16 @@ void AliAnalysisTaskFilteredTree::ProcessCosmics(AliESDEvent *const event, AliES if (TMath::Abs(TMath::Abs(track0->GetAlpha()-track1->GetAlpha())-TMath::Pi())>kMaxDelta[2]) isPair=kFALSE; //delta with correct sign /* - TCut cut0="abs(t1.fP[0]+t0.fP[0])<2" - TCut cut3="abs(t1.fP[3]+t0.fP[3])<0.02" - TCut cut4="abs(t1.fP[4]+t0.fP[4])<0.2" - */ + TCut cut0="abs(t1.fP[0]+t0.fP[0])<2" + TCut cut3="abs(t1.fP[3]+t0.fP[3])<0.02" + TCut cut4="abs(t1.fP[4]+t0.fP[4])<0.2" + */ if (TMath::Abs(par0[0]+par1[0])>kMaxDelta[0]) isPair=kFALSE; //delta y opposite sign if (TMath::Abs(par0[3]+par1[3])>kMaxDelta[3]) isPair=kFALSE; //delta tgl opposite sign if (TMath::Abs(AliTracker::GetBz())>1 && TMath::Abs(par0[4]+par1[4])>kMaxDelta[4]) isPair=kFALSE; //delta 1/pt opposite sign if (!isPair) continue; TString filename(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName()); Int_t eventNumber = event->GetEventNumberInFile(); - //Bool_t hasFriend = kFALSE; - //Bool_t hasITS=(track0->GetNcls(0)+track1->GetNcls(0)>4); - //printf("DUMPHPTCosmic:%s|%f|%d|%d|%d\n",filename.Data(),(TMath::Min(track0->Pt(),track1->Pt())), eventNumber,hasFriend,hasITS); - // const AliExternalTrackParam * trackIn1 = track1->GetInnerParam(); // // Int_t ntracksSPD = vertexSPD->GetNContributors(); @@ -453,56 +440,67 @@ void AliAnalysisTaskFilteredTree::ProcessCosmics(AliESDEvent *const event, AliES Float_t magField = event->GetMagneticField(); TObjString triggerClass = event->GetFiredTriggerClasses().Data(); - // - // Dump to the tree - // vertex - // TPC-ITS tracks - // - - //fCosmicPairsTree->Branch("fileName",&fileName,32000,0); - //fCosmicPairsTree->Branch("runNumber",&runNumber,"runNumber/I"); - //fCosmicPairsTree->Branch("timeStamp",&timeStamp,"timeStamp/I"); - //fCosmicPairsTree->Branch("eventNumber",&eventNumber,"eventNumber/I"); - //fCosmicPairsTree->Branch("triggerMask",&triggerMask,32000,0); - //fCosmicPairsTree->Branch("triggerClass",&triggerClass,32000,0); - //fCosmicPairsTree->Branch("magField",&magField,"magField/F"); - //fCosmicPairsTree->Branch("ntracksSPD",&ntracksSPD,"ntracksSPD/I"); - //fCosmicPairsTree->Branch("ntracksTPC",&ntracksTPC,"ntracksTPC/I"); - //fCosmicPairsTree->Branch("vertexSPD",vertexSPD,32000,0); - //fCosmicPairsTree->Branch("vertexTPC",vertexTPC,32000,0); - //fCosmicPairsTree->Branch("track0",track0,32000,0); - //fCosmicPairsTree->Branch("track1",track1,32000,0); - // - //fCosmicPairsTree->Fill(); + // Global event id calculation using orbitID, bunchCrossingID and periodID + ULong64_t orbitID = (ULong64_t)event->GetOrbitNumber(); + ULong64_t bunchCrossID = (ULong64_t)event->GetBunchCrossNumber(); + ULong64_t periodID = (ULong64_t)event->GetPeriodNumber(); + ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID); + - Bool_t newFriendTrack1=kFALSE; AliESDfriendTrack* friendTrack1=NULL; - if (esdFriend) {if (!esdFriend->TestSkipBit()) friendTrack1 = esdFriend->GetTrack(itrack1);} //this guy can be NULL - if (!friendTrack1) {friendTrack1=new AliESDfriendTrack(); newFriendTrack1=kTRUE;} + if (esdFriend &&!esdFriend->TestSkipBit()){ + if (itrack1GetTrack(itrack1); + } //this guy can be NULL + } + // + AliESDfriendTrack *friendTrackStore0=friendTrack0; // store friend track0 for later processing + AliESDfriendTrack *friendTrackStore1=friendTrack1; // store friend track1 for later processing + if (fFriendDownscaling>=1){ // downscaling number of friend tracks + if (gRandom->Rndm()>1./fFriendDownscaling){ + friendTrackStore0 = 0; + friendTrackStore1 = 0; + } + } + if (fFriendDownscaling<=0){ + if (((*fTreeSRedirector)<<"CosmicPairs").GetTree()){ + TTree * tree = ((*fTreeSRedirector)<<"CosmicPairs").GetTree(); + if (tree){ + Double_t sizeAll=tree->GetZipBytes(); + TBranch * br= tree->GetBranch("friendTrack0.fPoints"); + Double_t sizeFriend=(br!=NULL)?br->GetZipBytes():0; + br= tree->GetBranch("friendTrack0.fCalibContainer"); + if (br) sizeFriend+=br->GetZipBytes(); + if (sizeFriend*TMath::Abs(fFriendDownscaling)>sizeAll) { + friendTrackStore0=0; + friendTrackStore1=0; + } + } + } + } if(!fFillTree) return; if(!fTreeSRedirector) return; (*fTreeSRedirector)<<"CosmicPairs"<< - "fileName.="<<&fileName<< // file name - "runNumber="<GetInputEventHandler(); - if (!inputHandler) - { - Printf("ERROR: Could not receive input handler"); - return; - } - - // get file name - TTree *chain = (TChain*)GetInputData(0); - if(!chain) { - Printf("ERROR: Could not receive input chain"); - return; - } - TObjString fileName(chain->GetCurrentFile()->GetName()); - // trigger if(evtCuts->IsTriggerRequired()) { @@ -628,7 +608,7 @@ void AliAnalysisTaskFilteredTree::Process(AliESDEvent *const esdEvent, AliMCEven if(!vtxSPD) return; Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); - //printf("isEventOK %d, isEventTriggered %d, status %d, vz %f \n",isEventOK, isEventTriggered, vtxESD->GetStatus(), vtxESD->GetZv()); + //printf("isEventOK %d, isEventTriggered %d, status %d, vz %f \n",isEventOK, isEventTriggered, vtxESD->GetStatus(), vtxESD->GetZ()); //printf("GetAnalysisMode() %d \n",GetAnalysisMode()); Int_t ntracks = esdEvent->GetNumberOfTracks(); @@ -651,9 +631,9 @@ void AliAnalysisTaskFilteredTree::Process(AliESDEvent *const esdEvent, AliMCEven // //Double_t vert[3] = {0}; - //vert[0] = vtxESD->GetXv(); - //vert[1] = vtxESD->GetYv(); - //vert[2] = vtxESD->GetZv(); + //vert[0] = vtxESD->GetX(); + //vert[1] = vtxESD->GetY(); + //vert[2] = vtxESD->GetZ(); Int_t mult = vtxESD->GetNContributors(); Int_t multSPD = vtxSPD->GetNContributors(); Int_t multTPC = vtxTPC->GetNContributors(); @@ -662,6 +642,13 @@ void AliAnalysisTaskFilteredTree::Process(AliESDEvent *const esdEvent, AliMCEven Int_t runNumber = esdEvent->GetRunNumber(); Int_t evtTimeStamp = esdEvent->GetTimeStamp(); Int_t evtNumberInFile = esdEvent->GetEventNumberInFile(); + + // Global event id calculation using orbitID, bunchCrossingID and periodID + ULong64_t orbitID = (ULong64_t)esdEvent->GetOrbitNumber(); + ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber(); + ULong64_t periodID = (ULong64_t)esdEvent->GetPeriodNumber(); + ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID); + // high pT tracks for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++) @@ -692,29 +679,11 @@ void AliAnalysisTaskFilteredTree::Process(AliESDEvent *const esdEvent, AliMCEven // TPC-ITS tracks // TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data(); - - //fHighPtTree->Branch("fileName",&fileName,32000,0); - //fHighPtTree->Branch("runNumber",&runNumber,"runNumber/I"); - //fHighPtTree->Branch("evtTimeStamp",&evtTimeStamp,"evtTimeStamp/I"); - //fHighPtTree->Branch("evtNumberInFile",&evtNumberInFile,"evtNumberInFile/I"); - //fHighPtTree->Branch("triggerClass",&triggerClass,32000,0); - //fHighPtTree->Branch("Bz",&bz,"Bz/F"); - //fHighPtTree->Branch("vtxESD",vtxESD,32000,0); - //fHighPtTree->Branch("IRtot",&ir1,"IRtot/I"); - //fHighPtTree->Branch("IRint2",&ir2,"IRint2/I"); - //fHighPtTree->Branch("mult",&mult,"mult/I"); - //fHighPtTree->Branch("esdTrack",track,32000,0); - //fHighPtTree->Branch("centralityF",¢ralityF,"centralityF/F"); - - //fHighPtTree->Fill(); - - //Double_t vtxX=vtxESD->GetX(); - //Double_t vtxY=vtxESD->GetY(); - //Double_t vtxZ=vtxESD->GetZ(); if(!fFillTree) return; if(!fTreeSRedirector) return; (*fTreeSRedirector)<<"highPt"<< - "fileName.="<<&fileName<< + "gid="<GetCurrentFile()->GetName()); - + const Double_t kMinPt = 5; if(!fFillTree) return; if(!fTreeSRedirector) return; - - // laser events const AliESDHeader* esdHeader = esdEvent->GetHeader(); - if(esdHeader && esdHeader->GetEventSpecie()==AliRecoParam::kCalib) - { + if(esdHeader && esdHeader->GetEventSpecie()==AliRecoParam::kCalib) { Int_t countLaserTracks = 0; Int_t runNumber = esdEvent->GetRunNumber(); Int_t evtTimeStamp = esdEvent->GetTimeStamp(); Int_t evtNumberInFile = esdEvent->GetEventNumberInFile(); Float_t bz = esdEvent->GetMagneticField(); - TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data(); - - for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++) - { + TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data(); + // Global event id calculation using orbitID, bunchCrossingID and periodID + ULong64_t orbitID = (ULong64_t)esdEvent->GetOrbitNumber(); + ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber(); + ULong64_t periodID = (ULong64_t)esdEvent->GetPeriodNumber(); + ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID); + for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++){ AliESDtrack *track = esdEvent->GetTrack(iTrack); if(!track) continue; - - if(track->GetTPCInnerParam()) countLaserTracks++; - - Bool_t newFriendTrack=kFALSE; + if(track->GetTPCInnerParam()) countLaserTracks++; AliESDfriendTrack* friendTrack=NULL; - if (esdFriend) {if (!esdFriend->TestSkipBit()) friendTrack = esdFriend->GetTrack(iTrack);} //this guy can be NULL - if (!friendTrack) {friendTrack=new AliESDfriendTrack(); newFriendTrack=kTRUE;} - + // suppress beam background and CE random reacks + if (track->GetInnerParam()->Pt()Rndm()>1/(1+TMath::Abs(fFriendDownscaling)); + if (skipTrack) continue; + if (esdFriend) {if (!esdFriend->TestSkipBit()) friendTrack = esdFriend->GetTrack(iTrack);} //this guy can be NULL (*fTreeSRedirector)<<"Laser"<< - "fileName.="<<&fileName<< + "gid="<GetInputEventHandler(); - if (!inputHandler) - { - Printf("ERROR: Could not receive input handler"); - return; - } - // get file name - TTree *chain = (TChain*)GetInputData(0); - if(!chain) { - Printf("ERROR: Could not receive input chain"); - return; - } - TObjString fileName(chain->GetCurrentFile()->GetName()); // trigger if(evtCuts->IsTriggerRequired()) @@ -972,21 +908,24 @@ void AliAnalysisTaskFilteredTree::ProcessAll(AliESDEvent *const esdEvent, AliMCE // Double_t vert[3] = {0}; - vert[0] = vtxESD->GetXv(); - vert[1] = vtxESD->GetYv(); - vert[2] = vtxESD->GetZv(); + vert[0] = vtxESD->GetX(); + vert[1] = vtxESD->GetY(); + vert[2] = vtxESD->GetZ(); Int_t mult = vtxESD->GetNContributors(); Float_t bz = esdEvent->GetMagneticField(); Int_t runNumber = esdEvent->GetRunNumber(); Int_t evtTimeStamp = esdEvent->GetTimeStamp(); Int_t evtNumberInFile = esdEvent->GetEventNumberInFile(); + Int_t numberOfTracks=esdEvent->GetNumberOfTracks(); // high pT tracks - for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++) + for (Int_t iTrack = 0; iTrack < numberOfTracks; iTrack++) { AliESDtrack *track = esdEvent->GetTrack(iTrack); AliESDfriendTrack* friendTrack=NULL; - if (esdFriend) {if (!esdFriend->TestSkipBit()) friendTrack = esdFriend->GetTrack(iTrack);} //this guy can be NULL + Int_t numberOfFriendTracks=0; + if (esdFriend) numberOfFriendTracks=esdFriend->GetNumberOfTracks(); + if (esdFriend && iTrackTestSkipBit()) friendTrack = esdFriend->GetTrack(iTrack);} //this guy can be NULL if(!track) continue; if(track->Charge()==0) continue; if(!esdTrackCuts->AcceptTrack(track)) continue; @@ -1181,361 +1120,325 @@ void AliAnalysisTaskFilteredTree::ProcessAll(AliESDEvent *const esdEvent, AliMCE Bool_t isOKtrackInnerC3 = kFALSE; AliExternalTrackParam *trackInnerC3 = new AliExternalTrackParam(*(track->GetInnerParam())); - if(mcEvent) + if(mcEvent && stack) { - if(!stack) return; - multMCTrueTracks = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts); - // - // global track - // - Int_t label = TMath::Abs(track->GetLabel()); - if (label >= mcStackSize) continue; - particle = stack->Particle(label); - if (!particle) continue; - if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0.) - { - particleMother = GetMother(particle,stack); - mech = particle->GetUniqueID(); - isPrim = stack->IsPhysicalPrimary(label); - isFromStrangess = IsFromStrangeness(label,stack); - isFromConversion = IsFromConversion(label,stack); - isFromMaterial = IsFromMaterial(label,stack); - } - - // - // TPC track - // - Int_t labelTPC = TMath::Abs(track->GetTPCLabel()); - if (labelTPC >= mcStackSize) continue; - particleTPC = stack->Particle(labelTPC); - if (!particleTPC) continue; - if(particleTPC && particleTPC->GetPDG() && particleTPC->GetPDG()->Charge()!=0.) + do //artificial loop (once) to make the continue statements jump out of the MC part { - particleMotherTPC = GetMother(particleTPC,stack); - mechTPC = particleTPC->GetUniqueID(); - isPrimTPC = stack->IsPhysicalPrimary(labelTPC); - isFromStrangessTPC = IsFromStrangeness(labelTPC,stack); - isFromConversionTPC = IsFromConversion(labelTPC,stack); - isFromMaterialTPC = IsFromMaterial(labelTPC,stack); - } + multMCTrueTracks = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts); + // + // global track + // + Int_t label = TMath::Abs(track->GetLabel()); + if (label >= mcStackSize) continue; + particle = stack->Particle(label); + if (!particle) continue; + if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0.) + { + particleMother = GetMother(particle,stack); + mech = particle->GetUniqueID(); + isPrim = stack->IsPhysicalPrimary(label); + isFromStrangess = IsFromStrangeness(label,stack); + isFromConversion = IsFromConversion(label,stack); + isFromMaterial = IsFromMaterial(label,stack); + } - // - // store first track reference - // for TPC track - // - TParticle *part=0; - TClonesArray *trefs=0; - Int_t status = mcEvent->GetParticleAndTR(TMath::Abs(labelTPC), part, trefs); + // + // TPC track + // + Int_t labelTPC = TMath::Abs(track->GetTPCLabel()); + if (labelTPC >= mcStackSize) continue; + particleTPC = stack->Particle(labelTPC); + if (!particleTPC) continue; + if(particleTPC && particleTPC->GetPDG() && particleTPC->GetPDG()->Charge()!=0.) + { + particleMotherTPC = GetMother(particleTPC,stack); + mechTPC = particleTPC->GetUniqueID(); + isPrimTPC = stack->IsPhysicalPrimary(labelTPC); + isFromStrangessTPC = IsFromStrangeness(labelTPC,stack); + isFromConversionTPC = IsFromConversion(labelTPC,stack); + isFromMaterialTPC = IsFromMaterial(labelTPC,stack); + } - if(status>0 && part && trefs && part->GetPDG() && part->GetPDG()->Charge()!=0.) - { - Int_t nTrackRef = trefs->GetEntries(); - //printf("nTrackRef %d \n",nTrackRef); + // + // store first track reference + // for TPC track + // + TParticle *part=0; + TClonesArray *trefs=0; + Int_t status = mcEvent->GetParticleAndTR(TMath::Abs(labelTPC), part, trefs); - Int_t countITS = 0; - for (Int_t iref = 0; iref < nTrackRef; iref++) + if(status>0 && part && trefs && part->GetPDG() && part->GetPDG()->Charge()!=0.) { - AliTrackReference *ref = (AliTrackReference *)trefs->At(iref); + Int_t nTrackRef = trefs->GetEntries(); + //printf("nTrackRef %d \n",nTrackRef); - // Ref. in the middle ITS - if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kITS) + Int_t countITS = 0; + for (Int_t iref = 0; iref < nTrackRef; iref++) { - nrefITS++; - if(!refITS && countITS==2) { - refITS = ref; - //printf("refITS %p \n",refITS); + AliTrackReference *ref = (AliTrackReference *)trefs->At(iref); + + // Ref. in the middle ITS + if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kITS) + { + nrefITS++; + if(!refITS && countITS==2) { + refITS = ref; + //printf("refITS %p \n",refITS); + } + countITS++; } - countITS++; - } - // TPC - if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kTPC) - { - nrefTPC++; - refTPCOut=ref; - if(!refTPCIn) { - refTPCIn = ref; - //printf("refTPCIn %p \n",refTPCIn); - //break; + // TPC + if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kTPC) + { + nrefTPC++; + refTPCOut=ref; + if(!refTPCIn) { + refTPCIn = ref; + //printf("refTPCIn %p \n",refTPCIn); + //break; + } } - } - // TRD - if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kTRD) - { - nrefTRD++; - if(!refTRD) { - refTRD = ref; + // TRD + if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kTRD) + { + nrefTRD++; + if(!refTRD) { + refTRD = ref; + } } - } - // TOF - if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kTOF) - { - nrefTOF++; - if(!refTOF) { - refTOF = ref; + // TOF + if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kTOF) + { + nrefTOF++; + if(!refTOF) { + refTOF = ref; + } + } + // EMCAL + if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kEMCAL) + { + nrefEMCAL++; + if(!refEMCAL) { + refEMCAL = ref; + } } + // PHOS + // if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kPHOS) + // { + // nrefPHOS++; + // if(!refPHOS) { + // refPHOS = ref; + // } + // } } - // EMCAL - if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kEMCAL) + + // transform inner params to TrackRef + // reference frame + if(refTPCIn && trackInnerC3) { - nrefEMCAL++; - if(!refEMCAL) { - refEMCAL = ref; - } + Double_t kRefPhi = TMath::ATan2(refTPCIn->Y(),refTPCIn->X()); + isOKtrackInnerC3 = trackInnerC3->Rotate(kRefPhi); + isOKtrackInnerC3 = AliTracker::PropagateTrackToBxByBz(trackInnerC3,refTPCIn->R(),track->GetMass(),kStep,kFALSE); } - // PHOS - // if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kPHOS) - // { - // nrefPHOS++; - // if(!refPHOS) { - // refPHOS = ref; - // } - // } } - // transform inner params to TrackRef - // reference frame - if(refTPCIn && trackInnerC3) + // + // ITS track + // + Int_t labelITS = TMath::Abs(track->GetITSLabel()); + if (labelITS >= mcStackSize) continue; + particleITS = stack->Particle(labelITS); + if (!particleITS) continue; + if(particleITS && particleITS->GetPDG() && particleITS->GetPDG()->Charge()!=0.) { - Double_t kRefPhi = TMath::ATan2(refTPCIn->Y(),refTPCIn->X()); - isOKtrackInnerC3 = trackInnerC3->Rotate(kRefPhi); - isOKtrackInnerC3 = AliTracker::PropagateTrackToBxByBz(trackInnerC3,refTPCIn->R(),track->GetMass(),kStep,kFALSE); + particleMotherITS = GetMother(particleITS,stack); + mechITS = particleITS->GetUniqueID(); + isPrimITS = stack->IsPhysicalPrimary(labelITS); + isFromStrangessITS = IsFromStrangeness(labelITS,stack); + isFromConversionITS = IsFromConversion(labelITS,stack); + isFromMaterialITS = IsFromMaterial(labelITS,stack); } } - - // - // ITS track - // - Int_t labelITS = TMath::Abs(track->GetITSLabel()); - if (labelITS >= mcStackSize) continue; - particleITS = stack->Particle(labelITS); - if (!particleITS) continue; - if(particleITS && particleITS->GetPDG() && particleITS->GetPDG()->Charge()!=0.) - { - particleMotherITS = GetMother(particleITS,stack); - mechITS = particleITS->GetUniqueID(); - isPrimITS = stack->IsPhysicalPrimary(labelITS); - isFromStrangessITS = IsFromStrangeness(labelITS,stack); - isFromConversionITS = IsFromConversion(labelITS,stack); - isFromMaterialITS = IsFromMaterial(labelITS,stack); - } + while (0); } - // - Bool_t dumpToTree=kFALSE; - - if(isOKtpcInnerC && isOKtrackInnerC) dumpToTree = kTRUE; - //if(fUseESDfriends && isOKtrackInnerC2 && isOKouterITSc) dumpToTree = kTRUE; - if(isOKtrackInnerC2 && isOKouterITSc) dumpToTree = kTRUE; - if(mcEvent && isOKtrackInnerC3) dumpToTree = kTRUE; - TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data(); - if (fReducePileUp){ - // - // 18.03 - Reduce pile-up chunks, done outside of the ESDTrackCuts for 2012/2013 data pile-up about 95 % of tracks - // Done only in case no MC info // - Float_t dcaTPC[2]; - track->GetImpactParametersTPC(dcaTPC[0],dcaTPC[1]); - Bool_t isRoughPrimary = TMath::Abs(dcaTPC[1])<10; - Bool_t hasOuter=(track->IsOn(AliVTrack::kITSin))||(track->IsOn(AliVTrack::kTOFout))||(track->IsOn(AliVTrack::kTRDin)); - Bool_t keepPileUp=gRandom->Rndm()<0.05; - if ( (!hasOuter) && (!isRoughPrimary) && (!keepPileUp)){ - dumpToTree=kFALSE; + Bool_t dumpToTree=kFALSE; + + if(isOKtpcInnerC && isOKtrackInnerC) dumpToTree = kTRUE; + //if(fUseESDfriends && isOKtrackInnerC2 && isOKouterITSc) dumpToTree = kTRUE; + if(isOKtrackInnerC2 && isOKouterITSc) dumpToTree = kTRUE; + if(mcEvent && isOKtrackInnerC3) dumpToTree = kTRUE; + TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data(); + if (fReducePileUp){ + // + // 18.03 - Reduce pile-up chunks, done outside of the ESDTrackCuts for 2012/2013 data pile-up about 95 % of tracks + // Done only in case no MC info + // + Float_t dcaTPC[2]; + track->GetImpactParametersTPC(dcaTPC[0],dcaTPC[1]); + Bool_t isRoughPrimary = TMath::Abs(dcaTPC[1])<10; + Bool_t hasOuter=(track->IsOn(AliVTrack::kITSin))||(track->IsOn(AliVTrack::kTOFout))||(track->IsOn(AliVTrack::kTRDin)); + Bool_t keepPileUp=gRandom->Rndm()<0.05; + if ( (!hasOuter) && (!isRoughPrimary) && (!keepPileUp)){ + dumpToTree=kFALSE; + } } - } - ///////////////// - //book keeping of created dummy objects (to avoid NULL in trees) - Bool_t newvtxESD=kFALSE; - Bool_t newtrack=kFALSE; - Bool_t newtpcInnerC=kFALSE; - Bool_t newtrackInnerC=kFALSE; - Bool_t newtrackInnerC2=kFALSE; - Bool_t newouterITSc=kFALSE; - Bool_t newtrackInnerC3=kFALSE; - Bool_t newrefTPCIn=kFALSE; - Bool_t newrefITS=kFALSE; - Bool_t newparticle=kFALSE; - Bool_t newparticleMother=kFALSE; - Bool_t newparticleTPC=kFALSE; - Bool_t newparticleMotherTPC=kFALSE; - Bool_t newparticleITS=kFALSE; - Bool_t newparticleMotherITS=kFALSE; - Bool_t newFriendTrack=kFALSE; - - //check that the vertex is there and that it is OK, - //i.e. no null member arrays, otherwise a problem with merging - //later on. - //this is a very ugly hack! - if (!vtxESD) - { - AliInfo("fixing the ESD vertex for streaming"); - vtxESD=new AliESDVertex(); - //vtxESD->SetNContributors(1); - //UShort_t pindices[1]; pindices[0]=0; - //vtxESD->SetIndices(1,pindices); - //vtxESD->SetNContributors(0); - newvtxESD=kTRUE; - } - // - if (!track) {track=new AliESDtrack();newtrack=kTRUE;} - if (!friendTrack) {friendTrack=new AliESDfriendTrack(); newFriendTrack=kTRUE;} - if (!tpcInnerC) {tpcInnerC=new AliExternalTrackParam();newtpcInnerC=kTRUE;} - if (!trackInnerC) {trackInnerC=new AliExternalTrackParam();newtrackInnerC=kTRUE;} - if (!trackInnerC2) {trackInnerC2=new AliExternalTrackParam();newtrackInnerC2=kTRUE;} - if (!outerITSc) {outerITSc=new AliExternalTrackParam();newouterITSc=kTRUE;} - if (!trackInnerC3) {trackInnerC3=new AliExternalTrackParam();newtrackInnerC3=kTRUE;} - if (mcEvent) - { - if (!refTPCIn) {refTPCIn=new AliTrackReference(); newrefTPCIn=kTRUE;} - if (!refITS) {refITS=new AliTrackReference();newrefITS=kTRUE;} - if (!particle) {particle=new TParticle(); newparticle=kTRUE;} - if (!particleMother) {particleMother=new TParticle();newparticleMother=kTRUE;} - if (!particleTPC) {particleTPC=new TParticle();newparticleTPC=kTRUE;} - if (!particleMotherTPC) {particleMotherTPC=new TParticle();newparticleMotherTPC=kTRUE;} - if (!particleITS) {particleITS=new TParticle();newparticleITS=kTRUE;} - if (!particleMotherITS) {particleMotherITS=new TParticle();newparticleMotherITS=kTRUE;} - } - ///////////////////////// - - //Double_t vtxX=vtxESD->GetX(); - //Double_t vtxY=vtxESD->GetY(); - //Double_t vtxZ=vtxESD->GetZ(); - - //AliESDVertex* pvtxESD = (AliESDVertex*)vtxESD->Clone(); - //AliESDtrack* ptrack=(AliESDtrack*)track->Clone(); - //AliExternalTrackParam* ptpcInnerC = (AliExternalTrackParam*)tpcInnerC->Clone(); - //AliExternalTrackParam* ptrackInnerC = (AliExternalTrackParam*)trackInnerC->Clone(); - //AliExternalTrackParam* ptrackInnerC2 = (AliExternalTrackParam*)trackInnerC2->Clone(); - //AliExternalTrackParam* pouterITSc = (AliExternalTrackParam*)outerITSc->Clone(); - //AliExternalTrackParam* ptrackInnerC3 = (AliExternalTrackParam*)trackInnerC3->Clone(); - //AliESDVertex* pvtxESD = new AliESDVertex(*vtxESD); - //AliESDtrack* ptrack= new AliESDtrack(*track); - //AliExternalTrackParam* ptpcInnerC = new AliExternalTrackParam(*tpcInnerC); - //AliExternalTrackParam* ptrackInnerC = new AliExternalTrackParam(*trackInnerC); - //AliExternalTrackParam* ptrackInnerC2 = new AliExternalTrackParam(*trackInnerC2); - //AliExternalTrackParam* pouterITSc = new AliExternalTrackParam(*outerITSc); - //AliExternalTrackParam* ptrackInnerC3 = new AliExternalTrackParam(*trackInnerC3); - - Int_t ntracks = esdEvent->GetNumberOfTracks(); - - // fill histograms - FillHistograms(track, tpcInnerC, centralityF, (Double_t)chi2(0,0)); - - if(fTreeSRedirector && dumpToTree && fFillTree) - { - (*fTreeSRedirector)<<"highPt"<< - "fileName.="<<&fileName<< // name of the chunk file (hopefully full) - "runNumber="<GetZipBytes(); + TBranch * br= tree->GetBranch("friendTrack.fPoints"); + Double_t sizeFriend=(br!=NULL)?br->GetZipBytes():0; + br= tree->GetBranch("friendTrack.fCalibContainer"); + if (br) sizeFriend+=br->GetZipBytes(); + if (sizeFriend*TMath::Abs(fFriendDownscaling)>sizeAll) friendTrackStore=0; + } + } + } + + + // if (!friendTrackStore && fFriendDownscaling<=1) {friendTrack=fDummyFriendTrack;} + if (!vtxESD) {vtxESD=&dummyvtxESD;} if (mcEvent) { - AliTrackReference refDummy; - if (!refITS) refITS = &refDummy; - if (!refTRD) refTRD = &refDummy; - if (!refTOF) refTOF = &refDummy; - if (!refEMCAL) refEMCAL = &refDummy; - if (!refPHOS) refPHOS = &refDummy; - (*fTreeSRedirector)<<"highPt"<< - "multMCTrueTracks="<GetNumberOfTracks(); + + // Global event id calculation using orbitID, bunchCrossingID and periodID + ULong64_t orbitID = (ULong64_t)esdEvent->GetOrbitNumber(); + ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber(); + ULong64_t periodID = (ULong64_t)esdEvent->GetPeriodNumber(); + ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID); + + // fill histograms + FillHistograms(track, tpcInnerC, centralityF, (Double_t)chi2(0,0)); + + if(fTreeSRedirector && dumpToTree && fFillTree) { + (*fTreeSRedirector)<<"highPt"<< + "gid="<GetInputEventHandler(); - if (!inputHandler) - { - Printf("ERROR: Could not receive input handler"); - return; - } - // get file name - TTree *chain = (TChain*)GetInputData(0); - if(!chain) { - Printf("ERROR: Could not receive input chain"); - return; - } - TObjString fileName(chain->GetCurrentFile()->GetName()); // trigger if(evtCuts->IsTriggerRequired()) { // always MB isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB; - physicsSelection = static_cast (inputHandler->GetEventSelection()); if(!physicsSelection) return; @@ -1702,9 +1586,9 @@ void AliAnalysisTaskFilteredTree::ProcessMCEff(AliESDEvent *const esdEvent, AliM // reco event info Double_t vert[3] = {0}; - vert[0] = vtxESD->GetXv(); - vert[1] = vtxESD->GetYv(); - vert[2] = vtxESD->GetZv(); + vert[0] = vtxESD->GetX(); + vert[1] = vtxESD->GetY(); + vert[2] = vtxESD->GetZ(); Int_t mult = vtxESD->GetNContributors(); Double_t bz = esdEvent->GetMagneticField(); Double_t runNumber = esdEvent->GetRunNumber(); @@ -1780,7 +1664,7 @@ void AliAnalysisTaskFilteredTree::ProcessMCEff(AliESDEvent *const esdEvent, AliM } else if (trackIndex >-1) { recTrack = esdEvent->GetTrack(trackIndex); } else { - recTrack = new AliESDtrack(); + recTrack = fDummyTrack; } particleMother = GetMother(particle,stack); @@ -1798,7 +1682,7 @@ void AliAnalysisTaskFilteredTree::ProcessMCEff(AliESDEvent *const esdEvent, AliM // if(fTreeSRedirector && fFillTree) { (*fTreeSRedirector)<<"MCEffTree"<< - "fileName.="<<&fileName<< + "fileName.="<<&fCurrentFileName<< "triggerClass.="<<&triggerClass<< "runNumber="<GetInputEventHandler(); - if (!inputHandler) - { - Printf("ERROR: Could not receive input handler"); - return; - } - - // get file name - TTree *chain = (TChain*)GetInputData(0); - if(!chain) { - Printf("ERROR: Could not receive input chain"); - return; - } - TObjString fileName(chain->GetCurrentFile()->GetName()); // trigger if(evtCuts->IsTriggerRequired()) @@ -1952,6 +1823,12 @@ void AliAnalysisTaskFilteredTree::ProcessV0(AliESDEvent *const esdEvent, AliMCEv Int_t evNr=esdEvent->GetEventNumberInFile(); Double_t bz = esdEvent->GetMagneticField(); + // Global event id calculation using orbitID, bunchCrossingID and periodID + ULong64_t orbitID = (ULong64_t)esdEvent->GetOrbitNumber(); + ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber(); + ULong64_t periodID = (ULong64_t)esdEvent->GetPeriodNumber(); + ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID); + for (Int_t iv0=0; iv0TestSkipBit()) - { - friendTrack0 = esdFriend->GetTrack(v0->GetIndex(0)); //this guy can be NULL - friendTrack1 = esdFriend->GetTrack(v0->GetIndex(1)); //this guy can be NULL + if (esdFriend) { + if (!esdFriend->TestSkipBit()){ + Int_t ntracksFriend = esdFriend->GetNumberOfTracks(); + if (v0->GetIndex(0)GetTrack(v0->GetIndex(0)); //this guy can be NULL + } + if (v0->GetIndex(1)GetTrack(v0->GetIndex(1)); //this guy can be NULL + } } } - Bool_t newFriendTrack0=kFALSE; - Bool_t newFriendTrack1=kFALSE; - if (!friendTrack0) {friendTrack0=new AliESDfriendTrack(); newFriendTrack0=kTRUE;} - if (!friendTrack1) {friendTrack1=new AliESDfriendTrack(); newFriendTrack1=kTRUE;} if (track0->GetSign()<0) { track1 = esdEvent->GetTrack(v0->GetIndex(0)); track0 = esdEvent->GetTrack(v0->GetIndex(1)); } + + // + AliESDfriendTrack *friendTrackStore0=friendTrack0; // store friend track0 for later processing + AliESDfriendTrack *friendTrackStore1=friendTrack1; // store friend track1 for later processing + if (fFriendDownscaling>=1){ // downscaling number of friend tracks + if (gRandom->Rndm()>1./fFriendDownscaling){ + friendTrackStore0 = 0; + friendTrackStore1 = 0; + } + } + if (fFriendDownscaling<=0){ + if (((*fTreeSRedirector)<<"V0s").GetTree()){ + TTree * tree = ((*fTreeSRedirector)<<"V0s").GetTree(); + if (tree){ + Double_t sizeAll=tree->GetZipBytes(); + TBranch * br= tree->GetBranch("friendTrack0.fPoints"); + Double_t sizeFriend=(br!=NULL)?br->GetZipBytes():0; + br= tree->GetBranch("friendTrack0.fCalibContainer"); + if (br) sizeFriend+=br->GetZipBytes(); + if (sizeFriend*TMath::Abs(fFriendDownscaling)>sizeAll) { + friendTrackStore0=0; + friendTrackStore1=0; + } + } + } + } + // Bool_t isDownscaled = IsV0Downscaled(v0); if (isDownscaled) continue; @@ -1990,25 +1893,24 @@ void AliAnalysisTaskFilteredTree::ProcessV0(AliESDEvent *const esdEvent, AliMCEv if(!fFillTree) return; if(!fTreeSRedirector) return; (*fTreeSRedirector)<<"V0s"<< - "isDownscaled="<GetCurrentFile()->GetName()); // trigger Bool_t isEventTriggered = kTRUE; @@ -2049,11 +1944,6 @@ void AliAnalysisTaskFilteredTree::ProcessdEdx(AliESDEvent *const esdEvent, AliMC // AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler(); - if (!inputHandler) - { - Printf("ERROR: Could not receive input handler"); - return; - } if(evtCuts->IsTriggerRequired()) { @@ -2093,24 +1983,34 @@ void AliAnalysisTaskFilteredTree::ProcessdEdx(AliESDEvent *const esdEvent, AliMC if(isEventOK && isEventTriggered) { Double_t vert[3] = {0}; - vert[0] = vtxESD->GetXv(); - vert[1] = vtxESD->GetYv(); - vert[2] = vtxESD->GetZv(); + vert[0] = vtxESD->GetX(); + vert[1] = vtxESD->GetY(); + vert[2] = vtxESD->GetZ(); Int_t mult = vtxESD->GetNContributors(); Double_t bz = esdEvent->GetMagneticField(); Double_t runNumber = esdEvent->GetRunNumber(); Double_t evtTimeStamp = esdEvent->GetTimeStamp(); Int_t evtNumberInFile = esdEvent->GetEventNumberInFile(); - + + // Global event id calculation using orbitID, bunchCrossingID and periodID + ULong64_t orbitID = (ULong64_t)esdEvent->GetOrbitNumber(); + ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber(); + ULong64_t periodID = (ULong64_t)esdEvent->GetPeriodNumber(); + ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID); + // large dEdx for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++) { AliESDtrack *track = esdEvent->GetTrack(iTrack); if(!track) continue; - Bool_t newFriendTrack=kFALSE; AliESDfriendTrack* friendTrack=NULL; - if (esdFriend) {if (!esdFriend->TestSkipBit()) friendTrack = esdFriend->GetTrack(iTrack);} //this guy can be NULL - if (!friendTrack) {friendTrack=new AliESDfriendTrack(); newFriendTrack=kTRUE;} + if (esdFriend && !esdFriend->TestSkipBit()) { + Int_t ntracksFriend = esdFriend->GetNumberOfTracks(); + if (iTrackGetTrack(iTrack); + } //this guy can be NULL + } + if(track->Charge()==0) continue; if(!esdTrackCuts->AcceptTrack(track)) continue; if(!accCuts->AcceptTrack(track)) continue; @@ -2120,19 +2020,19 @@ void AliAnalysisTaskFilteredTree::ProcessdEdx(AliESDEvent *const esdEvent, AliMC if(!fFillTree) return; if(!fTreeSRedirector) return; - (*fTreeSRedirector)<<"dEdx"<< - "fileName.="<<&fileName<< + (*fTreeSRedirector)<<"dEdx"<< // high dEdx tree + "gid="<GetTgl())>1) return 0; if ((track0->GetTPCClusterInfo(2,1))<100) return 0; if ((track1->GetTPCClusterInfo(2,1))<100) return 0; - //if ((track0->GetITSclusters(0))<2) return 0; - //if ((track1->GetITSclusters(0))<2) return 0; Float_t pos0[2]={0}, cov0[3]={0}; Float_t pos1[2]={0}, cov1[3]={0}; track0->GetImpactParameters(pos0,cov0); @@ -2375,7 +2273,7 @@ TParticle *AliAnalysisTaskFilteredTree::GetMother(TParticle *const particle, Ali } //_____________________________________________________________________________ -Bool_t AliAnalysisTaskFilteredTree::IsFromConversion(const Int_t label, AliStack *const stack) +Bool_t AliAnalysisTaskFilteredTree::IsFromConversion(Int_t label, AliStack *const stack) { Bool_t isFromConversion = kFALSE; @@ -2407,7 +2305,7 @@ Bool_t AliAnalysisTaskFilteredTree::IsFromConversion(const Int_t label, AliStack } //_____________________________________________________________________________ -Bool_t AliAnalysisTaskFilteredTree::IsFromMaterial(const Int_t label, AliStack *const stack) +Bool_t AliAnalysisTaskFilteredTree::IsFromMaterial(Int_t label, AliStack *const stack) { Bool_t isFromMaterial = kFALSE; @@ -2437,7 +2335,7 @@ Bool_t AliAnalysisTaskFilteredTree::IsFromMaterial(const Int_t label, AliStack * } //_____________________________________________________________________________ -Bool_t AliAnalysisTaskFilteredTree::IsFromStrangeness(const Int_t label, AliStack *const stack) +Bool_t AliAnalysisTaskFilteredTree::IsFromStrangeness(Int_t label, AliStack *const stack) { Bool_t isFromStrangeness = kFALSE; @@ -2479,67 +2377,6 @@ void AliAnalysisTaskFilteredTree::FinishTaskOutput() // Called one at the end // locally on working node // - - //// must be deleted to store trees - //if(fTreeSRedirector) delete fTreeSRedirector; fTreeSRedirector=0; - - //// open temporary file and copy trees to the ouptut container - - //TChain* chain = 0; - //// - //chain = new TChain("highPt"); - //if(chain) { - // chain->Add("jotwinow_Temp_Trees.root"); - // fHighPtTree = chain->CopyTree("1"); - // delete chain; chain=0; - //} - //if(fHighPtTree) fHighPtTree->Print(); - - //// - //chain = new TChain("V0s"); - //if(chain) { - // chain->Add("jotwinow_Temp_Trees.root"); - // fV0Tree = chain->CopyTree("1"); - // delete chain; chain=0; - //} - //if(fV0Tree) fV0Tree->Print(); - - //// - //chain = new TChain("dEdx"); - //if(chain) { - // chain->Add("jotwinow_Temp_Trees.root"); - // fdEdxTree = chain->CopyTree("1"); - // delete chain; chain=0; - //} - //if(fdEdxTree) fdEdxTree->Print(); - - //// - //chain = new TChain("Laser"); - //if(chain) { - // chain->Add("jotwinow_Temp_Trees.root"); - // fLaserTree = chain->CopyTree("1"); - // delete chain; chain=0; - //} - //if(fLaserTree) fLaserTree->Print(); - - //// - //chain = new TChain("MCEffTree"); - //if(chain) { - // chain->Add("jotwinow_Temp_Trees.root"); - // fMCEffTree = chain->CopyTree("1"); - // delete chain; chain=0; - //} - //if(fMCEffTree) fMCEffTree->Print(); - - //// - //chain = new TChain("CosmicPairs"); - //if(chain) { - // chain->Add("jotwinow_Temp_Trees.root"); - // fCosmicPairsTree = chain->CopyTree("1"); - // delete chain; chain=0; - //} - //if(fCosmicPairsTree) fCosmicPairsTree->Print(); - Bool_t deleteTrees=kTRUE; if ((AliAnalysisManager::GetAnalysisManager())) { @@ -2549,40 +2386,14 @@ void AliAnalysisTaskFilteredTree::FinishTaskOutput() } if (deleteTrees) delete fTreeSRedirector; fTreeSRedirector=NULL; - - //OpenFile(1); - - // Post output data. - //PostData(1, fHighPtTree); - //PostData(2, fV0Tree); - //PostData(3, fdEdxTree); - //PostData(4, fLaserTree); - //PostData(5, fMCEffTree); - //PostData(6, fCosmicPairsTree); } //_____________________________________________________________________________ void AliAnalysisTaskFilteredTree::Terminate(Option_t *) { + // // Called one at the end - /* - fOutputSummary = dynamic_cast (GetOutputData(1)); - if(fOutputSummary) delete fOutputSummary; fOutputSummary=0; - TChain* chain = new TChain("highPt"); - if(!chain) return; - chain->Add("jotwinow_HighPt_TrackAndV0_Trees.root"); - TTree *tree = chain->CopyTree("1"); - if (chain) { delete chain; chain=0; } - if(!tree) return; - tree->Print(); - fOutputSummary = tree; - - if (!fOutputSummary) { - Printf("ERROR: AliAnalysisTaskFilteredTree::Terminate(): Output data not avaiable %p \n", GetOutputData(1)); - return; - } - */ - + // } //_____________________________________________________________________________ @@ -2636,7 +2447,7 @@ Int_t AliAnalysisTaskFilteredTree::GetMCTrueTrackMult(AliMCEvent *const mcEvent, } //_____________________________________________________________________________ -void AliAnalysisTaskFilteredTree::FillHistograms(AliESDtrack* const ptrack, AliExternalTrackParam* const ptpcInnerC, const Double_t centralityF, const Double_t chi2TPCInnerC) +void AliAnalysisTaskFilteredTree::FillHistograms(AliESDtrack* const ptrack, AliExternalTrackParam* const ptpcInnerC, Double_t centralityF, Double_t chi2TPCInnerC) { // // Fill pT relative resolution histograms for @@ -2706,3 +2517,205 @@ void AliAnalysisTaskFilteredTree::FillHistograms(AliESDtrack* const ptrack, AliE fPtResCentPtTPCc->Fill(ptpcInnerC->Pt(),centralityF,1./abs(ptpcInnerC->GetSigned1Pt())*TMath::Sqrt(ptpcInnerC->GetSigma1Pt2())); } } + + +void AliAnalysisTaskFilteredTree::ProcessITSTPCmatchOut(AliESDEvent *const esdEvent, AliESDfriend *const esdFriend){ + // + // Process ITS standalone tracks find match with closest TPC(or combined tracks) tracks + // marian.ivanov@cern.ch + // 0.) Init variables + // 1.) GetTrack parameters at TPC inner wall + // 2.) Match closest TPC track (STANDALONE/global) - chi2 match criteria + // + // Logic to be used in reco: + // 1.) Find matching ITSalone->TPCalone + // 2.) if (!TPCalone.FindClose(TPCother)) TPCalone.Addopt(ITSalone) + // 3.) ff ((ITSalone.FindClose(Global)==0) CreateGlobaltrack + const Double_t radiusMatch=84.; // redius to propagate + // + const Double_t dFastPhiCut=0.2; // 6 sigma (200 MeV) fast angular cut + const Double_t dFastThetaCut=0.12; // 6 sigma (200 MeV) fast angular cut + const Double_t dFastPosPhiCut=0.06; // 6 sigma (200 MeV) fast angular cut + const Double_t dFastZCut=6; // 6 sigma (200 MeV) fast z difference cut + const Double_t dFastPtCut=2.; // 6 sigma (200 MeV) fast 1/pt cut + const Double_t chi2Cut=100; // chi2 matching cut + // + if (!esdFriend) return; // not ITS standalone track + if (esdFriend->TestSkipBit()) return; // friends tracks not stored + Int_t ntracks=esdEvent->GetNumberOfTracks(); + Float_t bz = esdEvent->GetMagneticField(); + // + // 0.) Get parameters in reference radius TPC Inner wall + // + // + TMatrixD vecPosR0(ntracks,6); // possition and momentum estimate at reference radius + TMatrixD vecMomR0(ntracks,6); // + TMatrixD vecPosR1(ntracks,6); // possition and momentum estimate at reference radius TPC track + TMatrixD vecMomR1(ntracks,6); // + Double_t xyz[3], pxyz[3]; // + for (Int_t iTrack=0; iTrackGetTrack(iTrack); + if(!track) continue; + if (track->GetInnerParam()){ + const AliExternalTrackParam *trackTPC=track->GetInnerParam(); + trackTPC->GetXYZAt(radiusMatch,bz,xyz); + trackTPC->GetPxPyPzAt(radiusMatch,bz,pxyz); + for (Int_t i=0; i<3; i++){ + vecPosR1(iTrack,i)=xyz[i]; + vecMomR1(iTrack,i)=pxyz[i]; + } + vecPosR1(iTrack,3)= TMath::ATan2(xyz[1],xyz[0]); // phi pos angle + vecMomR1(iTrack,3)= TMath::ATan2(pxyz[1],pxyz[0]); // phi mom angle + vecMomR1(iTrack,4)= trackTPC->GetSigned1Pt();; + vecMomR1(iTrack,5)= trackTPC->GetTgl();; + } + AliESDfriendTrack* friendTrack=esdFriend->GetTrack(iTrack); + if(!friendTrack) continue; + if (friendTrack->GetITSOut()){ + const AliExternalTrackParam *trackITS=friendTrack->GetITSOut(); + trackITS->GetXYZAt(radiusMatch,bz,xyz); + trackITS->GetPxPyPzAt(radiusMatch,bz,pxyz); + for (Int_t i=0; i<3; i++){ + vecPosR0(iTrack,i)=xyz[i]; + vecMomR0(iTrack,i)=pxyz[i]; + } + vecPosR0(iTrack,3)= TMath::ATan2(xyz[1],xyz[0]); + vecMomR0(iTrack,3)= TMath::ATan2(pxyz[1],pxyz[0]); + vecMomR0(iTrack,4)= trackITS->GetSigned1Pt();; + vecMomR0(iTrack,5)= trackITS->GetTgl();; + } + } + // + // 1.) Find closest matching tracks, between the ITS standalone track + // and the all other tracks + // a.) caltegory - All + // b.) category - without ITS + // + // + Int_t ntracksPropagated=0; + AliExternalTrackParam extTrackDummy; + AliESDtrack esdTrackDummy; + AliExternalTrackParam itsAtTPC; + AliExternalTrackParam itsAtITSTPC; + for (Int_t iTrack0=0; iTrack0GetTrack(iTrack0); + if(!track0) continue; + if (track0->IsOn(AliVTrack::kTPCin)) continue; + AliESDfriendTrack* friendTrack0=esdFriend->GetTrack(iTrack0); + if (!friendTrack0) continue; + //if (!track0->IsOn(AliVTrack::kITSpureSA)) continue; + //if (!friendTrack0->GetITSOut()) continue; // is there flag for ITS standalone? + ntracksPropagated++; + // + // 2.) find clostest TPCtrack + // a.) all tracks + Double_t minChi2All=10000000; + Double_t minChi2TPC=10000000; + Double_t minChi2TPCITS=10000000; + Int_t indexAll=-1; + Int_t indexTPC=-1; + Int_t indexTPCITS=-1; + Int_t ncandidates0=0; // n candidates - rough cut + Int_t ncandidates1=0; // n candidates - rough + chi2 cut + itsAtTPC=*(friendTrack0->GetITSOut()); + itsAtITSTPC=*(friendTrack0->GetITSOut()); + for (Int_t iTrack1=0; iTrack1GetTrack(iTrack1); + if(!track1) continue; + if (!track1->IsOn(AliVTrack::kTPCin)) continue; + // fast checks + // + if (TMath::Abs(vecPosR1(iTrack1,2)-vecPosR0(iTrack0,2))>dFastZCut) continue; + if (TMath::Abs(vecPosR1(iTrack1,3)-vecPosR0(iTrack0,3))>dFastPosPhiCut) continue; + if (TMath::Abs(vecMomR1(iTrack1,3)-vecMomR0(iTrack0,3))>dFastPhiCut) continue; + if (TMath::Abs(vecMomR1(iTrack1,5)-vecMomR0(iTrack0,5))>dFastThetaCut) continue; + if (TMath::Abs(vecMomR1(iTrack1,4)-vecMomR0(iTrack0,4))>dFastPtCut) continue; + ncandidates0++; + // + const AliExternalTrackParam * param1= track1->GetInnerParam(); + if (!friendTrack0->GetITSOut()) continue; + AliExternalTrackParam outerITS = *(friendTrack0->GetITSOut()); + if (!outerITS.Rotate(param1->GetAlpha())) continue; + if (!outerITS.PropagateTo(param1->GetX(),bz)) continue; // assume track close to the TPC inner wall + Double_t chi2 = outerITS.GetPredictedChi2(param1); + if (chi2>chi2Cut) continue; + ncandidates1++; + if (chi2IsOn(AliVTrack::kITSin)==0){ + minChi2TPC=chi2; + indexTPC=iTrack1; + itsAtTPC=outerITS; + } + if (chi2IsOn(AliVTrack::kITSin)){ + minChi2TPCITS=chi2; + indexTPCITS=iTrack1; + itsAtITSTPC=outerITS; + } + } + // + AliESDtrack * trackAll= (indexAll>=0)? esdEvent->GetTrack(indexAll):&esdTrackDummy; + AliESDtrack * trackTPC= (indexTPC>=0)? esdEvent->GetTrack(indexTPC):&esdTrackDummy; + AliESDtrack * trackTPCITS= (indexTPCITS>=0)? esdEvent->GetTrack(indexTPCITS):&esdTrackDummy; + (*fTreeSRedirector)<<"itsTPC"<< + "indexAll="<