PWG/FLOW/Tasks/AliAnalysisTaskJetFlowMC.cxx

   1 // Simple class to generate toy events which can be fed into the jet finder
   2 //
   3 // Author: Redmer Alexander Bertens, Utrecht University, 2013
   4 // rbertens@cern.ch, rbertens@nikhef.nl, r.a.bertens@uu.nl
   5
   6 // root includes
   7 #include "TChain.h"
   8 #include "TList.h"
   9 #include "TClonesArray.h"
  10 #include "TArrayI.h"
  11 #include "TRandom3.h"
  12 // aliroot includes
  13 #include "AliAODEvent.h"
  14 #include "AliAnalysisManager.h"
  15 #include "AliAnalysisTaskJetFlowMC.h"
  16 #include "AliLog.h"
  17 #include "AliPicoTrack.h"
  18
  19 class AliAnalysisTaskJetFlowMC;
  20 using namespace std;
  21
  22 ClassImp(AliAnalysisTaskJetFlowMC)
  23
  24 //_____________________________________________________________________________
  25 AliAnalysisTaskJetFlowMC::AliAnalysisTaskJetFlowMC() : AliAnalysisTaskSE("AliAnalysisTaskJetFlowMC"), fTracksOutName("JetFlowMC"),  fTracksInName("PicoTrack"), fTracksIn(0),  fTracksOut(0), fReuseTracks(kFALSE), fMult(2200), fCenBin(-1), fCentralityClasses(0), fFuncVn(0), fOutputList(0), fTrackSpectrum(0), fRandomizeEta(kTRUE), fJetSpectrumSF(0), fNoOfSFJets(0), fHistIntV2(0), fHistIntV3(0), fFlowFluctuations(-10), fMaxNumberOfIterations(100), fPsi2(-10), fPsi3(-10), fPrecisionPhi(1e-10), fDetectorType(kVZEROC), fHistSFJetSpectrum(0), fHistSFJetEtaPhi(0) {
  26     // default constructor for root IO
  27     for(Int_t i(0); i < 10; i++) {
  28         fFuncDiffV2[i]                  = 0x0;
  29         fFuncDiffV3[i]                  = 0x0;
  30         fHistDiffV2[i]                  = 0x0;
  31         fHistDiffV3[i]                  = 0x0;
  32         fHistOriginalSpectrum[i]        = 0x0;
  33         fHistOriginalEtaPhi[i]          = 0x0;
  34         fHistToySpectrum[i]             = 0x0;
  35         fHistToyEtaPhi[i]               = 0x0;
  36         fHistOriginalDeltaPhi[i]        = 0x0;
  37         fHistToyDeltaPhi[i]             = 0x0;
  38         fHistToyVn[i]                   = 0x0;
  39     }
  40 }
  41 //_____________________________________________________________________________
  42 AliAnalysisTaskJetFlowMC::AliAnalysisTaskJetFlowMC(const char *name) : AliAnalysisTaskSE(name), fTracksOutName("JetFlowMC"), fTracksInName("PicoTrack"), fTracksIn(0), fTracksOut(0), fReuseTracks(kFALSE), fMult(2200), fCenBin(-1), fCentralityClasses(0), fFuncVn(0), fOutputList(0), fTrackSpectrum(0), fRandomizeEta(kTRUE), fJetSpectrumSF(0), fNoOfSFJets(0), fHistIntV2(0), fHistIntV3(0), fFlowFluctuations(-10), fMaxNumberOfIterations(100), fPsi2(-10), fPsi3(-10), fPrecisionPhi(1e-10), fDetectorType(kVZEROC), fHistSFJetSpectrum(0), fHistSFJetEtaPhi(0) {
  43     // constructor
  44     DefineInput(0, TChain::Class());
  45     DefineOutput(1, TList::Class());
  46     for(Int_t i(0); i < 10; i++) {
  47         fFuncDiffV2[i]                  = 0x0;
  48         fFuncDiffV3[i]                  = 0x0;
  49         fHistDiffV2[i]                  = 0x0;
  50         fHistDiffV3[i]                  = 0x0;
  51         fHistOriginalSpectrum[i]        = 0x0;
  52         fHistOriginalEtaPhi[i]          = 0x0;
  53         fHistToySpectrum[i]             = 0x0;
  54         fHistToyEtaPhi[i]               = 0x0;
  55         fHistOriginalDeltaPhi[i]        = 0x0;
  56         fHistToyDeltaPhi[i]             = 0x0;
  57         fHistToyVn[i]                   = 0x0;
  58     }
  59 }
  60 //_____________________________________________________________________________
  61 AliAnalysisTaskJetFlowMC::~AliAnalysisTaskJetFlowMC()
  62 {
  63     // desctructor, claim ownership of deleted objects by setting member pointers to NULL
  64     if(fOutputList) {
  65       delete fOutputList;
  66       fOutputList= NULL;
  67     }
  68 }
  69 //_____________________________________________________________________________
  70 void AliAnalysisTaskJetFlowMC::UserCreateOutputObjects()
  71 {
  72     // Create my user objects.
  73     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
  74     fTracksOut = new TClonesArray("AliPicoTrack");
  75     fTracksOut->SetName(fTracksOutName);
  76     fOutputList = new TList();
  77     fOutputList->SetOwner(kTRUE);
  78     if(!fCentralityClasses) { // classes must be defined at this point
  79         Int_t c[] = {0, 90};
  80         fCentralityClasses = new TArrayI(sizeof(c)/sizeof(c[0]), c);
  81     }
  82     if(fHistIntV2 && !fHistIntV3) {      // define function
  83         fFuncVn = new TF1("kV2", "[0]*([1]+[2]*[3]*TMath::Cos([2]*(x-[4])))", 0, TMath::TwoPi());
  84         fFuncVn->SetParameter(0, 1.);        // normalization
  85         fFuncVn->SetParameter(3, 0.2);       // v2
  86         fFuncVn->FixParameter(1, 1.);        // constant
  87         fFuncVn->FixParameter(2, 2.);        // constant
  88     } else if (!fHistIntV2 && fHistIntV3) {
  89         fFuncVn = new TF1("kV3", "[0]*([1]+[2]*[3]*TMath::Cos([2]*(x-[4])))", 0, TMath::TwoPi());
  90         fFuncVn->SetParameter(0, 1.);        // normalization
  91         fFuncVn->SetParameter(3, 0.2);       // v3
  92         fFuncVn->FixParameter(1, 1.);        // constant
  93         fFuncVn->FixParameter(2, 3.);        // constant
  94     } else if (fHistIntV2 && fHistIntV3) {
  95         fFuncVn = new TF1("kCombined", "[0]*([1]+[2]*([3]*TMath::Cos([2]*(x-[4]))+[7]*TMath::Cos([5]*(x-[6]))))", 0, TMath::TwoPi());
  96         fFuncVn->SetParameter(0, 1.);       // normalization
  97         fFuncVn->SetParameter(3, 0.2);      // v2
  98         fFuncVn->FixParameter(1, 1.);       // constant
  99         fFuncVn->FixParameter(2, 2.);       // constant
 100         fFuncVn->FixParameter(5, 3.);       // constant
 101         fFuncVn->SetParameter(7, 0.2);      // v3
 102     }
 103     // add the generator objects that have been added to the task
 104     if(fTrackSpectrum)  fOutputList->Add(fTrackSpectrum);
 105     if(fJetSpectrumSF)  fOutputList->Add(fJetSpectrumSF);
 106     if(fHistIntV2)      fOutputList->Add(fHistIntV2);
 107     if(fHistIntV3)      fOutputList->Add(fHistIntV3);
 108     // create the QA histos
 109     for(Int_t i(0); i < fCentralityClasses->GetSize()-1; i++) {
 110         fHistOriginalSpectrum[i] = BookTH1F("fHistOriginalSpectrum", "p_{t} [GeV/c]", 200, 0, 200, i);
 111         fHistOriginalEtaPhi[i] = BookTH2F("fHistOriginalEtaPhi", "#eta", "#varphi", 100, -1., 1., 100, 0, TMath::TwoPi(), i);
 112         fHistToySpectrum[i] = BookTH1F("fHistToySpectrum", "p_{t} [GeV/c]", 200, 0, 200, i);
 113         fHistToyEtaPhi[i] = BookTH2F("fHistToyEtaPhi", "#eta", "#varphi", 100, -1., 1., 100, 0, TMath::TwoPi(), i);
 114         fHistOriginalDeltaPhi[i] = BookTH1F("fHistOriginalDeltaPhi", "#varphi - #Psi", 100, 0., TMath::Pi(), i);
 115         fHistToyDeltaPhi[i] = BookTH1F("fHistToyDeltaPhi", "#varphi - #Psi", 100, 0., TMath::Pi(), i);
 116         fHistToyVn[i] = BookTH2F("fHistToyVn", "p_{t} [GeV/c]", "v_{n}", 1000, 0, 200, 200, 0, .2, i);
 117         // add to outputlist
 118         if(fFuncDiffV2[i]) fOutputList->Add(fFuncDiffV2[i]);
 119         if(fFuncDiffV3[i]) fOutputList->Add(fFuncDiffV3[i]);
 120         if(fHistDiffV2[i]) fOutputList->Add(fHistDiffV2[i]);
 121         if(fHistDiffV3[i]) fOutputList->Add(fHistDiffV3[i]);
 122     }
 123     if(fJetSpectrumSF) {
 124         fHistSFJetSpectrum = BookTH1F("fHistSFJetSpectrum", "p_{t} SF jets [GeV/c]", 100, 0, 200);
 125         fHistSFJetEtaPhi = BookTH2F("fHistSFJetEtaPhi", "#eta", "#varphi", 100, -1., 1., 100, 0, TMath::TwoPi());
 126     }
 127     // reset the random seed
 128     if(gRandom) {
 129         delete gRandom;
 130         gRandom = new TRandom3(0);
 131     }
 132
 133     fOutputList->Sort();
 134     PostData(1, fOutputList);
 135 }
 136 //_____________________________________________________________________________
 137 TH1F* AliAnalysisTaskJetFlowMC::BookTH1F(const char* name, const char* x, Int_t bins, Double_t min, Double_t max, Int_t c, Bool_t append)
 138 {
 139     // book a TH1F and connect it to the output container
 140     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 141     if(!fOutputList) return 0x0;
 142     TString title(name);
 143     if(c!=-1) { // format centrality dependent histograms accordingly
 144         name = Form("%s_%i", name, c);
 145         title += Form("_%i-%i", fCentralityClasses->At(c), fCentralityClasses->At(1+c));
 146     }
 147     title += Form(";%s;[counts]", x);
 148     TH1F* histogram = new TH1F(name, title.Data(), bins, min, max);
 149     histogram->Sumw2();
 150     if(append) fOutputList->Add(histogram);
 151     return histogram;
 152 }
 153 //_____________________________________________________________________________
 154 TH2F* AliAnalysisTaskJetFlowMC::BookTH2F(const char* name, const char* x, const char*y, Int_t binsx, Double_t minx, Double_t maxx, Int_t binsy, Double_t miny, Double_t maxy, Int_t c, Bool_t append)
 155 {
 156     // book a TH2F and connect it to the output container
 157     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 158     if(!fOutputList) return 0x0;
 159     TString title(name);
 160     if(c!=-1) { // format centrality dependent histograms accordingly
 161         name = Form("%s_%i", name, c);
 162         title += Form("_%i-%i", fCentralityClasses->At(c), fCentralityClasses->At(1+c));
 163     }
 164     title += Form(";%s;%s", x, y);
 165     TH2F* histogram = new TH2F(name, title.Data(), binsx, minx, maxx, binsy, miny, maxy);
 166     histogram->Sumw2();
 167     if(append) fOutputList->Add(histogram);
 168     return histogram;
 169 }
 170 //_____________________________________________________________________________
 171 void AliAnalysisTaskJetFlowMC::UserExec(Option_t *)
 172 {
 173     // user exec, called for each event.
 174     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 175     if(!AliAnalysisManager::GetAnalysisManager()) return;
 176     // retrieve tracks from input.
 177     if (!fTracksIn) {
 178         fTracksIn = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksInName));
 179         if (!fTracksIn) {
 180           AliError(Form("Could not retrieve tracks %s!", fTracksInName.Data()));
 181           return;
 182         }
 183     }
 184     // get the centrality bin for qa plots
 185     Double_t cent(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"));
 186     fCenBin = -1;
 187     for(Int_t i(0); i < fCentralityClasses->GetSize()-1; i++) {
 188         if(cent >= fCentralityClasses->At(i) && cent <= fCentralityClasses->At(1+i)) {
 189             fCenBin = i;
 190             break; }
 191     }
 192     if(fCenBin < 0) return;
 193     // add tracks to event if not yet there
 194     fTracksOut->Delete();
 195     if (!(InputEvent()->FindListObject(fTracksOutName))) InputEvent()->AddObject(fTracksOut);
 196     fTracksOut->Clear();
 197     // get the event plane
 198     CalculateEventPlane();
 199     Int_t nacc(0);
 200     if(fReuseTracks) {
 201         const Int_t Ntracks = fTracksIn->GetEntriesFast();
 202         for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
 203             AliPicoTrack* track = static_cast<AliPicoTrack*>(fTracksIn->At(iTracks));
 204             if(!track) continue;
 205             Double_t phi(track->Phi()), pt((fTrackSpectrum) ? GetTrackPt() : track->Pt()), eta(fRandomizeEta ? GetTrackEta() : track->Eta());
 206             // fill qa histo's before applying any (possible) afterburner
 207             FillHistogramsOriginalData(pt, eta, phi);
 208             if(fHistDiffV2[fCenBin] || fFuncDiffV2[fCenBin])        V2AfterBurner(phi, eta, pt);
 209             else if(fHistDiffV3[fCenBin] || fFuncDiffV3[fCenBin])   V3AfterBurner(phi, eta, pt);
 210             else if(fHistIntV2 || fHistIntV3)                       SampleVnFromTF1(phi);
 211             /*AliPicoTrack *picotrack =*/ new ((*fTracksOut)[nacc]) AliPicoTrack(pt, eta, phi, track->Charge(), track->GetLabel(), 4, track->GetTrackEtaOnEMCal(), track->GetTrackPhiOnEMCal(), track->GetTrackPtOnEMCal(), 1);
 212             nacc++;
 213         }
 214     } else {
 215         Double_t pt(0), eta(0), phi(0);
 216         for (Int_t iTracks = 0; iTracks < fMult; ++iTracks) {
 217             pt = GetTrackPt();
 218             eta = gRandom->Uniform(-.9, .9);
 219             phi = gRandom->Uniform(0., TMath::TwoPi());
 220             // fill qa histo's before applying any (possible) afterburner
 221             FillHistogramsOriginalData(pt, eta, phi);
 222             if(fHistDiffV2[fCenBin] || fFuncDiffV2[fCenBin])        V2AfterBurner(phi, eta, pt);
 223             else if(fHistDiffV3[fCenBin] || fFuncDiffV3[fCenBin])   V3AfterBurner(phi, eta, pt);
 224             else if(fHistIntV2 || fHistIntV3)                       SampleVnFromTF1(phi);
 225             /*AliPicoTrack *picotrack =*/ new ((*fTracksOut)[nacc]) AliPicoTrack(pt, eta, phi, 1, 0, 4, eta, phi, pt, 1);
 226             nacc++;
 227         }
 228     }
 229     if(fJetSpectrumSF && fNoOfSFJets > 0) InjectSingleFragmentationJetSpectrum(nacc);
 230     PostData(1, fOutputList);
 231     if(fDebug > 0) PrintInfo();
 232 }
 233 //_____________________________________________________________________________
 234 void AliAnalysisTaskJetFlowMC::V2AfterBurner(Double_t &phi, Double_t &eta, Double_t &pt) const
 235 {
 236     // similar to AliFlowTrackSimple::AddV2, except for the flow fluctuations
 237     if(fDebug > 1) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 238     phi = gRandom->Uniform(0, TMath::TwoPi());
 239     Double_t phi0(phi), v2(GetV2(pt)), f(0.), fp(0.), phiprev(0.);
 240     if(TMath::AreEqualAbs(v2, 0, 1e-5)) {
 241         FillHistogramsToyData(pt, eta, phi, v2);
 242         return;
 243     }
 244     // introduce flow fluctuations (gaussian)
 245     if(fFlowFluctuations > -10) GetFlowFluctuation(v2);
 246     for (Int_t i=0; i!=fMaxNumberOfIterations; ++i) {
 247         phiprev=phi; //store last value for comparison
 248         f =  phi-phi0+v2*TMath::Sin(2.*(phi-fPsi2));
 249         fp = 1.0+2.0*v2*TMath::Cos(2.*(phi-fPsi2)); //first derivative
 250         phi -= f/fp;
 251         if (TMath::AreEqualAbs(phiprev,phi,fPrecisionPhi)) break;
 252     }
 253     FillHistogramsToyData(pt, eta, phi, v2);
 254 }
 255 //_____________________________________________________________________________
 256 void AliAnalysisTaskJetFlowMC::V3AfterBurner(Double_t &phi, Double_t &eta, Double_t &pt) const
 257 {
 258     // similar to AliFlowTrackSimple::AddV3, except for the flow fluctuations
 259     if(fDebug > 1) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 260     phi = gRandom->Uniform(0, TMath::TwoPi());
 261     Double_t phi0(phi), v3(GetV3(pt)), f(0.), fp(0.), phiprev(0.);
 262     if(TMath::AreEqualAbs(v3, 0, 1e-5)) {
 263         FillHistogramsToyData(pt, eta, phi, v3);
 264         return;
 265     }
 266     // introduce flow fluctuations (gaussian)
 267     if(fFlowFluctuations > -10) GetFlowFluctuation(v3);
 268     for (Int_t i=0; i<fMaxNumberOfIterations; i++)  {
 269         phiprev=phi; //store last value for comparison
 270         f =  phi-phi0+2./3.*v3*TMath::Sin(3.*(phi-fPsi3));
 271         fp = 1.0+2.0*v3*TMath::Cos(3.*(phi-fPsi3)); //first derivative
 272         phi -= f/fp;
 273         if (TMath::AreEqualAbs(phiprev,phi,fPrecisionPhi)) break;
 274     }
 275     FillHistogramsToyData(pt, eta, phi, v3);
 276 }
 277 //_____________________________________________________________________________
 278 void AliAnalysisTaskJetFlowMC::InjectSingleFragmentationJetSpectrum(Int_t nacc)
 279 {
 280     // inject single fragmentation jet spectrum to the tclones array, note that emcal params
 281     // equal the barrel kinematics to pass the track and jet cuts later on
 282     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 283     for(Int_t i(nacc); i < (nacc + fNoOfSFJets); i++) {
 284         Double_t eta(gRandom->Uniform(-.5, .5)), phi(gRandom->Uniform(0, TMath::TwoPi())), pt(fJetSpectrumSF->GetRandom());
 285         /*AliPicoTrack *picotrack =*/ new ((*fTracksOut)[i]) AliPicoTrack(pt, eta, phi, +1, 0, 0, eta, phi, pt, 0);
 286         fHistSFJetSpectrum->Fill(pt);
 287         fHistSFJetEtaPhi->Fill(eta, phi);
 288         ++i;
 289     }
 290     nacc = 0;
 291 }
 292 //_____________________________________________________________________________
 293 void AliAnalysisTaskJetFlowMC::CalculateEventPlane() {
 294     // grab the event plane orientation from the AliVEvent header
 295     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 296     Double_t a(0), b(0), e(0), f(0);
 297     switch (fDetectorType) {
 298         case kVZEROA : {
 299             fPsi2 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 8, 2, e, f);
 300             fPsi3 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 8, 3, a, b);
 301             break;
 302         }
 303         case kVZEROC : {
 304             fPsi2 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 9, 2, e, f);
 305             fPsi3 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 9, 3, a, b);
 306             break;
 307                        }
 308         case kVZEROComb : {
 309             fPsi2 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 10, 2, e, f) ;
 310             fPsi3 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 10, 3, a, b);
 311             break;
 312         }
 313         default : break;
 314     }
 315     // if requested, pass the event plane to the phi distribution
 316     if(fHistIntV2 && fHistIntV3) {
 317         fFuncVn->SetParameter(4, fPsi2);        fFuncVn->SetParameter(6, fPsi3);
 318         Double_t v2(fHistIntV2->GetBinContent(fHistIntV2->GetXaxis()->FindBin(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"))));
 319         Double_t v3(fHistIntV2->GetBinContent(fHistIntV3->GetXaxis()->FindBin(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"))));
 320         if(fFlowFluctuations > -10) {
 321             GetFlowFluctuation(v2);             GetFlowFluctuation(v3);
 322         }
 323         fFuncVn->SetParameter(3, v2);           fFuncVn->SetParameter(7, v3);
 324     } else if (fHistIntV2) {
 325         fFuncVn->SetParameter(4, fPsi2);
 326         Double_t v2(fHistIntV2->GetBinContent(fHistIntV2->GetXaxis()->FindBin(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"))));
 327         if(fFlowFluctuations > -10) GetFlowFluctuation(v2);
 328         fFuncVn->SetParameter(3, v2);
 329     } else if (fHistIntV3) {
 330         fFuncVn->SetParameter(4, fPsi3);
 331         Double_t v3(fHistIntV3->GetBinContent(fHistIntV2->GetXaxis()->FindBin(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"))));
 332         if(fFlowFluctuations > -10) GetFlowFluctuation(v3);
 333         fFuncVn->SetParameter(3, v3);
 334     }
 335 }
 336 //_____________________________________________________________________________
 337 void AliAnalysisTaskJetFlowMC::Terminate(Option_t *)
 338 {
 339     // terminate
 340 }
 341 //_____________________________________________________________________________
 342 void AliAnalysisTaskJetFlowMC::PrintInfo() const
 343 {
 344     // print info
 345     printf(" > No of retrieved tracks from %s \n \t %i \n", fTracksInName.Data(), fTracksIn->GetEntriesFast());
 346     printf(" > No of created tracks in %s \n \t %i \n", fTracksOutName.Data(), fTracksOut->GetEntriesFast());
 347
 348 }
 349 //_____________________________________________________________________________