Float_t multV0C=esdV0->GetMTotV0C();
multV0 = multV0A+multV0C;
}
+
if (fIsMC) multV0 = 0.85871 * multV0;
if (multV0 < fMultMin) return kFALSE;
if (multV0 > fMultMax) return kFALSE;
}
}
-
-Float_t AliAnalysisMultPbCentralitySelector::GetCorrV0(const AliESDEvent* esd, float &v0CorrResc) const
+//________________________________________________________________________
+Float_t AliAnalysisMultPbCentralitySelector::GetCorrV0(const AliESDEvent* esd, float &v0CorrResc) const
{
// correct V0 non-linearity, prepare a version rescaled to SPD2 corr
- const Double_t par0[64] = { 6.71e-02 , 6.86e-02 , 7.06e-02 , 6.32e-02 ,
- 5.91e-02 , 6.07e-02 , 5.78e-02 , 5.73e-02 , 5.91e-02 , 6.22e-02 ,
- 5.90e-02 , 6.11e-02 , 5.55e-02 , 5.29e-02 , 5.19e-02 , 5.56e-02 ,
- 6.25e-02 , 7.03e-02 , 5.64e-02 , 5.81e-02 , 4.57e-02 , 5.30e-02 ,
- 5.13e-02 , 6.43e-02 , 6.27e-02 , 6.48e-02 , 6.07e-02 , 1.01e-01 ,
- 6.68e-02 , 7.16e-02 , 6.36e-02 , 5.95e-02 , 2.52e-02 , 2.82e-02 ,
- 2.56e-02 , 2.86e-02 , 2.82e-02 , 2.10e-02 , 2.13e-02 , 2.32e-02 ,
- 2.75e-02 , 4.34e-02 , 3.78e-02 , 4.52e-02 , 4.11e-02 , 3.89e-02 ,
- 4.10e-02 , 3.73e-02 , 4.51e-02 , 5.07e-02 , 5.42e-02 , 4.74e-02 ,
- 4.33e-02 , 4.44e-02 , 4.64e-02 , 3.01e-02 , 6.38e-02 , 5.26e-02 ,
- 4.99e-02 , 5.26e-02 , 5.47e-02 , 3.84e-02 , 5.00e-02 , 5.20e-02 };
- const Double_t par1[64] = { -6.68e-05 , -7.78e-05 , -6.88e-05 , -5.92e-05 ,
- -2.43e-05 , -3.54e-05 , -2.91e-05 , -1.99e-05 , -1.40e-05 , -4.01e-05 ,
- -2.29e-05 , -3.68e-05 , -2.53e-05 , -2.44e-06 , -9.22e-06 , -1.51e-05 ,
- -2.80e-05 , -2.34e-05 , -1.72e-05 , -1.81e-05 , -1.29e-05 , -2.65e-05 ,
- -1.61e-05 , -2.86e-05 , -1.74e-05 , -4.23e-05 , -3.41e-05 , -1.05e-04 ,
- -2.76e-05 , -4.71e-05 , -3.06e-05 , -2.32e-05 , -1.55e-06 , 2.15e-05 ,
- 1.40e-05 , 2.16e-05 , 1.21e-05 , 3.05e-06 , 1.67e-05 , -3.84e-06 ,
- 3.09e-06 , 1.50e-05 , 3.47e-06 , 4.87e-06 , -3.71e-07 , -1.75e-06 ,
- -1.80e-06 , 9.99e-06 , -6.46e-06 , -4.91e-06 , 1.33e-05 , -2.52e-07 ,
- -3.85e-06 , 4.94e-06 , -2.48e-07 , -1.20e-05 , 2.07e-06 , 6.12e-06 ,
- -1.18e-06 , 4.54e-06 , -1.54e-05 , -1.25e-05 , 1.46e-06 , -6.67e-06 };
- const Double_t par2[64] = { 1.29e-08 , 1.51e-08 , 1.43e-08 , 1.11e-08 ,
- 5.04e-09 , 6.99e-09 , 5.58e-09 , 4.15e-09 , 4.00e-09 , 8.22e-09 ,
- 4.97e-09 , 7.66e-09 , 4.91e-09 , 1.10e-09 , 2.64e-09 , 3.64e-09 ,
- 5.76e-09 , 5.46e-09 , 3.38e-09 , 3.47e-09 , 2.43e-09 , 4.13e-09 ,
- 2.80e-09 , 5.80e-09 , 3.86e-09 , 7.46e-09 , 5.98e-09 , 2.58e-08 ,
- 5.50e-09 , 8.72e-09 , 5.23e-09 , 4.37e-09 , 2.33e-09 , -6.01e-10 ,
- 3.99e-11 , -2.02e-10 , 7.67e-10 , 2.03e-09 , 1.17e-10 , 2.56e-09 ,
- 1.16e-09 , -4.75e-10 , 1.28e-09 , 1.23e-09 , 1.62e-09 , 1.61e-09 ,
- 1.93e-09 , 2.97e-10 , 2.21e-09 , 2.16e-09 , 5.22e-10 , 1.03e-09 ,
- 1.56e-09 , 5.00e-10 , 1.01e-09 , 2.93e-09 , 1.05e-09 , 9.96e-11 ,
- 1.21e-09 , 7.45e-10 , 3.07e-09 , 2.31e-09 , 6.70e-10 , 1.89e-09 };
+ Double_t *par0;
+ Double_t *par1;
+ Double_t *par2;
+
+ Double_t par0_137161[64] = { 6.71e-02 , 6.86e-02 , 7.06e-02 , 6.32e-02 ,
+ 5.91e-02 , 6.07e-02 , 5.78e-02 , 5.73e-02 , 5.91e-02 , 6.22e-02 ,
+ 5.90e-02 , 6.11e-02 , 5.55e-02 , 5.29e-02 , 5.19e-02 , 5.56e-02 ,
+ 6.25e-02 , 7.03e-02 , 5.64e-02 , 5.81e-02 , 4.57e-02 , 5.30e-02 ,
+ 5.13e-02 , 6.43e-02 , 6.27e-02 , 6.48e-02 , 6.07e-02 , 1.01e-01 ,
+ 6.68e-02 , 7.16e-02 , 6.36e-02 , 5.95e-02 , 2.52e-02 , 2.82e-02 ,
+ 2.56e-02 , 2.86e-02 , 2.82e-02 , 2.10e-02 , 2.13e-02 , 2.32e-02 ,
+ 2.75e-02 , 4.34e-02 , 3.78e-02 , 4.52e-02 , 4.11e-02 , 3.89e-02 ,
+ 4.10e-02 , 3.73e-02 , 4.51e-02 , 5.07e-02 , 5.42e-02 , 4.74e-02 ,
+ 4.33e-02 , 4.44e-02 , 4.64e-02 , 3.01e-02 , 6.38e-02 , 5.26e-02 ,
+ 4.99e-02 , 5.26e-02 , 5.47e-02 , 3.84e-02 , 5.00e-02 , 5.20e-02 };
+ Double_t par1_137161[64] = { -6.68e-05 , -7.78e-05 , -6.88e-05 , -5.92e-05 ,
+ -2.43e-05 , -3.54e-05 , -2.91e-05 , -1.99e-05 , -1.40e-05 , -4.01e-05 ,
+ -2.29e-05 , -3.68e-05 , -2.53e-05 , -2.44e-06 , -9.22e-06 , -1.51e-05 ,
+ -2.80e-05 , -2.34e-05 , -1.72e-05 , -1.81e-05 , -1.29e-05 , -2.65e-05 ,
+ -1.61e-05 , -2.86e-05 , -1.74e-05 , -4.23e-05 , -3.41e-05 , -1.05e-04 ,
+ -2.76e-05 , -4.71e-05 , -3.06e-05 , -2.32e-05 , -1.55e-06 , 2.15e-05 ,
+ 1.40e-05 , 2.16e-05 , 1.21e-05 , 3.05e-06 , 1.67e-05 , -3.84e-06 ,
+ 3.09e-06 , 1.50e-05 , 3.47e-06 , 4.87e-06 , -3.71e-07 , -1.75e-06 ,
+ -1.80e-06 , 9.99e-06 , -6.46e-06 , -4.91e-06 , 1.33e-05 , -2.52e-07 ,
+ -3.85e-06 , 4.94e-06 , -2.48e-07 , -1.20e-05 , 2.07e-06 , 6.12e-06 ,
+ -1.18e-06 , 4.54e-06 , -1.54e-05 , -1.25e-05 , 1.46e-06 , -6.67e-06 };
+ Double_t par2_137161[64] = { 1.29e-08 , 1.51e-08 , 1.43e-08 , 1.11e-08 ,
+ 5.04e-09 , 6.99e-09 , 5.58e-09 , 4.15e-09 , 4.00e-09 , 8.22e-09 ,
+ 4.97e-09 , 7.66e-09 , 4.91e-09 , 1.10e-09 , 2.64e-09 , 3.64e-09 ,
+ 5.76e-09 , 5.46e-09 , 3.38e-09 , 3.47e-09 , 2.43e-09 , 4.13e-09 ,
+ 2.80e-09 , 5.80e-09 , 3.86e-09 , 7.46e-09 , 5.98e-09 , 2.58e-08 ,
+ 5.50e-09 , 8.72e-09 , 5.23e-09 , 4.37e-09 , 2.33e-09 , -6.01e-10 ,
+ 3.99e-11 , -2.02e-10 , 7.67e-10 , 2.03e-09 , 1.17e-10 , 2.56e-09 ,
+ 1.16e-09 , -4.75e-10 , 1.28e-09 , 1.23e-09 , 1.62e-09 , 1.61e-09 ,
+ 1.93e-09 , 2.97e-10 , 2.21e-09 , 2.16e-09 , 5.22e-10 , 1.03e-09 ,
+ 1.56e-09 , 5.00e-10 , 1.01e-09 , 2.93e-09 , 1.05e-09 , 9.96e-11 ,
+ 1.21e-09 , 7.45e-10 , 3.07e-09 , 2.31e-09 , 6.70e-10 , 1.89e-09 };
+
+ Double_t par0_137366[64] = { 7.12e-02 , 7.34e-02 , 7.39e-02 , 6.54e-02 , 6.11e-02 , 6.31e-02 , 6.15e-02 ,
+ 6.00e-02 , 6.10e-02 , 6.49e-02 , 6.17e-02 , 6.33e-02 , 6.00e-02 , 5.48e-02 ,
+ 5.44e-02 , 5.81e-02 , 6.49e-02 , 7.07e-02 , 5.91e-02 , 6.18e-02 , 4.82e-02 ,
+ 5.67e-02 , 5.36e-02 , 6.60e-02 , 6.37e-02 , 6.78e-02 , 6.31e-02 , 1.04e-01 ,
+ 6.91e-02 , 7.32e-02 , 6.61e-02 , 6.16e-02 , 2.64e-02 , 2.81e-02 , 2.64e-02 ,
+ 2.85e-02 , 2.87e-02 , 2.18e-02 , 2.19e-02 , 2.43e-02 , 2.81e-02 , 4.37e-02 ,
+ 3.90e-02 , 4.66e-02 , 4.24e-02 , 4.09e-02 , 4.21e-02 , 3.88e-02 , 4.83e-02 ,
+ 5.23e-02 , 5.44e-02 , 4.85e-02 , 4.42e-02 , 4.58e-02 , 4.74e-02 , 3.14e-02 ,
+ 6.31e-02 , 5.30e-02 , 5.01e-02 , 5.33e-02 , 5.70e-02 , 3.95e-02 , 4.98e-02 , 5.31e-02 };
+ Double_t par1_137366[64] = { -6.99e-05 , -6.99e-05 , -6.94e-05 , -6.55e-05 , -3.55e-05 , -4.50e-05 ,
+ -3.10e-05 , -2.81e-05 , -2.29e-05 , -3.89e-05 , -2.53e-05 , -4.25e-05 ,
+ -1.87e-05 , -2.01e-05 , -1.53e-05 , -2.14e-05 , -2.86e-05 , -4.70e-05 ,
+ -2.23e-05 , -3.30e-05 ,-9.74e-06 , -2.62e-05 , -1.76e-05 , -2.38e-05 ,
+ -2.40e-05 , -3.43e-05 , -2.75e-05 , -6.86e-05 ,-2.35e-05 , -4.45e-05 ,
+ -2.51e-05 , -2.20e-05 , -1.25e-16 , -2.04e-17 , -2.06e-17 , -3.74e-19 ,
+ -1.18e-18 , -2.02e-15 , -3.78e-06 , -1.26e-06 , -2.71e-06 , -6.23e-17 ,
+ -7.39e-08 , -1.76e-16 , -8.98e-06 , -4.10e-18 , -1.34e-05 , -1.06e-16 ,
+ -3.34e-06 , -1.04e-05 , -5.28e-06 , -7.34e-06 , -1.05e-05 , -7.68e-06 ,
+ -1.78e-05 , -1.19e-05 , -1.78e-05 , -1.34e-06 , -9.23e-06 , -3.34e-06 ,
+ -8.02e-06 , -1.39e-05 , -1.38e-05 , -1.40e-05 };
+ Double_t par2_137366[64] = { 1.41e-08 , 1.47e-08 , 1.48e-08 , 1.24e-08 , 6.82e-09 , 8.73e-09 , 6.26e-09 ,
+ 5.53e-09 , 5.40e-09 , 7.93e-09 , 5.49e-09 , 8.77e-09 , 4.21e-09 , 3.93e-09 ,
+ 3.60e-09 , 4.67e-09 , 5.59e-09 , 8.81e-09 , 3.89e-09 , 6.19e-09 , 1.97e-09 ,
+ 4.38e-09 , 3.26e-09 , 5.00e-09 , 4.58e-09 , 6.39e-09 , 5.03e-09 , 1.30e-08 ,
+ 4.95e-09 , 8.26e-09 , 4.57e-09 , 4.10e-09 , 2.35e-09 , 2.30e-09 , 2.15e-09 ,
+ 2.27e-09 , 2.17e-09 , 2.27e-09 , 2.97e-09 , 2.25e-09 , 1.69e-09 , 1.44e-09 ,
+ 1.66e-09 , 1.75e-09 , 2.88e-09 , 1.82e-09 , 3.64e-09 , 1.80e-09 , 1.71e-09 ,
+ 2.66e-09 , 3.01e-09 , 1.95e-09 , 2.64e-09 , 2.42e-09 , 3.68e-09 , 2.66e-09 ,
+ 3.92e-09 , 1.18e-09 , 2.26e-09 , 1.57e-09 , 2.02e-09 , 2.71e-09 , 2.99e-09 , 3.04e-09 };
+
+
+ if (esd->GetRunNumber() <= 137165) {
+ par0=par0_137161;
+ par1=par1_137161;
+ par2=par2_137161;
+ } else {
+ par0=par0_137366;
+ par1=par1_137366;
+ par2=par2_137366;
+ }
//
Float_t multCorr = 0;
Float_t multCorr2 = 0;
for(Int_t i = 0; i < 64; ++i) {
Double_t b = (esdV0->GetMultiplicity(i)*par1[i]-par0[i]);
Double_t s = (b*b-4.*par2[i]*esdV0->GetMultiplicity(i)*esdV0->GetMultiplicity(i));
- Double_t n = 0;
- if (s<0) {
- printf("FPE %d %.2f %.2f %.2e\n",i,esdV0->GetMultiplicity(i),b,(b*b-4.*par2[i]*esdV0->GetMultiplicity(i)*esdV0->GetMultiplicity(i)));
- n = -b;
- }
- else {
- n = (-b + TMath::Sqrt(s));
- }
+ Double_t n = (s<0) ? -b : (-b + TMath::Sqrt(s));
multChCorr[i] = 2.*esdV0->GetMultiplicity(i)/n*par0[i];
multCorr += multChCorr[i];
multCorr2 += (multChCorr[i]/par0[i]/64.);
ClassImp(AliAnalysisMultPbTrackHistoManager)
-const char * AliAnalysisMultPbTrackHistoManager::kStatStepNames[] = { "All Events", "After centrality selection", "After physics Selection", "With Vertex", "After ZDC cut" };
+const char * AliAnalysisMultPbTrackHistoManager::kStatStepNames[] = { "All Events", "After centrality selection", "After physics Selection", "With Vertex (quality cuts)", "After ZDC cut", "After vertex range cut (10 cm)" };
const char * AliAnalysisMultPbTrackHistoManager::kHistoPtEtaVzNames[] = { "hGenPtEtaVz", "hRecPtEtaVz", "hRecPtEtaVzPrim",
"hRecPtEtaVzSecWeak", "hRecPtEtaVzSecMaterial", "hRecPtEtaVzFake"};
const char * AliAnalysisMultPbTrackHistoManager::kHistoDCANames[] = { "hGenDCA", "hRecDCA", "hRecDCAPrim", "hRecDCASecWeak","hRecDCASecMaterial", "hRecDCAFake"};
-const char * AliAnalysisMultPbTrackHistoManager::kHistoPrefix[] = { "hGen", "hRec", "hRecPrim", "hRecSecWeak","hRecSecMaterial", "hRecFake"};
+const char * AliAnalysisMultPbTrackHistoManager::kHistoPrefix[] = { "hGen", "hRec", "hRecPrim", "hRecSecWeak","hRecSecMaterial", "hRecFake", "hRecHighestMeanPt", "hRecLowestMeanPt"};
const char * AliAnalysisMultPbTrackHistoManager::kSpeciesName[] = { "pi+", "K+", "p", "l+", "pi-", "K-", "barp", "l-", "Others"};
}
+
+TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoMeanPt(Histo_t id) {
+ // mean pt computed event by event
+ TString name = TString(kHistoPrefix[id])+"_MeanPt";
+
+ TH1D * h = (TH1D*) GetHisto(name);
+ if (!h) {
+ name+=fHNameSuffix;
+ Bool_t oldStatus = TH1::AddDirectoryStatus();
+ TH1::AddDirectory(kFALSE);
+
+ AliInfo(Form("Booking histo %s",name.Data()));
+
+ h = new TH1D (name.Data(), Form("Pt Event by Event(%s)",kHistoPrefix[id]), 100, 0, 1);
+ h->Sumw2();
+ TH1::AddDirectory(oldStatus);
+ fList->Add(h);
+
+
+ }
+ return h;
+
+}
+
+TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoPtEvent(Histo_t id) {
+
+ // Returns of dNdpt, used to compute mean pt event by event
+
+ TString name = TString(kHistoPrefix[id])+"_PtEvent";
+
+ TH1D * h = (TH1D*) GetHisto(name);
+ if (!h) {
+ name+=fHNameSuffix;
+ Bool_t oldStatus = TH1::AddDirectoryStatus();
+ TH1::AddDirectory(kFALSE);
+
+ AliInfo(Form("Booking histo %s",name.Data()));
+ const Int_t nptbins = 68;
+ const Double_t binsPt[] = {0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.2,2.4,2.6,2.8,3.0,3.2,3.4,3.6,3.8,4.0,4.5,5.0,5.5,6.0,6.5,7.0,8.0,9.0,10.0,11.0,12.0,13.0,14.0,15.0,16.0,18.0,20.0,22.0,24.0,26.0,28.0,30.0,32.0,34.0,36.0,40.0,45.0,50.0};
+
+ h = new TH1D (name.Data(), Form("Pt Event by Event(%s)",kHistoPrefix[id]), nptbins,binsPt);
+ h->SetBit(kKeepMaxMean,1);
+ if(id == kHistoRecLowestMeanPt ) h->SetBit(kKeepMinMean,1);
+ h->Sumw2();
+ TH1::AddDirectory(oldStatus);
+ fList->Add(h);
+
+
+ }
+ return h;
+
+
+}
+
+
TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoVzEvent(Histo_t id) {
// Returns histogram with Vz of the event
+Long64_t AliAnalysisMultPbTrackHistoManager::Merge(TCollection* list)
+{
+ // Merge a list of AliHistoListWrapper objects with this.
+ // Returns the number of merged objects (including this).
+
+ // We have to make sure that all the list contain the same histos in
+ // the same order. We thus also have to sort the list (sorting is
+ // done by name in TList).
+
+ // AliInfo("Merging");
+
+ if (!list)
+ return 0;
+
+ if (list->IsEmpty())
+ return 1;
+
+ TIterator* iter = list->MakeIterator();
+ TObject* obj;
+
+ // collections of all histograms
+ TList collections;
+
+ Int_t count = 0;
+
+ while ((obj = iter->Next())) {
+ Bool_t foundDiffinThisIterStep = kFALSE;
+
+ // Printf("%d - %s",count, obj->GetName());
+ AliHistoListWrapper* entry = dynamic_cast<AliHistoListWrapper*> (obj);
+ if (entry == 0)
+ continue;
+
+ TList * hlist = entry->GetList();
+
+ // Check if all histos in this fList are also in the one from entry and viceversa
+ // Use getters to automatically book non defined histos
+
+ Bool_t areListsDifferent=kTRUE;
+ Int_t iloop = 0;
+ Int_t max_loops = hlist->GetSize() + fList->GetSize(); // In the worst case all of the histos will be different...
+ while(areListsDifferent) {
+ if(iloop>max_loops) AliFatal("Infinite Loop?");
+ iloop++;
+ // sort
+ hlist->Sort();
+ fList->Sort();
+ // loop over the largest
+ TObject * hist =0;
+ TIterator * iterlist = 0;
+ TList * thislist = 0; // the list over which I'm iterating
+ TList * otherlist = 0; // the other
+
+ if (hlist->GetSize() >= fList->GetSize()) {
+ thislist = hlist;
+ otherlist = fList;
+ }
+ else{
+ thislist = fList;
+ otherlist = hlist;
+ }
+ iterlist = thislist->MakeIterator();
+
+ while ((hist= iterlist->Next())){
+ if(!otherlist->FindObject(hist->GetName())){
+ AliInfo(Form("Adding object %s",hist->GetName()));
+ TH1 * hclone = (TH1*) hist->Clone();
+ if (!hclone->InheritsFrom("TH1")) AliFatal(Form("Found a %s. This class only supports objects inheriting from TH1",hclone->ClassName()));
+ hclone->Reset();
+ otherlist->Add(hclone);
+ foundDiffinThisIterStep=kTRUE;
+ }
+ }
+
+ // re-sort before checking
+ hlist->Sort();
+ fList->Sort();
+
+ // check if everything is fine
+ areListsDifferent=kFALSE;
+ if (hlist->GetSize() == fList->GetSize()) {
+ Int_t nhist = fList->GetSize();
+ for(Int_t ihist = 0; ihist < nhist; ihist++){
+ if(strcmp(fList->At(ihist)->GetName(),hlist->At(ihist)->GetName())) areListsDifferent = kTRUE;
+ }
+ } else {
+ areListsDifferent=kTRUE;
+ }
+ }
+
+ // last check: if something is not ok die loudly
+ if (hlist->GetSize() != fList->GetSize()) {
+ AliFatal("Mismatching size!");
+ }
+ Int_t nhist = fList->GetSize();
+ for(Int_t ihist = 0; ihist < nhist; ihist++){
+ if(strcmp(fList->At(ihist)->GetName(),hlist->At(ihist)->GetName())){
+ AliFatal(Form("Mismatching histos: %s -> %s", fList->At(ihist)->GetName(),hlist->At(ihist)->GetName()));
+ } else {
+ // Specific merging strategies, according to the bits...
+ if (fList->At(ihist)->TestBits(kKeepMinMean|kKeepMaxMean)){
+ TH1 * h1 = (TH1*) fList->At(ihist);
+ TH1 * h2 = (TH1*) hlist->At(ihist);
+ if (h1->GetEntries()>0 && h2->GetEntries()>0) {
+ if (h1->TestBit(kKeepMinMean)) {
+ AliInfo(Form("Keeping only the histo with the lowest mean [%s][%f][%f]",h1->GetName(), h1->GetMean(), h2->GetMean()));
+ if(h1->GetMean() > h2->GetMean()) h1->Reset();
+ else h2->Reset();
+ }
+ if (h1->TestBit(kKeepMaxMean)) {
+ AliInfo(Form("Keeping only the histo with the highest mean [%s][%f][%f]",h1->GetName(), h1->GetMean(), h2->GetMean()));
+ if(h2->GetMean() > h1->GetMean()) h1->Reset();
+ else h2->Reset();
+ }
+ }
+ }
+ }
+ }
+
+ if (foundDiffinThisIterStep){
+ iter->Reset(); // We found a difference: previous lists could
+ // also be affected... We start from scratch
+ collections.Clear();
+ count = 0;
+ }
+ else {
+
+ collections.Add(hlist);
+
+ count++;
+ }
+ }
+
+ fList->Merge(&collections);
+
+ delete iter;
+
+ return count+1;
+}
public:
- typedef enum {kHistoGen, kHistoRec, kHistoRecPrim, kHistoRecSecWeak, kHistoRecSecMat, kHistoRecFake, kNHistos} Histo_t;
- typedef enum {kStatAll, kStatCentr, kStatPhysSel, kStatVtx, kStatZDCCut, kNStatBins} Stat_t;
+ typedef enum {kHistoGen, kHistoRec, kHistoRecPrim, kHistoRecSecWeak, kHistoRecSecMat, kHistoRecFake, kHistoRecHighestMeanPt, kHistoRecLowestMeanPt, kNHistos} Histo_t;
+ typedef enum {kStatAll, kStatCentr, kStatPhysSel, kStatVtx, kStatZDCCut, kStatVtxRangeCut, kNStatBins} Stat_t;
typedef enum {kPartPiPlus, kPartKPlus, kPartP, kPartLPlus, kPartPiMinus, kPartKMinus, kPartPBar, kPartLMinus, kPartOther, kNPart} Part_t;
+ // these bits define the behaviour at merging
+ // if kKeepMax(Mib)Mean is set, the histogram with the highest (lowest) mean is kept
+ enum {kKeepMaxMean = BIT(22), kKeepMinMean = BIT(23)};
+
+
AliAnalysisMultPbTrackHistoManager();
AliAnalysisMultPbTrackHistoManager(const char * name,const char * title);
TH1D * GetHistoSpecies(Histo_t id);
TH1D * GetHistoProcess(Histo_t id);
TH1D * GetHistoVzEvent(Histo_t id);
-
+ TH1D * GetHistoPtEvent(Histo_t id);
+ TH1D * GetHistoMeanPt (Histo_t id);
// Misch utils
//
TH1 * GetHisto(const char * name);
+ virtual void Print(Option_t* option = "") const {fList->Print();}
+
+
+ Long64_t Merge(TCollection* list); // need to override this because of the mean pts
+
private:
static const char * kStatStepNames[]; // names of the step hist
#include "AliGenEventHeader.h"
#include "AliESDCentrality.h"
#include "AliMultiplicity.h"
+#include "TFile.h"
#include <iostream>
#include "AliAnalysisMultPbCentralitySelector.h"
#include "AliTriggerAnalysis.h"
PostData(2,fTrackCuts);
PostData(3,fCentrSelector);
- // Cache histogram pointers
+ // Cache (some) histogram pointers to avoid loop in the histo manager
static TH3D * hTracks [AliAnalysisMultPbTrackHistoManager::kNHistos];
static TH2D * hDCA [AliAnalysisMultPbTrackHistoManager::kNHistos];
static TH1D * hNTracks [AliAnalysisMultPbTrackHistoManager::kNHistos];
hNTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen] = fHistoManager->GetHistoMult(AliAnalysisMultPbTrackHistoManager::kHistoGen );
hNTracks[AliAnalysisMultPbTrackHistoManager::kHistoRec] = fHistoManager->GetHistoMult(AliAnalysisMultPbTrackHistoManager::kHistoRec );
+
+ fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec)->Reset();
+
fESD = dynamic_cast<AliESDEvent*>(fInputEvent);
if (strcmp(fESD->ClassName(),"AliESDEvent")) {
AliFatal("Not processing ESDs");
// Fill generated histo
hTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen]->Fill(mcPart->Pt(),mcPart->Eta(),zvGen);
Int_t partCode = fHistoManager->GetLocalParticleID(mcPart);
- cout << "Part " << partCode << endl;
+ // cout << "Part " << partCode << endl;
if (partCode == AliAnalysisMultPbTrackHistoManager::kPartPiPlus ||
partCode == AliAnalysisMultPbTrackHistoManager::kPartPiMinus ||
partCode == AliAnalysisMultPbTrackHistoManager::kPartKPlus ||
partCode == AliAnalysisMultPbTrackHistoManager::kPartP ||
partCode == AliAnalysisMultPbTrackHistoManager::kPartPBar
){
- cout << "Found " << partCode << endl;
+ // cout << "Found " << partCode << endl;
fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoGen, partCode);
}
if (!(zdcA && zdcC) && (!fIsMC)) return;
fHistoManager->GetHistoStats()->Fill(AliAnalysisMultPbTrackHistoManager::kStatZDCCut);
+
+ if(TMath::Abs(vtxESD->GetZ()) > 10) return;
+ fHistoManager->GetHistoStats()->Fill(AliAnalysisMultPbTrackHistoManager::kStatVtxRangeCut);
+
+ // FIXME
+ Float_t ntracksOk = fTrackCuts->CountAcceptedTracks(fESD);
+ const AliMultiplicity * mult = fESD->GetMultiplicity();
+ if (ntracksOk < 1000 && mult->GetNumberOfITSClusters(1) > 4500) {
+ Float_t dummy;
+ Printf("IEV: %d, Orbit: %x, Period: %d, BC: %d\n",fESD->GetEventNumberInFile(), fESD->GetOrbitNumber(),fESD->GetPeriodNumber(),fESD->GetBunchCrossNumber());
+ printf("%s ----> Processing event # %lld\n", CurrentFileName(), Entry());
+ cout << "File " << AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName() << endl;
+ cout << "Nt " << ntracksOk << ", "<< fESD->GetNumberOfTracks() << " V0 " << fCentrSelector->GetCorrV0(fESD,dummy) << " SPD1 " << mult->GetNumberOfITSClusters(1) << endl;
+ }
+
+
// Fill Vertex
}
}
// for each track
- if(accepted)
+ if(accepted) {
hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRec]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
+ if (TMath::Abs(esdTrack->Eta())<0.5 && TMath::Abs(vtxESD->GetZ())<7)
+ fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec)->Fill(esdTrack->Pt());
+ }
if(acceptedNoDCA)
hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRec]->Fill(weightedDCA,esdTrack->Pt());
hNTracks[AliAnalysisMultPbTrackHistoManager::kHistoRec] ->Fill(acceptedTracks);
Float_t v0;
fHistoManager->GetHistoV0vsNtracks(AliAnalysisMultPbTrackHistoManager::kHistoRec)->Fill(acceptedTracks, fCentrSelector->GetCorrV0(fESD,v0));
- // FIXME
- // hNTracks[AliAnalysisMultPbTrackHistoManager::kHistoRec] ->Fill(fESD->GetMultiplicity()->GetNumberOfTracklets());
-
+ // Fill <pt> analysis histos
+ fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec)->Scale(1,"width");// correct bin width
+ if (fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec)->GetEntries()>0)
+ fHistoManager->GetHistoMeanPt(AliAnalysisMultPbTrackHistoManager::kHistoRec)->Fill(fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec)->GetMean());
+ if (fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec)->GetMean() >
+ fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRecHighestMeanPt)->GetMean()) {
+ // Found a new highest pt: first reset than add it to the container for the highest
+ fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRecHighestMeanPt)->Reset();
+ fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRecHighestMeanPt)
+ ->Add(fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec));
+ }
+ if ((fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec)->GetMean() <
+ fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRecLowestMeanPt)->GetMean() &&
+ fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec)->GetEntries()>0) ||
+ !(fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRecLowestMeanPt)->GetEntries()>0)) {
+ // Found a new lowest pt: first reset than add it to the container for the lowest
+ fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRecLowestMeanPt)->Reset();
+ fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRecLowestMeanPt)
+ ->Add(fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec));
+ }
}
void AliAnalysisTaskMultPbTracks::Terminate(Option_t *){
#include "AliESDtrackCuts.h"
#include "AliESDVZERO.h"
#include "TH2F.h"
+#include "AliESDUtils.h"
using namespace std;
// Called once
fHistoList = new AliHistoListWrapper("histoList","histogram list for trigger studies");
fTriggerAnalysis = new AliTriggerAnalysis();
+ if (fIsMC) fTriggerAnalysis->SetAnalyzeMC(1);
}
AliESDVZERO* esdV0 = fESD->GetVZEROData();
Float_t multV0A=esdV0->GetMTotV0A();
Float_t multV0C=esdV0->GetMTotV0C();
- Float_t multV0 = multV0A+multV0C;
+ Float_t dummy = 0;
+ Float_t multV0 = AliESDUtils::GetCorrV0(fESD, dummy);
// Get number of clusters in layer 1
Float_t outerLayerSPD = mult->GetNumberOfITSClusters(1);
Float_t innerLayerSPD = mult->GetNumberOfITSClusters(0);
Float_t totalClusSPD = outerLayerSPD+innerLayerSPD;
+ if ( !fIsMC &&(!(fESD->IsTriggerClassFired("CMBS2A-B-NOPF-ALL")|| fESD->IsTriggerClassFired("CMBS2C-B-NOPF-ALL") || fESD->IsTriggerClassFired("CMBAC-B-NOPF-ALL")) )) return;
GetHistoSPD1 ("All", "All events before any selection")->Fill(outerLayerSPD);
+ GetHistoV0M ("All", "All events before any selection")->Fill(multV0);
// Physics selection
Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
if(!isSelected) return;
- // FIXME trigger classes
- // if ( !(fESD->IsTriggerClassFired("CMBS2A-B-NOPF-ALL")|| fESD->IsTriggerClassFired("CMBS2C-B-NOPF-ALL") || fESD->IsTriggerClassFired("CMBAC-B-NOPF-ALL")) ) return;
- GetHistoSPD1 ("AllPSNoPrino", "All events after physsel, before Francesco")->Fill(outerLayerSPD);
+
+ GetHistoSPD1 ("AllPSNoPrino", "All events after physsel (no ZDC time)")->Fill(outerLayerSPD);
+ GetHistoV0M ("AllPSNoPrino", "All events after physsel (no ZDC time)")->Fill(multV0);
// Francesco's cuts
- const AliESDVertex * vtxESDTPC= fESD->GetPrimaryVertexTPC();
- if(vtxESDTPC->GetNContributors()<1) return;
- if (vtxESDTPC->GetNContributors()<(-10.+0.25*fESD->GetMultiplicity()->GetNumberOfITSClusters(0))) return;
- const AliESDVertex * vtxESDSPD= fESD->GetPrimaryVertexSPD();
- Float_t tpcContr=vtxESDTPC->GetNContributors();
+ // const AliESDVertex * vtxESDTPC= fESD->GetPrimaryVertexTPC();
+ // if(vtxESDTPC->GetNContributors()<1) return;
+ // if (vtxESDTPC->GetNContributors()<(-10.+0.25*fESD->GetMultiplicity()->GetNumberOfITSClusters(0))) return;
+ // const AliESDVertex * vtxESDSPD= fESD->GetPrimaryVertexSPD();
+ // Float_t tpcContr=vtxESDTPC->GetNContributors();
Bool_t c0OM3 = h->IsTriggerInputFired("0OM3"); // thr >= 3 (input 20)
// ZDC triggers
- Bool_t zdcA = kFALSE;
- Bool_t zdcC = kFALSE;
- Bool_t zdcBar = kFALSE;
+ Bool_t zdcA = kFALSE;
+ Bool_t zdcC = kFALSE;
+ Bool_t zdcBar = kFALSE;
+ Bool_t zdcTime = fTriggerAnalysis->ZDCTimeTrigger(fESD);
if (!fIsMC) {
// If it's data, we use the TDCs
vdArray[kVDC0VBC] = c0v0C;
//vdArray[kVDC0OM2] = c0OM2;
- // Plots requested by Jurgen on 18/11/2010
+ // Plots requested by Jurgen on 18/11/2010 + later updates (including plots for the note)
// FIXME: will skip everything else
- GetHistoSPD1 ("PhysSel", "All events after physics selection and Francesco's cut")->Fill(outerLayerSPD);
- GetHistoTracks("PhysSel", "All events after physics selection and Francesco's cut")->Fill(ntracks);
- GetHistoV0M ("PhysSel", "All events after physics selection and Francesco's cut")->Fill(multV0);
+
+ if(zdcTime) {
+ GetHistoSPD1 ("ZDCTIME", "ZDC Time Cut")->Fill(outerLayerSPD);
+ GetHistoTracks("ZDCTIME", "ZDC Time Cut")->Fill(ntracks);
+ GetHistoV0M ("ZDCTIME", "ZDC Time Cut")->Fill(multV0);
+ }
+
+ if(zdcTime && ntracks > 1) {
+ GetHistoSPD1 ("ZDCTIME1TRACK", "ZDC Time Cut & 1 Track")->Fill(outerLayerSPD);
+ GetHistoTracks("ZDCTIME1TRACK", "ZDC Time Cut & 1 Track")->Fill(ntracks);
+ GetHistoV0M ("ZDCTIME1TRACK", "ZDC Time Cut & 1 Track")->Fill(multV0);
+ }
+
+ // GetHistoSPD1 ("PhysSel", "All events after physics selection and Francesco's cut")->Fill(outerLayerSPD);
+ // GetHistoTracks("PhysSel", "All events after physics selection and Francesco's cut")->Fill(ntracks);
+ // GetHistoV0M ("PhysSel", "All events after physics selection and Francesco's cut")->Fill(multV0);
if(c0v0A && c0v0C) {
GetHistoSPD1 ("V0AND", "V0A & V0C")->Fill(outerLayerSPD);
GetHistoTracks("V0AND", "V0A & V0C")->Fill(ntracks);
Bool_t oldStatus = TH1::AddDirectoryStatus();
TH1::AddDirectory(kFALSE);
// h = new TH1F (hname.Data(), title, 1000, -0.5, 10000-0.5);
- h = new TH1F (hname.Data(), title, 300, -0.5, 300-0.5);
+ h = new TH1F (hname.Data(), title, 500, -0.5, 500-0.5);
h->Sumw2();
h->SetXTitle("V0M");
fHistoList->GetList()->Add(h);
#include "AliESDtrackCuts.h"
#include "AliAnalysisMultPbCentralitySelector.h"
#include "TLegend.h"
+#include "AliBWTools.h"
using namespace std;
Float_t dNdetaE = TMath::Sqrt(f->IntegralError(0,0.1)*f->IntegralError(0,0.1) + errorData*errorData) / (etaMax-etaMin);
cout << "dN/deta " << dNdeta << " +- " << dNdetaE << endl;
cout << "(dN/deta)/0.5 Npart " << dNdeta/npart*2 << " +- " << dNdetaE/npart*2 << endl;
+ Double_t mean, meanerr;
+ AliBWTools::GetMeanDataAndExtrapolation(hCorrected, f, mean, meanerr,0,1000);
+ cout << "Mean Pt " << mean << "+-" << meanerr << endl;
+
TH1D * hDataGen = hManData->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoGen, etaMin,etaMax,zmin,zmax);
cout << "Generated dN/deta (data) = " << hDataGen->Integral("width") << endl;
hDataGen->Draw("same");
l->AddEntry(hMCPtGen, "Monte Carlo (generated)");
l->AddEntry(f, "Hagedorn Fit");
l->Draw();
+
+
+ Double_t errXCheck = 0, yieldXCheck=0;
+ AliBWTools::GetYield(hCorrected,f,yieldXCheck,errXCheck);
+ cout << "BWTools crosscheck: " << yieldXCheck << "+-" << errXCheck << endl;
+
+
}
void CheckSecondaries(Double_t &fracWeak, Double_t &fracMaterial) {
// Set the components which are used in HistoSum, the static
// function for GetFunctionHistoSum
// Project onti DCA axis
- // const Int_t ptBinsFit[] = {3,5,7,9,11,15,19,23,31,-1};
+ //const Int_t ptBinsFit[] = {3,5,7,9,11,15,19,23,31,-1};
const Int_t ptBinsFit[] = {3,20,-1};
Int_t ibinPt = -1;
while(ptBinsFit[(++ibinPt)+1]!=-1){
gROOT->ProcessLine(gSystem->ExpandPathName(".include $ALICE_ROOT/PWG0/multPbPb"));
gROOT->ProcessLine(gSystem->ExpandPathName(".include $ALICE_ROOT/PWG1/background"));
+ gROOT->ProcessLine(gSystem->ExpandPathName(".include $ALICE_ROOT/PWG2/SPECTRA/Fit"));
+
// Load helper classes
// TODO: replace this by a list of TOBJStrings
TString taskName("$ALICE_ROOT/PWG0/multPbPb/AliAnalysisTaskMultPbTracks.cxx+");
TF1 * GetFunctionHistoSum() {
TF1 * f = new TF1 ("fHistoSum",HistoSum,-1,1,kHistoFitCompoments);
- f->SetParNames("Primaries+Fakes", "Sec. Weak decays","Sec. Material");
+ // f->SetParNames("Primaries+Fakes", "Sec. Weak decays","Sec. Material");
+ f->SetParNames("Primaries+Sec. Material", "Sec. Weak decays");
f->SetNpx(1000);
return f;
gROOT->ProcessLine(".L $ALICE_ROOT/PWG0/multPbPb/AddTaskMultPbPbTracksAllCentrality.C");
centrSelector->SetUseV0Range(kTRUE);
Int_t ncentr = 11;
- const Float_t minCentr[] = {0, 79, 247,577, 1185,2155,3565,5527,8203, 12167,15073};
- const Float_t maxCentr[] = {79,247,577,1185,2155,3565,5527,8203,12167,15073,21000};
+
+ const Float_t minCentr[] = {0 ,79 ,239,559 ,1165,2135,3555,5525,8213 ,12191,15079};
+ const Float_t maxCentr[] = {79,239,559,1165,2135,3555,5525,8213,12191,15079,21000};
AliAnalysisTaskMultPbTracks ** tasks = AddTaskMultPbPbTracksAllCentrality("multPbPbtracks.root", cuts, centrSelector, ncentr,minCentr,maxCentr);
for(Int_t icentr = 0; icentr < ncentr; icentr++){
tasks[icentr]->Print();
runmode=1
dataset=/alice/sim/LHC10f8a_130844
ropt="-l"
-option="SAVE"
+option="DCA,SAVE"
workers=26
analysismode=9; #SPD + field on
centrBin=-1
npart=381.188
weakFactor=-1
useSingleBin=kTRUE
+OUTPATH=output.BAK2010
give_help() {
}
-while getopts "x:sr:c:gmd:o:w:n:e:b:t:k:vy:0:2:hz:a:l" opt; do
+while getopts "x:sr:c:gmd:o:w:n:e:b:t:k:vy:0:2:hz:a:lp:" opt; do
case $opt in
r)
run=yes
;;
c)
correct=yes
- dataDir="./output/${OPTARG%%,*}"
- mcDir="./output/${OPTARG##*,}"
+ dataDir="./$OUTPATH/${OPTARG%%,*}"
+ mcDir="./$OUTPATH/${OPTARG##*,}"
;;
z)
vzMin=${OPTARG%%,*}
// Add physics selection
gROOT->ProcessLine(".L $ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");
physicsSelectionTask = AddTaskPhysicsSelection(isMC,0);//FIXME
-
+ physicsSelectionTask->GetPhysicsSelection()->SetSkipZDCTime(1);// Skip ZDC - applyied later