]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGLF/totEt/AliAnalysisEtReconstructed.cxx
cleaning up streamers
[u/mrichter/AliRoot.git] / PWGLF / totEt / AliAnalysisEtReconstructed.cxx
index 18cc6fa44a8599b40f1d1379cc5a36739ea4a41d..4c024f975e8bbd0ff628ea1b124e55e58e408154 100644 (file)
@@ -73,6 +73,7 @@ AliAnalysisEtReconstructed::AliAnalysisEtReconstructed() :
        ,fHistMatchedTracksEvspTvsCent(0)
        ,fHistMatchedTracksEvspTvsCentEffCorr(0)
        ,fHistMatchedTracksEvspTvsCentEffTMCorr(0)
+       ,fHistPeripheralMatchedTracksEvspTvsCentEffTMCorr(0)
        ,fHistMatchedTracksEvspTvsCentEffTMCorr500MeV(0)
        ,fHistFoundHadronsvsCent(0)
        ,fHistNotFoundHadronsvsCent(0)
@@ -112,6 +113,7 @@ AliAnalysisEtReconstructed::AliAnalysisEtReconstructed() :
        ,fHistCentVsNchVsNclReco(0)
        ,fHistRawSignalReco(0)
        ,fHistEffCorrSignalReco(0)
+       ,fHistRecoRCorrVsPtVsCent(0)
 {
 
 }
@@ -141,6 +143,7 @@ AliAnalysisEtReconstructed::~AliAnalysisEtReconstructed()
     delete fHistMatchedTracksEvspTvsCent;
     delete fHistMatchedTracksEvspTvsCentEffCorr;
     delete fHistMatchedTracksEvspTvsCentEffTMCorr;
+    delete fHistPeripheralMatchedTracksEvspTvsCentEffTMCorr;
     delete fHistMatchedTracksEvspTvsCentEffTMCorr500MeV;
     delete fHistFoundHadronsvsCent;
     delete fHistNotFoundHadronsvsCent;
@@ -179,6 +182,7 @@ AliAnalysisEtReconstructed::~AliAnalysisEtReconstructed()
     delete fHistCentVsNchVsNclReco;
     delete fHistRawSignalReco;
     delete fHistEffCorrSignalReco;
+    delete fHistRecoRCorrVsPtVsCent;
 }
 
 Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
@@ -213,8 +217,8 @@ Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
 
 
   //for PID
-  AliESDpid *pID = new AliESDpid();
-  pID->MakePID(event);
+  //AliESDpid *pID = new AliESDpid();
+  //pID->MakePID(event);
   Float_t etPIDProtons = 0.0;
   Float_t etPIDAntiProtons = 0.0;
   Float_t etPiKPMatched = 0.0;
@@ -237,7 +241,7 @@ Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
          }
        else{
          totalPt +=track->Pt();
-         pID->MakeITSPID(track);
+         //pID->MakeITSPID(track);
 
 
        }
@@ -253,8 +257,8 @@ Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
     Float_t effHighRawEt = 0;
     Float_t effLowRawEt = 0;
     Float_t uncorrEt = 0;
-    Float_t rawSignal;
-    Float_t effCorrSignal;
+    Float_t rawSignal = 0;
+    Float_t effCorrSignal = 0;
 
     nChargedHadronsMeasured = 0.0;
     nChargedHadronsTotal = 0.0;
@@ -300,14 +304,14 @@ Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
         cluster->GetPosition(pos);
         TVector3 cp(pos);
        fClusterPositionAll->Fill(cp.Phi(), cp.PseudoRapidity());
-       fClusterPositionAllEnergy->Fill(cp.Phi(), cp.PseudoRapidity(),cluster->E());
+       Float_t fReconstructedE = cluster->E();
+       fClusterPositionAllEnergy->Fill(cp.Phi(), cp.PseudoRapidity(),GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE);
 
        //if(TMath::Abs(cp.Eta())> fCuts->fCuts->GetGeometryEmcalEtaAccCut() || cp.Phi() >  fCuts->GetGeometryEmcalPhiAccMaxCut()*TMath::Pi()/180. ||  cp.Phi() >  fCuts->GetGeometryEmcalPhiAccMinCut()*TMath::Pi()/180.) continue;//Do not accept if cluster is not in the acceptance
-       fTotAllRawEt += TMath::Sin(cp.Theta())*cluster->E();
-       fTotAllRawEtEffCorr += CorrectForReconstructionEfficiency(*cluster,cent);
-
-       fClusterEnergyCent->Fill(cluster->E(),cent);
+       fTotAllRawEt += TMath::Sin(cp.Theta())*GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE;
+       fTotAllRawEtEffCorr +=GetCorrectionModification(*cluster,0,0,cent)* CorrectForReconstructionEfficiency(*cluster,cent);
 
+       fClusterEnergyCent->Fill(GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE,cent);
         Bool_t matched = kTRUE;//default to no track matched
        Int_t trackMatchedIndex = cluster->GetTrackMatchedIndex();//find the index of the matched track
        matched = !(fSelector->PassTrackMatchingCut(*cluster));//PassTrackMatchingCut is false if there is a matched track
@@ -317,6 +321,15 @@ Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
            AliESDtrack *track = event->GetTrack(trackMatchedIndex);
            //if this is a good track, accept track will return true.  The track matched is good, so not track matched is false
            matched = fEsdtrackCutsTPC->AcceptTrack(track);//If the track is bad, don't count it
+           if(matched){//if it is still matched see if the track p was less than the energy
+             Float_t rcorr = TMath::Sin(cp.Theta())*GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE;
+             fHistRecoRCorrVsPtVsCent->Fill(rcorr,track->Pt(), fCentClass);
+             if(fSelector->PassMinEnergyCut( (fReconstructedE - fsub* track->P())*TMath::Sin(cp.Theta()) )){//if more energy was deposited than the momentum of the track  and more than one particle led to the cluster
+               //            if(fReconstructedE - fsub* track->P() > 0.0){
+               fReconstructedE = fReconstructedE - fsub* track->P();
+               matched = kFALSE;
+             }
+           }
          }
        }
 
@@ -332,7 +345,7 @@ Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
                 }
                 else {
                  totalMatchedPt +=track->Pt();
-                 fClusterEnergyCentMatched->Fill(cluster->E(),cent);
+                 fClusterEnergyCentMatched->Fill(GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE,cent);
                  fHistMatchedClusterSizeVsCent->Fill(cluster->GetNCells(),cent);
 
                  float eff = fTmCorrections->TrackMatchingEfficiency(track->Pt(),cent);
@@ -340,41 +353,49 @@ Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
                  //cout<<"pt "<<track->Pt()<<" eff "<<eff<<" total "<<nChargedHadronsTotal<<endl;
                  nChargedHadronsMeasured++;
                  nChargedHadronsTotal += 1/eff;
-                 Double_t effCorrEt = CorrectForReconstructionEfficiency(*cluster,cent);
-                 nChargedHadronsEtMeasured+= TMath::Sin(cp.Theta())*cluster->E();
+                 Double_t effCorrEt = GetCorrectionModification(*cluster,0,0,cent) * CorrectForReconstructionEfficiency(*cluster,fReconstructedE,cent);
+                 nChargedHadronsEtMeasured+= TMath::Sin(cp.Theta())*GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE;
                  //One efficiency is the gamma efficiency and the other is the track matching efficiency.
-                 nChargedHadronsEtTotal+= 1/eff *TMath::Sin(cp.Theta())*cluster->E();
-                 //cout<<"nFound "<<1<<" nFoundTotal "<<1/eff<<" etMeas "<<TMath::Sin(cp.Theta())*cluster->E()<<" ET total "<< 1/eff *TMath::Sin(cp.Theta())*cluster->E()<<endl;
+                 nChargedHadronsEtTotal+= 1/eff *TMath::Sin(cp.Theta())*GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE;
+                 //cout<<"nFound "<<1<<" nFoundTotal "<<1/eff<<" etMeas "<<TMath::Sin(cp.Theta())*fReconstructedE<<" ET total "<< 1/eff *TMath::Sin(cp.Theta())*fReconstructedE<<endl;
 
                  Float_t nSigmaPion = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kPion); 
                  Float_t nSigmaProton = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kProton); 
                  bool isProton = (nSigmaPion>3.0 && nSigmaProton<3.0 && track->Pt()<0.9);
                  //cout<<"NSigmaProton "<<nSigmaProton<<endl;
                  etPiKPMatched += effCorrEt;
-                 etPiKPMatchedNoEff  +=TMath::Sin(cp.Theta())*cluster->E();
+                 etPiKPMatchedNoEff  +=TMath::Sin(cp.Theta())*GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE;
                  if(isProton){
                    if(track->Charge()>0){
                      etPIDProtons += effCorrEt;
-                     etPIDProtonsNoEff +=TMath::Sin(cp.Theta())*cluster->E();
+                     etPIDProtonsNoEff +=TMath::Sin(cp.Theta())*GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE;
                    }
                    else{
-                     etPIDAntiProtonsNoEff +=TMath::Sin(cp.Theta())*cluster->E();
+                     etPIDAntiProtonsNoEff +=TMath::Sin(cp.Theta())*GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE;
                      etPIDAntiProtons += effCorrEt;
                    }
                  }
-                 if(TMath::Sin(cp.Theta())*cluster->E()>0.5){
+                 if(TMath::Sin(cp.Theta())*GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE>0.5){
                    nChargedHadronsMeasured500MeV++;
                    nChargedHadronsTotal500MeV += 1/eff;
-                   nChargedHadronsEtMeasured500MeV+= TMath::Sin(cp.Theta())*cluster->E();
-                   nChargedHadronsEtTotal500MeV+= 1/eff *TMath::Sin(cp.Theta())*cluster->E();
+                   nChargedHadronsEtMeasured500MeV+= TMath::Sin(cp.Theta())*GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE;
+                   nChargedHadronsEtTotal500MeV+= 1/eff *TMath::Sin(cp.Theta())*GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE;
                  }
-                 fHistMatchedTracksEvspTvsCent->Fill(track->P(),TMath::Sin(cp.Theta())*cluster->E(),cent);
+                 fHistMatchedTracksEvspTvsCent->Fill(track->P(),TMath::Sin(cp.Theta())*GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE,cent);
                  fHistMatchedTracksEvspTvsCentEffCorr->Fill(track->P(),effCorrEt,cent);
                  //Weighed by the number of tracks we didn't find
                  fHistMatchedTracksEvspTvsCentEffTMCorr->Fill(track->P(), effCorrEt,cent, (1/eff-1) );
+                 if(cent<16 && cent>11){//centralities 60-80% where false track matches are low
+                   for(int cbtest = 0; cbtest<20; cbtest++){//then we calculate the deposit matched to hadrons with different centrality bins' efficiencies
+                     float efftest = fTmCorrections->TrackMatchingEfficiency(track->Pt(),cbtest);
+                     if(TMath::Abs(efftest)<1e-5) efftest = 1.0;
+                     Double_t effCorrEttest = GetCorrectionModification(*cluster,0,0,cent)*CorrectForReconstructionEfficiency(*cluster,fReconstructedE,cbtest);
+                     fHistPeripheralMatchedTracksEvspTvsCentEffTMCorr->Fill(track->P(), effCorrEttest,cbtest, (1/efftest-1) );
+                   }
+                 }
                  cluster->GetPosition(pos);      
                  TVector3 p2(pos);
-                 uncorrEt += TMath::Sin(p2.Theta())*cluster->E();
+                 uncorrEt += TMath::Sin(p2.Theta())*GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE;
                  if(uncorrEt>=0.5) fHistMatchedTracksEvspTvsCentEffTMCorr500MeV->Fill(track->P(), effCorrEt,cent, (1/eff-1) );
                     const Double_t *pidWeights = track->PID();
 
@@ -393,7 +414,7 @@ Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
                         }
                         if (fCuts->GetHistMakeTreeDeposit() && fDepositTree)
                         {
-                            fEnergyDeposited = cluster->E();
+                            fEnergyDeposited =GetCorrectionModification(*cluster,0,0,cent)* fReconstructedE;
                             fMomentumTPC = track->P();
                             fCharge = track->Charge();
                             fParticlePid = maxpid;
@@ -414,30 +435,30 @@ Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
 
                             //Float_t theta = TMath::ATan(pos[2]/dist)+TMath::Pi()/2;
 
-                            //Float_t et = cluster->E() * TMath::Sin(theta);
+                            //Float_t et = fReconstructedE * TMath::Sin(theta);
                             if (maxpid == AliPID::kProton)
                             {
 
                                 if (track->Charge() == 1)
                                 {
-                                    fHistProtonEnergyDeposit->Fill(cluster->E(), track->E());
+                                    fHistProtonEnergyDeposit->Fill(GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE, track->E());
                                 }
                                 else if (track->Charge() == -1)
                                 {
-                                    fHistAntiProtonEnergyDeposit->Fill(cluster->E(), track->E());
+                                    fHistAntiProtonEnergyDeposit->Fill(GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE, track->E());
                                 }
                             }
                             else if (maxpid == AliPID::kPion)
                             {
-                                fHistChargedPionEnergyDeposit->Fill(cluster->E(), track->E());
+                                fHistChargedPionEnergyDeposit->Fill(GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE, track->E());
                             }
                             else if (maxpid == AliPID::kKaon)
                             {
-                                fHistChargedKaonEnergyDeposit->Fill(cluster->E(), track->E());
+                                fHistChargedKaonEnergyDeposit->Fill(GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE, track->E());
                             }
                             else if (maxpid == AliPID::kMuon)
                             {
-                                fHistMuonEnergyDeposit->Fill(cluster->E(), track->E());
+                                fHistMuonEnergyDeposit->Fill(GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE, track->E());
                             }
                         }
                     }
@@ -449,26 +470,26 @@ Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
          fCutFlow->Fill(x++);
          //std::cout << x++ << std::endl;
          
-         //if (cluster->E() >  fSingleCellEnergyCut && cluster->GetNCells() == fCuts->GetCommonSingleCell()) continue;
-         //if (cluster->E() < fClusterEnergyCut) continue;
+         //if (fReconstructedE >  fSingleCellEnergyCut && cluster->GetNCells() == fCuts->GetCommonSingleCell()) continue;
+         //if (fReconstructedE < fClusterEnergyCut) continue;
          cluster->GetPosition(pos);
          
            TVector3 p2(pos);
            
            fClusterPositionAccepted->Fill(p2.Phi(), p2.PseudoRapidity());
-           fClusterPositionAcceptedEnergy->Fill(p2.Phi(), p2.PseudoRapidity(),cluster->E());
-           fClusterEnergy->Fill(cluster->E());
-           fClusterEnergyCentNotMatched->Fill(cluster->E(),cent);
+           fClusterPositionAcceptedEnergy->Fill(p2.Phi(), p2.PseudoRapidity(),GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE);
+           fClusterEnergy->Fill(GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE);
+           fClusterEnergyCentNotMatched->Fill(GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE,cent);
            fHistClusterSizeVsCent->Fill(cluster->GetNCells(),cent);
-           fClusterEt->Fill(TMath::Sin(p2.Theta())*cluster->E());
-           uncorrEt += TMath::Sin(p2.Theta())*cluster->E();
-           float myuncorrEt = TMath::Sin(p2.Theta())*cluster->E();
+           fClusterEt->Fill(TMath::Sin(p2.Theta())*GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE);
+           uncorrEt += TMath::Sin(p2.Theta())*GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE;
+           float myuncorrEt = TMath::Sin(p2.Theta())*GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE;
            fTotRawEt += myuncorrEt;
 
-           Double_t effCorrEt = CorrectForReconstructionEfficiency(*cluster,cent);
+           Double_t effCorrEt = CorrectForReconstructionEfficiency(*cluster,fReconstructedE,cent)*GetCorrectionModification(*cluster,0,0,cent);
            rawSignal += myuncorrEt;
            effCorrSignal +=effCorrEt;
-           //cout<<"cluster energy "<<cluster->E()<<" eff corr Et "<<effCorrEt<<endl;
+           //cout<<"cluster energy "<<fReconstructedE<<" eff corr Et "<<effCorrEt<<endl;
            fTotRawEtEffCorr += effCorrEt;
            fTotNeutralEt += effCorrEt;
            nominalRawEt += effCorrEt;
@@ -555,7 +576,7 @@ Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
     fHistCentVsNchVsNclReco->Fill(cent,multiplicity,nCluster);
     fHistPiKPTrackMatchedDepositedVsNch->Fill(etPiKPMatched,multiplicity);
     fHistPiKPTrackMatchedDepositedVsNchNoEff->Fill(etPiKPMatchedNoEff,multiplicity);
-    delete pID;
+    //delete pID;
     return 0;
 }
 
@@ -641,6 +662,7 @@ void AliAnalysisEtReconstructed::FillOutputList(TList* list)
     list->Add(fHistMatchedTracksEvspTvsCent);
     list->Add(fHistMatchedTracksEvspTvsCentEffCorr);
     list->Add(fHistMatchedTracksEvspTvsCentEffTMCorr);
+    list->Add(fHistPeripheralMatchedTracksEvspTvsCentEffTMCorr);
     list->Add(fHistMatchedTracksEvspTvsCentEffTMCorr500MeV);
     list->Add(fHistFoundHadronsvsCent);
     list->Add(fHistNotFoundHadronsvsCent);
@@ -679,6 +701,7 @@ void AliAnalysisEtReconstructed::FillOutputList(TList* list)
     list->Add(fHistCentVsNchVsNclReco);
     list->Add(fHistRawSignalReco);
     list->Add(fHistEffCorrSignalReco);
+    list->Add(fHistRecoRCorrVsPtVsCent);
 }
 
 void AliAnalysisEtReconstructed::CreateHistograms()
@@ -787,6 +810,7 @@ void AliAnalysisEtReconstructed::CreateHistograms()
     fHistMatchedTracksEvspTvsCent = new TH3F("fHistMatchedTracksEvspTvsCent", "fHistMatchedTracksEvspTvsCent",100, 0, 3,100,0,3,20,-0.5,19.5);
     fHistMatchedTracksEvspTvsCentEffCorr = new TH3F("fHistMatchedTracksEvspTvsCentEffCorr", "fHistMatchedTracksEvspTvsCentEffCorr",100, 0, 3,100,0,3,20,-0.5,19.5);
     fHistMatchedTracksEvspTvsCentEffTMCorr = new TH3F("fHistMatchedTracksEvspTvsCentEffTMCorr", "fHistMatchedTracksEvspTvsCentEffTMCorr",100, 0, 3,100,0,3,20,-0.5,19.5);
+    fHistPeripheralMatchedTracksEvspTvsCentEffTMCorr = new TH3F("fHistPeripheralMatchedTracksEvspTvsCentEffTMCorr", "fHistPeripheralMatchedTracksEvspTvsCentEffTMCorr",100, 0, 3,100,0,3,20,-0.5,19.5);
     fHistMatchedTracksEvspTvsCentEffTMCorr500MeV = new TH3F("fHistMatchedTracksEvspTvsCentEffTMCorr500MeV", "fHistMatchedTracksEvspTvsCentEffTMCorr500MeV",100, 0, 3,100,0,3,20,-0.5,19.5);
 
     float max = 200;
@@ -850,6 +874,7 @@ void AliAnalysisEtReconstructed::CreateHistograms()
 
    fHistRawSignalReco = new TH1F("fHistRawSignalReco","fHistRawSignalReco",20,-0.5,19.5);
    fHistEffCorrSignalReco = new TH1F("fHistEffCorrSignalReco","fHistEffCorrSignalReco",20,-0.5,19.5);
+   fHistRecoRCorrVsPtVsCent = new TH3F("fHistRecoRCorrVsPtVsCent","fHistRecoRCorrVsPtVsCent",72,0,2,50,0,10,20,-0.5,19.5);
 
 }
 Double_t AliAnalysisEtReconstructed::ApplyModifiedCorrections(const AliESDCaloCluster& cluster,Int_t nonLinCorr, Int_t effCorr, Int_t cent)