Update to newest tree classes.
authorloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 21 Sep 2010 06:33:11 +0000 (06:33 +0000)
committerloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 21 Sep 2010 06:33:11 +0000 (06:33 +0000)
PWG4/TwoPartCorr/AutoCorr.C
PWG4/TwoPartCorr/TreeClasses.h
PWG4/TwoPartCorr/anaCorr.C
PWG4/TwoPartCorr/runAutoCorr.C

index dd2e3e7..257df0e 100644 (file)
@@ -68,8 +68,8 @@ Double_t AutoCorr::DeltaPhi(const MyPart &t1, const MyPart &t2,
 {
   Double_t dphi = -999;
   Double_t pi = TMath::Pi();
-  Double_t phia = t1.Phi();  
-  Double_t phib = t2.Phi();  
+  Double_t phia = t1.fPhi;  
+  Double_t phib = t2.fPhi;  
   
   if (phia < 0)         phia += 2*pi;
   else if (phia > 2*pi) phia -= 2*pi;
@@ -84,7 +84,7 @@ Double_t AutoCorr::DeltaPhi(const MyPart &t1, const MyPart &t2,
 
 Double_t AutoCorr::DeltaEta(const MyPart &t1, const MyPart &t2) const
 {
-  return t1.Eta() - t2.Eta();
+  return t1.fEta - t2.fEta;
 }
 
 Bool_t AutoCorr::InBounds(Double_t val, Double_t min, Double_t max) const
@@ -116,14 +116,14 @@ Bool_t AutoCorr::IsEventOk(const MyHeader &ev, Int_t minVc,
 
 Bool_t AutoCorr::IsTrackOk(const MyPart &t, Double_t etaMin, Double_t etaMax) const
 {
-  return InBounds(t.Eta(), etaMin, etaMax);    
+  return InBounds(t.fEta, etaMin, etaMax);    
 }
 
 Bool_t AutoCorr::IsTrackOk(const MyPart &t, Double_t etaMin, Double_t etaMax,
                           Double_t ptMin, Double_t ptMax) const
 {
-  Bool_t etaOk = InBounds(t.Eta(), etaMin, etaMax);
-  Bool_t ptOk  = InBounds(t.Pt(), ptMin, ptMax);    
+  Bool_t etaOk = InBounds(t.fEta, etaMin, etaMax);
+  Bool_t ptOk  = InBounds(t.fPt, ptMin, ptMax);    
   return  etaOk && ptOk;
 }
 
index 357afbc..49195b3 100644 (file)
 class MyHeader
 {
 public:
-  MyHeader() : fRun(0), fOrbit(0), fTime(0), fPeriod(0), fBx(0),
+  MyHeader() : fRun(0), fOrbit(0), fTime(0), fPeriod(0), fBx(0), fL0(0), fL1(0), fL2(0),
                fNChips1(0), fNChips2(0), fNTracks(0), fNSelTracks(0), fNTracklets(0), 
-               fVz(0), fVc(-1), fIsPileupSPD(0), fNPileupSPD(0), fNPileup(0),
+               fVx(0), fVy(0), fVz(0), fVc(-1), fIsPileupSPD(0), fNPileupSPD(0), fNPileup(0),
                fTrClassMask(0), fTrCluster(0), fEvNumberInFile(-1), fFileId(-1),
-               fVzSPD(0), fVcSPD(-1) {;}
+               fVxSPD(0), fVySPD(0), fVzSPD(0), fVcSPD(-1) {;}
   virtual ~MyHeader() {;}
   ULong64_t     GetEventId() const {
                   return ((ULong64_t)fBx+
@@ -27,18 +27,24 @@ public:
                           (ULong64_t)fPeriod*16777215*3564);
                 }
 
+public:
   Int_t         fRun;
   UInt_t        fOrbit; 
   UInt_t        fTime;   
   UInt_t        fPeriod;
   UShort_t      fBx;
+  UInt_t        fL0;
+  UInt_t        fL1;
+  UShort_t      fL2;
   Short_t       fNChips1;
   Short_t       fNChips2;
   Short_t       fNTracks;
   Short_t       fNSelTracks;
   Short_t       fNTracklets;
-  Double_t      fVz; //[0,0,16]
-  Double_t      fVc; //[0,0,16]
+  Double_t      fVx;             //[0,0,16]
+  Double_t      fVy;             //[0,0,16]
+  Double_t      fVz;             //[0,0,16]
+  Double_t      fVc;             //[0,0,16]
   Bool_t        fIsPileupSPD; 
   Char_t        fNPileupSPD;
   Char_t        fNPileup;
@@ -46,32 +52,73 @@ public:
   UChar_t       fTrCluster;
   Int_t         fEvNumberInFile;
   Short_t       fFileId;
-  Double_t      fVzSPD; //[0,0,16]
-  Double_t      fVcSPD; //[0,0,16]
+  Double_t      fVxSPD;          //[0,0,16]
+  Double_t      fVySPD;          //[0,0,16]
+  Double_t      fVzSPD;          //[0,0,16]
+  Double_t      fVcSPD;          //[0,0,16]
+  Double_t      fVxTPC;          //[0,0,16]
+  Double_t      fVyTPC;          //[0,0,16]
+  Double_t      fVzTPC;          //[0,0,16]
+  Double_t      fVcTPC;          //[0,0,16]
 
-  ClassDef(MyHeader,3) // My header class
+  ClassDef(MyHeader,4) // My header class
 };
 
 class MyPart : public TObject
 {
 public:
-  MyPart(ULong_t st=0, Char_t c=0, Double_t pt=0, Double_t eta=0, Double_t phi=0) :
-    TObject(), fSt(st), fC(c), fPt(pt), fEta(eta), fPhi(phi) {;}
-  Char_t      Charge() const { return fC;   }
-  Double_t    Eta()    const { return fEta; }
-  Double_t    Phi()    const { return fPhi; }
-  Double_t    Pt()     const { return fPt;  }
-  Double_t    Px()     const { return fPt*TMath::Cos(fPhi);  }
-  Double_t    Py()     const { return fPt*TMath::Sin(fPhi);  }
-  Double_t    Pz()     const { return fPt*TMath::SinH(fEta); }
-  ULong_t     Status() const { return fSt;  }
-protected:
+  enum { /*from AliESDtrack.h */
+    kITSin=0x0001,kITSout=0x0002,kITSrefit=0x0004,kITSpid=0x0008,
+    kTPCin=0x0010,kTPCout=0x0020,kTPCrefit=0x0040,kTPCpid=0x0080,
+    kTRDin=0x0100,kTRDout=0x0200,kTRDrefit=0x0400,kTRDpid=0x0800,
+    kTOFin=0x1000,kTOFout=0x2000,kTOFrefit=0x4000,kTOFpid=0x8000,
+    kHMPIDout=0x10000,kHMPIDpid=0x20000,
+    kEMCALmatch=0x40000,
+    kPHOSmatch=0x200000,
+    kTRDbackup =0x80000,
+    kTRDStop=0x20000000,
+    kESDpid=0x40000000,
+    kTIME=0x80000000,
+    kGlobalMerge=0x08000000,
+    kITSpureSA=0x10000000,
+    kMultInV0=0x2000000,    //BIT(25): assumed to be belong to V0 in multiplicity estimates
+    kMultSec=0x4000000     //BIT(26): assumed to be secondary (due to the DCA) in multiplicity estimates
+  }; 
+
+  MyPart(ULong_t st=0, Char_t c=0, Double_t pt=0, Double_t eta=0, Double_t phi=0,
+         Short_t ncltpc=0, Short_t ncltpc1=0, Short_t ncltpcs=0,Char_t nclits=0, Double_t chi2tpc=0,
+         Double_t chi2tpc1=0, Double_t chi2its=0, Double_t d=0, Double_t z=0, Double_t dtpc=0, Double_t ztpc=0) :
+    TObject(), fSt(st), fC(c), fPt(pt), fEta(eta), fPhi(phi), fNClTPC(ncltpc), fNClTPC1(ncltpc1), 
+    fNClTPCShared(ncltpcs),fNClITS(nclits), fChi2TPC(chi2tpc), fChi2TPC1(chi2tpc1), fChi2ITS(chi2its), 
+    fD(d), fZ(z), fDTPC(dtpc), fZTPC(ztpc) {;}
+
+
+  Double_t    Px()         const { return fPt*TMath::Cos(fPhi);  }
+  Double_t    Py()         const { return fPt*TMath::Sin(fPhi);  }
+  Double_t    Pz()         const { return fPt*TMath::SinH(fEta); }
+  Bool_t      IsITSRefit() const { return (fSt&kITSrefit)==0;    }
+  Bool_t      IsTPCIn()    const { return (fSt&kTPCin)==0;       }
+  Bool_t      IsTPCRefit() const { return (fSt&kITSrefit)==0;    }
+
+public:
   ULong_t     fSt;
   Char_t      fC;
-  Double32_t  fPt;  //[0,0,14]
-  Double32_t  fEta; //[0,0,10]
-  Double32_t  fPhi; //[0,0,10]
-  ClassDef(MyPart,1) // My particle class in cylindrical coordinates
+  Double32_t  fPt;           //[0,0,16]
+  Double32_t  fEta;          //[0,0,10]
+  Double32_t  fPhi;          //[0,0,10]
+  Short_t     fNClTPC;
+  Short_t     fNClTPC1;
+  Short_t     fNClTPCShared;
+  Char_t      fNClITS;
+  Double32_t  fChi2TPC;      //[0,0,10]
+  Double32_t  fChi2TPC1;      //[0,0,10]
+  Double32_t  fChi2ITS;      //[0,0,10]
+  Double32_t  fD;            //[0,0,16]
+  Double32_t  fZ;            //[0,0,16]
+  Double32_t  fDTPC;         //[0,0,16]
+  Double32_t  fZTPC;         //[0,0,16]
+
+  ClassDef(MyPart,2) // My particle class in cylindrical coordinates
 };
 
 class MyTracklet : public TObject
@@ -79,16 +126,13 @@ class MyTracklet : public TObject
 public:
   MyTracklet(Double_t dphi=0, Double_t dth=0, Double_t eta=0, Double_t phi=0) :
     TObject(), fDPhi(dphi), fDTh(dth), fEta(eta), fPhi(phi) {;}
-  Double_t    DPhi()   const { return fDPhi; }
-  Double_t    DTh()    const { return fDTh;  }
-  Double_t    Eta()    const { return fEta;  }
-  Double_t    Phi()    const { return fPhi;  }
-protected:
+
+public:
   Double32_t  fDPhi; //[0,0,10]
   Double32_t  fDTh;  //[0,0,10]
   Double32_t  fEta;  //[0,0,10]
   Double32_t  fPhi;  //[0,0,10]
+
   ClassDef(MyTracklet,1) // My tracklet class
 };
-
 #endif
index c059f7a..a404493 100644 (file)
@@ -104,15 +104,15 @@ void anaCorr(Int_t nEvents,
       Int_t ntracks = tracks->GetEntries();
       for (Int_t t1=0;t1<ntracks;++t1) {
         MyPart *part1 = (MyPart*)tracks->At(t1);
-        if (!InBounds(part1->Pt(),ptmin1,ptmax1))
+        if (!InBounds(part1->fPt,ptmin1,ptmax1))
           continue;
-        if (!InBounds(part1->Eta(),etamin1,etamax1))
+        if (!InBounds(part1->fEta,etamin1,etamax1))
           continue;
         for (Int_t t2=t1+1;t2<ntracks;++t2) {
           MyPart *part2 = (MyPart*)tracks->At(t2);
-          if (!InBounds(part2->Pt(),ptmin2,ptmax2))
+          if (!InBounds(part2->fPt,ptmin2,ptmax2))
             continue;
-          if (!InBounds(part2->Eta(),etamin2,etamax2))
+          if (!InBounds(part2->fEta,etamin2,etamax2))
             continue;
           //if (InvMass(*part1,*part2)<0.001)
           //  continue;
@@ -125,9 +125,9 @@ void anaCorr(Int_t nEvents,
           Int_t mtracks = nmix;
           for (Int_t t3=0;t3<mtracks;) {
             MyPart *part2 = pool->GetRandomTrack();
-            if (!InBounds(part2->Pt(),ptmin2,ptmax2))
+            if (!InBounds(part2->fPt,ptmin2,ptmax2))
               continue;
-            if (!InBounds(part2->Eta(),etamin2,etamax2))
+            if (!InBounds(part2->fEta,etamin2,etamax2))
               continue;
             //if (InvMass(*part1,*part2)<0.001)
             //  continue;
@@ -159,7 +159,7 @@ void anaCorr(Int_t nEvents,
   TCanvas *c4 = new TCanvas("cMix2");
   hMix2->Draw("surf1");
   TTimeStamp t;
-  TString fname(Form("histout-%d.root",t.GetSec()));
+  TString fname(Form("histout-%d.root",(Int_t)t.GetSec()));
   TFile outfile(fname,"recreate");
   hSig->Write();
   hBkg->Write();
@@ -183,8 +183,8 @@ Double_t DeltaPhi(const MyPart &t1, const MyPart &t2,
 {
   Double_t dphi = -999;
   Double_t pi = TMath::Pi();
-    Double_t phia = t1.Phi();  
-  Double_t phib = t2.Phi();  
+    Double_t phia = t1.fPhi;  
+  Double_t phib = t2.fPhi;  
   
   if (phia < 0)         phia += 2*pi;
   else if (phia > 2*pi) phia -= 2*pi;
@@ -199,7 +199,7 @@ Double_t DeltaPhi(const MyPart &t1, const MyPart &t2,
 
 Double_t DeltaEta(const MyPart &t1, const MyPart &t2)
 {
-  return t1.Eta() - t2.Eta();
+  return t1.fEta - t2.fEta;
 }
 
 Double_t InvMass(const MyPart &p1, const MyPart &p2)
index 221ebf8..70dfdd1 100644 (file)
@@ -23,38 +23,38 @@ void runAutoCorr(const char* dataFile =
   const double PI = TMath::Pi();
   Int_t poolsize = 50;
   Int_t nMix = 5;
-  const Int_t nMultBins = 2;
-  Double_t multbin[nMultBins+1] = {100, 200, 400}; 
-  const Int_t nZvtxBins = 2;
-  Double_t zvtxbin[nZvtxBins+1] = {-15, 0, 15};
-  const Int_t nPtBins = 3;
-  Double_t ptbin[nPtBins+1] = {0.15, 0.50, 2.0, 10.};
+  const Int_t nMultBins = 1;
+  Double_t multbin[nMultBins+1] = {100, 400}; 
+  const Int_t nZvtxBins = 3;
+  Double_t zvtxbin[nZvtxBins+1] = {-10, -3, 4, 11};
+  const Int_t nPtBins = 1;
+  Double_t ptbin[nPtBins+1] = {1.0, 2.0};
   Double_t etaMin = -1.4; // for track cuts
   Double_t etaMax = 1.4;
-  Double_t ptMin = 0.150; // for track cuts
+  Double_t ptMin = 0.40; // for track cuts
   Double_t ptMax = 10.0;
-  Double_t dEtaMax = 2.5; // for histogram limits
+  Double_t dEtaMax = 3.0; // for histogram limits
   TString sDataSet = "res_LHC10e_09122010";
 
   // Event cuts
   Int_t minVc = 25;
   Int_t maxNTracklets = 200;
-  Double_t zMin = -15.0;
-  Double_t zMax = 15.0;
+  Double_t zMin = zvtxbin[0];
+  Double_t zMax = zvtxbin[nZvtxBins];
   Int_t runNumber = -1;
   Int_t firstRunNumber = -1;
 
   TString sMultBins(""), sZvtxBins(""), sPtBins("");
-  for (int i=0; i<nMultBins+1; i++) sMultBins.Append(Form("%d ", (Int_t)multbin[i]));
-  for (int i=0; i<nZvtxBins+1; i++) sZvtxBins.Append(Form("%d ", (Int_t)zvtxbin[i]));
-  for (int i=0; i<nPtBins+1;   i++)   sPtBins.Append(Form("%d ", (Int_t)ptbin[i]));
+  for (int i=0; i<nMultBins+1; i++) sMultBins.Append(Form("%.0f ", multbin[i]));
+  for (int i=0; i<nZvtxBins+1; i++) sZvtxBins.Append(Form("%.1f ", zvtxbin[i]));
+  for (int i=0; i<nPtBins+1;   i++)   sPtBins.Append(Form("%.2f ", ptbin[i]));
 
   // Document the applied cuts, get quick look with cuts->Print()
   TList* cuts = new TList();
   cuts->Add(new TNamed("data_sample", sDataSet.Data()));
   cuts->Add(new TNamed("mult_bins", sMultBins.Data()));
   cuts->Add(new TNamed("zvtx_bins", sZvtxBins.Data()));
-  cuts->Add(new TNamed("mult_bins", sPtBins.Data()));
+  cuts->Add(new TNamed("pt_bins", sPtBins.Data()));
   cuts->Add(new TNamed("mass_cut", "No mass cut"));
   cuts->Add(new TNamed("eta_range", Form("%.1f < #eta < %.1f", etaMin, etaMax)));
   cuts->Add(new TNamed("z_vtx_range", Form("%.1f < z-vtx < %.1f", zMin, zMax)));
@@ -118,31 +118,51 @@ void runAutoCorr(const char* dataFile =
   TH1F* hZvtx = new TH1F("hZvtx", "Event Z-vertex", 60, -30, 30);
   TH1F* hMass = new TH1F("hMass", "hMass", 1000, 0, 10);
   TH1F* hMassBg = new TH1F("hMassBg", "hMassBg", 1000, 0, 10);
+  TH1F* hNevts = new TH1F("hNevts", "# events: available, sampled, accepted", 3, 0, 3);
+
+  //----------------------
+ TString sMultBin[nMultBins];
+  for (int i=0; i<nMultBins; i++) {
+    double m1 = hMultBins->GetBinLowEdge(i+1);
+    double m2 = m1 + hMultBins->GetBinWidth(i+1);
+    sMultBin[i] = Form("%.0f - %.0f tracks", m1, m2);
+  }
+  TString sPtBin[nPtBins];
+  for (int i=0; i<nPtBins; i++) {
+    double pt1 = hPtBins->GetBinLowEdge(i+1);
+    double pt2 = pt1 + hPtBins->GetBinWidth(i+1);
+    sPtBin[i] = Form("%.2f - %.2f GeV/c", pt1, pt2);
+  }
+  //---------------------
 
   for (int iM=0; iM<nMultBins; iM++) {
     char* name;
+    char* title;
     name = Form("hMult_%i", iM);
-    hMult[iM] = new TH1F(name, name, 1000, 0, 1000);
+    title = Form("%s, %.2f - %.2f GeV/c", sMultBin[iM].Data(), ptMin, ptMax);
+
+    hMult[iM] = new TH1F(name, title, 1000, 0, 1000);
 
     name = Form("hSig_%i", iM);
-    hSig[iM] = new TH2F(name, name, 60, -PI/2, 3*PI/2, 60, -dEtaMax, dEtaMax);
+    hSig[iM] = new TH2F(name, title, 60, -PI/2, 3*PI/2, 60, -dEtaMax, dEtaMax);
     hSig[iM]->Sumw2();
     name = Form("hBkg_%i", iM);
-    hBkg[iM] = new TH2F(name, name, 60, -PI/2, 3*PI/2, 60, -dEtaMax, dEtaMax);
+    hBkg[iM] = new TH2F(name, title, 60, -PI/2, 3*PI/2, 60, -dEtaMax, dEtaMax);
     hBkg[iM]->Sumw2();
     name = Form("hSB_%i", iM);
-    hSB[iM] = new TH2F(name, name, 60, -PI/2, 3*PI/2, 60, -dEtaMax, dEtaMax);
+    hSB[iM] = new TH2F(name, title, 60, -PI/2, 3*PI/2, 60, -dEtaMax, dEtaMax);
     hSB[iM]->Sumw2();
 
     for (int iPt=0; iPt<nPtBins; iPt++) {
       name = Form("hSigPt_%i_%i", iM, iPt);
-      hSigPt[iM][iPt] = new TH2F(name, name, 60, -PI/2, 3*PI/2, 60, -dEtaMax, dEtaMax);
+      title = Form("%s, %s", sMultBin[iM].Data(), sPtBin[iPt].Data());
+      hSigPt[iM][iPt] = new TH2F(name, title, 60, -PI/2, 3*PI/2, 60, -dEtaMax, dEtaMax);
       hSigPt[iM][iPt]->Sumw2();
       name = Form("hBkgPt_%i_%i", iM, iPt);
-      hBkgPt[iM][iPt] = new TH2F(name, name, 60, -PI/2, 3*PI/2, 60, -dEtaMax, dEtaMax);
+      hBkgPt[iM][iPt] = new TH2F(name, title, 60, -PI/2, 3*PI/2, 60, -dEtaMax, dEtaMax);
       hBkgPt[iM][iPt]->Sumw2();
       name = Form("hSBPt_%i_%i", iM, iPt);
-      hSBPt[iM][iPt] = new TH2F(name, name, 60, -PI/2, 3*PI/2, 60, -dEtaMax, dEtaMax);
+      hSBPt[iM][iPt] = new TH2F(name, title, 60, -PI/2, 3*PI/2, 60, -dEtaMax, dEtaMax);
       hSBPt[iM][iPt]->Sumw2();
     }
     
@@ -160,7 +180,7 @@ void runAutoCorr(const char* dataFile =
                    "accepted/sampled tracks");
   // ---------------------------------
 
-  Int_t nEvents = tree->GetEntries();
+  Int_t nEvents = 400;//tree->GetEntries();
   Int_t nAcceptedEvents = 0;
 
   TStopwatch* watch = new TStopwatch();
@@ -202,13 +222,13 @@ void runAutoCorr(const char* dataFile =
     if (pool->IsPoolReady()) {
       for (int i=0; i<ntracks; i++) {
        MyPart* t1 = (MyPart*)trk->At(i);
-       hPtAll->Fill(t1->Pt());
+       hPtAll->Fill(t1->fPt);
        if (!ac->IsTrackOk(*t1, etaMin, etaMax, ptMin, ptMax))
          continue;
        
-       hEta->Fill(t1->Eta());
-       hPtSel->Fill(t1->Pt());
-       Int_t ptBin1 = hPtBins->FindBin(t1->Pt()) - 1;
+       hEta->Fill(t1->fEta);
+       hPtSel->Fill(t1->fPt);
+       Int_t ptBin1 = hPtBins->FindBin(t1->fPt) - 1;
 
        // --- Same-event correlation loop ---
        for (int j=i+1; j<ntracks; j++) {
@@ -228,7 +248,7 @@ void runAutoCorr(const char* dataFile =
            hDRcut->Fill(dphi, deta); // Show elliptical hole from cut 
            hMass->Fill(ac->InvMass(*t1, *t2));
 
-           Int_t ptBin2 = hPtBins->FindBin(t2->Pt()) - 1;
+           Int_t ptBin2 = hPtBins->FindBin(t2->fPt) - 1;
            if (ptBin1==ptBin2)
              hSigPt[multBin][ptBin1]->Fill(dphi, deta);
          }
@@ -251,7 +271,7 @@ void runAutoCorr(const char* dataFile =
              Double_t dphi_bg = ac->DeltaPhi(*t1, *tbg);
              hBkg[multBin]->Fill(dphi_bg, deta_bg);
 
-             Int_t ptBin3 = hPtBins->FindBin(tbg->Pt()) - 1;
+             Int_t ptBin3 = hPtBins->FindBin(tbg->fPt) - 1;
              if (ptBin1==ptBin3)
                hBkgPt[multBin][ptBin3]->Fill(dphi_bg, deta_bg);
              
@@ -293,7 +313,10 @@ void runAutoCorr(const char* dataFile =
       hSBPt[iM][iPt]->Divide(hSigPt[iM][iPt], hBkgPt[iM][iPt]);  
     }
   }
-  
+  hNevts->SetBinContent(1, tree->GetEntries());
+  hNevts->SetBinContent(2, nEvents);
+  hNevts->SetBinContent(3, nAcceptedEvents);
+
   outFile->Write();
   cuts->Write("cuts", TObject::kSingleKey);
   outFile->Close();