/************************************************************************* * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. * * * * Author: The ALICE Off-line Project. * * Contributors are mentioned in the code where appropriate. * * * * Permission to use, copy, modify and distribute this software and its * * documentation strictly for non-commercial purposes is hereby granted * * without fee, provided that the above copyright notice appears in all * * copies and that both the copyright notice and this permission notice * * appear in the supporting documentation. The authors make no claims * * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ /////////////////////////////////////////////////////////////////////////// // Class AliTrackletTaskMultipp // // Analysis task to produce data and MC histos needed for tracklets // // dNdEta extraction in multiple bins in one go // // Author: ruben.shahoyan@cern.ch // /////////////////////////////////////////////////////////////////////////// /* Important parameters to set: 1) make sure to initialize correct geometry in UserCreateOutputObjects 2) The cut on signal selection variable (delta, dphi ...) should be decided beforehand ... */ #include "TChain.h" #include "TTree.h" #include "TRandom.h" #include "TH1F.h" #include "TH2F.h" #include "TH3F.h" #include "THnSparse.h" #include "TList.h" #include "TNtuple.h" #include "TObjArray.h" #include "TGeoGlobalMagField.h" #include "AliAnalysisManager.h" #include "AliMultiplicity.h" #include "AliESDEvent.h" #include "AliESDInputHandler.h" #include "AliESDInputHandlerRP.h" #include "../ANALYSIS/EventMixing/AliMixEventInputHandler.h" #include "AliCDBPath.h" #include "AliCDBManager.h" #include "AliCDBEntry.h" #include "AliCDBStorage.h" #include "AliGeomManager.h" #include "AliMagF.h" #include "AliESDVZERO.h" #include "AliESDZDC.h" #include "AliRunLoader.h" #include "AliMCEventHandler.h" #include "AliMCEvent.h" #include "AliMCParticle.h" #include "AliStack.h" #include "AliGenEventHeader.h" #include "AliCentrality.h" #include "../ITS/AliITSRecPoint.h" #include "../ITS/AliITSgeomTGeo.h" #include "../ITS/AliITSMultReconstructor.h" #include "../ITS/AliITSsegmentationSPD.h" #include "AliLog.h" #include "AliPhysicsSelection.h" #include "AliTrackletTaskMultipp.h" #include "AliITSMultRecBg.h" #include "AliGenEventHeader.h" #include "AliGenHijingEventHeader.h" #include "AliGenDPMjetEventHeader.h" #include "AliESDtrackCuts.h" ClassImp(AliTrackletTaskMultipp) // centrality percentile (inverted: 100% - most central) const Float_t AliTrackletTaskMultipp::fgkCentPerc[] = {0,100};//{0,5,10,20,30}; //const Float_t AliTrackletTaskMultipp::fgkCentPerc[] = {0,5,10,20,30,40}; //{0,10,20,30,40,50,60,70,80,90,95,101}; const char* AliTrackletTaskMultipp::fgkCentSelName[] = {"V0M","FMD","TRK","TKL","CL0","CL1","V0MvsFMD","TKLvsV0M","ZENvsZDC"}; const char* AliTrackletTaskMultipp::fgkPDGNames[] = { "#pi^{+}", "p", "K^{+}", "K^{*+}", "e^{-}", "#mu^{-}", "#rho^{+}", "D^{+}", "D^{*+}", "D_{s}^{+}", "D_{s}^{*+}", "#Delta^{-}", "#Delta^{+}", "#Delta^{++}", "#Sigma^{-}", "#Sigma^{+}", "#Sigma^{*-}", "#Sigma^{*+}", "#Sigma^{*+}_{c}", "#Sigma^{*++}_{c}", "#Xi^{-}", "#Xi^{*-}", "#Lambda^{+}_{c}", "n", "#Delta^{0}", "#gamma", "K^{0}_{S}", "K^{0}_{L}", "K^{0}", "K^{*}", "#eta", "#pi^{0}", "#rho^{0}", "#varphi", "#eta'", "#omega", "#Lambda", "#Sigma^{0}", "#Sigma^{*0}_{c}", "#Sigma^{*0}", "D^{0}", "D^{*0}", "#Xi_{0}", "#Xi^{*0}", "#Xi^{0}_{c}", "#Xi^{*0}_{c}", "Nuclei", "Others" }; const int AliTrackletTaskMultipp::fgkPDGCodes[] = { 211, 2212, 321, 323, 11, 13, 213, 411, 413, 431, 433, 1114, 2214, 2224, 3112, 3222, 3114, 3224, 4214, 4224, 3312, 3314, 4122, 2112, 2114, 22, 310, 130, 311, 313, 221, 111, 113, 333, 331, 223, 3122, 3212, 4114, 3214, 421, 423, 3322, 3324, 4132, 4314 // nuclei // unknown }; //________________________________________________________________________ /*//Default constructor AliTrackletTaskMultipp::AliTrackletTaskMultipp(const char *name) : AliAnalysisTaskSE(name), */ //________________________________________________________________________ AliTrackletTaskMultipp::AliTrackletTaskMultipp(const char *name) : AliAnalysisTaskSE(name), // fOutput(0), // fDoNormalReco(kFALSE), fDoInjection(kFALSE), fDoRotation(kFALSE), fDoMixing(kFALSE), // fUseMC(kFALSE), fCheckReconstructables(kFALSE), // fHistosTrData(0), fHistosTrInj(0), fHistosTrRot(0), fHistosTrMix(0), // fHistosTrPrim(0), fHistosTrSec(0), fHistosTrComb(0), fHistosTrCombU(0), // fHistosTrRcblPrim(0), fHistosTrRcblSec(0), fHistosCustom(0), // fEtaMin(-3.0), fEtaMax(3.0), fZVertexMin(-20), fZVertexMax( 20), // fScaleDTBySin2T(kFALSE), fCutOnDThetaX(kFALSE), fNStdDev(1.), fDPhiWindow(0.08), fDThetaWindow(0.025), fDPhiShift(0.0045), fPhiOverlapCut(0.005), fZetaOverlap(0.05), fPhiRot(0.), fInjScale(1.), fRemoveOverlaps(kFALSE), // fDPhiSCut(0.06), fNStdCut(1.), fMCV0Scale(0.7520), // fMultReco(0), fRPTree(0), fRPTreeMix(0), fStack(0), fMCevent(0), // fNPart(0), fNBColl(0), fCurrCentBin(-1), fNCentBins(0), fUseCentralityVar(kCentV0M) { // Constructor DefineOutput(1, TList::Class()); // SetScaleDThetaBySin2T(); SetNStdDev(); SetPhiWindow(); SetThetaWindow(); SetPhiShift(); SetPhiOverlapCut(); SetZetaOverlapCut(); SetPhiRot(); SetRemoveOverlaps(); // } //________________________________________________________________________ AliTrackletTaskMultipp::~AliTrackletTaskMultipp() { // Destructor // histograms are in the output list and deleted when the output // list is deleted by the TSelector dtor if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) { //RRR printf("Deleteing output\n"); delete fOutput; } // delete fMultReco; // delete fHistosTrData; delete fHistosTrPrim; delete fHistosTrSec; delete fHistosTrComb; delete fHistosTrCombU; delete fHistosTrInj; delete fHistosTrRot; delete fHistosTrMix; delete fHistosTrRcblPrim; delete fHistosTrRcblSec; delete fHistosCustom; // } //________________________________________________________________________ void AliTrackletTaskMultipp::UserCreateOutputObjects() { // create output fOutput = new TList(); fOutput->SetOwner(); // // Bool_t needGeom = GetDoNormalReco() || GetDoInjection() || GetDoRotation() || GetDoMixing(); if (needGeom) { AliCDBManager *man = AliCDBManager::Instance(); if (fUseMC) { Bool_t newGeom = kTRUE; man->SetDefaultStorage("alien://Folder=/alice/simulation/2008/v4-15-Release/Residual"); if (newGeom) { // new geom AliCDBEntry* obj = man->Get("GRP/Geometry/Data",130844,8); AliGeomManager::SetGeometry((TGeoManager*) obj->GetObject()); if (!AliGeomManager::ApplyAlignObjsToGeom("ITS",130844,6,-1)) AliFatal("Failed to misalign geometry"); } else { // old geom AliCDBEntry* obj = man->Get("GRP/Geometry/Data",130845,7); AliGeomManager::SetGeometry((TGeoManager*) obj->GetObject()); if (!AliGeomManager::ApplyAlignObjsToGeom("ITS",130845,5,-1)) AliFatal("Failed to misalign geometry"); } } else { man->SetDefaultStorage("alien://Folder=/alice/data/2010/OCDB"); //man->SetRun(137366); AliCDBEntry* obj = man->Get("GRP/Geometry/Data",137366); AliGeomManager::SetGeometry((TGeoManager*) obj->GetObject()); if (!AliGeomManager::ApplyAlignObjsToGeom("ITS",137366,-1,-1)) AliFatal("Failed to misalign geometry"); } } // // Create histograms fNCentBins = sizeof(fgkCentPerc)/sizeof(Float_t)-1; //---------------------------------------------Standard histos per tracklet type--->> UInt_t hPattern = 0xffffffff; fHistosTrData = BookHistosSet("TrData",hPattern); hPattern &= ~(BIT(kHEtaZvSPD1)); // fill single clusters for "data" only if (GetDoInjection()) fHistosTrInj = BookHistosSet("TrInj",hPattern); if (GetDoRotation()) fHistosTrRot = BookHistosSet("TrRot",hPattern); if (GetDoMixing()) fHistosTrMix = BookHistosSet("TrMix",hPattern); if (fUseMC) { fHistosTrPrim = BookHistosSet("TrPrim",hPattern); fHistosTrSec = BookHistosSet("TrSec",hPattern); fHistosTrComb = BookHistosSet("TrComb",hPattern); fHistosTrCombU = BookHistosSet("TrCombU",hPattern); if (fCheckReconstructables) { fHistosTrRcblPrim = BookHistosSet("TrRcblPrim",hPattern); fHistosTrRcblSec = BookHistosSet("TrRcblSec",hPattern); } } //---------------------------------------------Standard histos per tracklet type---<< // //---------------------------------------------Custom Histos----------------------->> // put here any non standard histos fHistosCustom = BookCustomHistos(); // //---------------------------------------------Custom Histos-----------------------<< int nhist = fOutput->GetEntries(); for (int i=0;iAt(i); if (!hst || !(hst->InheritsFrom(TH1::Class()))) continue; ((TH1*)hst)->Sumw2(); } // if (fUseCentralityVar<0||fUseCentralityVar>kNCentTypes) { printf("Wrong centrality type %d\n",fUseCentralityVar); exit(1); } AliInfo(Form("Centrality type selected: %s\n",fgkCentSelName[fUseCentralityVar])); PostData(1, fOutput); // } //________________________________________________________________________ void AliTrackletTaskMultipp::UserExec(Option_t *) { // Main loop // Bool_t needRecPoints = GetDoNormalReco() || GetDoInjection() || GetDoRotation() || GetDoMixing(); // AliAnalysisManager* anMan = AliAnalysisManager::GetAnalysisManager(); fRPTree = fRPTreeMix = 0; AliESDInputHandler *handler = (AliESDInputHandler*)anMan->GetInputEventHandler(); AliESDInputHandlerRP *handRP = 0; if (needRecPoints) { handRP = (AliESDInputHandlerRP*)handler; if (!handRP) { printf("No RP handler\n"); return; } } AliESDEvent *esd = handler->GetEvent(); if (!esd) { printf("No AliESDEvent\n"); return; } // // do we need to initialize the field? AliMagF* field = (AliMagF*)TGeoGlobalMagField::Instance()->GetField(); if (!field && !esd->InitMagneticField()) {printf("Failed to initialize the B field\n");return;} // // TH1* hstat = (TH1*)fHistosCustom->UncheckedAt(kHStat); hstat->Fill(kEvTot0); // RS const AliESDVertex* vtxESD = esd->GetPrimaryVertexSPD(); const AliMultiplicity* multESD = esd->GetMultiplicity(); // if (vtxESD->GetNContributors()<1) return; TString vtxTyp = vtxESD->GetTitle(); if (vtxTyp.Contains("vertexer: Z")) { if (vtxESD->GetDispersion()>0.04) return; if (vtxESD->GetZRes()>0.25) return; } // AliCentrality *centrality = esd->GetCentrality(); // hstat->Fill(kEvTot); // RS // Double_t esdvtx[3]; vtxESD->GetXYZ(esdvtx); for (int i=3;i--;) fESDVtx[i] = esdvtx[i]; // float vtxf[3] = {vtxESD->GetX(),vtxESD->GetY(),vtxESD->GetZ()}; // // registed Ntracklets and ZVertex of the event ((TH1*)fHistosCustom->UncheckedAt(kHZVtxNoSel))->Fill(esdvtx[2]); // if(vtxf[2] < fZVertexMin || vtxf[2] > fZVertexMax) return; // // centrality->Dump(); Float_t multV0=0; AliESDVZERO* esdV0 = esd->GetVZEROData(); if (esdV0) { multV0 = esdV0->GetMTotV0A()+esdV0->GetMTotV0C(); if (fUseMC) multV0 *= fMCV0Scale; } ((TH1*)fHistosCustom->UncheckedAt(kHV0NoSel))->Fill(multV0); // float nSPD2 = multESD->GetNumberOfITSClusters(1); ((TH1*)fHistosCustom->UncheckedAt(kHNClSPD2NoSel))->Fill(nSPD2); // //------------------------------------------------------ AliESDZDC *esdZDC = esd->GetESDZDC(); float zdcEnergy=0,zemEnergy=0; if (esdZDC) { zdcEnergy = (esdZDC->GetZDCN1Energy() + esdZDC->GetZDCP1Energy() + esdZDC->GetZDCN2Energy()+ esdZDC->GetZDCP2Energy()); zemEnergy = (esdZDC->GetZDCEMEnergy(0)+esdZDC->GetZDCEMEnergy(1))/8.; } ((TH2*)fHistosCustom->UncheckedAt(kHZDCZEMNoSel))->Fill(zemEnergy,zdcEnergy); // Float_t centPercentile = centrality->GetCentralityPercentileUnchecked(fgkCentSelName[fUseCentralityVar]); // temporary >>>>>>>>>>>>>>>>>>>>>>>> if (fUseCentralityVar==kCentZEMvsZDC) { float zdcEn = zdcEnergy; float zemEn = zemEnergy; Float_t slope; Float_t zdcPercentile; if (zemEn > 295.) { slope = (zdcEn + 15000.)/(zemEn - 295.); slope += 2.23660e+02; zdcPercentile = (TMath::ATan(slope) - 1.56664)/8.99571e-05; if (zdcPercentile<0) zdcPercentile = 0; } else zdcPercentile = 100; centPercentile = zdcPercentile; } // temporary >>>>>>>>>>>>>>>>>>>>>>>> fCurrCentBin = GetCentralityBin(centPercentile); // // printf("CentPerc: %f : Bin %d\n",centPercentile, fCurrCentBin); if (fCurrCentBin<0) { //printf("Reject: %.1f : V0:%.1f V0Cor:%.1f V0CR:%.1f SPD2c:%.1f\n",mltTst, multV0,multV0Corr,multV0CorrResc,nSPD2Corr); return; } // ((TH1*)fHistosCustom->UncheckedAt(kHStatCentBin))->Fill(fCurrCentBin); ((TH1*)fHistosCustom->UncheckedAt(kHStatCent))->Fill(centPercentile); // AliMCEventHandler* eventHandler = 0; fMCevent = 0; fStack = 0; // if (fUseMC) { eventHandler = (AliMCEventHandler*)anMan->GetMCtruthEventHandler(); if (!eventHandler) { printf("ERROR: Could not retrieve MC event handler\n"); return; } fMCevent = eventHandler->MCEvent(); if (!fMCevent) { printf("ERROR: Could not retrieve MC event\n"); return; } fStack = fMCevent->Stack(); if (!fStack) { printf("Stack not available\n"); return; } } // if (needRecPoints) { fRPTree = handRP->GetTreeR("ITS"); if (!fRPTree) { AliError(" Invalid ITS cluster tree !\n"); return; } } // // =============================================================================>>> // MC Generator info AliGenEventHeader* mcGenH = 0; fNPart = 0; fNBColl = 0; if (fUseMC) { mcGenH = fMCevent->GenEventHeader(); if (mcGenH->InheritsFrom(AliGenHijingEventHeader::Class())) { AliGenHijingEventHeader* hHijing = (AliGenHijingEventHeader*)mcGenH; fNPart = (hHijing->ProjectileParticipants()+hHijing->TargetParticipants())/2.; fNBColl = hHijing->NN()+hHijing->NNw()+hHijing->NwN()+hHijing->NwNw(); } else if (mcGenH->InheritsFrom(AliGenDPMjetEventHeader::Class())) { AliGenDPMjetEventHeader* hDpmJet = (AliGenDPMjetEventHeader*)mcGenH; fNPart = (hDpmJet->ProjectileParticipants()+hDpmJet->TargetParticipants())/2.; fNBColl = hDpmJet->NN()+hDpmJet->NNw()+hDpmJet->NwN()+hDpmJet->NwNw(); } else {} // unknown generator } // // register Ntracklets and ZVertex of the event ((TH2*)fHistosCustom->UncheckedAt(kHZVtx))->Fill(esdvtx[2],fCurrCentBin); ((TH2*)fHistosCustom->UncheckedAt(kHV0))->Fill(multV0,fCurrCentBin); ((TH2*)fHistosCustom->UncheckedAt(kHNClSPD2))->Fill(nSPD2,fCurrCentBin); ((TH3*)fHistosCustom->UncheckedAt(kHZDCZEM))->Fill(zemEnergy,zdcEnergy,fCurrCentBin); // if (fUseMC) FillMCPrimaries(); // // normal reconstruction hstat->Fill(kBinEntries+kEvProcData + kEntriesPerBin*fCurrCentBin); // if (GetDoNormalReco() || GetDoInjection()) { // for the injection the normal reco should be done InitMultReco(); fMultReco->Run(fRPTree, vtxf); printf("Multiplicity Reconstructed:\n"); AliMultiplicity* mlt = fMultReco->GetMultiplicity(); if (mlt) mlt->Print(); if (GetDoNormalReco()) FillHistos(kData,mlt); FillClusterInfo(); // } if (!GetDoNormalReco()) { FillHistos(kData,multESD); // fill data histos from ESD FillClusterInfoFromMult(multESD, vtxf[2] ); FillClusterAutoCorrelationFromMult(multESD, vtxf[2]); } // // Injection: it must come right after the normal reco since needs its results if (GetDoInjection()) { if (!fMultReco) InitMultReco(); // in principle, not needed, the reco is created above fMultReco->SetRecType(AliITSMultRecBg::kBgInj); fMultReco->Run(fRPTree, vtxf); printf("Multiplicity from Injection:\n"); AliMultiplicity* mlt = fMultReco->GetMultiplicity(); if (mlt) mlt->Print(); hstat->Fill(kBinEntries + kEvProcInj + kEntriesPerBin*fCurrCentBin); FillHistos(kBgInj,mlt); } // // Rotation if (GetDoRotation()) { InitMultReco(); fMultReco->SetRecType(AliITSMultRecBg::kBgRot); fMultReco->SetPhiRotationAngle(fPhiRot); fMultReco->Run(fRPTree, vtxf); printf("Multiplicity from Rotation:\n"); AliMultiplicity* mlt = fMultReco->GetMultiplicity(); if (mlt) mlt->Print(); hstat->Fill(kBinEntries + kEvProcRot + kEntriesPerBin*fCurrCentBin); FillHistos(kBgRot,mlt); } // if (GetDoMixing()) { /* AliMixEventInputHandler* handToMix = (AliMixEventInputHandler*)handRP->MixingHandler(); if (!handToMix) { printf("No Mixing handler\n"); return; } handToMix->GetEntry(); if(handToMix->MixedEventNumber()<1) {printf("Mixing: No enough events in pool\n"); return;} AliESDInputHandlerRP* handRPMix = (AliESDInputHandlerRP*) handToMix->InputEventHandler(0); if (!handRPMix) { printf("No Mixing RP handler\n"); return; } fRPTreeMix = handRPMix->GetTreeR("ITS"); if (!fRPTreeMix) { AliError(" Invalid ITS cluster tree of the 2nd event!\n"); return; } // AliESDEvent *esdMix = handRPMix->GetEvent(); const AliESDVertex* vtxESDmix = esdMix->GetVertex(); ((TH2*)fHistosCustom->UncheckedAt(kHZVtxMixDiff))->Fill(vtxESDmix->GetZ()-esdvtx[2],fCurrCentBin); ((TH2*)fHistosCustom->UncheckedAt(kHNTrMixDiff) )-> Fill(esdMix->GetMultiplicity()->GetNumberOfTracklets() - multESD->GetNumberOfTracklets(),fCurrCentBin); // InitMultReco(); fMultReco->SetRecType(AliITSMultRecBg::kBgMix); fMultReco->Run(fRPTree, vtxf,fRPTreeMix); printf("Multiplicity from Mixing:\n"); AliMultiplicity* mlt = fMultReco->GetMultiplicity(); if (mlt) mlt->Print(); hstat->Fill(kBinEntries + kEvProcMix + kEntriesPerBin*fCurrCentBin); FillHistos(kBgMix,mlt); */ // } // =============================================================================<<< // if (fMultReco) delete fMultReco; fMultReco = 0; // } //________________________________________________________________________ void AliTrackletTaskMultipp::Terminate(Option_t *) { // finish processing Printf("Terminating..."); TH1* hstat; TList *lst = dynamic_cast(GetOutputData(1)); printf("Term: %p %p %p\n",fOutput,lst,fHistosCustom); if (lst && (hstat=(TH1*)lst->FindObject("hStat"))) { Info("Terminate","Registering used settings"); // fill used settings hstat->Fill(kOneUnit,1.); hstat->Fill(kCentVar,fUseCentralityVar); hstat->Fill(kDPhi,fDPhiWindow); hstat->Fill(kDTht,fDThetaWindow); hstat->Fill(kNStd,fNStdDev); hstat->Fill(kPhiShift,fDPhiShift); hstat->Fill(kThtS2,fScaleDTBySin2T); hstat->Fill(kThtCW,fCutOnDThetaX); hstat->Fill(kPhiOvl,fPhiOverlapCut); hstat->Fill(kZEtaOvl,fZetaOverlap); hstat->Fill(kNoOvl,fRemoveOverlaps); // hstat->Fill(kPhiRot,fPhiRot); hstat->Fill(kInjScl,fInjScale); hstat->Fill(kEtaMin,fEtaMin); hstat->Fill(kEtaMax,fEtaMax); hstat->Fill(kZVMin,fZVertexMin); hstat->Fill(kZVMax,fZVertexMax); // hstat->Fill(kDPiSCut,fDPhiSCut); hstat->Fill(kNStdCut,fNStdCut); hstat->Fill(kMCV0Scale, fMCV0Scale); // } // // AliAnalysisTaskSE::Terminate(); } //_________________________________________________________________________ void AliTrackletTaskMultipp::InitMultReco() { // create mult reconstructor if (fMultReco) delete fMultReco; fMultReco = new AliITSMultRecBg(); fMultReco->SetCreateClustersCopy(kTRUE); fMultReco->SetScaleDThetaBySin2T(fScaleDTBySin2T); fMultReco->SetNStdDev(fNStdDev); fMultReco->SetPhiWindow( fDPhiWindow ); fMultReco->SetThetaWindow( fDThetaWindow ); fMultReco->SetPhiShift( fDPhiShift ); fMultReco->SetRemoveClustersFromOverlaps(fRemoveOverlaps); fMultReco->SetPhiOverlapCut(fPhiOverlapCut); fMultReco->SetZetaOverlapCut(fZetaOverlap); fMultReco->SetHistOn(kFALSE); fMultReco->SetRecType( AliITSMultRecBg::kData ); } //_________________________________________________________________________ TObjArray* AliTrackletTaskMultipp::BookCustomHistos() { // book custom histos, not related to specific tracklet type TObjArray* histos = new TObjArray(); TH1F* hstat; // // ------------ job parameters, statistics ------------------------------>>> int nbs = kBinEntries + fNCentBins*kEntriesPerBin; hstat = new TH1F("hStat","Run statistics",nbs,0.5,nbs+0.5); // hstat->GetXaxis()->SetBinLabel(kEvTot, "Ev.Tot0"); hstat->GetXaxis()->SetBinLabel(kEvTot, "Ev.Tot"); hstat->GetXaxis()->SetBinLabel(kOneUnit,"ScaleMerge"); hstat->GetXaxis()->SetBinLabel(kNWorkers,"Workers"); // hstat->GetXaxis()->SetBinLabel(kCentVar, fgkCentSelName[fUseCentralityVar]); hstat->GetXaxis()->SetBinLabel(kDPhi, "#Delta#varphi"); hstat->GetXaxis()->SetBinLabel(kDTht, "#Delta#theta"); hstat->GetXaxis()->SetBinLabel(kNStd, "N.std"); hstat->GetXaxis()->SetBinLabel(kPhiShift,"#delta#varphi"); hstat->GetXaxis()->SetBinLabel(kThtS2,"scale #Delta#theta"); hstat->GetXaxis()->SetBinLabel(kPhiOvl,"#varpho_{Ovl}"); hstat->GetXaxis()->SetBinLabel(kZEtaOvl,"#z_{Ovl}"); hstat->GetXaxis()->SetBinLabel(kNoOvl, "rem.ovl"); // hstat->GetXaxis()->SetBinLabel(kPhiRot,"#varphi_{rot}"); hstat->GetXaxis()->SetBinLabel(kInjScl,"inj"); hstat->GetXaxis()->SetBinLabel(kEtaMin,"#eta_{min}"); hstat->GetXaxis()->SetBinLabel(kEtaMax,"#eta_{max}"); hstat->GetXaxis()->SetBinLabel(kZVMin,"ZV_{min} cut"); hstat->GetXaxis()->SetBinLabel(kZVMax,"ZV_{max} cut"); // hstat->GetXaxis()->SetBinLabel(kDPiSCut,"#Delta#varphi-#delta_{#phi} cut"); hstat->GetXaxis()->SetBinLabel(kNStdCut,"#Delta cut"); // hstat->GetXaxis()->SetBinLabel(kMCV0Scale,"MC V0 scale"); // for (int i=0;iGetXaxis()->SetBinLabel(offs + kEvProcData, bnt+" Ev.ProcData"); hstat->GetXaxis()->SetBinLabel(offs + kEvProcInj, bnt+" Ev.ProcInj"); hstat->GetXaxis()->SetBinLabel(offs + kEvProcRot, bnt+" Ev.ProcRot"); hstat->GetXaxis()->SetBinLabel(offs + kEvProcMix, bnt+" Ev.ProcMix"); // } // hstat->Fill(kNWorkers); // AddHisto(histos,hstat,kHStat); // // ------------------------ events per centrality bin ---------------------- TH1D* hCentAx = new TH1D("EvCentr","Events per centrality",fNCentBins,fgkCentPerc); hCentAx->GetXaxis()->SetTitle("Centrality parameter"); AddHisto(histos,hCentAx,kHStatCent); // TH1D* hCentBin = new TH1D("EvCentrBin","Events per centrality bin",fNCentBins,-0.5,fNCentBins-0.5); hCentBin->GetXaxis()->SetTitle("Centrality Bin"); AddHisto(histos,hCentBin,kHStatCentBin); // // ------------ job parameters, statistics ------------------------------<<< // double etaMn=-3,etaMx=3; double zMn=-30, zMx=30; int nEtaBins = int((etaMx-etaMn)/0.1); if (nEtaBins<1) nEtaBins = 1; // int nZVBins = int(zMx-zMn); if (nZVBins<1) nZVBins = 1; // // Z vertex distribution for events before selection TH1F* hzvns = new TH1F("zvNoSel","Z vertex before selection",nZVBins,zMn,zMx); hzvns->GetXaxis()->SetTitle("Zvertex"); AddHisto(histos,hzvns,kHZVtxNoSel); // // V0 for events before selection int nbmltV0 = 250; double maxmltV0 = 25000; // TH1F* hnV0ns = new TH1F("V0NoSel","V0 signal Before Cent. Selection",nbmltV0,0,maxmltV0); hnV0ns->GetXaxis()->SetTitle("V0 signal"); AddHisto(histos,hnV0ns,kHV0NoSel); // // N SPD2 clusters int nbmltSPD2 = 175; double maxmltSPD2 = 7000; TH1F* hncl2ns = new TH1F("NClustersSPD2NoSel","N Clusters on SPD2 Before Cent Selection",nbmltSPD2,0,maxmltSPD2); hncl2ns->GetXaxis()->SetTitle("N Clus SPD2"); AddHisto(histos,hncl2ns,kHNClSPD2NoSel); // int nbzdc=50; double maxZDC=6000., maxZEM=2500.; TH2F* hzdczemns = new TH2F("ZDCZEMNoSel","ZDC vs ZEM Before Cent Selection", nbzdc,0,maxZEM,nbzdc,0,maxZDC); hzdczemns->GetXaxis()->SetTitle("ZEM"); hzdczemns->GetXaxis()->SetTitle("ZDC"); AddHisto(histos,hzdczemns,kHZDCZEMNoSel); // TH2F* hzv = new TH2F("zv","Z vertex after Selection per Cent.Bin",nZVBins,zMn,zMx, fNCentBins, -0.5,fNCentBins-0.5); hzv->GetXaxis()->SetTitle("Zvertex"); hzv->GetYaxis()->SetTitle("Cent.Bin ID"); AddHisto(histos,hzv,kHZVtx); // // V0 TH2F* hnV0 = new TH2F("V0","V0 signal per Cent.Bin ",nbmltV0,0,maxmltV0, fNCentBins, -0.5,fNCentBins-0.5); hnV0->GetXaxis()->SetTitle("V0 signal"); hnV0->GetYaxis()->SetTitle("Cent.Bin ID"); AddHisto(histos,hnV0,kHV0); // // N SPD2 clusters TH2F* hncl2 = new TH2F("NClustersSPD2","N Clusters on SPD2 per Cent.Bin",nbmltSPD2,0,maxmltSPD2, fNCentBins, -0.5,fNCentBins-0.5); hncl2->GetXaxis()->SetTitle("N Clus SPD2"); hncl2->GetYaxis()->SetTitle("Cent.Bin ID"); AddHisto(histos,hncl2,kHNClSPD2); // // ZDCZEM TH3F* hzdczem = new TH3F("ZDCZEM","ZDC vs ZEM per Cent.Bin",nbzdc,0,maxZEM,nbzdc,0,maxZDC, fNCentBins, -0.5,fNCentBins-0.5); hzdczem->GetXaxis()->SetTitle("ZEM"); hzdczem->GetYaxis()->SetTitle("ZDC"); AddHisto(histos,hzdczem,kHZDCZEM); // //---------------------------------------------------------------------- int nEtaBinsS = int((fEtaMax-fEtaMin)/0.1); if (nEtaBinsS<1) nEtaBins = 1; // int nZVBinsS = int(fZVertexMax-fZVertexMin); if (nZVBinsS<1) nZVBinsS = 1; if (fUseMC) { // Z vertex vs Eta distribution for primaries char buffn[100],bufft[500]; for (int ib=0;ibGetXaxis()->SetTitle("#eta"); hzvetap->GetYaxis()->SetTitle("Zvertex"); AddHisto(histos,hzvetap,kHZVEtaPrimMC+ib); } // // primaries according to MC generator TH1F* hnprimM = new TH1F("nPrimMean"," primaries",fNCentBins, -0.5,fNCentBins-0.5); hnprimM->GetXaxis()->SetTitle("Cent.Bin ID"); AddHisto(histos,hnprimM,kHNPrimMeanMC); // // primaries per part.pair according to MC generator TH1F* hnprim2partM = new TH1F("nPrim2Part"," primaries per part.pair",fNCentBins, -0.5,fNCentBins-0.5); hnprim2partM->GetXaxis()->SetTitle("Cent.Bin ID"); AddHisto(histos,hnprim2partM,kHNPrim2PartMC); // // primaries per part.pair vs npart.pair according to MC generator TH2F* hnprim2partNp = new TH2F("nPrim2Part_vs_NPart"," primaries per part.pair vs N part.pairs",105,0,210,200,0,40); hnprim2partNp->GetXaxis()->SetTitle("N.part.pairs"); hnprim2partNp->GetYaxis()->SetTitle("N.prim/N.part.pairs"); AddHisto(histos,hnprim2partNp,kHNPrim2PartNpMC); // // primaries per b.coll vs npart.pair according to MC generator TH2F* hnprim2BCollNp = new TH2F("nPrim2BColl_vs_NPart"," primaries per bin.coll vs N part.pairs",105,0,210,200,0,40); hnprim2BCollNp->GetXaxis()->SetTitle("N.part.pairs"); hnprim2BCollNp->GetYaxis()->SetTitle("N.prim/N.bin.coll."); AddHisto(histos,hnprim2BCollNp,kHNPrim2BCollNpMC); // // primaries per bin.coll. according to MC generator TH1F* hnprim2BCollM = new TH1F("nPrim2BColl"," primaries per bin.coll",fNCentBins, -0.5,fNCentBins-0.5); hnprim2BCollM->GetXaxis()->SetTitle("Cent.Bin ID"); AddHisto(histos,hnprim2BCollM,kHNPrim2BCollMC); // // n participants according to MC generator TH2F* hnpart = new TH2F("nPart","N participant pairs",210,0,210,fNCentBins, -0.5,fNCentBins-0.5); hnpart->GetXaxis()->SetTitle("N part. pairs"); hnpart->GetYaxis()->SetTitle("Cent.Bin ID"); AddHisto(histos,hnpart,kHNPartMC); // // participants according to MC generator TH1F* hnpartM = new TH1F("nPartMean"," participant pairs",fNCentBins, -0.5,fNCentBins-0.5); hnpartM->GetXaxis()->SetTitle("Cent.Bin ID"); AddHisto(histos,hnpartM,kHNPartMeanMC); // // n bin coll. according to MC generator TH2F* hnbcoll = new TH2F("nBColl","N bin. coll",2000,0,2000,fNCentBins, -0.5,fNCentBins-0.5); hnbcoll->GetXaxis()->SetTitle("N bin. coll"); hnbcoll->GetYaxis()->SetTitle("Cent.Bin ID"); AddHisto(histos,hnbcoll,kHNBCollMC); // // bin col according to MC generator TH1F* hnbcollM = new TH1F("nBCollMean"," bin.colls",fNCentBins, -0.5,fNCentBins-0.5); hnbcollM->GetXaxis()->SetTitle("Cent.Bin ID"); AddHisto(histos,hnbcollM,kHNBCollMeanMC); // } // if (GetDoMixing()) { // // Difference in Z vertex for mixed events TH2F* hzdiff = new TH2F("MixSPDVertexDiff","SPD #Delta Z Vertex distribution per mult bin ",100,-5,5, fNCentBins, -0.5,fNCentBins-0.5); hzdiff->GetXaxis()->SetTitle("#Delta Z Vertex [cm]"); hzdiff->GetYaxis()->SetTitle(Form("Entries / %1.2f [cm] per mult bin",10./100.)); AddHisto(histos,hzdiff,kHZVtxMixDiff); // // Difference in N tracklets for mixed events TH2F* hntdiff = new TH2F("MixNTrackletsDiff"," SPD tracklets Diff ",200,-1000,1000, fNCentBins, -0.5,fNCentBins-0.5); hntdiff->GetXaxis()->SetTitle("# tracklet diff"); AddHisto(histos,hntdiff,kHNTrMixDiff); } // // -------------------------------------------------- if (fUseMC) { int npdg = sizeof(fgkPDGNames)/sizeof(char*); TH2F* hpdgP = new TH2F("pdgPrim","primary PDG",npdg,0,npdg,fNCentBins, -0.5,fNCentBins-0.5); AddHisto(histos,hpdgP,kHPrimPDG); TH2F* hpdgS = new TH2F("pdgSec","secondary PDG",npdg,0,npdg,fNCentBins, -0.5,fNCentBins-0.5); AddHisto(histos,hpdgS,kHSecPDG); TH2F* hpdgPP = new TH2F("pdgPrimPar","primary parent PDG ",npdg,0,npdg,fNCentBins, -0.5,fNCentBins-0.5); AddHisto(histos,hpdgPP,kHPrimParPDG); TH2F* hpdgSP = new TH2F("pdgSecPar","secondary parent PDG",npdg,0,npdg,fNCentBins, -0.5,fNCentBins-0.5); AddHisto(histos,hpdgSP,kHSecParPDG); for (int i=0;iGetXaxis()->SetBinLabel(i+1,fgkPDGNames[i]); hpdgS->GetXaxis()->SetBinLabel(i+1,fgkPDGNames[i]); hpdgPP->GetXaxis()->SetBinLabel(i+1,fgkPDGNames[i]); hpdgSP->GetXaxis()->SetBinLabel(i+1,fgkPDGNames[i]); } } // // ------------------------------------------------- TH2F* hclinf=0; hclinf = new TH2F("cl0InfoUsed","#phi vs Z of used clusters, Lr0",60,-15,15, 80,0,2*TMath::Pi()); AddHisto(histos,hclinf,kHClUsedInfoL0); hclinf = new TH2F("cl1InfoUsed","#phi vs Z of used clusters, Lr1",60,-15,15, 2*80,0,2*TMath::Pi()); AddHisto(histos,hclinf,kHClUsedInfoL1); hclinf = new TH2F("cl0InfoAll","#phi vs Z of all clusters, Lr0",60,-15,15, 80,0,2*TMath::Pi()); AddHisto(histos,hclinf,kHClAllInfoL0); hclinf = new TH2F("cl1InfoAll","#phi vs Z of all clusters, Lr1",60,-15,15, 2*80,0,2*TMath::Pi()); AddHisto(histos,hclinf,kHClAllInfoL1); // // ------------------------------------------------- // correlations between SPD1 clusters hclinf = new TH2F("clDstVsZall","Distance between SPD1 clusters vs Z, all cl",30,-15,15, 50,0.,5.); AddHisto(histos,hclinf,kHclDstZAll); hclinf = new TH2F("clDstVsPhiall","Distance between SPD1 clusters vs #phi, all cl",30,0,2*TMath::Pi(),50,0.,5.); AddHisto(histos,hclinf,kHclDstPhiAll); // hclinf = new TH2F("clDstVsZused","Distance between SPD1 clusters vs Z, used-unused cl",30,-15,15, 50,0.,5.); AddHisto(histos,hclinf,kHclDstZUsed); hclinf = new TH2F("clDstVsPhiused","Distance between SPD1 clusters vs #phi, used-unused cl",30,0,2*TMath::Pi(),50,0.,5.); AddHisto(histos,hclinf,kHclDstPhiUsed); // histos->SetOwner(kFALSE); // return histos; } //_________________________________________________________________________ TObjArray* AliTrackletTaskMultipp::BookHistosSet(const char* pref, UInt_t selHistos) { // book standard set of histos attaching the pref in front of the name/title // const int kNDPhiBins = 100; const int kNDThtBins = 100; int nDistBins = int(fNStdDev)*5; // int nEtaBins = int((fEtaMax-fEtaMin)/0.1); if (nEtaBins<1) nEtaBins = 1; // int nZVBins = int(fZVertexMax-fZVertexMin); if (nZVBins<1) nZVBins = 1; float dphir = fDPhiWindow*TMath::Sqrt(fNStdDev); float dthtr = fDThetaWindow*TMath::Sqrt(fNStdDev); // TObjArray* histos = new TObjArray(); TH2F* h2; TH1F* h1; char buffn[100],bufft[500]; // for (int ib=0;ibGetXaxis()->SetTitle("#eta"); h2->GetYaxis()->SetTitle("Zv"); AddHisto(histos,h2,offs+kHEtaZvCut); } // if (selHistos & (0x1<GetXaxis()->SetTitle("#Delta#varphi [rad]"); h2->GetYaxis()->SetTitle("#Delta#theta [rad]"); AddHisto(histos,h2,offs+kHDPhiDTheta); } // if (selHistos & (0x1<GetXaxis()->SetTitle("#Delta#varphi-#delta_{#varphi} [rad]"); sprintf(bufft,"#Delta#theta%s",fScaleDTBySin2T ? "/sin^{2}(#theta)":""); h2->GetYaxis()->SetTitle(bufft); AddHisto(histos,h2,offs+kHDPhiSDThetaX); } // if (selHistos & (0x1<GetXaxis()->SetTitle(bufft); AddHisto(histos,h1,offs+kHWDist); } // if (selHistos & (0x1<GetXaxis()->SetTitle("#eta"); h2->GetYaxis()->SetTitle("Zv"); AddHisto(histos,h2,offs+kHEtaZvSPD1); } // } // histos->SetOwner(kFALSE); return histos; } //_________________________________________________________________________ void AliTrackletTaskMultipp::AddHisto(TObjArray* histos, TObject* h, Int_t at) { // add single histo to the set if (at>=0) histos->AddAtAndExpand(h,at); else histos->Add(h); fOutput->Add(h); } //_________________________________________________________________________ void AliTrackletTaskMultipp::FillHistos(Int_t type, const AliMultiplicity* mlt) { // fill histos of given type if (!mlt) return; // TObjArray* histos = 0; if (type == kData) histos = fHistosTrData; else if (type == kBgInj) histos = fHistosTrInj; else if (type == kBgRot) histos = fHistosTrRot; else if (type == kBgMix) histos = fHistosTrMix; // Bool_t fillMC = (type==kData) && fUseMC && fStack; // // //---------------------------------------- CHECK ------------------------------>>> TArrayF vtxMC; AliGenHijingEventHeader* pyHeader = 0; // if (fUseMC) { pyHeader = (AliGenHijingEventHeader*) fMCevent->GenEventHeader();//header->GenEventHeader(); pyHeader->PrimaryVertex(vtxMC); } //---------------------------------------- CHECK ------------------------------<<< // if (!histos) return; int ntr = mlt->GetNumberOfTracklets(); for (int itr=ntr;itr--;) { // //---------------------------------------- CHECK ------------------------------>>> /* if (fUseMC) { Bool_t reject = kFALSE; while(1) { int lab0 = mlt->GetLabel(itr,0); int lab1 = mlt->GetLabel(itr,1); if (lab0!=lab1) break; if (!fStack->IsPhysicalPrimary(lab0)) break; // TParticle* part = fStack->Particle(lab0); Float_t dz = part->Vz() - vtxMC[2]; if (TMath::Abs(dz)<1e-6) break; reject = kTRUE; break; } if (reject) continue; } */ //---------------------------------------- CHECK ------------------------------<<< // double theta = mlt->GetTheta(itr); double eta = -TMath::Log(TMath::Tan(theta/2)); if (etafEtaMax) continue; // double dtheta = mlt->GetDeltaTheta(itr); double dThetaX = dtheta; if (fScaleDTBySin2T) { double sint = TMath::Sin(theta); dThetaX /= (sint*sint); } if (fCutOnDThetaX && TMath::Abs(dThetaX)>fDThetaWindow) continue; // double phi = mlt->GetPhi(itr); double dphi = mlt->GetDeltaPhi(itr); double dist = mlt->CalcDist(itr); // FillHistosSet(histos,eta,/*phi,theta,*/dphi,dtheta,dThetaX,dist); // special handling for mc info if (fillMC && fStack) { int lab0 = mlt->GetLabel(itr,0); int lab1 = mlt->GetLabel(itr,1); int typeMC = 2; // comb.bg. if (lab0 == lab1) typeMC = fStack->IsPhysicalPrimary(lab0) ? 0:1; // prim or sec if (typeMC==0) FillHistosSet(fHistosTrPrim,eta,/*phi,theta,*/dphi,dtheta,dThetaX,dist); // primary else if (typeMC==1) FillHistosSet(fHistosTrSec, eta,/*phi,theta,*/dphi,dtheta,dThetaX,dist); // secondary else { FillHistosSet(fHistosTrComb,eta,/*phi,theta,*/dphi,dtheta,dThetaX,dist); // comb // for combinatorals fill also the uncorrelated part if (fMultReco) { float *trl = fMultReco->GetTracklet(itr); int clId0 = (int)trl[AliITSMultReconstructor::kClID1]; int clId1 = (int)trl[AliITSMultReconstructor::kClID2]; float *clLabs0 = fMultReco->GetClusterOfLayer(0,clId0) + AliITSMultReconstructor::kClMC0; float *clLabs1 = fMultReco->GetClusterOfLayer(1,clId1) + AliITSMultReconstructor::kClMC0; if (!HaveCommonParent(clLabs0,clLabs1)) FillHistosSet(fHistosTrCombU,eta,/*phi,theta,*/dphi,dtheta,dThetaX,dist); } } // combinatorials if (dist>> int offsH = fCurrCentBin*kNStandardH; TH2* hSingles = (TH2*)histos->UncheckedAt(offsH+kHEtaZvSPD1); if (hSingles) { int nclS = mlt->GetNumberOfSingleClusters(); double *thtS = mlt->GetThetaSingle(); for (int ics=nclS;ics--;) { double etaS = -TMath::Log(TMath::Tan(thtS[ics]/2)); if (etaSfEtaMax) continue; hSingles->Fill(etaS,fESDVtx[2]); } } //-------------------------------------------------------------TMP RS - singles -------<<< // } //_________________________________________________________________________ void AliTrackletTaskMultipp::FillMCPrimaries() { // fill all MC primaries Zv vs Eta if (!fStack || !fMCevent) return; //---------------------------------------- CHECK ------------------------------>>> TArrayF vtxMC; AliGenHijingEventHeader* pyHeader = 0; // if (fUseMC) { pyHeader = (AliGenHijingEventHeader*) fMCevent->GenEventHeader();//header->GenEventHeader(); pyHeader->PrimaryVertex(vtxMC); } //---------------------------------------- CHECK ------------------------------<<< // int ntr = fStack->GetNtrack(); TH2* hprimEtaZ = (TH2F*)fHistosCustom->UncheckedAt(kHZVEtaPrimMC+fCurrCentBin); int nprim = 0; for (int itr=ntr;itr--;) { if (!fStack->IsPhysicalPrimary(itr)) continue; AliMCParticle *part = (AliMCParticle*)fMCevent->GetTrack(itr); if (!part->Charge()) continue; // //---------------------------------------- CHECK ------------------------------>>> /* Float_t dz = part->Zv() - vtxMC[2]; if (TMath::Abs(dz)>1e-6) continue; // reject */ //---------------------------------------- CHECK ------------------------------<<< // Float_t theta = part->Theta(); if (theta<1e-6 || theta>TMath::Pi()-1e-6) continue; Float_t eta = part->Eta(); if (etafEtaMax) continue; hprimEtaZ->Fill(eta, fESDVtx[2]); nprim++; } // ((TH1*)fHistosCustom->UncheckedAt(kHNPrimMeanMC))->Fill(fCurrCentBin,nprim); if (fNPart>0) { ((TH1*)fHistosCustom->UncheckedAt(kHNPrim2PartMC))->Fill(fCurrCentBin,nprim/fNPart); ((TH2*)fHistosCustom->UncheckedAt(kHNPrim2PartNpMC))->Fill(fNPart,nprim/fNPart); ((TH2*)fHistosCustom->UncheckedAt(kHNPartMC))->Fill(fNPart,fCurrCentBin); ((TH1*)fHistosCustom->UncheckedAt(kHNPartMeanMC))->Fill(fCurrCentBin,fNPart); } if (fNBColl>0) { ((TH1*)fHistosCustom->UncheckedAt(kHNPrim2BCollMC))->Fill(fCurrCentBin,nprim/fNBColl); ((TH2*)fHistosCustom->UncheckedAt(kHNPrim2BCollNpMC))->Fill(fNPart,nprim/fNBColl); ((TH2*)fHistosCustom->UncheckedAt(kHNBCollMC))->Fill(fNBColl,fCurrCentBin); ((TH1*)fHistosCustom->UncheckedAt(kHNBCollMeanMC))->Fill(fCurrCentBin,fNBColl); } // } //_________________________________________________________________________ void AliTrackletTaskMultipp::FillHistosSet(TObjArray* histos, double eta, //double /*phi*/,double /*theta*/, double dphi,double dtheta,double dThetaX, double dist) { // fill standard set of histos if (dist>fNStdDev) return; // double dphiS = TMath::Abs(dphi) - fDPhiShift; if (dphi<0) dphiS = -dphiS; // int offs = fCurrCentBin*kNStandardH; // if (histos->UncheckedAt(offs+kHDPhiSDThetaX)) ((TH2*)histos->UncheckedAt(offs+kHDPhiSDThetaX))->Fill(dphiS,dThetaX); // if (histos->UncheckedAt(offs+kHDPhiDTheta)) ((TH2*)histos->UncheckedAt(offs+kHDPhiDTheta))->Fill(dphi,dtheta); // if (histos->UncheckedAt(kHWDist)) ((TH2*)histos->UncheckedAt(offs+kHWDist))->Fill(dist); // if (distUncheckedAt(offs+kHEtaZvCut)) ((TH2*)histos->UncheckedAt(offs+kHEtaZvCut))->Fill(eta,fESDVtx[2]); // } //__________________________________________________________________ void AliTrackletTaskMultipp::FillSpecies(Int_t primsec, Int_t id) { // fill PDGcode TH1 *hPart=0,*hParent=0; if (primsec==0) { hPart = (TH1*)fHistosCustom->UncheckedAt(kHPrimPDG); hParent = (TH1*)fHistosCustom->UncheckedAt(kHPrimParPDG); } else if (primsec==1) { hPart = (TH1*)fHistosCustom->UncheckedAt(kHSecPDG); hParent = (TH1*)fHistosCustom->UncheckedAt(kHSecParPDG); } else return; int ntr = fStack->GetNtrack(); TParticle* part = fStack->Particle(id); int pdgCode = TMath::Abs(part->GetPdgCode()); int pdgBin = GetPdgBin(pdgCode); int parID = part->GetFirstMother(); int pdgCodePar = -1; int pdgBinPar = -1; while (parID>=0 && parIDParticle(parID); pdgCodePar = TMath::Abs(part->GetPdgCode()); parID = part->GetFirstMother(); } if (pdgCodePar>0) pdgBinPar = GetPdgBin(pdgCodePar); // hPart->Fill(pdgBin,fCurrCentBin); hParent->Fill(pdgBinPar,fCurrCentBin); // } //_________________________________________________________________________ Int_t AliTrackletTaskMultipp::GetCentralityBin(Float_t perc) const { // calculate centrality bin for (int i=0;i=fgkCentPerc[i] && perc<=fgkCentPerc[i+1]) return i; return -1; } //_________________________________________________________________________ Int_t AliTrackletTaskMultipp::GetPdgBin(Int_t pdgCode) { // return my pdg bin int ncodes = sizeof(fgkPDGCodes)/sizeof(int); int pdgBin=0; for (pdgBin=0;pdgBin=ncodes) { if (float(pdgCode)>1e9) pdgBin = ncodes; // nuclei else pdgBin = ncodes+1; // unknown } return pdgBin; } //_________________________________________________________________________ Bool_t AliTrackletTaskMultipp::HaveCommonParent(const float* clLabs0,const float* clLabs1) { // do 2 clusters have common parrent const int kMaxPar = 50; static int pars[2][50]; int npars[2]={0,0}; const float *labs[2] = {clLabs0,clLabs1}; int ntr = fStack->GetNtrack(); // printf("\nHaveCommonParent \n"); for (int il=0;il<2;il++) { for (int ilb=0;ilb<3;ilb++) { int lbl = (int)labs[il][ilb]; if (lbl<2 || lbl>=ntr) continue; // while (npars[il]Particle(lbl); if (!part) break; int code = TMath::Abs(part->GetPdgCode()); int q = (int)TMath::Abs(part->GetPDG()->Charge()); if (code==21 || code<10 || q==1 || q==2 || q==4 ) break; //printf("%d/%d/%d: %4d (%d)%s|",il,ilb,npars[il],lbl,part->GetStatusCode(),part->GetName()); pars[il][ npars[il]++ ] = lbl; lbl = part->GetFirstMother(); if (lbl<1 || lbl>=ntr) break; } // printf("\n"); } } // compare array of parents for (int i0=npars[0];i0--;) for (int i1=npars[1];i1--;) if (pars[0][i0]==pars[1][i1]) return kTRUE; return kFALSE; } //_________________________________________________________________________ void AliTrackletTaskMultipp::CheckReconstructables() { // fill reconstructable tracklets hitsos static TArrayI trInd; static TBits isPrimArr; // if (!fMultReco || !fMultReco->IsRecoDone()) {AliInfo("To check reconstructables the reco had to be requested"); return;} if (!fStack) {AliInfo("MC Stack is not availalble"); return;} const double kPtMin = 0.05; // TClonesArray *clArr[2]; for (int ilr=0;ilr<2;ilr++) { clArr[ilr] = fMultReco->GetClustersOfLayer(ilr); if (!clArr[ilr]) {AliInfo("Clusters are not available"); return;} } // int ntr = fStack->GetNtrack(); if (!ntr) return; trInd.Reset(); if (trInd.GetSize()Particle(it); if (TMath::Abs(part->Eta())>2.2) continue; if (TMath::Abs(part->Pt())IsPhysicalPrimary(it)) { isPrimArr.SetBitNumber(it); ntrStorePrim++; } else { // check if secondary is worth cheking TParticlePDG* pdgPart = part->GetPDG(); if (TMath::Abs(pdgPart->Charge())!=3) continue; if (part->R()>5.) continue; } trIndArr[it] = ++ntrStore; } // AliInfo(Form("Selected %d MC particles (%d prim) out of %d in the stack\n",ntrStore,ntrStorePrim,ntr)); // const int kMaxCl=3; AliITSRecPoint **clIndL[2]; clIndL[0] = new AliITSRecPoint*[kMaxCl*ntrStore]; // max 2 clusters per layer clIndL[1] = new AliITSRecPoint*[kMaxCl*ntrStore]; // max 2 clusters per layer memset(clIndL[0],0,kMaxCl*ntrStore*sizeof(AliITSRecPoint*)); memset(clIndL[1],0,kMaxCl*ntrStore*sizeof(AliITSRecPoint*)); // for (int ilr=0;ilr<2;ilr++) { TClonesArray *clusters = clArr[ilr]; int ncl = clusters->GetEntriesFast(); for (int icl=ncl;icl--;) { AliITSRecPoint *cl = (AliITSRecPoint*)clusters->UncheckedAt(icl); for (int ilb=3;ilb--;) { int lbl = cl->GetLabel(ilb); if (lbl<0 || lbl>=ntr) continue; int lblI = trIndArr[lbl]; if (--lblI<0) continue; // not kept for (int icc=0;iccGetGlobalXYZ( clusterLay[0] ); fMultReco->ClusterPos2Angles(clusterLay[0], fESDVtx); for (int icl1=0;icl1GetGlobalXYZ( clusterLay[1] ); fMultReco->ClusterPos2Angles(clusterLay[1], fESDVtx); trComp[AliITSMultReconstructor::kTrPhi][ntrCand] = clusterLay[0][AliITSMultReconstructor::kClPh]; trComp[AliITSMultReconstructor::kTrTheta][ntrCand] = clusterLay[0][AliITSMultReconstructor::kClTh]; trComp[AliITSMultReconstructor::kTrDTheta][ntrCand] = clusterLay[0][AliITSMultReconstructor::kClTh] - clusterLay[1][AliITSMultReconstructor::kClTh]; trComp[AliITSMultReconstructor::kTrDPhi][ntrCand] = clusterLay[0][AliITSMultReconstructor::kClPh] - clusterLay[1][AliITSMultReconstructor::kClPh]; trComp[AliITSMultReconstructor::kTrLab1][ntrCand] = icl1*10 + icl0; double &dphi = trComp[ntrCand][3]; if (dphi>TMath::Pi()) dphi=2.*TMath::Pi()-dphi; // take into account boundary condition trComp[5][ntrCand] = fMultReco->CalcDist(trComp[AliITSMultReconstructor::kTrDPhi][ntrCand], trComp[AliITSMultReconstructor::kTrDTheta][ntrCand], trComp[AliITSMultReconstructor::kTrTheta][ntrCand]); ntrCand++; } } if (!ntrCand) continue; // no tracklets if (ntrCand>1) TMath::Sort(ntrCand,trComp[5],indQual,kFALSE); else indQual[0] = 0; // sort in weighted distance if (fRemoveOverlaps) ntrCand = 1; // select the best // // disable worst tracklet with shared cluster for (int itc=0;itcfEtaMax) continue; double dThetaX = trComp[AliITSMultReconstructor::kTrTheta][ind]; if (fScaleDTBySin2T) { double sint = TMath::Sin(trComp[AliITSMultReconstructor::kTrTheta][ind]); dThetaX /= (sint*sint); } FillHistosSet(histos,eta, //trComp[AliITSMultReconstructor::kTrPhi][ind],trComp[AliITSMultReconstructor::kTrTheta][ind], trComp[AliITSMultReconstructor::kTrDPhi][ind],trComp[AliITSMultReconstructor::kTrDTheta][ind], dThetaX,trComp[5][ind]); } } // delete[] clIndL[0]; delete[] clIndL[1]; } //_________________________________________________________________________ void AliTrackletTaskMultipp::FillClusterInfo() { // fill info on clusters associated to good tracklets if (!fMultReco) return; int ntr = fMultReco->GetNTracklets(); int clID[2]; TH2F *hclU[2] = {(TH2F*)fHistosCustom->UncheckedAt(kHClUsedInfoL0),(TH2F*)fHistosCustom->UncheckedAt(kHClUsedInfoL1)}; TH2F *hclA[2] = {(TH2F*)fHistosCustom->UncheckedAt(kHClAllInfoL0),(TH2F*)fHistosCustom->UncheckedAt(kHClAllInfoL1)}; for (int itr=ntr;itr--;) { Float_t *trc = fMultReco->GetTracklet(itr); if (TMath::Abs(TMath::Abs(trc[AliITSMultReconstructor::kTrDPhi])-fDPhiShift)>fDPhiSCut) continue; if (fMultReco->CalcDist(trc[AliITSMultReconstructor::kTrDPhi], trc[AliITSMultReconstructor::kTrDTheta], trc[AliITSMultReconstructor::kTrTheta]) > fNStdCut) continue; clID[0] = (int)trc[AliITSMultReconstructor::kClID1]; clID[1] = (int)trc[AliITSMultReconstructor::kClID2]; for (int il=0;il<2;il++) { Float_t *clinf = fMultReco->GetClusterOfLayer(il,clID[il]); hclU[il]->Fill( clinf[AliITSMultReconstructor::kClZ], clinf[AliITSMultReconstructor::kClPh]); } } // for (int il=0;il<2;il++) for (int ic=fMultReco->GetNClustersLayer(il);ic--;) { Float_t *clinf = fMultReco->GetClusterOfLayer(il,ic); hclA[il]->Fill( clinf[AliITSMultReconstructor::kClZ], clinf[AliITSMultReconstructor::kClPh]); } // } //_________________________________________________________________________ void AliTrackletTaskMultipp::FillClusterInfoFromMult(const AliMultiplicity* mlt, double zVertex) { // fill info on clusters taking them from Multiplicity object const double kRSPD1 = 3.9; TH2F *hclU = (TH2F*)fHistosCustom->UncheckedAt(kHClUsedInfoL0); TH2F *hclA = (TH2F*)fHistosCustom->UncheckedAt(kHClAllInfoL0); int ntr = mlt->GetNumberOfTracklets(); for (int itr=ntr;itr--;) { Bool_t goodTracklet = kTRUE; if (TMath::Abs( mlt->GetDeltaPhi(itr)-fDPhiShift)>fDPhiSCut) goodTracklet = kFALSE; if (mlt->CalcDist(itr) > fNStdCut) goodTracklet = kFALSE; double phi = mlt->GetPhi(itr); double z = kRSPD1/TMath::Tan(mlt->GetTheta(itr)) + zVertex; if (goodTracklet) hclU->Fill(z,phi); hclA->Fill(z,phi); } // int ncl = mlt->GetNumberOfSingleClusters(); for (int icl=ncl;icl--;) { double phi = mlt->GetPhiSingle(icl); double z = kRSPD1/TMath::Tan(mlt->GetThetaSingle(icl)) + zVertex; hclA->Fill(z,phi); } // } //_________________________________________________________________________ void AliTrackletTaskMultipp::FillClusterAutoCorrelationFromMult(const AliMultiplicity* mlt, double zVertex) { // fill mutual distance between SPD1 clusters const double kRSPD1 = 3.9; enum {kX=0,kY,kZ,kPhi,kNV}; // TH2F *hclDstZAll = (TH2F*)fHistosCustom->UncheckedAt(kHclDstZAll); TH2F *hclDstPhiAll = (TH2F*)fHistosCustom->UncheckedAt(kHclDstPhiAll); TH2F *hclDstZUsed = (TH2F*)fHistosCustom->UncheckedAt(kHclDstZUsed); TH2F *hclDstPhiUsed = (TH2F*)fHistosCustom->UncheckedAt(kHclDstPhiUsed); if (!hclDstZAll && !hclDstPhiAll && !hclDstZUsed && !hclDstPhiUsed) return; // nothing to fill int nCl = mlt->GetNumberOfTracklets() + mlt->GetNumberOfSingleClusters(); if (nCl<2) return; double *clXYZPhi = new Double_t[kNV*nCl]; int cnt = 0; for (int itr=mlt->GetNumberOfTracklets();itr--;) { // extract coordinates of used SPD1 clusters double phi = mlt->GetPhi(itr); if (phi<0) phi += 2*TMath::Pi(); int offs = cnt*kNV; clXYZPhi[offs + kX] = TMath::Cos(phi)*kRSPD1; clXYZPhi[offs + kY] = TMath::Sin(phi)*kRSPD1; clXYZPhi[offs + kZ] = kRSPD1/TMath::Tan(mlt->GetTheta(itr)) + zVertex; clXYZPhi[offs + kPhi] = phi; cnt++; } // for (int icl=mlt->GetNumberOfSingleClusters();icl--;) { // extract coordinates of unused SPD1 clusters double phi = mlt->GetPhiSingle(icl); if (phi<0) phi += 2*TMath::Pi(); int offs = cnt*kNV; clXYZPhi[offs + kX] = TMath::Cos(phi)*kRSPD1; clXYZPhi[offs + kY] = TMath::Sin(phi)*kRSPD1; clXYZPhi[offs + kZ] = kRSPD1/TMath::Tan(mlt->GetThetaSingle(icl)) + zVertex; clXYZPhi[offs + kPhi] = phi; cnt++; } // for (int icl=0;iclFill(clXYZPhi[offs+kZ],-1,2); // fill in underflow, to count the clusters if (hclDstPhiAll) hclDstPhiAll->Fill(clXYZPhi[offs+kPhi],-1,2); // fill in underflow, to count the clusters if (iclGetNumberOfTracklets()) { if (hclDstZUsed) hclDstZUsed->Fill(clXYZPhi[offs+kZ],-1); // fill in underflow, to count the clusters if (hclDstPhiUsed) hclDstPhiUsed->Fill(clXYZPhi[offs+kPhi],-1); // fill in underflow, to count the clusters } for (int icl1=icl+1;icl1Fill(clXYZPhi[offs+kZ],dst); hclDstZAll->Fill(clXYZPhi[offs1+kZ],dst); } if (hclDstPhiAll) { hclDstPhiAll->Fill(clXYZPhi[offs+kPhi],dst); hclDstPhiAll->Fill(clXYZPhi[offs1+kPhi],dst); } // // tracklets with singles only if (icl>=mlt->GetNumberOfTracklets() || icl1GetNumberOfTracklets()) continue; if (hclDstZUsed) hclDstZUsed->Fill(clXYZPhi[offs+kZ],dst); if (hclDstPhiUsed) hclDstPhiUsed->Fill(clXYZPhi[offs+kPhi],dst); // } } }