PWG1/ITS/AliAnalysisTaskSDDRP.cxx

   1 #include "AliAnalysisTask.h"
   2 #include "AliAnalysisManager.h"
   3 #include "AliAnalysisDataContainer.h"
   4 #include "AliITSRecPoint.h"
   5 #include "AliESDEvent.h"
   6 #include "AliTrackPointArray.h"
   7 #include "AliITSgeomTGeo.h"
   8 #include "AliESDfriend.h"
   9 #include "AliCDBManager.h"
  10 #include "AliCDBEntry.h"
  11 #include "AliITSCalibrationSDD.h"
  12 #include "AliITSresponseSDD.h"
  13 #include "AliGeomManager.h"
  14 #include <TSystem.h>
  15 #include <TTree.h>
  16 #include <TH1F.h>
  17 #include <TChain.h>
  18 #include <TGeoGlobalMagField.h>
  19 #include "AliESDInputHandlerRP.h"
  20 #include "AliAnalysisTaskSDDRP.h"
  21
  22 ClassImp(AliAnalysisTaskSDDRP)
  23 //______________________________________________________________________________
  24 AliAnalysisTaskSDDRP::AliAnalysisTaskSDDRP() : AliAnalysisTaskSE("SDD RecPoints"),
  25   fOutput(0),
  26   fHistNEvents(0),
  27   fRecPMod(0),
  28   fTrackPMod(0),
  29   fGoodAnMod(0),
  30   fRecPLadLay3(0),
  31   fRecPLadLay4(0),
  32   fTrackPLadLay3(0),
  33   fTrackPLadLay4(0),
  34   fGoodAnLadLay3(0),
  35   fGoodAnLadLay4(0),
  36   fDriftTimeRP(0),
  37   fDriftTimeTP(0),
  38   fESD(0),
  39   fESDfriend(0),
  40   fResp(0),
  41   fRunNumber(0),
  42   fMinITSpts(3),
  43   fMinPfordEdx(1.5),
  44   fOnlyCINT1BTrig(0),
  45   fInitialised(0)
  46 {
  47   //
  48   DefineOutput(1, TList::Class());
  49 }
  50
  51
  52 //___________________________________________________________________________
  53 AliAnalysisTaskSDDRP::~AliAnalysisTaskSDDRP(){
  54   //
  55   if (fOutput) {
  56     delete fOutput;
  57     fOutput = 0;
  58   }
  59 }
  60
  61
  62 //___________________________________________________________________________
  63
  64 void AliAnalysisTaskSDDRP::UserCreateOutputObjects() {
  65
  66   fOutput = new TList();
  67   fOutput->SetOwner();
  68   fOutput->SetName("OutputHistos");
  69
  70   fHistNEvents = new TH1F("hNEvents", "Number of processed events",3,-1.5,1.5);
  71   fHistNEvents->Sumw2();
  72   fHistNEvents->SetMinimum(0);
  73   fOutput->Add(fHistNEvents);
  74
  75   // -- Module histos
  76
  77   fRecPMod = new TH1F("hRPMod","Rec Points per Module",260,239.5,499.5);
  78   fRecPMod->Sumw2();
  79   fRecPMod->SetMinimum(0);
  80   fOutput->Add(fRecPMod);
  81
  82   fTrackPMod = new TH1F("hTPMod","Track Points per Module",260,239.5,499.5);
  83   fTrackPMod->Sumw2();
  84   fTrackPMod->SetMinimum(0);
  85   fOutput->Add(fTrackPMod);
  86
  87   fGoodAnMod = new TH1F("hGAMod","Good Anodes per Module",260,239.5,499.5);
  88   fOutput->Add(fGoodAnMod);
  89
  90   // -- Ladder histos
  91
  92   fRecPLadLay3 = new TH1F("hRPLad3","Rec Points per Ladder Layer 3",14,-0.5,13.5);
  93   fRecPLadLay3->Sumw2();
  94   fRecPLadLay3->SetMinimum(0);
  95   fOutput->Add(fRecPLadLay3);
  96
  97   fRecPLadLay4 = new TH1F("hRPLad4","Rec Points per Ladder Layer 4",22,-0.5,21.5);
  98   fRecPLadLay4->Sumw2();
  99   fRecPLadLay4->SetMinimum(0);
 100   fOutput->Add(fRecPLadLay4);
 101
 102   fTrackPLadLay3 = new TH1F("hTPLad3","Track Points per Ladder Layer 3",14,-0.5,13.5);
 103   fTrackPLadLay3->Sumw2();
 104   fTrackPLadLay3->SetMinimum(0);
 105   fOutput->Add(fTrackPLadLay3);
 106
 107   fTrackPLadLay4 = new TH1F("hTPLad4","Track Points per Ladder Layer 4",22,-0.5,21.5);
 108   fTrackPLadLay4->Sumw2();
 109   fTrackPLadLay4->SetMinimum(0);
 110   fOutput->Add(fTrackPLadLay4);
 111
 112   fGoodAnLadLay3 = new TH1F("hGALad3","Good Anodes per Ladder Layer 3",14,-0.5,13.5);
 113   fOutput->Add(fGoodAnLadLay3);
 114
 115   fGoodAnLadLay4 = new TH1F("hGALad4","Good Anodes per Ladder Layer 4",22,-0.5,21.5);
 116   fOutput->Add(fGoodAnLadLay4);
 117
 118   fDriftTimeRP=new TH1F("hDrTimRP","Drift Time from Rec Points (ns)",100,0.,6400.);
 119   fDriftTimeRP->Sumw2();
 120   fDriftTimeRP->SetMinimum(0.);
 121   fOutput->Add(fDriftTimeRP);
 122
 123   fDriftTimeTP=new TH1F("hDrTimTP","Drift Time from Track Points (ns)",100,0.,6400.);
 124   fDriftTimeTP->Sumw2();
 125   fDriftTimeTP->SetMinimum(0.);
 126   fOutput->Add(fDriftTimeTP);
 127
 128   for(Int_t it=0; it<8; it++){
 129     fSignalTime[it]=new TH1F(Form("hSigTimeInt%d",it),Form("hSigTimeInt%d",it),100,0.,300.);
 130     fSignalTime[it]->Sumw2();
 131     fSignalTime[it]->SetMinimum(0);
 132     fOutput->Add(fSignalTime[it]);
 133   }
 134 }
 135 //______________________________________________________________________________
 136 void AliAnalysisTaskSDDRP::UserExec(Option_t *)
 137 {
 138   //
 139   fESD = (AliESDEvent*) (InputEvent());
 140
 141   if(!fESD) {
 142     printf("AliAnalysisTaskSDDRP::Exec(): bad ESD\n");
 143     return;
 144   }
 145
 146   fESDfriend = static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
 147
 148
 149
 150   if(!fESDfriend) {
 151     printf("AliAnalysisTaskSDDRP::Exec(): bad ESDfriend\n");
 152     return;
 153   }
 154   PostData(1, fOutput);
 155   fESD->SetESDfriend(fESDfriend);
 156
 157   if (!fInitialised) {
 158     fInitialised = 1;
 159
 160     AliCDBManager* man = AliCDBManager::Instance();
 161     man->SetDefaultStorage("raw://");
 162     man->SetRun(fESD->GetRunNumber());
 163
 164
 165     AliCDBEntry* eR=(AliCDBEntry*)man->Get("ITS/Calib/RespSDD");
 166     fResp=(AliITSresponseSDD*)eR->GetObject();
 167
 168     AliCDBEntry* eC=(AliCDBEntry*)man->Get("ITS/Calib/CalibSDD");
 169     TObjArray* calsdd=(TObjArray*)eC->GetObject();
 170     Int_t countGood3[14];
 171     Int_t countGood4[22];
 172     Int_t countGoodMod[260];
 173     for(Int_t ilad=0;ilad<14;ilad++) countGood3[ilad]=0;
 174     for(Int_t ilad=0;ilad<22;ilad++) countGood4[ilad]=0;
 175     for(Int_t imod=0;imod<260;imod++) countGoodMod[imod]=0;
 176     for(Int_t imod=0;imod<260;imod++){
 177       AliITSCalibrationSDD* cal=(AliITSCalibrationSDD*)calsdd->At(imod);
 178       if(cal->IsBad()) continue;
 179       Int_t modId=imod+AliITSgeomTGeo::GetModuleIndex(3,1,1);
 180       Int_t lay,lad,det;
 181       AliITSgeomTGeo::GetModuleId(modId,lay,lad,det);
 182       if(!CheckModule(lay,lad,det)) continue;
 183       for(Int_t ian=0; ian<512; ian++){
 184         if(cal->IsBadChannel(ian)) continue;
 185         countGoodMod[imod]++;
 186         if(lay==3) countGood3[lad-1]++;
 187         else if(lay==4) countGood4[lad-1]++;
 188       }
 189     }
 190     for(Int_t imod=0;imod<260;imod++) fGoodAnMod->SetBinContent(imod+1,countGoodMod[imod]);
 191     fGoodAnMod->SetMinimum(0);
 192     for(Int_t ilad=0;ilad<14;ilad++) fGoodAnLadLay3->SetBinContent(ilad+1,countGood3[ilad]);
 193     fGoodAnLadLay3->SetMinimum(0);
 194     for(Int_t ilad=0;ilad<22;ilad++) fGoodAnLadLay4->SetBinContent(ilad+1,countGood4[ilad]);
 195     fGoodAnLadLay4->SetMinimum(0);
 196   }
 197
 198
 199   fHistNEvents->Fill(0);
 200   if(fOnlyCINT1BTrig){
 201     if(!fESD->IsTriggerClassFired("CINT1B-ABCE-NOPF-ALL")) return;
 202     fHistNEvents->Fill(1);
 203   }
 204   const AliTrackPointArray *array = 0;
 205   Int_t ntracks = fESD->GetNumberOfTracks();
 206   for (Int_t itrack=0; itrack < ntracks; itrack++) {
 207     AliESDtrack * track = fESD->GetTrack(itrack);
 208     if (!track) continue;
 209     if(track->GetNcls(1)>0) continue;
 210     if(track->GetNcls(0) < fMinITSpts) continue;
 211     Double_t dedx[4];
 212     track->GetITSdEdxSamples(dedx);
 213     array = track->GetTrackPointArray();
 214     if(!array) continue;
 215     for(Int_t ipt=0; ipt<array->GetNPoints(); ipt++) {
 216       AliTrackPoint point;
 217       Int_t modId;
 218       array->GetPoint(point,ipt);
 219       Int_t volId = point.GetVolumeID();
 220       Int_t layerId = AliGeomManager::VolUIDToLayer(volId,modId);
 221       if(layerId<3 || layerId>4) continue;
 222       modId+=AliITSgeomTGeo::GetModuleIndex(layerId,1,1);
 223       Int_t lay,lad,det;
 224       AliITSgeomTGeo::GetModuleId(modId,lay,lad,det);
 225       if(!CheckModule(lay,lad,det)) continue;
 226       fTrackPMod->Fill(modId);
 227       fDriftTimeTP->Fill(point.GetDriftTime());
 228       Float_t dtime=point.GetDriftTime()-fResp->GetTimeZero(modId);
 229       Int_t theBin=int(dtime/6500.*8.);
 230       if(layerId==3){
 231         fTrackPLadLay3->Fill(lad-1);
 232         if(dedx[0]>0. && track->P()>fMinPfordEdx) fSignalTime[theBin]->Fill(dedx[0]);
 233       }
 234       if(layerId==4){
 235         fTrackPLadLay4->Fill(lad-1);
 236         if(dedx[0]>0.&& track->P()>fMinPfordEdx) fSignalTime[theBin]->Fill(dedx[1]);
 237       }
 238     }
 239   }
 240
 241   AliESDInputHandlerRP *hand = dynamic_cast<AliESDInputHandlerRP*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 242   TTree* tR = hand->GetTreeR("ITS");
 243   TClonesArray *ITSrec= new TClonesArray("AliITSRecPoint");
 244   TBranch *branch =tR->GetBranch("ITSRecPoints");
 245   branch->SetAddress(&ITSrec);
 246   for (Int_t modId=240; modId<500; modId++){
 247     Int_t lay,lad,det;
 248     AliITSgeomTGeo::GetModuleId(modId,lay,lad,det);
 249     if(!CheckModule(lay,lad,det)) continue;
 250     branch->GetEvent(modId);
 251     Int_t nrecp = ITSrec->GetEntries();
 252     fRecPMod->Fill(modId,nrecp);
 253     if(lay==3) fRecPLadLay3->Fill(lad-1,nrecp);
 254     if(lay==4) fRecPLadLay4->Fill(lad-1,nrecp);
 255     for(Int_t irec=0;irec<nrecp;irec++) {
 256       AliITSRecPoint *recp = (AliITSRecPoint*)ITSrec->At(irec);
 257       fDriftTimeRP->Fill(recp->GetDriftTime());
 258     }
 259   }
 260   ITSrec->Delete();
 261   delete ITSrec;
 262
 263   PostData(1,fOutput);
 264
 265 }
 266 //______________________________________________________________________________
 267 Bool_t AliAnalysisTaskSDDRP::CheckModule(Int_t lay, Int_t lad, Int_t det) const{
 268   //
 269   if(lay==4){
 270     if(lad==3 && det==5) return kFALSE; // 1500 V
 271     if(lad==3 && det==6) return kFALSE; // 0 V
 272     if(lad==3 && det==7) return kFALSE; // 1500 V
 273     if(lad==4 && det==1) return kFALSE; // 0 V
 274     if(lad==4 && det==2) return kFALSE; // 1500 V
 275     if(lad==7 && det==5) return kFALSE; // 0 MV
 276     if(lad==9 && det==3) return kFALSE; // 1500 V
 277     if(lad==9 && det==4) return kFALSE; // 0 V
 278     if(lad==9 && det==5) return kFALSE; // 1500 V
 279     if(lad==11 && det==6) return kFALSE; // 1500 V
 280     if(lad==11 && det==7) return kFALSE; // 0 V
 281     if(lad==11 && det==8) return kFALSE; // 1500 V
 282     if(lad==18 && det==5) return kFALSE; // 1500 V
 283     if(lad==18 && det==6) return kFALSE; // 0 V
 284     if(lad==18 && det==7) return kFALSE; // 1500 V
 285     if(lad==22 && det==1) return kFALSE; // 0 V
 286     if(lad==22 && det==2) return kFALSE; // 1500 V
 287   }
 288   if(lay==3){
 289     if(lad==4 && det==4) return kFALSE; // 1500 V
 290     if(lad==3) return kFALSE;  // swapped in geometry
 291   }
 292   return kTRUE;
 293 }
 294
 295 //______________________________________________________________________________
 296 void AliAnalysisTaskSDDRP::Terminate(Option_t */*option*/)
 297 {
 298   // Terminate analysis
 299   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 300   if (!fOutput) {
 301     printf("ERROR: fOutput not available\n");
 302     return;
 303   }
 304   fHistNEvents= dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
 305   fRecPMod= dynamic_cast<TH1F*>(fOutput->FindObject("hRPMod"));
 306   fTrackPMod= dynamic_cast<TH1F*>(fOutput->FindObject("hTPMod"));
 307   fGoodAnMod= dynamic_cast<TH1F*>(fOutput->FindObject("hGAMod"));
 308
 309   fRecPLadLay3= dynamic_cast<TH1F*>(fOutput->FindObject("hRPLad3"));
 310   fRecPLadLay4= dynamic_cast<TH1F*>(fOutput->FindObject("hRPLad4"));
 311   fTrackPLadLay3= dynamic_cast<TH1F*>(fOutput->FindObject("hTPLad3"));
 312   fTrackPLadLay4= dynamic_cast<TH1F*>(fOutput->FindObject("hTPLad4"));
 313   fGoodAnLadLay3= dynamic_cast<TH1F*>(fOutput->FindObject("hGALad3"));
 314   fGoodAnLadLay4= dynamic_cast<TH1F*>(fOutput->FindObject("hGALad4"));
 315
 316   fDriftTimeRP= dynamic_cast<TH1F*>(fOutput->FindObject("hDrTimRP"));
 317   fDriftTimeTP= dynamic_cast<TH1F*>(fOutput->FindObject("hDrTimTP"));
 318
 319   for(Int_t it=0; it<8; it++){
 320     fSignalTime[it]= dynamic_cast<TH1F*>(fOutput->FindObject(Form("hSigTimeInt%d",it)));
 321   }
 322
 323   return;
 324 }
 325
 326
 327
 328
 329