ITS/ITSSPDVertexDiamondda.cxx

   1 /*
   2 Contact: cvetan.cheshkov@cern.ch
   3 Link: missing
   4 Run Type: PHYSICS
   5 DA Type: MON
   6 Number of events needed: 10000
   7 Input Files:
   8 Output Files:
   9 Trigger types used: PHYSICS
  10 */
  11
  12 #define X_LIMIT 2.0
  13 #define Y_LIMIT 2.0
  14 #define Z_LIMIT 50.0
  15 #define X_DELTA 0.02
  16 #define Y_DELTA 0.02
  17 #define Z_DELTA 0.2
  18 #define OUTPUT_FILE "SPDVertexDiamondDA.root"
  19 #define CDB_STORAGE "local://$ALICE_ROOT"
  20 #define N_EVENTS_AUTOSAVE 50
  21 #define NFITPARAMS 5
  22
  23 extern "C" {
  24 #include "daqDA.h"
  25 }
  26
  27 #include "event.h"
  28 #include "monitor.h"
  29
  30 #ifdef ALI_AMORE
  31 #include <AmoreDA.h>
  32 //int amore::da::Updated(char const*) {}
  33 #endif
  34
  35 #include <TTree.h>
  36 #include <TH1.h>
  37 #include <TH2.h>
  38 #include <TF2.h>
  39 #include <TFile.h>
  40 #include <TPluginManager.h>
  41 #include <TROOT.h>
  42 #include <TFitter.h>
  43
  44 #include "AliRawReaderDate.h"
  45 #include "AliGeomManager.h"
  46 #include "AliCDBManager.h"
  47 #include "AliESDVertex.h"
  48 #include "AliITSDetTypeRec.h"
  49 #include "AliITSInitGeometry.h"
  50 #include "AliITSVertexer3DTapan.h"
  51
  52 TF2 *fitFcn = 0x0;
  53
  54 AliESDVertex* FitVertexDiamond(TH2F *hXY, TH1F *hZ);
  55
  56 int main(int argc, char **argv) {
  57
  58   gROOT->GetPluginManager()->AddHandler("TVirtualStreamerInfo",
  59                                         "*",
  60                                         "TStreamerInfo",
  61                                         "RIO",
  62                                         "TStreamerInfo()");
  63   TFitter *minuitFit = new TFitter(NFITPARAMS);
  64   TVirtualFitter::SetFitter(minuitFit);
  65
  66   int status;
  67   if (argc<2) {
  68     printf("Wrong number of arguments\n");
  69     return -1;
  70   }
  71
  72   /* define data source : this is argument 1 */
  73   status=monitorSetDataSource( argv[1] );
  74   if (status!=0) {
  75     printf("monitorSetDataSource() failed : %s\n",monitorDecodeError(status));
  76     return -1;
  77   }
  78
  79   /* declare monitoring program */
  80   status=monitorDeclareMp( __FILE__ );
  81   if (status!=0) {
  82     printf("monitorDeclareMp() failed : %s\n",monitorDecodeError(status));
  83     return -1;
  84   }
  85
  86   /* define wait event timeout - 1s max */
  87   monitorSetNowait();
  88   monitorSetNoWaitNetworkTimeout(1000);
  89
  90   /* log start of process */
  91   printf("Vertex-Diamond SPD DA started\n");
  92
  93   /* init some counters */
  94   int nevents_with_vertex = 0;
  95   int nevents_physics=0;
  96   int nevents_total=0;
  97
  98   struct eventHeaderStruct *event;
  99   eventTypeType eventT;
 100
 101   // Get run number
 102   if (getenv("DATE_RUN_NUMBER")==0) {
 103     printf("DATE_RUN_NUMBER not properly set.\n");
 104     return -1;
 105   }
 106   int runNr = atoi(getenv("DATE_RUN_NUMBER"));
 107
 108   // Histograms initialization
 109   TH2F *hXY = new TH2F("hXY","Vertex Diamond (Y vs X)",
 110                        2*(Int_t)(X_LIMIT/X_DELTA),-X_LIMIT,X_LIMIT,
 111                        2*(Int_t)(Y_LIMIT/Y_DELTA),-Y_LIMIT,Y_LIMIT);
 112   TH1F *hZ  = new TH1F("hZ"," Longitudinal Vertex Profile",
 113                        2*(Int_t)(Z_LIMIT/Z_DELTA),-Z_LIMIT,Z_LIMIT);
 114
 115   // Global initializations
 116   AliCDBManager *man = AliCDBManager::Instance();
 117   man->SetDefaultStorage(CDB_STORAGE);
 118   man->SetRun(runNr);
 119   AliGeomManager::LoadGeometry("geometry.root");
 120   AliGeomManager::ApplyAlignObjsFromCDB("ITS");
 121
 122   // ITS initializations
 123   AliITSInitGeometry initgeom;
 124   AliITSgeom *geom = initgeom.CreateAliITSgeom();
 125   printf("Geometry name: %s\n",(initgeom.GetGeometryName()).Data());
 126
 127   AliITSDetTypeRec *detTypeRec = new AliITSDetTypeRec();
 128   detTypeRec->SetITSgeom(geom);
 129   detTypeRec->SetDefaults();
 130   detTypeRec->SetDefaultClusterFindersV2(kTRUE);
 131
 132   // Initialization of AMORE sender
 133 #ifdef ALI_AMORE
 134   amore::da::AmoreDA vtxAmore(amore::da::AmoreDA::kSender);
 135 #endif
 136   /* main loop (infinite) */
 137   for(;;) {
 138
 139     /* check shutdown condition */
 140     if (daqDA_checkShutdown()) {break;}
 141
 142     /* get next event (blocking call until timeout) */
 143     status=monitorGetEventDynamic((void **)&event);
 144     if (status==MON_ERR_EOF) {
 145       printf ("End of File detected\n");
 146       break; /* end of monitoring file has been reached */
 147     }
 148
 149     if (status!=0) {
 150       printf("monitorGetEventDynamic() failed : %s\n",monitorDecodeError(status));
 151       break;
 152     }
 153
 154     /* retry if got no event */
 155     if (event==NULL) {
 156       continue;
 157     }
 158
 159     nevents_total++;
 160     eventT=event->eventType;
 161     switch (event->eventType){
 162
 163       /* START OF RUN */
 164     case START_OF_RUN:
 165       break;
 166       /* END START OF RUN */
 167
 168     /* END OF RUN */
 169     case END_OF_RUN:
 170       break;
 171
 172     case PHYSICS_EVENT:
 173       nevents_physics++;
 174       AliRawReader *rawReader = new AliRawReaderDate((void*)event);
 175
 176       // Run SPD cluster finder
 177       TTree* clustersTree = new TTree("TreeR", "Reconstructed Points Container"); //make a tree
 178       detTypeRec->DigitsToRecPoints(rawReader,clustersTree,"SPD");
 179
 180       // Run vertex-finder
 181       AliITSVertexer3DTapan *vertexer = new AliITSVertexer3DTapan(geom,1000);
 182       AliESDVertex *vtx = vertexer->FindVertexForCurrentEvent(clustersTree);
 183
 184       if (TMath::Abs(vtx->GetChi2()) < 0.1) {
 185         // Fill the vertex into the histos
 186         nevents_with_vertex++;
 187         hXY->Fill(vtx->GetXv(),vtx->GetYv());
 188         hZ->Fill(vtx->GetZv());
 189
 190         // Auto save
 191         if ((nevents_with_vertex%N_EVENTS_AUTOSAVE) == 0) {
 192           TFile outFile(OUTPUT_FILE, "update");
 193           AliESDVertex *fitVtx = FitVertexDiamond(hXY,hZ);
 194           if (fitVtx) {
 195             fitVtx->Write(fitVtx->GetName(),TObject::kOverwrite);
 196             TH1 *fithXY = fitFcn->CreateHistogram();
 197             fithXY->Write(fithXY->GetName(),TObject::kOverwrite);
 198             delete fithXY;
 199           }
 200           hXY->Write(hXY->GetName(),TObject::kOverwrite);
 201           hZ->Write(hZ->GetName(),TObject::kOverwrite);
 202           outFile.Close();
 203           delete fitVtx;
 204         }
 205       }
 206
 207       delete vtx;
 208       delete vertexer;
 209
 210       delete clustersTree;
 211       delete rawReader;
 212
 213     }
 214
 215     /* free resources */
 216     free(event);
 217
 218     /* exit when last event received, no need to wait for TERM signal */
 219     if (eventT==END_OF_RUN) {
 220       printf("EOR event detected\n");
 221       break;
 222     }
 223   }
 224
 225   if (detTypeRec) delete detTypeRec;
 226
 227   // Store the final histograms
 228   TFile outFile(OUTPUT_FILE, "update");
 229   if (nevents_with_vertex > N_EVENTS_AUTOSAVE) {
 230     // Fit XY & Z histograms
 231     AliESDVertex *fitVtx = FitVertexDiamond(hXY,hZ);
 232     if (fitVtx) {
 233       fitVtx->Write(fitVtx->GetName(),TObject::kOverwrite);
 234       TH1 *fithXY = fitFcn->CreateHistogram();
 235       fithXY->Write(fithXY->GetName(),TObject::kOverwrite);
 236       delete fithXY;
 237     }
 238     delete fitVtx;
 239   }
 240   hXY->Write(hXY->GetName(),TObject::kOverwrite);
 241   hZ->Write(hZ->GetName(),TObject::kOverwrite);
 242   outFile.Close();
 243
 244 #ifdef ALI_AMORE
 245   // send the histos to AMORE pool
 246   printf("AMORE send status: %d",vtxAmore.Send(hXY->GetName(),hXY));
 247 #endif
 248
 249   delete minuitFit;
 250   TVirtualFitter::SetFitter(0);
 251
 252   /* write report */
 253   printf("Run #%s, received %d events with vertex, out of %d physics and out of %d total events\n",getenv("DATE_RUN_NUMBER"),nevents_with_vertex,nevents_physics,nevents_total);
 254
 255   status=0;
 256
 257   /* export file to FXS */
 258   if (daqDA_FES_storeFile(OUTPUT_FILE, "VertexDiamond")) {
 259     status=-2;
 260   }
 261
 262   return status;
 263 }
 264
 265 Double_t fitFunction(Double_t *x, Double_t *par) {
 266
 267   Double_t t1 =   x[0] - par[1];
 268   Double_t t2 =   x[1] - par[2];
 269
 270   return par[0]*TMath::Exp(-0.5*(t1*t1/(par[3]*par[3])+t2*t2/(par[4]*par[4])));
 271 }
 272
 273 AliESDVertex* FitVertexDiamond(TH2F *hXY, TH1F *hZ)
 274 {
 275
 276   if (!fitFcn) {
 277     fitFcn = new TF2("fitFcn",fitFunction,
 278                      -X_LIMIT,X_LIMIT,
 279                      -Y_LIMIT,Y_LIMIT,NFITPARAMS);
 280     fitFcn->SetNpx(2*(Int_t)(X_LIMIT/X_DELTA));
 281     fitFcn->SetNpy(2*(Int_t)(Y_LIMIT/Y_DELTA));
 282     fitFcn->SetParameters(hXY->GetMaximum(),0,0,hXY->GetRMS(1),hXY->GetRMS(2));
 283     fitFcn->Update();
 284   }
 285
 286   if (hXY->Fit("fitFcn","L0V+") != 0) {
 287     printf("XY fit failed!");
 288     return 0x0;
 289   }
 290
 291   Double_t pos[3],poserr[3];
 292   pos[0] = fitFcn->GetParameter(1);
 293   pos[1] = fitFcn->GetParameter(2);
 294   poserr[0] = fitFcn->GetParameter(3);
 295   poserr[1] = fitFcn->GetParameter(4);
 296
 297   // Could be replaced by something more robust...
 298   pos[2] = hZ->GetMean();
 299   poserr[2] = hZ->GetRMS();
 300
 301   return new AliESDVertex(pos,poserr);
 302 }