]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
correct.C
authormfloris <mfloris@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 24 Nov 2010 09:34:47 +0000 (09:34 +0000)
committermfloris <mfloris@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 24 Nov 2010 09:34:47 +0000 (09:34 +0000)
- Fitting with Hagedorn as Default
- DCA fit vs Pt, with only 2 components (prim+material) and (weak)
run.C
- added "_" to suffixes
AliAnalysisMultPbTrackHistoManager
- 2D DCA histos (DCA vs pt)
AliAnalysisTaskMultPbTracks
- 2D DCA histos
- Skipping displaced primaries in buggy hijing (up to 10h6 production)
- Looping over all generated (not just primaries)
AliAnalysisTaskTriggerStudy
- Default set of plots changed.

PWG0/multPbPb/AliAnalysisMultPbTrackHistoManager.cxx
PWG0/multPbPb/AliAnalysisMultPbTrackHistoManager.h
PWG0/multPbPb/AliAnalysisTaskMultPbTracks.cxx
PWG0/multPbPb/AliAnalysisTaskTriggerStudy.cxx
PWG0/multPbPb/AliAnalysisTaskTriggerStudy.h
PWG0/multPbPb/correct.C
PWG0/multPbPb/run.C
PWG0/multPbPb/runTriggerStudy.C

index e0f1e00371abf895f8a6a086fb0b46f99e35cde1..9849abc865ae10993a7e351c8267c3748cc78fed 100644 (file)
@@ -57,10 +57,10 @@ TH3D * AliAnalysisMultPbTrackHistoManager::GetHistoPtEtaVz(Histo_t id, Int_t par
 
 }
 
-TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoDCA(Histo_t id) {
+TH2D * AliAnalysisMultPbTrackHistoManager::GetHistoDCA(Histo_t id) {
   // Returns a 3D histo of Pt/eta/vtx. It it does not exist, books it.
 
-  TH1D * h = (TH1D*) GetHisto(kHistoDCANames[id]);
+  TH2D * h = (TH2D*) GetHisto(kHistoDCANames[id]);
   if (!h) {
     h = BookHistoDCA(kHistoDCANames[id], Form("Pt Eta Vz distribution (%s)",kHistoDCANames[id]));
   }
@@ -97,14 +97,14 @@ TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoPt (Histo_t id,
   // Get range in terms of bin numners.  If the float range is
   // less than -11111 take the range from the first to the last bin (i.e. no
   // under/over-flows)
-  // FIXME: UNDERFLOWS
+
   Int_t min1 = minEta  < -11111 ? -1 : h3D ->GetYaxis()->FindBin(minEta);
   Int_t min2  = minVz  < -11111 ? -1 : h3D ->GetZaxis()->FindBin(minVz) ;
 
-  // Int_t max1 = maxEta  < -11111 ? h3D->GetNbinsY() : h3D ->GetYaxis()->FindBin(maxEta-0.00001);
-  // Int_t max2  = maxVz  < -11111 ? h3D->GetNbinsZ() : h3D ->GetZaxis()->FindBin(maxVz -0.00001);
-  Int_t max1 = maxEta  < -11111 ? -1 : h3D ->GetYaxis()->FindBin(maxEta-0.00001);
-  Int_t max2  = maxVz  < -11111 ? -1 : h3D ->GetZaxis()->FindBin(maxVz -0.00001);
+  Int_t max1 = maxEta  < -11111 ? h3D->GetNbinsY() : h3D ->GetYaxis()->FindBin(maxEta-0.00001);
+  Int_t max2  = maxVz  < -11111 ? h3D->GetNbinsZ() : h3D ->GetZaxis()->FindBin(maxVz -0.00001);
+  // Int_t max1 = maxEta  < -11111 ? -1 : h3D ->GetYaxis()->FindBin(maxEta-0.00001);
+  // Int_t max2  = maxVz  < -11111 ? -1 : h3D ->GetZaxis()->FindBin(maxVz -0.00001);
 
 
   TString hname = h3D->GetName();
@@ -412,9 +412,25 @@ TH3D * AliAnalysisMultPbTrackHistoManager::BookHistoPtEtaVz(const char * name, c
   return h;
 }
 
-TH1D * AliAnalysisMultPbTrackHistoManager::BookHistoDCA(const char * name, const char * title) {
+TH2D * AliAnalysisMultPbTrackHistoManager::BookHistoDCA(const char * name, const char * title) {
   // Books a DCA histo 
 
+  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};
+  const Int_t ndcabins=100;
+  Double_t binsDca[ndcabins+1];
+  Float_t minDca = -3;
+  Float_t maxDca =  3;
+  //  const Double_t binsDCA[] = {-3,-2.8,-2.6,-2.4,-2.2,-2.0,-1.9,};
+  Float_t dcaStep = (maxDca-minDca)/ndcabins;
+  for(Int_t ibin = 0; ibin < ndcabins; ibin++){    
+    binsDca[ibin]   = minDca + ibin*dcaStep;
+    binsDca[ibin+1] = minDca + ibin*dcaStep + dcaStep;
+  }
+
+
+
   Bool_t oldStatus = TH1::AddDirectoryStatus();
   TH1::AddDirectory(kFALSE);
 
@@ -424,9 +440,14 @@ TH1D * AliAnalysisMultPbTrackHistoManager::BookHistoDCA(const char * name, const
   AliInfo(Form("Booking %s",hname.Data()));
   
 
+#if defined WEIGHTED_DCA
   TH1D * h = new TH1D (hname,title, 200,0,200);
-
   h->SetXTitle("#Delta DCA");
+#elif defined TRANSVERSE_DCA
+  TH2D * h = new TH2D (hname,title, ndcabins,binsDca,nptbins,binsPt);
+  h->SetYTitle("p_{T} (GeV)");
+  h->SetXTitle("d_{0} r#phi (cm)");
+#endif
   h->Sumw2();
   
   fList->Add(h);
index b130d960a1a8d12da664f0c170c642aa90ce1319..9aa03ae6d959d96fc8c396ec00ed9ece458d2d56 100644 (file)
@@ -9,6 +9,10 @@ class TH1D;
 class TH1I;
 class AliMCParticle;
 
+//#define WEIGHTED_DCA
+#define TRANSVERSE_DCA
+
+
 //-------------------------------------------------------------------------
 //                      AliAnalysisMultPbTrackHistoManager
 // 
@@ -44,7 +48,7 @@ public:
   TH2D * GetHistoPtVz (Histo_t id, Float_t minEta = -22222, Float_t maxEta = -22222, Bool_t scaleWidth = kFALSE);
 
   TH1I * GetHistoStats();
-  TH1D * GetHistoDCA(Histo_t id);
+  TH2D * GetHistoDCA(Histo_t id);
   TH1D * GetHistoMult(Histo_t id);
 
   TH1D * GetHistoSpecies(Histo_t id);
@@ -60,7 +64,7 @@ public:
 
   // Histo bookers
   TH3D * BookHistoPtEtaVz(const char * name, const char * title);
-  TH1D * BookHistoDCA(const char * name, const char * title);
+  TH2D * BookHistoDCA(const char * name, const char * title);
   TH1I * BookHistoStats();
   TH1D * BookHistoMult(const char * name, const char * title);
 
index 2df37a6efd4dbe32fd5b311c44fba7eb46d93e11..bb61d888d34a3cc7e36b501bac9cb925c6c7e723 100644 (file)
@@ -13,6 +13,7 @@
 #include "AliStack.h"
 #include "TH1I.h"
 #include "TH3D.h"
+#include "TH2D.h"
 #include "AliMCParticle.h"
 #include "AliGenEventHeader.h"
 #include "AliESDCentrality.h"
@@ -108,7 +109,7 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *)
 
   // Cache histogram pointers
   static TH3D * hTracks  [AliAnalysisMultPbTrackHistoManager::kNHistos];
-  static TH1D * hDCA     [AliAnalysisMultPbTrackHistoManager::kNHistos];
+  static TH2D * hDCA     [AliAnalysisMultPbTrackHistoManager::kNHistos];
   static TH1D * hNTracks [AliAnalysisMultPbTrackHistoManager::kNHistos];
   hTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen]        = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoGen       );
   hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRec]        = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec       );
@@ -132,7 +133,6 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *)
     AliFatal("Not processing ESDs");
   }
   
-  // FIXME: use physics selection here to keep track of events lost?
   fHistoManager->GetHistoStats()->Fill(AliAnalysisMultPbTrackHistoManager::kStatAll);
 
 
@@ -171,7 +171,9 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *)
       
       //loop on the MC event, only over primaries, which are always
       //      the first in stack.
-      Int_t nMCTracks = fMCEvent->GetNumberOfPrimaries();
+      // WRONG: D&B decays are produced in the transport... Need To Loop over all particles
+      //      Int_t nMCTracks = fMCEvent->GetNumberOfPrimaries();
+      Int_t nMCTracks = fMCEvent->GetNumberOfTracks();
       Int_t nPhysicalPrimaries = 0;
       for (Int_t ipart=0; ipart<nMCTracks; ipart++) { 
        
@@ -182,7 +184,13 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *)
 
        //check if current particle is a physical primary
        if(!IsPhysicalPrimaryAndTransportBit(ipart)) continue;
+       if(TMath::Abs(mcPart->Zv()-zvGen)>1e-6) {
+         // This cures a bug in Hijing
+         // A little hack here: I put those in the underflow bin of the process type to keep track of them
+         fHistoManager->GetHistoProcess(AliAnalysisMultPbTrackHistoManager::kHistoGen)->Fill(-1);
+         continue;
+       }
+
        nPhysicalPrimaries++;
        // Fill species histo and particle species
        fHistoManager->GetHistoProcess(AliAnalysisMultPbTrackHistoManager::kHistoGen)->Fill(mcPart->Particle()->GetUniqueID());
@@ -192,10 +200,18 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *)
        // Fill generated histo
        hTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen]->Fill(mcPart->Pt(),mcPart->Eta(),zvGen);
        Int_t partCode = fHistoManager->GetLocalParticleID(mcPart);
+       cout << "Part " << partCode << endl;
        if (partCode == AliAnalysisMultPbTrackHistoManager::kPartPiPlus  || 
-           partCode == AliAnalysisMultPbTrackHistoManager::kPartPiMinus)
+           partCode == AliAnalysisMultPbTrackHistoManager::kPartPiMinus ||
+           partCode == AliAnalysisMultPbTrackHistoManager::kPartKPlus   || 
+           partCode == AliAnalysisMultPbTrackHistoManager::kPartKMinus  ||
+           partCode == AliAnalysisMultPbTrackHistoManager::kPartP       || 
+           partCode == AliAnalysisMultPbTrackHistoManager::kPartPBar  
+           ){
+         cout << "Found " << partCode << endl;
+         
          fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoGen, partCode);
-       
+       }
       }
       hNTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen]->Fill(nPhysicalPrimaries);
       fHistoManager->GetHistoVzEvent(AliAnalysisMultPbTrackHistoManager::kHistoGen)->Fill(zvGen);
@@ -205,21 +221,10 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *)
   }
   
 
-  // FIXME: shall I take the primary vertex?
   const AliESDVertex* vtxESD = fESD->GetPrimaryVertex();
   if(!vtxESD) return;
-  // FIXME: leave the cuts here or find a better place?)
-  // FIXME: implement these latest cuts, 
-// const AliESDVertex* vtxESD = esd->GetPrimaryVertexSPD();
-// if (vtxESD->GetNContributors()<1) return;
-// if (vtxESD->GetDispersion()>0.04) return;
-// if (vtxESD->GetZRes()>0.25) return;
-// const AliMultiplicity* multESD = esd->GetMultiplicity();
-// const AliESDVertex* vtxESDTPC = esd->GetPrimaryVertexTPC();
-// if(vtxESDTPC->GetNContributors()<1) return;
-// if(vtxESDTPC->GetNContributors()<(-10.+0.25*multESD->GetNumberOfITSClusters(0))) 
- // Quality cut on vertexer Z, as suggested by Francesco Prino
+  // TODO: leave the cuts here or find a better place?)
+  // Quality cut on vertexer Z, as suggested by Francesco Prino
   if(vtxESD->IsFromVertexerZ()) {
     if (vtxESD->GetNContributors() <= 0) return;
     if (vtxESD->GetDispersion() >= 0.04) return;
@@ -227,7 +232,10 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *)
   }
   // "Beam gas" vertex cut
   const AliESDVertex * vtxESDTPC= fESD->GetPrimaryVertexTPC(); 
+  if(vtxESDTPC->GetNContributors()<1) return;
   if (vtxESDTPC->GetNContributors()<(-10.+0.25*fESD->GetMultiplicity()->GetNumberOfITSClusters(0)))     return;
+
+  // Fill vertex and statistics
   fHistoManager->GetHistoStats()->Fill(AliAnalysisMultPbTrackHistoManager::kStatVtx);
   fHistoManager->GetHistoVzEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec)->Fill(vtxESD->GetZ());
 
@@ -247,13 +255,13 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *)
 
     // Compute weighted offset
     // TODO: other choiches for the DCA?
-    Double_t b = fESD->GetMagneticField();
-    Double_t dca[2];
-    Double_t cov[3];
     Double_t weightedDCA = 10;
     
 
-
+#if defined WEIGHTED_DCA
+    Double_t b = fESD->GetMagneticField();
+    Double_t dca[2];
+    Double_t cov[3];
     if (esdTrack->PropagateToDCA(vtxESD, b,3., dca, cov)) {
       Double_t det = cov[0]*cov[2]-cov[1]*cov[1]; 
       if (det<=0) {
@@ -266,7 +274,11 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *)
       //      printf("dR:%e dZ%e  sigR2:%e sigRZ:%e sigZ2:%e\n",dca[0],dca[1],cov[0],cov[1],cov[2]);
     }
     
-
+#elif defined TRANSVERSE_DCA
+    Float_t xz[2]; 
+    esdTrack->GetDZ(vtxESD->GetX(),vtxESD->GetY(),vtxESD->GetZ(), fESD->GetMagneticField(), xz); 
+    weightedDCA = xz[0];
+#endif
     
 
     // Alternative: p*DCA
@@ -274,11 +286,23 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *)
     // esdTrack->GetDZ(vtxESD->GetX(),vtxESD->GetY(),vtxESD->GetZ(), fESD->GetMagneticField(), xz); 
     // Float_t dca = xz[0]*esdTrack->P();
 
-    // for each track
-    if(accepted) 
-      hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRec]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
-    if(acceptedNoDCA)
-      hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRec]->Fill(weightedDCA);
+  // This cures a bug in Hijing (displaced primaries)
+  if (fIsMC) {
+    Int_t label = TMath::Abs(esdTrack->GetLabel()); // no fakes!!!    
+    if (IsPhysicalPrimaryAndTransportBit(label)) {
+      AliMCParticle *mcPart  =  (AliMCParticle*)fMCEvent->GetTrack(label);
+      if(!mcPart) AliFatal("No particle");
+      TArrayF   vertex;
+      fMCEvent->GenEventHeader()->PrimaryVertex(vertex);
+      Float_t zvGen = vertex[2];    
+      if(TMath::Abs(mcPart->Zv()-zvGen)>1e-6) continue;    
+    }
+  }
+  // for each track
+  if(accepted) 
+    hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRec]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
+  if(acceptedNoDCA)
+    hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRec]->Fill(weightedDCA,esdTrack->Pt());
 
     // Fill separately primaries and secondaries
     // FIXME: fakes? Is this the correct way to account for them?
@@ -291,7 +315,7 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *)
        if(accepted)
          hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecFake]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
        if(acceptedNoDCA)
-         hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecFake]->Fill(weightedDCA);
+         hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecFake]->Fill(weightedDCA,esdTrack->Pt());
       }
       else {
        if(IsPhysicalPrimaryAndTransportBit(label)) {
@@ -301,11 +325,16 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *)
            hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecPrim]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
            Int_t partCode = fHistoManager->GetLocalParticleID(mcPart);
            if (partCode == AliAnalysisMultPbTrackHistoManager::kPartPiPlus  || 
-               partCode == AliAnalysisMultPbTrackHistoManager::kPartPiMinus)
+               partCode == AliAnalysisMultPbTrackHistoManager::kPartPiMinus ||
+               partCode == AliAnalysisMultPbTrackHistoManager::kPartKPlus   || 
+               partCode == AliAnalysisMultPbTrackHistoManager::kPartKMinus  ||
+               partCode == AliAnalysisMultPbTrackHistoManager::kPartP       || 
+               partCode == AliAnalysisMultPbTrackHistoManager::kPartPBar  
+               )
              fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim, partCode);
          }
          if(acceptedNoDCA)
-           hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecPrim]->Fill(weightedDCA);
+           hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecPrim]->Fill(weightedDCA,esdTrack->Pt());
        } 
        else {
          Int_t mfl=0;
@@ -320,13 +349,13 @@ void AliAnalysisTaskMultPbTracks::UserExec(Option_t *)
            if(accepted)
              hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
            if(acceptedNoDCA)
-             hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak]->Fill(weightedDCA);      
+             hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak]->Fill(weightedDCA,esdTrack->Pt());       
          }else{ // material
            fHistoManager->GetHistoProcess(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat)->Fill(mcPart->Particle()->GetUniqueID());
            if(accepted)
              hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
            if(acceptedNoDCA)
-             hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat]->Fill(weightedDCA);       
+             hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat]->Fill(weightedDCA,esdTrack->Pt());        
 
          }
 
@@ -353,7 +382,7 @@ void   AliAnalysisTaskMultPbTracks::Terminate(Option_t *){
 
 
 Bool_t AliAnalysisTaskMultPbTracks::IsPhysicalPrimaryAndTransportBit(Int_t ipart) {
-
+  // Check if it is a primary and if it was transported
   Bool_t physprim=fMCEvent->IsPhysicalPrimary(ipart);
   if (!physprim) return kFALSE;
   Bool_t transported = ((AliMCParticle*)fMCEvent->GetTrack(ipart))->Particle()->TestBit(kTransportBit);
index 54866d7d297d7285f539a1110f0825334d3ac752..f76f7d1720104e04dd0c621b96c7824e024d67d5 100644 (file)
@@ -23,6 +23,8 @@
 #include "TFile.h"
 #include "AliLog.h"
 #include "AliESDtrackCuts.h"
+#include "AliESDVZERO.h"
+#include "TH2F.h"
 
 using namespace std;
 
@@ -93,6 +95,8 @@ void AliAnalysisTaskTriggerStudy::UserExec(Option_t *)
 {
   // User code
 
+  // FIXME: make sure you have the right cuts here
+
   /* PostData(0) is taken care of by AliAnalysisTaskSE */
   PostData(1,fHistoList);
 
@@ -101,10 +105,60 @@ void AliAnalysisTaskTriggerStudy::UserExec(Option_t *)
     AliFatal("Not processing ESDs");
   }
 
-  
   // get the multiplicity object
   const AliMultiplicity* mult = fESD->GetMultiplicity();
   Int_t ntracklets = mult->GetNumberOfTracklets();
+  // Get Number of tracks
+  Int_t ntracks    = AliESDtrackCuts::GetReferenceMultiplicity(fESD,kTRUE); // tpc only
+  
+  // Get V0 Multiplicity
+  AliESDVZERO* esdV0 = fESD->GetVZEROData();
+  Float_t multV0A=esdV0->GetMTotV0A();
+  Float_t multV0C=esdV0->GetMTotV0C();
+  Float_t multV0 = multV0A+multV0C;
+
+  // Get number of clusters in layer 1
+  Float_t outerLayerSPD = mult->GetNumberOfITSClusters(1);  
+  Float_t innerLayerSPD = mult->GetNumberOfITSClusters(0);  
+  Float_t totalClusSPD = outerLayerSPD+innerLayerSPD;
+
+  GetHistoSPD1  ("All", "All events before any selection")->Fill(outerLayerSPD);
+
+
+  // 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);
+
+
+  // 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();
+
+  
+
+  // GetT0 Stuff
+  const Double32_t *meanT0   = fESD->GetT0TOF();
+  const Double32_t  meanT0A  = 0.001* meanT0[1];
+  const Double32_t  meanT0C  = 0.001* meanT0[2];
+  const Double32_t  meanT0AC = 0.001* meanT0[0];
+  Double32_t  T0Vertex = fESD->GetT0zVertex();
+  //  Double32_t  *ampT0 =Esdevent ->GetT0amplitude();
+
+//   cut1yesF = ( (meanC < 95. && meanA < 95.) && (meanC < -2.) ) && francescoscut
+// cut1notF = ( (meanC < 95. && meanA < 95.) && (meanC < -2.) ) && ! francescoscut
+// cut2 = ( (meanC < 95. && meanA < 95.) && ( (meanC-meanA) <=-0.7) && meanC > -2) )
+// cut3 = ( (meanC < 95. && meanA < 95.) && ( (meanC-meanA) < 0.7 && (meanC-meanA) > -0.7 ) )
+//   cut4 = ( (meanC < 95. && meanA < 95.) && (meanA < -2.)
+
+  Bool_t cut1T0 =  ( (meanT0C < 95. && meanT0A < 95.) && (meanT0C < -2.) );
+  Bool_t cut2T0 = ( (meanT0C < 95. && meanT0A < 95.) && ( (meanT0C-meanT0A) <=-0.7) && meanT0C > -2) ;
 
   if(ntracklets > fNTrackletsCut) return;
 
@@ -139,9 +193,27 @@ void AliAnalysisTaskTriggerStudy::UserExec(Option_t *)
   Bool_t c0OM3 = h->IsTriggerInputFired("0OM3"); // thr >= 3 (input 20)
 
   // ZDC triggers
-  Bool_t zdcA   = fTriggerAnalysis->ZDCTrigger(fESD, AliTriggerAnalysis::kASide) ;
-  Bool_t zdcC   = fTriggerAnalysis->ZDCTrigger(fESD, AliTriggerAnalysis::kCSide) ;
-  Bool_t zdcBar = fTriggerAnalysis->ZDCTrigger(fESD, AliTriggerAnalysis::kCentralBarrel) ;
+  Bool_t zdcA   = kFALSE;
+  Bool_t zdcC   = kFALSE;
+  Bool_t zdcBar = kFALSE;
+
+  if (!fIsMC) {
+    // If it's data, we use the TDCs
+    zdcA   = fTriggerAnalysis->ZDCTDCTrigger(fESD, AliTriggerAnalysis::kASide, kTRUE, kFALSE) ; 
+    zdcC   = fTriggerAnalysis->ZDCTDCTrigger(fESD, AliTriggerAnalysis::kCSide, kTRUE, kFALSE) ;                        
+    zdcBar = fTriggerAnalysis->ZDCTDCTrigger(fESD, AliTriggerAnalysis::kCentralBarrel) ;                               
+  } else {
+    // If it's MC, we use the energy
+    Double_t minEnergy = 0;
+    AliESDZDC *esdZDC = fESD->GetESDZDC();    
+    Double_t  zNCEnergy = esdZDC->GetZDCN1Energy();
+    Double_t  zPCEnergy = esdZDC->GetZDCP1Energy();
+    Double_t  zNAEnergy = esdZDC->GetZDCN2Energy();
+    Double_t  zPAEnergy = esdZDC->GetZDCP2Energy();
+    zdcA = (zNAEnergy>minEnergy || zPAEnergy>minEnergy);
+    zdcC = (zNCEnergy>minEnergy || zPCEnergy>minEnergy);
+  }
+
 
   // Some macros for the online triggers
   Bool_t cMBS2A = c0sm2 && c0v0A;
@@ -156,6 +228,72 @@ void AliAnalysisTaskTriggerStudy::UserExec(Option_t *)
   vdArray[kVDC0VBC]  = c0v0C;
   //vdArray[kVDC0OM2]  = c0OM2;
 
+  // Plots requested by Jurgen on 18/11/2010
+  // 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(c0v0A && c0v0C) {
+    GetHistoSPD1  ("V0AND", "V0A & V0C")->Fill(outerLayerSPD);
+    GetHistoTracks("V0AND", "V0A & V0C")->Fill(ntracks);
+    GetHistoV0M   ("V0AND", "V0A & V0C")->Fill(multV0);
+  }
+  if(zdcA && zdcC) {
+    GetHistoSPD1  ("ZDCAND", "ZDCA & ZDCC")->Fill(outerLayerSPD);
+    GetHistoTracks("ZDCAND", "ZDCA & ZDCC")->Fill(ntracks);
+    GetHistoV0M   ("ZDCAND", "ZDCA & ZDCC")->Fill(multV0);
+  }
+  if((c0v0A && c0v0C) && !(zdcA && zdcC)) {
+    GetHistoSPD1  ("V0ANDNOTZDCAND", "(V0A & V0C) & !(ZDCA & ZDCC)")->Fill(outerLayerSPD);
+    GetHistoTracks("V0ANDNOTZDCAND", "(V0A & V0C) & !(ZDCA & ZDCC)")->Fill(ntracks);
+    GetHistoV0M   ("V0ANDNOTZDCAND", "(V0A & V0C) & !(ZDCA & ZDCC)")->Fill(multV0);
+  }
+  if((c0v0A && c0v0C) && !zdcA && !zdcC) {
+    GetHistoSPD1  ("V0ANDNOTZDC", "(V0A & V0C) & !ZDCA & !ZDCC)")->Fill(outerLayerSPD);
+    GetHistoTracks("V0ANDNOTZDC", "(V0A & V0C) & !ZDCA & !ZDCC)")->Fill(ntracks);
+    GetHistoV0M   ("V0ANDNOTZDC", "(V0A & V0C) & !ZDCA & !ZDCC)")->Fill(multV0);
+}
+  if((c0v0A && c0v0C) && ((!zdcA && zdcC) || (zdcA && !zdcC))) {
+    GetHistoSPD1  ("V0ANDONEZDC", "(V0A & V0C) &  ((!ZDCA && ZDCC) || (ZDCA && !ZDCC)))")->Fill(outerLayerSPD);
+    GetHistoTracks("V0ANDONEZDC", "(V0A & V0C) &  ((!ZDCA && ZDCC) || (ZDCA && !ZDCC)))")->Fill(ntracks);
+    GetHistoV0M   ("V0ANDONEZDC", "(V0A & V0C) &  ((!ZDCA && ZDCC) || (ZDCA && !ZDCC)))")->Fill(multV0);
+  }
+
+  if(((c0v0A && !c0v0C) || (!c0v0A && c0v0C)) && (zdcA && zdcC)) {
+    GetHistoSPD1  ("V0ONEZDCBOTH", "((V0A && !V0C) || (!V0A && V0C)) && (ZDCA && ZDCC)")->Fill(outerLayerSPD);
+    GetHistoTracks("V0ONEZDCBOTH", "((V0A && !V0C) || (!V0A && V0C)) && (ZDCA && ZDCC)")->Fill(ntracks);
+    GetHistoV0M   ("V0ONEZDCBOTH", "((V0A && !V0C) || (!V0A && V0C)) && (ZDCA && ZDCC)")->Fill(multV0);
+  }
+
+  if((c0v0A && c0v0C) && (zdcA && zdcC)) {
+    GetHistoSPD1  ("V0ANDZDCAND", "(V0A & V0C) & (ZDCA & ZDCC)")->Fill(outerLayerSPD);
+    GetHistoTracks("V0ANDZDCAND", "(V0A & V0C) & (ZDCA & ZDCC)")->Fill(ntracks);
+    GetHistoV0M   ("V0ANDZDCAND", "(V0A & V0C) & (ZDCA & ZDCC)")->Fill(multV0);
+  }
+
+  if((c0v0A && zdcA && !zdcC && !c0v0C) || (c0v0C && zdcC && !zdcA && !c0v0A)) {
+    GetHistoSPD1  ("OneSided", "(V0A & ZDCA & !ZDCC & !V0C) || (V0C & ZDCC & !ZDCA & !V0A)")->Fill(outerLayerSPD);
+    GetHistoTracks("OneSided", "(V0A & ZDCA & !ZDCC & !V0C) || (V0C & ZDCC & !ZDCA & !V0A)")->Fill(ntracks);
+    GetHistoV0M   ("OneSided", "(V0A & ZDCA & !ZDCC & !V0C) || (V0C & ZDCC & !ZDCA & !V0A)")->Fill(multV0);
+  }
+
+  // GetHistoCorrelationSPDTPCVz("All", "After physics selection and Francesco's cut")->Fill(vtxESDSPD->GetZ(),vtxESDTPC->GetZ());
+  // if(cut1T0)   GetHistoCorrelationSPDTPCVz("Cut1T0", "T0 Cut 1")->Fill(vtxESDSPD->GetZ(),vtxESDTPC->GetZ());
+  // if(cut2T0)   GetHistoCorrelationSPDTPCVz("Cut2T0", "T0 Cut 2")->Fill(vtxESDSPD->GetZ(),vtxESDTPC->GetZ());
+
+  // GetHistoCorrelationContrTPCSPDCls("All", "After physics selection and Francesco's cut")->Fill(totalClusSPD,tpcContr);
+  // if(cut1T0)   GetHistoCorrelationContrTPCSPDCls("Cut1T0", "T0 Cut 1")->Fill(totalClusSPD,tpcContr);
+  // if(cut2T0)   GetHistoCorrelationContrTPCSPDCls("Cut2T0", "T0 Cut 2")->Fill(totalClusSPD,tpcContr);
+
+  // GetHistoCorrelationTrackletsSPDCls("All", "After physics selection and Francesco's cut")->Fill(totalClusSPD,ntracklets);
+  // if(cut1T0)   GetHistoCorrelationTrackletsSPDCls("Cut1T0", "T0 Cut 1")->Fill(totalClusSPD,ntracklets);
+  // if(cut2T0)   GetHistoCorrelationTrackletsSPDCls("Cut2T0", "T0 Cut 2")->Fill(totalClusSPD,ntracklets);
+
+
+  return; // FIXME
+
+
   // Reject background
   if (v0BG && fRejectBGWithV0) {
     cout << "Rejection BG" << endl;
@@ -355,6 +493,89 @@ TH1 *   AliAnalysisTaskTriggerStudy::GetHistoTracklets(const char * name, const
   return h;
 }
 
+
+TH1 *   AliAnalysisTaskTriggerStudy::GetHistoSPD1(const char * name, const char * title){
+  // Book histo of events vs ntracklets, if needed
+
+  TString hname = "hSPD1_";
+  hname+=name;  
+  hname+=fHistoSuffix;
+  TH1 * h = (TH1*) fHistoList->GetList()->FindObject(hname.Data());
+  
+  // const Int_t nbin = 118;
+  // Float_t bins[nbin] = {0};
+  // bins[0] = 0;
+  // for(Int_t ibin = 1; ibin <= nbin; ibin++){
+  //   if (ibin < 100) 
+  //     bins[ibin] = bins [ibin-1]+1; 
+  //   else if (ibin < 1000) 
+  //     bins[ibin] = bins [ibin-1]+100; 
+  //   else if (ibin < 10000) 
+  //     bins[ibin] = bins [ibin-1]+1000; 
+  // }
+  
+
+  if(!h) {
+    AliInfo(Form("Booking histo %s",hname.Data()));
+    Bool_t oldStatus = TH1::AddDirectoryStatus();
+    TH1::AddDirectory(kFALSE);
+    //    h = new TH1F (hname.Data(), title, nbin, bins);
+    //h = new TH1F (hname.Data(), title, 100, -0.5, 100-0.5);
+    h = new TH1F (hname.Data(), title, 10000, -0.5, 10000-0.5);
+    h->Sumw2();
+    h->SetXTitle("SPD1");
+    fHistoList->GetList()->Add(h);
+    TH1::AddDirectory(oldStatus);
+  }
+  return h;
+}
+
+TH1 *   AliAnalysisTaskTriggerStudy::GetHistoV0M(const char * name, const char * title){
+  // Book histo of events vs ntracklets, if needed
+
+  TString hname = "hV0M_";
+  hname+=name;  
+  hname+=fHistoSuffix;
+  TH1 * h = (TH1*) fHistoList->GetList()->FindObject(hname.Data());
+  
+  if(!h) {
+    AliInfo(Form("Booking histo %s",hname.Data()));
+    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->Sumw2();
+    h->SetXTitle("V0M");
+    fHistoList->GetList()->Add(h);
+    TH1::AddDirectory(oldStatus);
+  }
+  return h;
+}
+
+
+TH1 *   AliAnalysisTaskTriggerStudy::GetHistoTracks(const char * name, const char * title){
+  // Book histo of events vs ntracklets, if needed
+
+  TString hname = "hTPCTracks_";
+  hname+=name;  
+  hname+=fHistoSuffix;
+  TH1 * h = (TH1*) fHistoList->GetList()->FindObject(hname.Data());
+  
+  if(!h) {
+    AliInfo(Form("Booking histo %s",hname.Data()));
+    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, 50, -0.5, 50-0.5);
+    h->Sumw2();
+    h->SetXTitle("ntracks");
+    fHistoList->GetList()->Add(h);
+    TH1::AddDirectory(oldStatus);
+  }
+  return h;
+}
+
+
 TH1 *   AliAnalysisTaskTriggerStudy::GetHistoEta(const char * name, const char * title){
   // Book histo of events vs ntracklets, if needed
 
@@ -399,6 +620,75 @@ TH1 *   AliAnalysisTaskTriggerStudy::GetHistoPt(const char * name, const char *
   return h;
 }
 
+TH1 *   AliAnalysisTaskTriggerStudy::GetHistoCorrelationSPDTPCVz(const char * name, const char * title){
+  // Book histo of events vs ntracklets, if needed
+
+  TString hname = "hTPCvsSPD_";
+  hname+=name;  
+  hname+=fHistoSuffix;
+  TH1 * h = (TH1*) fHistoList->GetList()->FindObject(hname.Data());
+  
+  if(!h) {
+    AliInfo(Form("Booking histo %s",hname.Data()));
+    Bool_t oldStatus = TH1::AddDirectoryStatus();
+    TH1::AddDirectory(kFALSE);
+    // h = new TH1F (hname.Data(), title, 1000, -0.5, 10000-0.5);
+    h = new TH2F (hname.Data(), title, 80, -20, 20,  80, -20, 20);
+    h->Sumw2();
+    h->SetXTitle("SPD Vz");
+    h->SetYTitle("TPC Vz");
+    fHistoList->GetList()->Add(h);
+    TH1::AddDirectory(oldStatus);
+  }
+  return h;
+}
+
+TH1 *   AliAnalysisTaskTriggerStudy::GetHistoCorrelationContrTPCSPDCls(const char * name, const char * title){
+  // Book histo of events vs ntracklets, if needed
+
+  TString hname = "hContrTPCvsSPDCls_";
+  hname+=name;  
+  hname+=fHistoSuffix;
+  TH1 * h = (TH1*) fHistoList->GetList()->FindObject(hname.Data());
+  
+  if(!h) {
+    AliInfo(Form("Booking histo %s",hname.Data()));
+    Bool_t oldStatus = TH1::AddDirectoryStatus();
+    TH1::AddDirectory(kFALSE);
+    // h = new TH1F (hname.Data(), title, 1000, -0.5, 10000-0.5);
+    h = new TH2F (hname.Data(), title, 1000, 0, 9000,  1000, 0, 5000);
+    h->Sumw2();
+    h->SetXTitle("SPD Clusters");
+    h->SetYTitle("TPC V Contributors");
+    fHistoList->GetList()->Add(h);
+    TH1::AddDirectory(oldStatus);
+  }
+  return h;
+}
+TH1 *   AliAnalysisTaskTriggerStudy::GetHistoCorrelationTrackletsSPDCls(const char * name, const char * title){
+  // Book histo of events vs ntracklets, if needed
+
+  TString hname = "hTrackletsvsSPDCls_";
+  hname+=name;  
+  hname+=fHistoSuffix;
+  TH1 * h = (TH1*) fHistoList->GetList()->FindObject(hname.Data());
+  
+  if(!h) {
+    AliInfo(Form("Booking histo %s",hname.Data()));
+    Bool_t oldStatus = TH1::AddDirectoryStatus();
+    TH1::AddDirectory(kFALSE);
+    // h = new TH1F (hname.Data(), title, 1000, -0.5, 10000-0.5);
+    h = new TH2F (hname.Data(), title, 1000, 0, 9000,  1000, 0, 5000);
+    h->Sumw2();
+    h->SetXTitle("SPD Clusters");
+    h->SetYTitle("N tracklets");
+    fHistoList->GetList()->Add(h);
+    TH1::AddDirectory(oldStatus);
+  }
+  return h;
+}
+
+
 void AliAnalysisTaskTriggerStudy::FillTriggerOverlaps (const char * name, const char * title, Bool_t * vdArray){
   //Fills a histo with the different trigger statistics in a venn like diagramm. Books it if needed.
 
index 3aede4ccd7c93f466919c68a44ee667f598057e3..232e89145d74dbee156077eb9b6e7dee1b98ab62 100644 (file)
@@ -39,6 +39,12 @@ public:
   TH1 * GetHistoTracklets   (const char * name, const char * title);
   TH1 * GetHistoPt(const char * name, const char * title);
   TH1 * GetHistoEta(const char * name, const char * title);
+  TH1 * GetHistoV0M(const char * name, const char * title);
+  TH1 * GetHistoSPD1(const char * name, const char * title);
+  TH1 * GetHistoTracks(const char * name, const char * title);
+  TH1 * GetHistoCorrelationSPDTPCVz(const char * name, const char * title);
+  TH1 * GetHistoCorrelationContrTPCSPDCls  (const char * name, const char * title);
+  TH1 * GetHistoCorrelationTrackletsSPDCls (const char * name, const char * title);
   void  FillTriggerOverlaps (const char * name, const char * title, Bool_t * vdArray) ;
 
   void SetNTrackletsCut(Int_t cut){ fNTrackletsCut = cut;}
index 35f9ec39a7259340de63c262001d4c99e72c70e0..7a24c1333f9c17a23a9ca76b36216af3b509671b 100644 (file)
@@ -20,7 +20,7 @@ AliAnalysisMultPbTrackHistoManager * hManCorr = 0;
 TH2D * hEvStatData = 0;
 TH2D * hEvStatCorr = 0;
 
-const Int_t kHistoFitCompoments = 3;
+const Int_t kHistoFitCompoments = 2;
 TH1D * gHistoCompoments[kHistoFitCompoments];
 
 void LoadLibs(  Bool_t debug=0);
@@ -69,7 +69,7 @@ void correct(TString dataFolder = "./output/LHC10g2d_130844_V0M_bin_10/", TStrin
   CheckSanity();
 
   Double_t fractionWeak = 1, fractionMaterial=1; 
-  //  CheckSecondaries(fractionWeak, fractionMaterial);
+  CheckSecondaries(fractionWeak, fractionMaterial);
   cout << "Rescaling secondaries correction, weak: " << fractionWeak << ", material: " << fractionMaterial <<endl;
   if(scaleWeakByHand >= 0) {
     fractionWeak = scaleWeakByHand;
@@ -210,8 +210,8 @@ void correct(TString dataFolder = "./output/LHC10g2d_130844_V0M_bin_10/", TStrin
   hCorrected->SetMarkerColor(kBlack);
   hCorrected->Draw("same");
   hMCPtGen->DrawClone("same");
-  TF1 * f = GetLevy();
-  //TF1 * f = GetHagedorn();
+  //  TF1 * f = GetLevy();
+  TF1 * f = GetHagedorn();
   hCorrected->Fit(f,"", "same");
   hCorrected->Fit(f,"IME", "same",0,2);
   cout << "dN/deta (function)  = " << f->Integral(0,100) << " +- " << f->IntegralError(0,100) << endl;
@@ -225,8 +225,8 @@ void correct(TString dataFolder = "./output/LHC10g2d_130844_V0M_bin_10/", TStrin
   Double_t errorData = 0;
   Float_t dNdeta  =  (f->Integral(0,0.1) + hCorrected->IntegralAndError(1,-1,errorData, "width")) / (etaMax-etaMin);
   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)/Npart  " << dNdeta/npart << " +- " << dNdetaE/npart << endl;
+  cout << "dN/deta              " << dNdeta << " +- " << dNdetaE << endl;
+  cout << "(dN/deta)/0.5 Npart  " << dNdeta/npart*2 << " +- " << dNdetaE/npart*2 << endl;
   
   TH1D * hDataGen  = hManData->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoGen,        etaMin,etaMax,zmin,zmax);
   cout << "Generated dN/deta (data) =       " << hDataGen->Integral("width") << endl;
@@ -235,7 +235,7 @@ void correct(TString dataFolder = "./output/LHC10g2d_130844_V0M_bin_10/", TStrin
   l->AddEntry(hDataPt, "Raw data");
   l->AddEntry(hCorrected, "Corrected data");
   l->AddEntry(hMCPtGen, "Monte Carlo (generated)");
-  l->AddEntry(f, "Levy Fit");
+  l->AddEntry(f, "Hagedorn Fit");
   l->Draw();
 }
 
@@ -243,24 +243,25 @@ void CheckSecondaries(Double_t &fracWeak, Double_t &fracMaterial) {
   // Returns the fraction you need to rescale the secondaries from weak decays for.
 
   // Some shorthands
-  TH1D * hDataDCA   = hManData->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec       );
+  TH2D * hDataDCAPt   = hManData->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec       );
   //  TH1D * hMCDCAGen  = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoGen       );
-  TH1D * hMCDCARec  = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec       );
-  TH1D * hMCDCAPri  = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim   );
-  TH1D * hMCDCASW  = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak    );
-  TH1D * hMCDCASM  = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat    );
-  TH1D * hMCDCAFak  = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecFake   );
+  TH2D * hMCDCARecPt  = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec       );
+  TH2D * hMCDCAPriPt  = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim   );
+  TH2D * hMCDCASWPt   = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak    );
+  TH2D * hMCDCASMPt   = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat    );
+  TH2D * hMCDCAFakPt  = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecFake   );
  
 
   TCanvas * cCumulative = new TCanvas("cDCAculumative","DCA cumulative distributions");
   cCumulative->cd();
-  GetCumulativeHisto(hMCDCAPri )->Draw();
-  GetRatioIntegratedFractions(hMCDCASW, hMCDCARec  )->Draw("same");
-  GetRatioIntegratedFractions(hMCDCASM, hMCDCARec  )->Draw("same");
-  GetRatioIntegratedFractions(hMCDCAPri,hMCDCARec  )->Draw("same");
+  GetCumulativeHisto(hMCDCAPriPt )->Draw();
+  GetRatioIntegratedFractions(hMCDCASWPt, hMCDCARecPt  )->Draw("same");
+  GetRatioIntegratedFractions(hMCDCASMPt, hMCDCARecPt  )->Draw("same");
+  GetRatioIntegratedFractions(hMCDCAPriPt,hMCDCARecPt  )->Draw("same");
 
 
   TCanvas * c1 = new TCanvas("cDCAFit", "Fit to the DCA distributions");  
+  c1->Divide(4,2);
   c1->SetLogy();
   // Draw all
   //  hDataDCA->Draw();
@@ -272,34 +273,62 @@ void CheckSecondaries(Double_t &fracWeak, Double_t &fracMaterial) {
   // return;
   
   // Fit the DCA distribution. Uses a TF1 made by summing histograms
-  TH1D * hMCPrimSMFak = (TH1D*) hMCDCAPri->Clone("hMCPrimSMFak");
-  //  hMCPrimSMFak->Add(hMCDCASM);
-  hMCPrimSMFak->Add(hMCDCAFak);
+  TH2D * hMCPrimSMFakPt = (TH2D*) hMCDCAPriPt->Clone("hMCPrimSMFak");
+  hMCPrimSMFakPt->Add(hMCDCASMPt);
+  hMCPrimSMFakPt->Add(hMCDCAFakPt);
 
   // Set the components which are used in HistoSum, the static
   // function for GetFunctionHistoSum
-  gHistoCompoments[0] = (TH1D*) hMCPrimSMFak->Clone();
-  gHistoCompoments[1] = (TH1D*) hMCDCASW->Clone();
-  gHistoCompoments[2] = (TH1D*) hMCDCASM->Clone();
-  TF1 * fHistos = GetFunctionHistoSum();
-  fHistos->SetParameters(1,1,1);
-  //fHistos->FixParameter(2,1);
-  fHistos->SetLineColor(kRed);
-  hDataDCA->Draw();
-  // Fit!
-  //  hDataDCA->Fit(fHistos,"","",0,200);
-  // Rescale the components and draw to see how it looks like
-  hMCPrimSMFak->Scale(fHistos->GetParameter(0));
-  hMCDCASW    ->Scale(fHistos->GetParameter(1));
-  hMCDCASM    ->Scale(fHistos->GetParameter(2));
-  hMCPrimSMFak->Draw("same");
-  hMCDCASW    ->Draw("same");
-  hMCDCASM    ->Draw("same");
-  fHistos->Draw("same");
-  // compute scaling factors
-  fracWeak     = fHistos->GetParameter(1)/fHistos->GetParameter(0);
-  fracMaterial = fHistos->GetParameter(2)/fHistos->GetParameter(0);
-
+  // Project onti DCA axis
+  const Int_t ptBinsFit[] = {3,5,7,9,11,15,19,23,31,-1};
+  Int_t ibinPt = -1;
+  while(ptBinsFit[(++ibinPt)+1]!=-1){
+    c1->cd(ibinPt+1);
+    gPad->SetLogy();
+    Int_t minPtBin=ptBinsFit[ibinPt];
+    Int_t maxPtBin=ptBinsFit[ibinPt+1]-1;
+  
+    TH1D * hDataDCA     = (TH1D*) hDataDCAPt    ->ProjectionX(TString(hDataDCAPt    ->GetName())+long(ibinPt)+"_px",minPtBin,maxPtBin,"E")->Clone();
+    TH1D * hMCPrimSMFak = (TH1D*) hMCPrimSMFakPt->ProjectionX(TString(hMCPrimSMFakPt->GetName())+long(ibinPt)+"_px",minPtBin,maxPtBin,"E")->Clone();
+    TH1D * hMCDCASW     = (TH1D*) hMCDCASWPt    ->ProjectionX(TString(hMCDCASWPt    ->GetName())+long(ibinPt)+"_px",minPtBin,maxPtBin,"E")->Clone();
+    TH1D * hMCDCASM     = (TH1D*) hMCDCASMPt    ->ProjectionX(TString(hMCDCASMPt    ->GetName())+long(ibinPt)+"_px",minPtBin,maxPtBin,"E")->Clone();
+    gHistoCompoments[0] = (TH1D*) hMCPrimSMFakPt->ProjectionX(TString(hMCPrimSMFakPt->GetName())+long(ibinPt)+"_px",minPtBin,maxPtBin,"E")->Clone();
+    gHistoCompoments[1] = (TH1D*) hMCDCASWPt    ->ProjectionX(TString(hMCDCASWPt    ->GetName())+long(ibinPt)+"_px",minPtBin,maxPtBin,"E")->Clone();
+    //  gHistoCompoments[2] = (TH1D*) hMCDCASMPt    ->ProjectionX("_px",minPtBin,maxPtBin,"E")->Clone();
+    TF1 * fHistos = GetFunctionHistoSum();
+    fHistos->SetParameters(1,1);
+    //  fHistos->FixParameter(2,1);
+    fHistos->SetLineColor(kRed);
+    hDataDCA= (TH1D*) hDataDCAPt->ProjectionX(TString(hDataDCAPt->GetName())+long(ibinPt)+"_px",minPtBin,maxPtBin,"E")->Clone();
+    hDataDCA->Draw();
+    // Fit!
+    hDataDCA->Fit(fHistos,"0Q","",-2,2);
+    // Rescale the components and draw to see how it looks like
+    hMCPrimSMFak->Scale(fHistos->GetParameter(0));
+    hMCDCASW    ->Scale(fHistos->GetParameter(1));
+    hMCDCASM    ->Scale(fHistos->GetParameter(2));
+    hMCPrimSMFak->Draw("same");
+    hMCDCASW    ->Draw("same");
+    hMCDCASM    ->Draw("same");
+    TH1D* hMCSum = (TH1D*) hMCPrimSMFak->Clone();
+    hMCSum->Add(hMCDCASW);
+    //    hMCSum->Add(hMCDCASM);
+    hMCSum->SetLineColor(kRed);
+    hMCSum->SetMarkerStyle(1);
+    hMCSum->Draw("lhist,same");
+    //    fHistos->Draw("same");
+    // compute scaling factors
+    //  fracWeak     = fHistos->GetParameter(1)/(fHistos->GetParameter(0)+fHistos->GetParameter(1));
+    cout << hDataDCAPt->GetYaxis()->GetBinLowEdge(minPtBin) << " - ";
+    cout << hDataDCAPt->GetYaxis()->GetBinLowEdge(maxPtBin+1) << " = ";
+    cout << 2*fHistos->GetParameter(1)/(fHistos->GetParameter(0)+fHistos->GetParameter(1));
+    cout << " ("<< fHistos->GetChisquare() <<")";
+    cout << " " << fHistos->GetParameter(0) << "," << fHistos->GetParameter(1) << endl;
+    
+    
+   fracWeak     = 2*fHistos->GetParameter(1)/(fHistos->GetParameter(0)+fHistos->GetParameter(1));
+     //  fracMaterial = fHistos->GetParameter(2)/(fHistos->GetParameter(0);
+  }
 
 }
 
@@ -379,7 +408,7 @@ void CheckSanity() {
 
 void CheckCorrections() {
 
-  TH1D * hDataPt   = (TH1D*) hManData->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRec, etaMin,etaMax,zmin,zmax)->Clone("hDataPt");
+  //  TH1D * hDataPt   = (TH1D*) hManData->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRec, etaMin,etaMax,zmin,zmax)->Clone("hDataPt");
   TH1D * hMCPtGen  = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoGen,         etaMin,etaMax,zmin,zmax);
   TH1D * hMCPtRec  = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRec,         etaMin,etaMax,zmin,zmax);
   TH1D * hMCPtPri  = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim,     etaMin,etaMax,zmin,zmax);
@@ -456,8 +485,8 @@ void SetStyle() {
   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetMarkerColor(kRed  );
   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->SetMarkerColor(kGreen);
   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat) ->SetMarkerColor(kBlue );
-  hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetMarkerColor(kBlue );
-  hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetMarkerColor(kCyan );
+  hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetMarkerColor(kCyan );
+  hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetMarkerColor(kYellow );
 
   hManData->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetMarkerStyle(kFullCircle);    
   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoGen)    ->SetMarkerStyle(kFullSquare);
@@ -472,16 +501,16 @@ void SetStyle() {
   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetLineColor(kRed  );
   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->SetLineColor(kGreen);
   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat) ->SetLineColor(kBlue );
-  hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetLineColor(kBlue-7 );
-  hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetLineColor(kCyan );
+  hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetLineColor(kCyan);
+  hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetLineColor(kYellow );
 
   hManData->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetMarkerColor(kBlack);    
   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoGen)    ->SetMarkerColor(kRed  );
   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetMarkerColor(kRed  );
   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->SetMarkerColor(kGreen);
   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat) ->SetMarkerColor(kBlue );
-  hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetMarkerColor(kBlue-7 );
-  hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetMarkerColor(kCyan );
+  hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetMarkerColor(kCyan );
+  hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetMarkerColor(kYellow );
 
   hManData->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetMarkerStyle(kFullCircle);    
   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoGen)    ->SetMarkerStyle(kFullSquare);
@@ -659,7 +688,7 @@ void ShowAcceptanceInVzSlices() {
 
 TF1 * GetFunctionHistoSum() {
 
-  TF1 * f = new TF1 ("fHistoSum",HistoSum,0,50,kHistoFitCompoments);
+  TF1 * f = new TF1 ("fHistoSum",HistoSum,-1,1,kHistoFitCompoments);
   f->SetParNames("Primaries+Fakes", "Sec. Weak decays","Sec. Material");
   f->SetNpx(1000);
   return f;
index c9dcec9da899b3b02454f3e82336fa32ce57ac77..ab2946e3612046fe319641a46cc16b3df4b0f355 100644 (file)
@@ -109,38 +109,35 @@ void run(Char_t* data, Long64_t nev = -1, Long64_t offset = 0, Bool_t debug = kF
     doSave = kTRUE;
   }
 
-  AliESDtrackCuts * cuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kFALSE); //FIXME
-  //  AliESDtrackCuts * cuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
-  TString pathsuffix = "";
-  // cuts->SetPtRange(0.15,0.2);// FIXME pt cut
-  // const char * pathsuffix = "_pt_015_020_nofakes";
+  AliESDtrackCuts * cuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kFALSE);  TString pathsuffix = "";
 
   if (optionStr.Contains("DCA")) {
     delete cuts;
     cuts = AliESDtrackCuts::GetStandardITSPureSATrackCuts2009();
     cout << ">>>> USING DCA cut" << endl;
     cuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
-    pathsuffix="DCAcut";
+    pathsuffix+="_DCAcut";
   }
 
   if (optionStr.Contains("ITSsa")) {
     delete cuts;
     cuts = AliESDtrackCuts::GetStandardITSPureSATrackCuts2009();
     cout << ">>>> USING ITS sa tracks" << endl;
-    pathsuffix="ITSsa";
+    pathsuffix+="_ITSsa";
   }
 
   if (optionStr.Contains("TPC")) {
     delete cuts;
     cuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
     cout << ">>>> USING TPC only tracks" << endl;
-    pathsuffix="TPC";
+    pathsuffix+="_TPC";
   }
 
   Bool_t useMCKinematics = isMC;
   if (optionStr.Contains("NOMCKIN")) {
     cout << ">>>> Ignoring MC kinematics" << endl;
     useMCKinematics=kFALSE;
+    pathsuffix+="_NOMCKIN";
   }
   
   
index f12c8904de7fb50aa63d9b1a721be3a52ea20b0e..94fa8f60217a5384673036cd6baf62d938b7ae64 100644 (file)
@@ -41,6 +41,13 @@ void runTriggerStudy(Char_t* data, Long64_t nev = -1, Long64_t offset = 0, Bool_
     mgr->SetGridHandler(alienHandler);  
   }
 
+  // Add physics selection
+  gROOT->ProcessLine(".L $ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");
+  physicsSelectionTask = AddTaskPhysicsSelection(isMC,0);//FIXME
+
+
+
+
   // Parse option strings
   TString optionStr(option);
   
@@ -149,7 +156,8 @@ void InitAndLoadLibs(Int_t runMode=kMyRunModeLocal, Int_t workers=0,Bool_t debug
     cout << "Init in CAF mode" << endl;
     
     gEnv->SetValue("XSec.GSI.DelegProxy", "2");
-    TProof::Open("alice-caf.cern.ch", workers>0 ? Form("workers=%d",workers) : "");
+    //    TProof::Open("alice-caf.cern.ch", workers>0 ? Form("workers=%d",workers) : "");
+    TProof::Open("skaf.saske.sk", workers>0 ? Form("workers=%d",workers) : "");
     
     // Enable the needed package
     gProof->UploadPackage("$ALICE_ROOT/STEERBase");