PWGCF/EBYE/macros/runBalanceFunctionOnHydro.C

   1 //========================================================//
   2 //Macro that reads the output of the hydro calculations
   3 //(Therminator) using the format that can be found at
   4 //arXiv:1102.0273
   5 //Author: Panos.Christakoglou@nikhef.nl
   6 //========================================================//
   7
   8 //========================================================//
   9 //Event structure
  10 struct StructEvent {
  11   UInt_t eventID;//unique identifier of the event (random number)
  12   UInt_t entries;//number of particles of each event
  13   UInt_t entriesprev;//set to 0 by default
  14 };
  15 //========================================================//
  16
  17 //========================================================//
  18 //Particle structure
  19 //Primordial particles: fathereid==-1 ==>
  20 //pid==fatherpid==rootpid
  21 struct StructParticle {
  22   Float_t mass;//the mass of the particle (in GeV)
  23   Float_t t;//the time coordinate of the particle in fm/c
  24   Float_t x;//the spacial coordinate x in fm/c
  25   Float_t y;//the spacial coordinate y in fm/c
  26   Float_t z;//the spacial coordinate z in fm/c
  27   Float_t e;//the energy in GeV
  28   Float_t px;//the x-coordinate of the particle's momentum in GeV
  29   Float_t py;//the y-coordinate of the particle's momentum in GeV
  30   Float_t pz;//the z-coordinate of the particle's momentum in GeV
  31   Int_t decayed;//decay flag (1: decayed, 0: no decay)
  32   Int_t pid;//PDG of particle
  33   Int_t fatherpid;//PDG of parent
  34   Int_t rootpid;//root (primordial) particle PDG number
  35   Int_t eid;//sequence number in the event
  36   Int_t fathereid;//parent sequence number in the event
  37   UInt_t eventid;//unique identifier of the event (random number)
  38 };
  39 //========================================================//
  40
  41 //========================================================//
  42 //Balance function analysis variables
  43 Bool_t gRunShuffling=kFALSE;
  44 Bool_t gRunMixing=kTRUE;
  45 Bool_t gRunMixingWithEventPlane=kFALSE;
  46
  47 Double_t gEtaMin = -0.8;
  48 Double_t gEtaMax = 0.8;
  49 Double_t gPtMin = 0.2;
  50 Double_t gPtMax = 20.0;
  51 //========================================================//
  52
  53
  54 void runBalanceFunctionOnHydro(TString aEventDir = "/glusterfs/alice1/alice2/pchrist/Therminator/pA/ChargeBalancing/Centrality1/lhyquid3v-LHCpPb5020s3.1Ti319t0.20Tf150e001/", Int_t aEventFiles = 2) {
  55   //Macro that reads the themrinator events
  56   //Author: Panos.Christakoglou@nikhef.nl
  57   // Time:
  58   TStopwatch timer;
  59   timer.Start();
  60
  61   //========================================================//
  62   //Load the aliroot libraries
  63   gSystem->Load("libSTEERBase.so");
  64   gSystem->Load("libESD.so");
  65   gSystem->Load("libAOD.so");
  66   gSystem->Load("libANALYSIS.so");
  67   gSystem->Load("libANALYSISalice.so");
  68   gSystem->Load("libEventMixing.so");
  69   gSystem->Load("libCORRFW.so");
  70   gSystem->Load("libPWGTools.so");
  71   gSystem->Load("libPWGCFebye.so");
  72   //========================================================//
  73
  74   //========================================================//
  75   //Configure the bf objects
  76   AliBalancePsi *bf = new AliBalancePsi();
  77   bf->SetAnalysisLevel("MC");
  78   bf->SetShuffle(gRunShuffling);
  79   bf->SetEventClass("EventPlane");
  80   bf->SetDeltaEtaMax(TMath::Abs(gEtaMax-gEtaMin));
  81   bf->InitHistograms();
  82
  83   AliBalancePsi *bfs = 0x0;
  84   if(gRunShuffling) {
  85     bfs = new AliBalancePsi();
  86     bfs->SetAnalysisLevel("MC");
  87     bfs->SetEventClass("EventPlane");
  88     bfs->SetDeltaEtaMax(TMath::Abs(gEtaMax-gEtaMin));
  89     bfs->InitHistograms();
  90   }
  91
  92   AliBalancePsi *bfm = 0x0;
  93   if(gRunMixing) {
  94     bfm = new AliBalancePsi();
  95     bfm->SetAnalysisLevel("MC");
  96     bfm->SetShuffle(gRunShuffling);
  97     bfm->SetEventClass("EventPlane");
  98     bfm->SetDeltaEtaMax(TMath::Abs(gEtaMax-gEtaMin));
  99     bfm->InitHistograms();
 100   }
 101   //========================================================//
 102
 103   //========================================================//
 104   //Output objects
 105   //QA list
 106   fList = new TList();
 107   fList->SetName("listQA");
 108   fList->SetOwner();
 109
 110   //Balance Function list
 111   TList *fListBF = new TList();
 112   fListBF->SetName("listBF");
 113   fListBF->SetOwner();
 114   fListBF->Add(bf->GetHistNp());
 115   fListBF->Add(bf->GetHistNn());
 116   fListBF->Add(bf->GetHistNpn());
 117   fListBF->Add(bf->GetHistNnn());
 118   fListBF->Add(bf->GetHistNpp());
 119   fListBF->Add(bf->GetHistNnp());
 120
 121   //Balance function list: shuffling
 122   TList *fListBFS = 0x0;
 123   if(gRunShuffling) {
 124     fListBFS = new TList();
 125     fListBFS->SetName("listBFShuffled");
 126     fListBFS->SetOwner();
 127     fListBFS->Add(bfs->GetHistNp());
 128     fListBFS->Add(bfs->GetHistNn());
 129     fListBFS->Add(bfs->GetHistNpn());
 130     fListBFS->Add(bfs->GetHistNnn());
 131     fListBFS->Add(bfs->GetHistNpp());
 132     fListBFS->Add(bfs->GetHistNnp());
 133   }
 134
 135   //Balance function list: event mixing
 136   TList *fListBFM = 0x0;
 137   if(gRunMixing) {
 138     fListBFM = new TList();
 139     fListBFM->SetName("listBFMixed");
 140     fListBFM->SetOwner();
 141     fListBFM->Add(bfm->GetHistNp());
 142     fListBFM->Add(bfm->GetHistNn());
 143     fListBFM->Add(bfm->GetHistNpn());
 144     fListBFM->Add(bfm->GetHistNnn());
 145     fListBFM->Add(bfm->GetHistNpp());
 146     fListBFM->Add(bfm->GetHistNnp());
 147   }
 148   //========================================================//
 149
 150   //========================================================//
 151   //Event Mixing
 152   if(gRunMixing){
 153     Int_t trackDepth = 50000;
 154     Int_t poolsize   = 1000;  // Maximum number of events, ignored in the present implemented of AliEventPoolManager
 155
 156     // centrality bins
 157     Double_t centralityBins[] = {0.,100.};
 158     Double_t* centbins        = centralityBins;
 159     Int_t nCentralityBins     = sizeof(centralityBins) / sizeof(Double_t) - 1;
 160
 161     // Zvtx bins
 162     Double_t vertexBins[] = {-10., 10.};
 163     Double_t* vtxbins     = vertexBins;
 164     Int_t nVertexBins     = sizeof(vertexBins) / sizeof(Double_t) - 1;
 165
 166     AliEventPoolManager *fPoolMgr = 0x0;
 167     fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins);
 168   }
 169   //========================================================//
 170
 171   //========================================================//
 172   //Create the TChain object: events
 173   TChain *eventChain = new TChain("events");
 174   StructEvent tStructEvents;
 175   eventChain->SetBranchAddress("event",&tStructEvents);
 176   //========================================================//
 177
 178   //========================================================//
 179   //Create the TChain object: particles
 180   TChain *particleChain = new TChain("particles");
 181   StructParticle tStructParticles;
 182   particleChain->SetBranchAddress("particle",&tStructParticles);
 183   //========================================================//
 184
 185   //========================================================//
 186   //Fill the TChain with the files in the directory
 187   for(Int_t iFile = 1; iFile <= 9; iFile++) {
 188     TString filename = aEventDir.Data();
 189     filename += "/event00"; filename += iFile;
 190     filename += ".root";
 191     cout<<"Adding file "<<filename.Data()<<" to the chain..."<<endl;
 192     eventChain->Add(filename.Data());
 193     particleChain->Add(filename.Data());
 194   }
 195   //========================================================//
 196
 197   //========================================================//
 198   //QA histograms
 199   //Event stats.
 200   TString gCutName[5] = {"Total","Offline trigger",
 201                          "Vertex","Analyzed","sel. Centrality"};
 202   TH1F *fHistEventStats = new TH1F("fHistEventStats",
 203                                    "Event statistics;;N_{events}",
 204                                    5,0.5,5.5);
 205   for(Int_t i = 1; i <= 5; i++)
 206     fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
 207   fList->Add(fHistEventStats);
 208
 209   //Number of accepted particles
 210   TH1F *fHistNumberOfAcceptedTracks = new TH1F("fHistNumberOfAcceptedTracks",";N_{acc.};Entries",4001,-0.5,4000.5);
 211   fList->Add(fHistNumberOfAcceptedTracks);
 212
 213   //Particle level: eta-phi
 214   TH2F *fHistEtaPhiPositive  = new TH2F("fHistEtaPhiPositive",";#eta;#varphi (rad)",80,gEtaMin,gEtaMax,72,0.0,2.*TMath::Pi());
 215   fList->Add(fHistEtaPhiPositive);
 216   TH2F *fHistEtaPhiNegative  = new TH2F("fHistEtaPhiNegative",";#eta;#varphi (rad)",80,gEtaMin,gEtaMax,72,0.0,2.*TMath::Pi());
 217   fList->Add(fHistEtaPhiNegative);
 218
 219   //Particle level: pt
 220   TH1F *fHistPtPositive  = new TH1F("fHistPtPositive",";p_{T} (GeV/c)",100,0.01,30.01);
 221   fList->Add(fHistPtPositive);
 222   TH1F *fHistPtNegative  = new TH1F("fHistPtNegative",";p_{T} (GeV/c)",100,0.01,30.01);
 223   fList->Add(fHistPtNegative);
 224   //========================================================//
 225
 226   //========================================================//
 227   //loop over the events
 228   Int_t nEvents = eventChain->GetEntries();
 229   cout<<"========================================="<<endl;
 230   cout<<"Number of events in the chain: "<<nEvents<<endl;
 231   cout<<"========================================="<<endl;
 232   Int_t iParticleCounter = 0;
 233   Int_t nTotalParticles = 0;
 234
 235   for(Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
 236     //for(Int_t iEvent = 0; iEvent < 1; iEvent++) {
 237     eventChain->GetEntry(iEvent);
 238
 239     //========================================================//
 240     //Create the TObjArray object to store the particles
 241     TObjArray* tracksAccepted = new TObjArray;
 242     tracksAccepted->SetOwner(kTRUE);
 243     //========================================================//
 244
 245     Int_t nParticles = tStructEvents.entries;
 246     if(((iEvent+1)%100)==0)
 247       cout<<"Event: "<<iEvent+1<<" - ID: "<<tStructEvents.eventID<<" - Entries: "<<tStructEvents.entries<<endl;
 248
 249     //========================================================//
 250     //Filling event level histos
 251     fHistEventStats->Fill(1);
 252     fHistEventStats->Fill(2);
 253     fHistEventStats->Fill(3);
 254     fHistEventStats->Fill(4);
 255     fHistEventStats->Fill(5);
 256
 257     Int_t gNumberOfAcceptedParticles = 0;
 258     Double_t gReactionPlane = 0.;
 259     Double_t gCharge = 0.;
 260     //========================================================//
 261     //loop over particles
 262     for(Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
 263       particleChain->GetEntry(nTotalParticles+iParticle);
 264
 265       //========================================================//
 266       //consider only primordial particles
 267       if(!IsPhysicalPrimary(tStructParticles.pid,
 268                             tStructParticles.fathereid,
 269                             gCharge)) continue;
 270
 271       //========================================================//
 272       //Calculate kinematic variables
 273       Double_t gPt = TMath::Sqrt(TMath::Power(tStructParticles.px,2) +
 274                                  TMath::Power(tStructParticles.py,2));
 275       Double_t gP = TMath::Sqrt(TMath::Power(tStructParticles.px,2) +
 276                                 TMath::Power(tStructParticles.py,2) +
 277                                 TMath::Power(tStructParticles.pz,2) );
 278       Double_t gEta = -100.;
 279       if(gP != tStructParticles.pz)
 280         gEta = 0.5*TMath::Log((gP + tStructParticles.pz)/(gP - tStructParticles.pz));
 281       Double_t gPhi = TMath::Pi()+TMath::ATan2(-tStructParticles.py,-tStructParticles.px);
 282
 283       //========================================================//
 284       //Fill QA
 285       if(gCharge > 0) {
 286         fHistEtaPhiPositive->Fill(gEta,gPhi);
 287         fHistPtPositive->Fill(gPt);
 288       }
 289       else if(gCharge < 0) {
 290         fHistEtaPhiNegative->Fill(gEta,gPhi);
 291         fHistPtNegative->Fill(gPt);
 292       }
 293       //========================================================//
 294
 295       //========================================================//
 296       //Apply cuts
 297       if((gEta > gEtaMax)||(gEta < gEtaMin)) continue;
 298       if((gPt > gPtMax)||(gPt < gPtMin)) continue;
 299
 300       tracksAccepted->Add(new AliBFBasicParticle(gEta,gPhi,gPt,gCharge, 1.));
 301       gNumberOfAcceptedParticles += 1;
 302
 303       iParticleCounter += 1;
 304       //cout<<"\t Particle counter: "<<iParticleCounter<<" - Particle in the event: "<<iParticle+1<<" - eventID: "<<tStructParticles.eventid<<" - pid: "<<tStructParticles.pid<<" - fatherpid: "<<tStructParticles.fatherpid<<" - rootpid: "<<tStructParticles.rootpid<<" - fathereid: "<<tStructParticles.fathereid<<" - eid: "<<tStructParticles.eid<<endl;
 305     }//particle loop
 306
 307     //========================================================//
 308     // Event mixing (borrowed code from the task)
 309     if (gRunMixing) {
 310       Int_t fMixingTracks = 50000;
 311       AliEventPool* pool = fPoolMgr->GetEventPool(1.,0.,0.);
 312
 313       if (!pool) {
 314         AliFatal(Form("No pool found"));
 315       }
 316       else {
 317         //pool->SetDebug(1);
 318         if (pool->IsReady() || pool->NTracksInPool() > fMixingTracks / 10 || pool->GetCurrentNEvents() >= 5){
 319
 320           Int_t nMix = pool->GetCurrentNEvents();
 321           //cout << "nMix = " << nMix << " tracks in pool = " << pool->NTracksInPool() << endl;
 322
 323           // Fill mixed-event histos here
 324           for (Int_t jMix=0; jMix<nMix; jMix++) {
 325             TObjArray* tracksMixed = pool->GetEvent(jMix);
 326             bfm->CalculateBalance(gReactionPlane,tracksAccepted,tracksMixed,1.,0.);
 327           }
 328         }
 329
 330         // Update the Event pool
 331         pool->UpdatePool(tracksAccepted);
 332         //pool->PrintInfo();
 333
 334       }//pool NULL check
 335     }//run mixing
 336       //========================================================//
 337
 338     //========================================================//
 339     // calculate balance function
 340     bf->CalculateBalance(gReactionPlane,tracksAccepted,NULL,1.,0.);
 341
 342     fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedParticles);
 343     nTotalParticles += nParticles;
 344   }
 345
 346   //========================================================//
 347   //Output file
 348   TFile *f = TFile::Open("AnalysisResults.root","recreate");
 349   fList->Write("listQA",TObject::kSingleKey);
 350   fListBF->Write("listBF",TObject::kSingleKey);
 351   if(gRunMixing) fListBFM->Write("listBFMixed",TObject::kSingleKey);
 352   f->Close();
 353   //========================================================//
 354
 355   // Print real and CPU time used for analysis:
 356   timer.Stop();
 357   timer.Print();
 358 }
 359
 360 //______________________________________________________________//
 361 Bool_t IsPhysicalPrimary(Int_t gPDGCode,
 362                          Int_t gFathereid,
 363                          Double_t &gCharge) {
 364   //Check whether the primordial particle belongs to the list
 365   //of known particles
 366   Bool_t kStatus = kFALSE;
 367   if(gFathereid != -1)
 368     return kStatus;
 369
 370   //List of stable particles
 371   const Int_t kNstable = 3;
 372   Int_t pdgStable[kNstable] = {
 373     211,         // Pion
 374     321,         // Kaon
 375     2212,        // Proton
 376   };
 377
 378   if(gPDGCode < 0) gCharge = -1;
 379   else if(gPDGCode > 0) gCharge = 1;
 380
 381   for(Int_t iParticle = 0; iParticle < kNstable; iParticle++) {
 382     if((TMath::Abs(gPDGCode) == pdgStable[iParticle])) {
 383       kStatus = kTRUE;
 384       return kStatus;
 385     }
 386   }
 387
 388   return kFALSE;
 389 }