FMD/scripts/AliFMDAnaFlow.C

   1 #include <TH1D.h>
   2 #include <TH2D.h>
   3 #include <TBrowser.h>
   4 #include <TArrayF.h>
   5 #include <TCanvas.h>
   6 #include <TStyle.h>
   7 #include <TLegend.h>
   8 #include <TProfile.h>
   9 #include <TProfile2D.h>
  10 #include <TROOT.h>
  11 #include "AliFMDAnaESD.h"
  12 #include "AliFMDAnaFlowRing.h"
  13 #include <TFile.h>
  14 #if 1
  15 # include "background.C"
  16 // # include "background2.C"
  17 #else
  18 TH2* bgFMD1i() { return 0; }
  19 TH2* bgFMD2i() { return 0; }
  20 TH2* bgFMD2o() { return 0; }
  21 TH2* bgFMD3o() { return 0; }
  22 TH2* bgFMD3o() { return 0; }
  23 #endif
  24 #if 1
  25 # include "reactionplane.C"
  26 #else
  27 Float_t reactionplane(int ev) { return 0; }
  28 #endif
  29 #include <iostream>
  30
  31
  32 //====================================================================
  33 class AliFMDAnaFlow : public AliFMDAnaESD
  34 {
  35 public:
  36   //__________________________________________________________________
  37   /** Constructor */
  38   AliFMDAnaFlow(Int_t n=2, Bool_t bg=false, Float_t cut0=.68)
  39     : fAxis(n, -3.7, 5.3),
  40       fV2("all", "Full FMD", 2, 1, fAxis)
  41   {
  42     for (int i = 0; i < 5; i++) {
  43       UShort_t det = (i == 0 ? 1 : (i <= 2 ? 2 : 3));
  44       Char_t   rng = (i == 0 ? 'I' :  (i % 2 == 1 ? 'I' : 'O'));
  45       AddRing(new AliFMDAnaFlowRing(det,rng,n,EtaMin(i)+.1,EtaMax(i)-.1,
  46                            Bg(bg,i),cut0,Cut1(i)));
  47     }
  48   }
  49   //__________________________________________________________________
  50   /** Constructor */
  51   AliFMDAnaFlow(Float_t deta, Bool_t bg=false, Float_t cut0=.68)
  52     : fAxis(Int_t((5.3- -3.7)/deta), -3.7, 5.3),
  53       fV2("all", "Full FMD", 2, 1, fAxis)
  54   {
  55     for (int i = 0; i < 5; i++) {
  56       UShort_t det = (i == 0 ? 1 : i <= 2 ? 2 : 3);
  57       Char_t   rng = (i % 2 == 0 ? 'I' : 'O');
  58       AddRing(new AliFMDAnaFlowRing(det,rng,fAxis,Bg(bg,i),cut0,Cut1(i)));
  59     }
  60   }
  61   //__________________________________________________________________
  62   /** Initialize */
  63   Bool_t Init()
  64   {
  65     fPhis.Set(51200);
  66     fEtas.Set(51200);
  67     fWeights.Set(51200);
  68     fNObserv = 0;
  69
  70     Float_t pi2 = 2 * TMath::Pi();
  71
  72     fPsi  = new TH2D("all_psi", "Event planes",
  73                      fAxis.N(), fAxis.Bins(), 40, 0, pi2);
  74     fPsi->SetXTitle("#eta");
  75     fPsi->SetYTitle("#Psi");
  76     fPsi->SetZTitle("Frequency");
  77     fPsi->SetMarkerStyle(20);
  78
  79     fDPsi = new TProfile("all_dpsi", "Distance to real #Psi_{R}",
  80                          fAxis.N(), fAxis.Bins(), -pi2, pi2);
  81     fDPsi->SetXTitle("#eta");
  82     fDPsi->SetYTitle("#Delta#Psi=|#Psi-#Psi_{R}|");
  83
  84     fDPhi = new TProfile2D("all_dphi","Distribution of phi signals",
  85                           fAxis.N(), -3.7, 5.3, 40, 0, pi2, 0, 20);
  86     fDPhi->SetXTitle("#eta");
  87     fDPhi->SetYTitle("#varphi");
  88     fDPhi->SetZTitle("signal");
  89     fDPhi->SetMarkerStyle(20);
  90
  91     return AliFMDAnaESD::Init();
  92   }
  93   //__________________________________________________________________
  94   /** Called at the beginning of an event */
  95   Bool_t Begin(Int_t ev)
  96   {
  97     fNObserv = 0;
  98     Float_t psir = reactionplane(fNEvents);
  99     AliFMDAnaFlowRing* fring = 0;
 100     for (size_t i = 0; i < 5; i++) {
 101       if (!fRing[i]) continue;
 102       fring = static_cast<AliFMDAnaFlowRing*>(fRing[i]);
 103       fring->fTrue.SetPsi(psir);
 104     }
 105     return AliFMDAnaESD::Begin(ev);
 106   }
 107   //__________________________________________________________________
 108   /** Fill in one strip */
 109   void Fill(Float_t phi, Float_t eta, Float_t mult)
 110   {
 111     Int_t imult = Int_t(mult+.33);
 112     if (fNObserv + imult >= fEtas.fN) {
 113       ULong_t n = Int_t(fEtas.fN * 1.5)+imult;
 114       fEtas.Set(n);
 115       fPhis.Set(n);
 116       fWeights.Set(n);
 117     }
 118     for (Int_t i = 0; i < imult; i++) {
 119       fEtas[fNObserv]    = eta;
 120       fPhis[fNObserv]    = phi;
 121       fWeights[fNObserv] = 1;
 122       fNObserv++;
 123     }
 124     fDPhi->Fill(eta, phi);
 125   }
 126   //__________________________________________________________________
 127   /** Called at the end of an event */
 128   Bool_t End()
 129   {
 130     float       psir    = reactionplane(fNEvents-1);
 131     // std::cout << "PsiR=" << psir << std::endl;
 132     fV2.Event(fNObserv, fPhis.fArray, fEtas.fArray,
 133               fWeights.fArray, fWeights.fArray);
 134     for (UShort_t i = 0; i < fAxis.N(); i++) {
 135       AliFMDFlowBin* bin = fV2.GetBin(i);
 136       if (bin && bin->Counts() > 0) {
 137         fPsi->Fill(fAxis.BinCenter(i),  bin->Psi() - psir);
 138         fDPsi->Fill(fAxis.BinCenter(i), bin->Psi() - psir);
 139       }
 140     }
 141     return AliFMDAnaESD::End();
 142   }
 143   //__________________________________________________________________
 144   /** Make a frame */
 145   TH1* MakeFrame(const char* xname, const char* yname,
 146                  Float_t ymin, Float_t ymax)
 147   {
 148     static int id = 0;
 149     id++;
 150     TH1* frame = new TH1D(Form("frame%d", id), "Frame",
 151                           fAxis.N(), fAxis.Bins()[0], fAxis.Bins()[fAxis.N()]);
 152     frame->SetMinimum(ymin);
 153     frame->SetMaximum(ymax);
 154     frame->SetXTitle(xname);
 155     frame->SetYTitle(yname);
 156     frame->SetLabelFont(132, "XYZ");
 157     frame->SetTitleFont(132, "XYZ");
 158     frame->SetLabelSize(.08, "XYZ");
 159     frame->SetTitleSize(.08, "XYZ");
 160     frame->SetTitleOffset(.5, "XYZ");
 161     frame->SetNdivisions(8, "XYZ");
 162     frame->Draw();
 163     return frame;
 164   }
 165   //__________________________________________________________________
 166   /** called at the end of a run */
 167   Bool_t Finish()
 168   {
 169     gStyle->SetPalette(1);
 170     gStyle->SetOptTitle(0);
 171     gStyle->SetOptStat(0);
 172     AliFMDAnaFlowRing* fring[5];
 173     for (size_t i = 0; i < 5; i++)
 174       fring[i] = static_cast<AliFMDAnaFlowRing*>(fRing[i]);
 175
 176     TCanvas* c = new TCanvas("c1", "Result", 600, 750);
 177     c->cd();
 178     c->SetGridy();
 179     c->SetTopMargin(0.02);
 180     c->SetRightMargin(0.02);
 181     c->SetFillColor(0);
 182     c->SetBorderSize(0);
 183     c->SetBorderMode(0);
 184     c->SetHighLightColor(0);
 185     c->Divide(1,3,0,0);
 186
 187     TVirtualPad* p = 0;
 188     TH1*         f = 0;
 189
 190     p = c->cd(1);
 191     p->SetGridx();
 192     p->SetGridy();
 193     p->SetFillColor(0);
 194     f = MakeFrame("#eta", "|#Psi_{2}-#Psi_{R}|", -TMath::Pi(), TMath::Pi());
 195     for (int i = 0; i < 5; i++) fring[i]->fDPsi.Draw("same");
 196     fDPsi->Draw("same");
 197
 198     p = c->cd(2);
 199     p->SetGridx();
 200     p->SetGridy();
 201     p->SetFillColor(0);
 202     f = MakeFrame("#eta", "R_{k}", .7, 1);
 203     for (int i = 0; i < 5; i++) fring[i]->fV2.Draw("tr:same");
 204     fV2.Draw("tr:same");
 205
 206     p = c->cd(3);
 207     p->SetGridx();
 208     p->SetGridy();
 209     p->SetFillColor(0);
 210     f = MakeFrame("#eta", "v_{2}", 0, .1);
 211     for (int i = 0; i < 5; i++) {
 212       fring[i]->fV2.Draw("th:same");
 213       fring[i]->fTrue.Draw("th:same");
 214     }
 215     fV2.Draw("th:same");
 216
 217
 218     for (int i = 0; i < 5; i++) {
 219       fring[i]->fV2.Print("s");
 220       fring[i]->fTrue.Print("s");
 221     }
 222     fV2.Print("s");
 223     gROOT->GetListOfSpecials()->Add(this);
 224     return kTRUE;
 225   }
 226   void ToFile(const char* file) {
 227     TFile* out = TFile::Open(file, "RECREATE");
 228     for (int i = 0; i < 5; i++) fRing[i]->Write(fRing[i]->Name());
 229     fV2.Write();
 230     if (fPsi) fPsi->Write();
 231     if (fDPsi) fDPsi->Write();
 232     if (fDPhi) fDPhi->Write();
 233     out->Close();
 234   }
 235
 236
 237
 238   //__________________________________________________________________
 239   /** Get min pseudo rapidity of ring idx */
 240   Float_t EtaMin(Int_t idx) const {
 241     switch (idx) {
 242     case 0: return  3.61728; break; // FMD1I
 243     case 1: return  2.18122; break; // FMD2I
 244     case 2: return  1.79821; break; // FMD2O
 245     case 3: return -3.39913; break; // FMD3I
 246     case 4: return -2.28925; break; // FMD3O
 247     }
 248     return 0;
 249   }
 250   //__________________________________________________________________
 251   /** Get max pseudo rapidity of ring idx */
 252   Float_t EtaMax(Int_t idx) const {
 253     switch (idx) {
 254     case 0: return  5.02643; break; // FMD1I
 255     case 1: return  3.57899; break; // FMD2I
 256     case 2: return  2.39085; break; // FMD2O
 257     case 3: return -2.00644; break; // FMD3I
 258     case 4: return -1.70080; break; // FMD3O
 259     }
 260     return 0;
 261   }
 262   //__________________________________________________________________
 263   /** Get cut 1 of ring idx */
 264   Float_t Cut1(Int_t idx) const {
 265     switch (idx) {
 266     case 0: return 1.33857; break; // FMD1I
 267     case 1: return 1.30464; break; // FMD2I
 268     case 2: return 1.27263; break; // FMD2O
 269     case 3: return 1.28052; break; // FMD3I
 270     case 4: return 1.23041; break; // FMD3O
 271     }
 272     return 0;
 273   }
 274   //__________________________________________________________________
 275   /** Get background correction of ring idx */
 276   TH2* Bg(Bool_t use, Int_t idx) const {
 277     if (!use) return 0;
 278     switch (idx) {
 279     case 0: return bgFMD1i(); break;
 280     case 1: return bgFMD2i(); break;
 281     case 2: return bgFMD2o(); break;
 282     case 3: return bgFMD3i(); break;
 283     case 4: return bgFMD3o(); break;
 284     }
 285     return 0;
 286   }
 287   /** Flow axis */
 288   AliFMDFlowAxis fAxis;
 289   /** Flow object */
 290   AliFMDFlowBinned1D fV2;
 291   /** Cache of phi */
 292   TArrayD            fPhis;
 293   /** Cache of eta */
 294   TArrayD            fEtas;
 295   /** Cache of weights */
 296   TArrayD            fWeights;
 297   /** Number of observables */
 298   Int_t              fNObserv;
 299   /** Histogram of Psi */
 300   TH2D*              fPsi;       // Event planes
 301   /** Histogram of dPsi */
 302   TProfile*          fDPsi;   // Mean distance to Psi
 303   /** Histogram of dPhi */
 304   TProfile2D*        fDPhi;
 305   ClassDef(AliFMDAnaFlow,1)
 306 };
 307
 308 //____________________________________________________________________
 309 //
 310 // EOF
 311 //