T0/AliT0HIanalysisTask.cxx

   1 #include "TChain.h"
   2 #include "TTree.h"
   3 #include "TList.h"
   4 #include "TH1F.h"
   5 #include "TH2F.h"
   6 #include "TCanvas.h"
   7 #include "TMath.h"
   8 #include "TObjString.h"
   9
  10 #include "AliAnalysisTask.h"
  11 #include "AliPhysicsSelection.h"
  12 #include "AliAnalysisManager.h"
  13
  14 #include "AliESDEvent.h"
  15 #include "AliESDfriend.h"
  16 #include "AliESDInputHandler.h"
  17 #include "AliESDTZEROfriend.h"
  18 #include "AliT0HIanalysisTask.h"
  19 #include "AliESDVZERO.h"
  20 #include "AliMultiplicity.h"
  21 #include "AliTriggerAnalysis.h"
  22 #include "AliESDpid.h"
  23 #include "AliESDtrackCuts.h"
  24 #include "AliCentrality.h"
  25 #include "AliVEvent.h"
  26
  27 ClassImp(AliT0HIanalysisTask)
  28
  29 //________________________________________________________________________
  30 AliT0HIanalysisTask::AliT0HIanalysisTask(const char *name)
  31   : AliAnalysisTaskSE(name), fESD(0), fOutputList(0),
  32   fT0OutTree(0x0)
  33                                         //,    fESDpid(new AliESDpid())
  34 {
  35   // Constructor
  36
  37   // Define input and output slots here
  38   // Input slot #0 works with a TChain
  39   DefineInput(0, TChain::Class());
  40   // Output slot #0 id reserved by the base class for AOD
  41   // Output slot #1 writes into a TH1 container
  42   DefineOutput(1, TList::Class());
  43 }
  44
  45 //________________________________________________________________________
  46 AliT0HIanalysisTask::~AliT0HIanalysisTask()
  47 {
  48 // Destructor
  49    if (fT0OutTree) delete fT0OutTree ;
  50 }
  51
  52 //_____________________________________________________________________________
  53 Bool_t AliT0HIanalysisTask::UserNotify()
  54 {
  55 //
  56 // Calls the mother class Notify()
  57 //
  58
  59 //  AliDebug(AliLog::kDebug,"<-");
  60   // AliDebug(AliLog::kDebug,"->");
  61
  62   return AliAnalysisTaskSE::UserNotify();
  63 }
  64
  65 //________________________________________________________________________
  66 void AliT0HIanalysisTask::UserCreateOutputObjects()
  67 {
  68   // Create histograms
  69   // Called once
  70    fOutputList = new TList();
  71   fOutputList->SetOwner(); // Will delete the histos on cleanup
  72
  73   fT0OutTree=new TTree("t0tree","None here");
  74   fT0OutTree->Branch("fNevent", &fEvent);
  75   fT0OutTree->Branch("fOrbit", &fOrbit);
  76   fT0OutTree->Branch("fBC", &fBC);
  77   fT0OutTree->Branch("fNcont", &fNcont, "fNcont/I");
  78   fT0OutTree->Branch("fNcontTPC", &fNcontTPC);
  79   fT0OutTree->Branch("fVertexTPC", &fVertexTPC);
  80   fT0OutTree->Branch("fVertexPrim", &fVertexPrim);
  81   fT0OutTree->Branch("vertexT0", &fVertex, "fVertex/F");
  82   fT0OutTree->Branch("vertexSPD", &fVertexSPD, "vertexSPD/F");
  83   fT0OutTree->Branch("meanAC", &fMeanAC,"meanAC/F");
  84   fT0OutTree->Branch("meanA", &fMeanA,"meanA/F");
  85   fT0OutTree->Branch("meanC", &fMeanC,"meanC/F");
  86   fT0OutTree->Branch("trackletSPD", &fTrackletSPD,"trackletSPD/I");
  87   fT0OutTree->Branch("clustersSPD", &fClustersSPD,"clustersSPD/I");
  88   fT0OutTree->Branch("TOFtracks", &fTOFtracks);
  89   fT0OutTree->Branch("ESDtracks", &fESDtracks);
  90   fT0OutTree->Branch("multV0A", &fMultV0A,"multV0A/F");
  91   fT0OutTree->Branch("multV0C", &fMultV0C,"multV0C/F");
  92   fT0OutTree->Branch("timeV0A", &fTimeV0A,"timeV0A/F");
  93   fT0OutTree->Branch("timeV0C", &fTimeV0C,"timeV0C/F");
  94   fT0OutTree->Branch("sumampA", &fSumampA);
  95   fT0OutTree->Branch("sumampC", &fSumampC);
  96   fT0OutTree->Branch("pileup", &fpileup);
  97   fT0OutTree->Branch("trigger", &fTrigger);
  98   fT0OutTree->Branch("triggerinput", &fTriggerinput);
  99   fT0OutTree->Branch("T0trigger", &fT0Trigger);
 100   fT0OutTree->Branch("t0tofTrack",  &t0tofTrack);
 101   fT0OutTree->Branch("ZDCcut", &fZDCcut);
 102   fT0OutTree->Branch("ftimestamp", &ftimestamp);
 103   fT0OutTree->Branch("centralityV0M", &fcentralityV0M);
 104   fT0OutTree->Branch("centralityZDC", &fcentralityZDC);
 105   fT0OutTree->Branch("centralityTRK", &fcentralityTRK);
 106   fT0OutTree->Branch("centralityCLA", &fcentralityCLA);
 107   for ( Int_t i=0; i<24; i++) {
 108     fT0OutTree->Branch(Form("amp%i", i+1), &famp[i]);
 109     fT0OutTree->Branch(Form("time%i", i+1), &ftime[i]);
 110   }
 111   for ( Int_t i=0; i<5; i++) {
 112     fT0OutTree->Branch(Form("fOrA%i", i+1), &fOrA[i]);
 113     fT0OutTree->Branch(Form("fOrC%i", i+1), &fOrC[i]);
 114     fT0OutTree->Branch(Form("fTVDC%i", i+1), &fTVDC[i]);
 115     for ( Int_t ii=0; ii<24; ii++)
 116       fT0OutTree->Branch(Form("fRawTime%i_%i", ii+1, i+1), &fRawTime[ii][i]);
 117   }
 118   TString pilename[3] = {"T0pileup", "T0background", "T0satellite"};
 119   for ( Int_t i=0; i<3; i++)
 120     fT0OutTree->Branch(pilename[i].Data(), &fT0pileup[i]);
 121   fT0OutTree->Branch("meanACbest", &fMeanACcalc,"meanACcalc/F");
 122   fT0OutTree->Branch("meanAbest", &fMeanAcalc,"meanAcalc/F");
 123   fT0OutTree->Branch("meanCbest", &fMeanCcalc,"meanCcalc/F");
 124   fT0OutTree->Branch("multEstimator", &fMultiplicity);
 125
 126
 127   fOutputList->Add(fT0OutTree);
 128
 129    PostData(1, fOutputList);
 130
 131 }
 132
 133 //________________________________________________________________________
 134 void AliT0HIanalysisTask::UserExec(Option_t *)
 135 {
 136   // Main loop
 137   // Called for each event
 138   /*
 139   AliAnalysisManager* anMan = AliAnalysisManager::GetAnalysisManager();
 140   AliESDInputHandler *handler =
 141     (AliESDInputHandler*)anMan->GetInputEventHandler();
 142   UInt_t selFlag = handler->IsEventSelected();
 143   Bool_t fIsSelected = selFlag & AliVEvent::kCINT7 ;
 144   cout<< selFlag << " SelFlag "
 145       << (selFlag&AliVEvent::kCINT5) <<" CINT5 "
 146       << (selFlag&AliVEvent::kINT7)  << " kCINT7 "
 147       << (selFlag&AliVEvent::kINT8)  << " kINT8"<<endl;
 148   */
 149   // printf (" @@@@@@@@@@@@@@@AliT0HIanalysisTask::UserExec() \n");
 150    fVertex= fVertexSPD =  fMeanAC = fMeanA =  fMeanC =  fTrackletSPD=-99999;
 151   fMeanACcalc = fMeanAcalc =  fMeanCcalc = fNcont = fNcontTPC = -99999;
 152   fMultV0A=fMultV0C=fTimeV0A=fTimeV0C=fMultiplicity = -99999;
 153
 154   for (Int_t i=0; i<24; i++) {
 155     famp[i] = ftime[i] = -9999;
 156     for ( Int_t i0=0; i0<5; i0++) {
 157       fRawTime[i][i0] = -9999;
 158       if(i0==0) {
 159         fTVDC[i0]=-9999;
 160         fOrC[i0]=-9999;
 161         fOrA[i0]=-9999;
 162       }
 163     }
 164   }
 165   Float_t orA , orC, tvdc;
 166   orA = orC = tvdc = -9999;
 167   fBC = -9999;
 168   // Post output data.
 169   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 170   if (!fESD) {
 171     printf("ERROR: fESD not available\n");
 172     return;
 173   }
 174
 175   if(fESD)
 176     {
 177       //    if(fIsSelected)
 178         {
 179           //  AliESDHeader* header = fESD-> GetHeader();
 180         AliESDVZERO* esdV0 = fESD->GetVZEROData();
 181         const AliMultiplicity *mult = fESD->GetMultiplicity();
 182         // printf(" VZERO MULT \n");
 183         AliESDTZERO* tz= (AliESDTZERO*) fESD->GetESDTZERO();
 184         //   printf(" TZERO \n");
 185         const Double32_t *amp, *time, *mean ;
 186         const Double32_t *meanbest;
 187         fVertex= fVertexSPD =  fMeanAC = fMeanA =  fMeanC =  fTrackletSPD=-99999;
 188         fMeanACcalc = fMeanAcalc =  fMeanCcalc = fNcont = fNcontTPC -99999;
 189         fMultV0A=fMultV0C=fTimeV0A=fTimeV0C=fMultiplicity = -99999;
 190
 191         fSumampA=0;
 192         fSumampC=0;
 193         fNcont=0;
 194           Float_t shift=0;
 195
 196           TString triggers = fESD->GetFiredTriggerClasses();
 197           fTrigger.SetString(triggers.Data());
 198           //      if ( !triggers.Contains("CINT7-B") )
 199           //  cout<<triggers<<endl;
 200           TString inputtriggers =fESD-> GetHeader()->GetFiredTriggerInputs();
 201           fTriggerinput.SetString(inputtriggers.Data());
 202           //      cout<<inputtriggers<<endl;
 203
 204           fEvent=fESD->GetEventNumberInFile();
 205           fT0Trigger=fESD->GetT0Trig();
 206
 207           fESDtracks=fESD->GetNumberOfTracks();
 208           AliTriggerAnalysis *trigAna = new AliTriggerAnalysis;
 209           ftimestamp=fESD->GetTimeStamp(); // - 1301217882;
 210           //      printf (" timestamp %d \n", ftimestamp);
 211
 212           //  fMultiplicity =  AliESDtrackCuts::GetReferenceMultiplicity(fESD);
 213           // cout<<" fMultiplicity "<<fMultiplicity<<endl;
 214           fOrbit=fESD->GetOrbitNumber();
 215           fBC=fESD->GetBunchCrossNumber();
 216           fNcont = fESD->GetPrimaryVertexSPD()->GetNContributors();
 217           fNcontTPC = fESD->GetPrimaryVertexTPC()->GetNContributors();
 218           const AliESDVertex *vtxESD = fESD->GetPrimaryVertex();
 219           TString vtxTyp = vtxESD->GetTitle();
 220           //      printf(" %s \n", vtxTyp.Data() );
 221           Bool_t fVtxOK = kFALSE;
 222           if ( !vtxTyp.Contains("vertexer: Z") ||
 223                (vtxESD->GetDispersion()<0.04 && vtxESD->GetZRes()<0.25)) {
 224             fVtxOK = kTRUE;
 225             fVertexPrim = vtxESD->GetZ();
 226             //   printf(" GetPrimaryVertex()  %f \n", fVertexPrim);
 227           }
 228
 229
 230           if(fNcontTPC>=2 && fVtxOK )
 231             fVertexTPC = fESD->GetPrimaryVertexTPC() ->GetZ();
 232           if(fNcont>=0 && fVtxOK) {
 233             fVertexSPD = fESD->GetPrimaryVertex()->GetZ();
 234           }
 235           //      printf(" SPD Vertex()  %f  TPC %f \n", fVertexSPD, fVertexTPC);
 236           fTrackletSPD = mult->GetNumberOfTracklets();
 237           fClustersSPD = mult->GetNumberOfITSClusters(0);
 238           fMultV0A =  esdV0->GetMTotV0A();
 239           fMultV0C =  esdV0->GetMTotV0C();
 240           //TOF hit
 241           Int_t ntracksMatchedToTOF = 0;
 242           Int_t ntracks = fESD->GetNumberOfTracks();
 243           for(Int_t itrk=0;itrk<ntracks;itrk++){
 244             AliESDtrack* track = fESD->GetTrack(itrk);
 245             if (!track) {
 246               Printf("ERROR: Could not receive track %d", itrk);
 247               continue;
 248             }
 249             //no track selection just TOF hit
 250             if (track->IsOn(AliESDtrack::kTOFout)) ntracksMatchedToTOF++;
 251           }
 252           fTOFtracks = ntracksMatchedToTOF;
 253              if(fESDpid){ //get T0_TOF
 254              fESDpid->SetTOFResponse(fESD,AliESDpid::kTOF_T0);
 255              t0tofTrack =(Float_t) (fESDpid->GetTOFResponse().GetStartTime(10.0)); //Get start time from all tracks
 256              }
 257
 258
 259           fpileup = fESD->IsPileupFromSPD();
 260           fTimeV0C = esdV0->GetV0CTime();
 261           fTimeV0A = esdV0->GetV0ATime();
 262           //      printf(" Vo info \n");
 263
 264
 265           mean = fESD->GetT0TOF();
 266           fVertex=fESD->GetT0zVertex();
 267           fMeanA = mean[1];
 268           fMeanC = mean[2];
 269           if (TMath::Abs(fMeanA)<9999 && TMath::Abs(fMeanC)<9999)
 270             fMeanAC = mean[0] ;
 271           fSumampC = tz->GetMultC();
 272           fSumampA = tz->GetMultA();
 273           //    printf("!!!!!!!!fSumampC %f fSumampA %f \n",fSumampC, fSumampA);
 274
 275           amp=fESD->GetT0amplitude();
 276           time=fESD->GetT0time();
 277           for (Int_t i=0; i<24; i++){
 278             if( time[i]>100 && amp[i]>0.1){
 279           //  if(i<12) fSumampC += Int_t (amp[i]+0.5);
 280           // if(i>11) fSumampA += Int_t (amp[i]+0.5);
 281               famp[i] = amp[i];
 282               ftime[i] = time[i];
 283             }
 284       }
 285
 286           //new OrA OrC TVDC
 287       // AliESDTZERO* tz= (AliESDTZERO*) fESD->GetESDTZERO();
 288       meanbest = tz->GetT0TOFbest();
 289       fMeanAcalc = meanbest[1];
 290       fMeanCcalc = meanbest[2];
 291       if (TMath::Abs(fMeanAcalc)<9999 && TMath::Abs(fMeanCcalc)<9999)
 292         fMeanACcalc = meanbest[0] ;
 293       for (Int_t ii=0; ii<5; ii++){
 294         orA   = tz->GetOrA(ii);
 295         orC   = tz->GetOrC(ii);
 296         tvdc  = tz->GetTVDC(ii);
 297         //      printf("  new signasl %d  %f %f %f \n",ii, orA, orC, tvdc);
 298         if(ii==0) {
 299           fOrA[ii] = orA;
 300           fOrC[ii] = orC;
 301           fTVDC[ii] = tvdc;
 302           //    printf("  new signasl %d  %f %f %f \n",ii, orA, orC, tvdc);
 303           //  cout<<tvdc<<" "<<fT0Trigger<<endl;
 304           if (tvdc>-5 && tvdc<5 && tvdc!=0) cout<<tvdc<<" "<<fT0Trigger<<endl;
 305         } else {
 306           if ( ii>0 && fOrA[ii-1] !=orA)        fOrA[ii] = orA;
 307           else fOrA[ii]=-9999;
 308           if ( ii>0 && fOrC[ii-1] !=orC)        fOrC[ii] = orC;
 309           else fOrC[ii]=-9999;
 310           if ( ii>0 && fTVDC[ii-1] != tvdc)     fTVDC[ii] = tvdc;
 311           else fTVDC[ii]=-9999;
 312         }
 313         for (Int_t iii=0; iii<24; iii++)
 314           fRawTime[iii][ii] =  tz->GetTimeFull(iii,ii);
 315       }
 316       fT0pileup[0] = tz->GetPileupFlag();
 317       fT0pileup[1] = tz->GetBackgroundFlag();
 318       fT0pileup[2] = tz->GetSatellite();
 319
 320        AliCentrality *centrality = fESD->GetCentrality();
 321
 322       //The centrality function you can use
 323
 324       fcentralityV0M = centrality->GetCentralityPercentile("V0M"); // returns the centrality
 325       fcentralityTRK = centrality->GetCentralityPercentile("TRK"); // returns the centrality
 326       fcentralityZDC = centrality->GetCentralityPercentile("ZEMvsZDC"); // returns the
 327 //      fcentralityTRK = centrality->GetCentralityClass10("TRK"); // returns centrality class
 328
 329   //    fcentralityZDC = centrality->GetCentralityClass10("ZEMvsZDC"); // returns centrality
 330
 331         } //selected
 332     } //ESD
 333
 334   fT0OutTree->Fill();
 335
 336   PostData(1, fOutputList);
 337
 338 }
 339
 340 //________________________________________________________________________
 341 void AliT0HIanalysisTask::Terminate(Option_t *)
 342 {
 343   // Draw result to the screen
 344   // Called once at the end of the query
 345
 346   /* fOutputList = dynamic_cast<TList*> (GetOutputData(1));
 347   if (!fOutputList) {
 348     printf("ERROR: Output list not available\n");
 349     return;
 350   }
 351     */
 352
 353 }