PWG4/JetCorrel/AliJetCorrelWriter.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15 /* $Id: $ */
  16
  17 //______________________________________________________
  18 // Class for output (histograms) definition and filling.
  19 // Contains also the methods for calculation of correlation parameters.
  20 //-- Author: Paul Constantin
  21
  22 #include "AliJetCorrelWriter.h"
  23
  24 using namespace std;
  25
  26 ClassImp(AliJetCorrelWriter)
  27
  28 AliJetCorrelWriter::AliJetCorrelWriter() :
  29   fSelector(NULL), fMaker(NULL), fHname(""), fHtit(""), fRecoTrigg(kFALSE), fRndm(6251093), fNtuParent(NULL),
  30   fHBinsCentr(NULL), fHBinsZVert(NULL), fHBinsTrigg(NULL), fHBinsAssoc(NULL), fHCentr(NULL), fHZVert(NULL) {
  31   // constructor
  32   for(UInt_t i=0; i<2; i++){
  33     fHTrkITSQA[i] = NULL; fHTrkTPCQA[i] = NULL; fHTrkVTXQA[i] = NULL;
  34     for(UInt_t ic=0; ic<AliJetCorrelSelector::kMaxCent; ic++)  fHTrkProx[i][ic] = NULL;
  35   }
  36   for(UInt_t ik=0; ik<AliJetCorrelSelector::kMaxCorrel; ik++){
  37     fHTriggAcc[ik] = NULL; fHAssocAcc[ik] = NULL;
  38     for(UInt_t ic=0; ic<AliJetCorrelSelector::kMaxCent; ic++){
  39       fHTriggPt[ik][ic] = NULL; fHAssocPt[ik][ic] = NULL;
  40       for(UInt_t iv=0; iv<AliJetCorrelSelector::kMaxVert; iv++)
  41         for(UInt_t it=0; it<AliJetCorrelSelector::kMaxTrig; it++)
  42           for(UInt_t ia=0; ia<AliJetCorrelSelector::kMaxAsso; ia++){
  43             fHReal[ik][ic][iv][it][ia] = NULL;
  44             fHMix[ik][ic][iv][it][ia]  = NULL;
  45           }
  46     }
  47   }
  48 }
  49
  50 AliJetCorrelWriter::~AliJetCorrelWriter(){
  51   // destructor
  52 }
  53
  54 void AliJetCorrelWriter::Init(AliJetCorrelSelector * const s, AliJetCorrelMaker * const m){
  55   // initialization method
  56   fSelector = s;
  57   fMaker = m;
  58   fRecoTrigg = m->RecoTrigger();
  59 }
  60
  61 ////////////////////////////////////////////
  62 // METHODS FOR CREATION OF OUTPUT HISTOGRAMS
  63 ////////////////////////////////////////////
  64
  65 void AliJetCorrelWriter::CreateGeneric(TList *histosContainer){
  66   // books generic histograms
  67   UInt_t nTypeTrigg = fMaker->NoOfTrigg();
  68   UInt_t nTypeAssoc = fMaker->NoOfAssoc();
  69   UInt_t nBinsCentr = fSelector->NoOfBins(centr);
  70   UInt_t nBinsZVert = fSelector->NoOfBins(zvert);
  71   UInt_t nBinsTrigg = fSelector->NoOfBins(trigg);
  72   UInt_t nBinsAssoc = fSelector->NoOfBins(assoc);
  73   Float_t minTrigg  = fSelector->MinLowBin(trigg);
  74   Float_t maxTrigg  = fSelector->MaxHighBin(trigg);
  75   Float_t minAssoc  = fSelector->MinLowBin(assoc);
  76   Float_t maxAssoc  = fSelector->MaxHighBin(assoc);
  77
  78   for(UInt_t ic=0; ic<nBinsCentr; ic++){ // loop over centrality bins
  79     for(UInt_t tt=0; tt<nTypeTrigg; tt++){ // loop over trigger types
  80       fHtit="type:"; fHtit+=tt; fHtit+=" cent:"; fHtit+=ic;
  81       fHname = "fHTriggPt"; fHname+=tt; fHname+=ic;
  82       fHTriggPt[tt][ic] = new TH1F(fHname,  fHtit, 10*nBinsTrigg, minTrigg, maxTrigg); // for <pT> each bin
  83       histosContainer->AddLast(fHTriggPt[tt][ic]);
  84     } // loop over trigger types
  85     for(UInt_t at=0; at<nTypeAssoc; at++){ // loop over trigger types
  86       fHtit="type:"; fHtit+=at; fHtit+=" cent:"; fHtit+=ic;
  87       fHname = "fHAssocPt"; fHname+=at; fHname+=ic;
  88       fHAssocPt[at][ic] = new TH1F(fHname,  fHtit, 10*nBinsAssoc, minAssoc, maxAssoc); // for <pT> each bin
  89       histosContainer->AddLast(fHAssocPt[at][ic]);
  90     } // loop over trigger types
  91   } // centrality binning loop
  92
  93   fHCentr = new TH1F("fHCentr","centrality distribution",100, 0., 100.);
  94   histosContainer->AddLast(fHCentr);
  95   fHZVert = new TH1F("fHZVert","vertex distribution",100, -50., 50.);
  96   histosContainer->AddLast(fHZVert);
  97
  98   fHBinsCentr = new TH1F("fHBinsCentr","centrality binning", nBinsCentr+1, 0, 1);
  99   histosContainer->AddLast(fHBinsCentr);
 100   for(UInt_t i=1;i<=nBinsCentr+1; i++)
 101     fHBinsCentr->SetBinContent(i,fSelector->BinBorder(centr,i-1));
 102   fHBinsZVert = new TH1F("fHBinsZVert","centrality binning", nBinsZVert+1, 0, 1);
 103   histosContainer->AddLast(fHBinsZVert);
 104   for(UInt_t i=1;i<=nBinsZVert+1; i++)
 105     fHBinsZVert->SetBinContent(i,fSelector->BinBorder(zvert,i-1));
 106   fHBinsTrigg = new TH1F("fHBinsTrigg","trigger binning", nBinsTrigg+1, 0, 1);
 107   histosContainer->AddLast(fHBinsTrigg);
 108   for(UInt_t i=1;i<=nBinsTrigg+1; i++)
 109     fHBinsTrigg->SetBinContent(i,fSelector->BinBorder(trigg,i-1));
 110   fHBinsAssoc = new TH1F("fHBinsAssoc","associated binning", nBinsAssoc+1, 0, 1);
 111   histosContainer->AddLast(fHBinsAssoc);
 112   for(UInt_t i=1;i<=nBinsAssoc+1; i++)
 113     fHBinsAssoc->SetBinContent(i,fSelector->BinBorder(assoc,i-1));
 114 }
 115
 116 void AliJetCorrelWriter::CreateQA(TList *histosContainer){
 117   // books QA histograms
 118   TString when[2] = {"before cuts","after cuts"};
 119   for(UInt_t i=0; i<2; i++){
 120     fHname = "fHTrkITSQA"; fHname += i;
 121     fHtit = "ITS nClust vs Chi2/nClust "; fHtit += when[i];
 122     fHTrkITSQA[i] = new TH2F(fHname,fHtit,20,0.,20.,50,0.,10.);
 123     histosContainer->AddLast(fHTrkITSQA[i]);
 124     fHname = "fHTrkTPCQA"; fHname += i;
 125     fHtit = "TPC nClust vs Chi2/nClust "; fHtit += when[i];
 126     fHTrkTPCQA[i] = new TH2F(fHname,fHtit,30,0.,150.,50,0.,10.);
 127     histosContainer->AddLast(fHTrkTPCQA[i]);
 128     fHname = "fHTrkVTXQA"; fHname += i;
 129     fHtit = "VTX KinkIndex vs nSigma "; fHtit += when[i];
 130     fHTrkVTXQA[i] = new TH2F(fHname,fHtit,21,-10.,10,50,0.,10.);
 131     histosContainer->AddLast(fHTrkVTXQA[i]);
 132   }
 133
 134   UInt_t nTypeTrigg = fMaker->NoOfTrigg();
 135   UInt_t nTypeAssoc = fMaker->NoOfAssoc();
 136   UInt_t nBinsCentr = fSelector->NoOfBins(centr);
 137   UInt_t nBinsTrigg = fSelector->NoOfBins(trigg);
 138   UInt_t nBinsAssoc = fSelector->NoOfBins(assoc);
 139   Float_t minTrigg  = fSelector->MinLowBin(trigg);
 140   Float_t maxTrigg  = fSelector->MaxHighBin(trigg);
 141   Float_t minAssoc  = fSelector->MinLowBin(assoc);
 142   Float_t maxAssoc  = fSelector->MaxHighBin(assoc);
 143   for(UInt_t ic=0; ic<nBinsCentr; ic++){ // centrality loop
 144     for(UInt_t it=0; it<2; it++){ // mixing type loop (real/mixed)
 145       fHname="fHTrkProx"; fHname+=it; fHname+=ic;
 146       fHtit="fHTrkProx"; fHtit+=it; fHtit+=ic;
 147       fHTrkProx[it][ic] = new TH3F(fHname,fHtit,100,0.,10.,
 148                                   nBinsTrigg,minTrigg,maxTrigg,nBinsAssoc,minAssoc,maxAssoc);
 149       histosContainer->AddLast(fHTrkProx[it][ic]);
 150     } // loop over mixing type
 151   } // loop over centrality
 152   for(UInt_t tt=0; tt<nTypeTrigg; tt++){ // loop over trigger types
 153     fHname = "fHTriggAcc"; fHname+=tt;
 154     fHTriggAcc[tt] = new TH3F(fHname,fHname,62,0.,2*TMath::Pi(),20,-1.,1.,nBinsTrigg,minTrigg,maxTrigg);
 155     histosContainer->AddLast(fHTriggAcc[tt]);
 156   } // loop over trigger types
 157   for(UInt_t ta=0; ta<nTypeAssoc; ta++){ // loop over associated types
 158     fHname="fHAssocAcc"; fHname+=ta;
 159     fHAssocAcc[ta] = new TH3F(fHname,fHname,62,0.,2*TMath::Pi(),20,-1.,1.,nBinsAssoc,minAssoc,maxAssoc);
 160     histosContainer->AddLast(fHAssocAcc[ta]);
 161   } // loop over associated types
 162 }
 163
 164 void AliJetCorrelWriter::CreateCorrelations(TList* histosContainer){
 165   // books correlation histograms
 166   const Float_t lr=-TMath::Pi()/3, ur=5*TMath::Pi()/3, bwPout=0.2;
 167   UInt_t nTypeCorrel = fMaker->NoOfCorrel();
 168   UInt_t nBinsCentr = fSelector->NoOfBins(centr);
 169   UInt_t nBinsZVert = fSelector->NoOfBins(zvert);
 170   UInt_t nBinsTrigg = fSelector->NoOfBins(trigg);
 171   UInt_t nBinsAssoc = fSelector->NoOfBins(assoc);
 172   Float_t maxAssoc  = fSelector->MaxHighBin(assoc);
 173   UInt_t kDPhiNumBins = fSelector->DPhiNumBins();
 174   UInt_t kDEtaNumBins = fSelector->DEtaNumBins();
 175   UInt_t nPoutBins = UInt_t(TMath::Ceil(2.2*maxAssoc/bwPout)); // since |p_out|<p_Ta
 176   if(fRecoTrigg) {  // if any correlation has reconstructed trigger, define ntuple; use id to differentiate
 177     fNtuParent = new TNtuple("fNtuParent","Reconstructed Parent Ntuple","id:q:m:pT:phi:eta:assym:open");
 178     histosContainer->AddLast(fNtuParent);
 179   }
 180   for(UInt_t htc=0; htc<nTypeCorrel; htc++){ // loop over correlation types
 181     for(UInt_t hic=0; hic<nBinsCentr; hic++){ // centrality loop
 182       for(UInt_t hiv=0; hiv<nBinsZVert; hiv++){ // vertex loop
 183         for(UInt_t hit=0; hit<nBinsTrigg; hit++){ // trigger loop
 184           for(UInt_t hia=0; hia<nBinsAssoc; hia++){ // associated loop
 185             fHtit = fMaker->Descriptor(htc); fHtit+=":"; fHtit+=hic; fHtit+=hiv; fHtit+=hit; fHtit+=hia;
 186             fHname="fHReal"; fHname+=htc; fHname+=hic; fHname+=hiv; fHname+=hit; fHname+=hia;
 187             fHReal[htc][hic][hiv][hit][hia] = new TH3F(fHname,fHtit, kDPhiNumBins,lr,ur,
 188                                                       kDEtaNumBins,-2.,2.,
 189                                                       nPoutBins,-1.1*maxAssoc,1.1*maxAssoc);
 190             histosContainer->AddLast(fHReal[htc][hic][hiv][hit][hia]);
 191             fHname="fHMix"; fHname+=htc; fHname+=hic; fHname+=hiv; fHname+=hit; fHname+=hia;
 192             fHMix[htc][hic][hiv][hit][hia] = new TH3F(fHname,fHtit, kDPhiNumBins,lr,ur,
 193                                                       kDEtaNumBins,-2.,2.,
 194                                                       nPoutBins,-1.1*maxAssoc,1.1*maxAssoc);
 195             histosContainer->AddLast(fHMix[htc][hic][hiv][hit][hia]);
 196           } // loop over associated bins
 197         } // loop over trigger bins
 198       } // loop over vertex bins
 199     } // loop over centrality bins
 200   } // loop over correlation types
 201 }
 202
 203 /////////////////////////////////////////////////////////
 204 // METHODS FOR FILLING THE OUTPUT HISTOGRAMS
 205 /////////////////////////////////////////////////////////
 206
 207 void AliJetCorrelWriter::FillGlobal(Float_t cent, Float_t zvert){
 208   // some global event histos
 209   fHCentr->Fill(cent);
 210   fHZVert->Fill(zvert);
 211 }
 212
 213 void AliJetCorrelWriter::FillSingleHistos(UInt_t cBin, CorrelList_t * const TriggList, UInt_t tIdx,
 214                                           CorrelList_t * const AssocList, UInt_t aIdx){
 215   // fills single-particle histograms
 216   if(TriggList->Size()>0){
 217     CorrelListIter_t trigIter = TriggList->Head();
 218     while(!trigIter.HasEnded()){
 219       CorrelParticle_t *particle = trigIter.Data();
 220       fHTriggPt[tIdx][cBin]->Fill(particle->Pt());
 221       if(fSelector->GenQA()) fHTriggAcc[tIdx]->Fill(particle->Phi(),particle->Eta(),particle->Pt());
 222       trigIter.Move();
 223     }
 224   }
 225   if(AssocList->Size()>0){
 226     CorrelListIter_t assoIter = AssocList->Head();
 227     while(!assoIter.HasEnded()){
 228       CorrelParticle_t *particle = assoIter.Data();
 229       fHAssocPt[aIdx][cBin]->Fill(particle->Pt());
 230       if(fSelector->GenQA()) fHAssocAcc[aIdx]->Fill(particle->Phi(),particle->Eta(),particle->Pt());
 231       assoIter.Move();
 232     }
 233   }
 234 }
 235
 236 void AliJetCorrelWriter::FillTrackQA(AliESDtrack * const track, UInt_t idx){
 237   // fills single-particle QA
 238   if(idx>1){std::cerr<<"AliJetCorrelWriter::FillTrackQA: wrong idx!"<<std::endl; exit(-1);}
 239
 240   UInt_t nClusITS = track->GetITSclusters(0);
 241   UInt_t nClusTPC = track->GetTPCclusters(0); // or track->GetTPCNcls() ?
 242   Float_t chi2ITS=-1., chi2TPC=-1.;
 243   if(nClusITS!=0) chi2ITS = track->GetITSchi2()/Float_t(nClusITS);
 244   if(nClusTPC!=0) chi2TPC = track->GetTPCchi2()/Float_t(nClusTPC);
 245   UInt_t kinkIndex = track->GetKinkIndex(0);
 246   Float_t sigVtx = fSelector->GetSigmaToVertex(track);
 247
 248   fHTrkITSQA[idx]->Fill(Float_t(nClusITS),chi2ITS);
 249   fHTrkTPCQA[idx]->Fill(Float_t(nClusTPC),chi2TPC);
 250   fHTrkVTXQA[idx]->Fill(Float_t(kinkIndex),sigVtx);
 251 }
 252
 253 void AliJetCorrelWriter::FillParentNtuple(CorrelList_t * const ParentList){
 254   // fills ntuple of triggers when they are reconstructed parents (pi0,Z0,etc.)
 255   if(!fRecoTrigg)
 256     {std::cerr<<"AliJetCorrelWriter::FillParentNtuple: you shouldn't be here!"<<std::endl; exit(-1);}
 257   if(ParentList->Size()<1) return;
 258   Float_t parVar[8]={0.};
 259   CorrelListIter_t parIter = ParentList->Head();
 260   while(!parIter.HasEnded()){
 261     CorrelRecoParent_t *parent = dynamic_cast<CorrelRecoParent_t*>(parIter.Data());
 262     if(!parent)
 263       {std::cerr<<"AliJetCorrelWriter::FillParentNtuple: failed casting!"<<std::endl; exit(-1);}
 264     parVar[0] = parent->ID();
 265     parVar[1] = parent->Q();
 266     parVar[2] = parent->M();
 267     parVar[3] = parent->Pt();
 268     parVar[4] = parent->Phi();
 269     parVar[5] = parent->Eta();
 270     parVar[6] = parent->Assym();
 271     parVar[7] = parent->OpenAng();
 272     fNtuParent->Fill(parVar);
 273     parIter.Move();
 274   }
 275 }
 276
 277 void AliJetCorrelWriter::FillCorrelations(FillType_t fTyp, UInt_t iCorr, UInt_t cBin, UInt_t vBin,
 278                                           CorrelParticle_t * const Trigg, CorrelParticle_t * const Assoc){
 279   // fills the correlation (two-particle) histograms
 280   // trigger information (this is why the first particle has to be the trigger):
 281   Float_t ptt  = Trigg->Pt();
 282   Float_t phit = Trigg->Phi();
 283   Float_t etat = Trigg->Eta();
 284   Int_t   tBin = fSelector->GetBin(trigg,ptt);
 285   // associated information:
 286   Float_t pta  = Assoc->Pt();
 287   Float_t phia = Assoc->Phi();
 288   Float_t etaa = Assoc->Eta();
 289   Int_t   aBin = fSelector->GetBin(assoc,pta);
 290   //  Short_t qprod= Trigg->Q()*Assoc->Q();
 291
 292   if(tBin<0 || aBin<0) return;  // one of them is not in the required pT range
 293   if(pta>=ptt) return; // use only associated particles below the trigger
 294   if(fabs(ptt-pta)<1.e-6 && fabs(phit-phia)<1.e-6 && fabs(etat-etaa)<1.e-6) return; // don't auto-correlate
 295
 296   // store track pair proximity
 297   if(Trigg->ID()==hadron && Assoc->ID()==hadron){
 298     CorrelTrack_t* trk1 = dynamic_cast<CorrelTrack_t*>(Trigg);
 299     CorrelTrack_t* trk2 = dynamic_cast<CorrelTrack_t*>(Assoc);
 300     if(!trk1 || !trk2)
 301       {std::cerr<<"AliJetCorrelWriter::FillCorrelations: failed casting!"<<std::endl; exit(-1);}
 302     Float_t pairDist = trk1->Dist(trk2);
 303     if(fSelector->CloseTrackPair(pairDist)) return; // proximity cut
 304     if(fSelector->GenQA()) fHTrkProx[fTyp][cBin]->Fill(pairDist,ptt,pta);
 305   }
 306
 307   // Fill correlation histograms:
 308   Float_t dphi = DeltaPhi(phit,phia);
 309   Float_t deta = etat-etaa;
 310   Float_t pout = pta*TMath::Sin(dphi);
 311   if(fTyp==real)
 312     fHReal[iCorr][cBin][vBin][tBin][aBin]->Fill(dphi,deta,pout);
 313   else
 314     fHMix[iCorr][cBin][vBin][tBin][aBin]->Fill(dphi,deta,pout);
 315 }
 316
 317 Float_t AliJetCorrelWriter::DeltaPhi(Float_t  phi1, Float_t phi2) {
 318 // wrapps dphi in (-pi/3,5pi/3)
 319   Float_t kPi = TMath::Pi();
 320   Float_t res = phi2-phi1;
 321   if(fRndm.Uniform()<0.5) res = phi1-phi2;
 322   if(5.*kPi/3.<res && res<=2.*kPi) res-=2*kPi;
 323   else if(-2.*kPi<=res && res<-kPi/3.) res+=2*kPi;
 324 //   if(3.*kPi/2.<res && res<=2*kPi) res-=2*kPi;
 325 //   else if(-2.*kPi<=res && res<-kPi/2.) res+=2*kPi;
 326   return res;
 327 }
 328